# MODULE:   r.plane for GRASS 5.7; based on r.plane for GRASS 5

#@atexit.register

def report():

    gs.debug("raster report...")

    for rast in tmp_rast:

        info = gs.raster_info(rast)

        gs.debug('%s: %s (%f-%f)' % (rast, info['datatype'],

            info['min'], info['max']))

    # set number of aggregation levels and neighborhood size

    NUM_LEVELS = 4

    NUM_CELLS = 3

    # level, by adding cells to all sides (with one extra on top/right

    # if an odd number of cells needs to be added)

        w = region['w']-math.floor(addx)*region['ewres'],

        e = region['e']+math.ceil(addx)*region['ewres'],

        s = region['s']-math.floor(addy)*region['nsres'],

        n = region['n']+math.ceil(addy)*region['nsres'])

    # arbitrary value used to fill null variance cells

        #TODO: is this the same as circle with radius 2 in arc?

        #TODO: what's happening at the edges here???

        #TODO: is order of ops faithful to aml w.r.t resolution changes?

    if (1 <= gs.debug_level):

        report()

# run in GRID, not ARC

# &run .aml ingrid sd prob bbox

# sd 0.0001 prob 0.05

# 9a - limit to 4 steps, see if that has any significant deterioration of smoothing performance. Should fix

# the problem with islands and headlands - turns out also need to remove water except

# for one cell adjacent to land, and give that a higher uncertainty. See cluster_multiscalesmooth9a_clean_216

#

# Version 9:

# Focus on implementing identical algorithm to directsmooth2 using multiscale method i.e. aggregating by factor of 3

# from already aggregated data, rather than from original resolution each time.

#

# Version 8:

# WARNING: have not yet checked that the additional weighting of the gaussian smoothing is not messing with

# the calculations of variance etc.

# Replaced simple 3x3 aggregation with gaussian smoothing

# Kernel is chosen to give appropriate level of smoothing and produce a fairly accurate approximation of the smoothed

# surface by interpolation of the coarser scale smoothed values

# Details in gaussian.xls on john's laptop

#

# Version 7:

# Further reworking of how the final values are selected - a mean is a candidate if its associated variance

# is less than the mean sample uncertainty, and the mean with the lowest variance among the candidates is the chosen one.

# Implement that in the nested sense by taking the lowest group variance divided by the chi^2 value, and its associated mean variance,

# and if that is lower than the data point variance the

#

# approximate critical value of chi^2/N with N degrees of freedom at 5% level as 1 + 2.45/sqrt(N) + 0.55/N

# for 1% level use 1 + 3.4/sqrt(N) + 2.4/N

#

# Version 6:

# Done from scratch after careful working through of theory.

# ingrid is the (potentially sparse) grid of data to be smoothed and

# interpolated, which can be of different extent to the output

# (resolution is assumed to be the same, could adapt this to relax that)

#

# var can be a constant or a grid

#

# bbox can be either a grid name or the 'xmin ymin xmax ymax' parameters

# for setwindow

# REGETZ NOTES

# - it looks like Ming used sd=0.001 (based on the Arc 'log' file in the

#   topo experimental directory)

#&type NB - using standard deviation as noise specification now, not variance!

multiscalesmooth <- function(ingrid, sd=0.0001 , prob=0.05, bbox) {

    # ingrid: grid to smooth

    # sd: stdev

    # prob: prob

    # set number of levels of aggregation

    # set aggregation factor

    AGG.FACTOR <- 3

    # (with one extra on top/right if an odd number of cells needs to be

    # added)

    max.size <- AGG.FACTOR^NUM.LEVELS

        xmin(ingrid)-floor(addx)*xres(ingrid),

        xmax(ingrid)+ceiling(addx)*xres(ingrid),

        ymin(ingrid)-floor(addy)*yres(ingrid),

        ymax(ingrid)+ceiling(addy)*yres(ingrid)

    # NB - only calculating sample variances here, not variances of estimated means.

    # Also note that v0_bg is an uncertainty, not a sample variance

    # and v1_bg is total variances, but both are labelled as "between-group" to simplify the smoothing

    # weighting for aggregation, based on total variance

    bigvar <- cellStats(v[[1]], max)

    #setcell minof

        hdiff <- z[[i-1]] - disaggregate(z[[i]], 3)

        v.bg <- aggregate(w.prev * hdiff^2, 3, sum) / w

            v.wg <- n - n # zero, but with window and cell size set for us

        chisq <- calc(n.eff, function(n) qchisq(0.05, n-1,

    bigvar  <- bigvar * 10

    #copy z[NUM.LEVELS] hs[NUM.LEVELS]

    #copy v[NUM.LEVELS] vs[NUM.LEVELS]

    #kill z[NUM.LEVELS]

    #kill v[NUM.LEVELS]

    #setcell hs[NUM.LEVELS]

    #setwindow hs[NUM.LEVELS]

    # smooth, refine and combine each layer in turn

    #hs <- z[[NUM.LEVELS+1]]

    #vs <- v[[NUM.LEVELS+1]]

    hs <- z

    vs <- v

    #todo: is this the same as circle with radius 2 in arc???

    #circle2 <- matrix(c(0,1,0, 1,1,1, 0,1,0), nrow=3)

    circle2 <- matrix(c(0,0,1,0,0, 0,1,1,1,0, 1,1,1,1,1, 0,1,1,1,0, 0,0,1,0,0), nrow=5)

        message("Refine from ", j, " to ", j-1)

        # for the first stage where the coarser grid is refined and smoothed, set window to the coarse grid

        #setcell z[j-1]

        #setwindow maxof

        # create smoothed higher resolution versions of z and v_bg, hopefully with no nulls!

        hs.tmp <- focal(disaggregate(hs[[j]], 3), w=circle2, fun=mean,

        vs.tmp <- focal(disaggregate(vs[[j]], 3), w=circle2, fun=mean,

           pad=TRUE, padValue=NA, na.rm=TRUE)

        h_c <- calc(z[[j-1]], function(x) ifelse(is.na(x), 0, x))

        # combine two values using least variance

        hs[[j-1]] <- (h_c/v_c + hs.tmp/vs.tmp ) / (1/v_c + 1/vs.tmp)

        vs[[j-1]] <- 1 / (1/v_c + 1/vs.tmp)

#todo: mimic all the AML setwindow/setcell stuff???

#        hs[[j-1]] <- expand(hs[[j-1]], 4)

#        vs[[j-1]] <- expand(vs[[j-1]], 4)

    result <- crop(stack(hs[[1]], vs[[1]]), extent(ingrid))

    layerNames(result) <- c("hs", "vs")