&args ingrid sd prob bbox:rest

/* 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

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

/* set up chisq parameters

&sv chisqa = [calc 2.807 - 0.6422 * [log10 %prob% ] - 3.410 * %prob% ** 0.3411 ]

&sv chisqb = [calc -5.871 - 3.675 * [log10 %prob% ] + 4.690 * %prob% ** 0.3377 ]

&type chisq parameters %chisqa% %chisqb%

setcell %ingrid%

/* work out maximum of ingrid and bbox extents

setwindow [unquote %bbox%] %ingrid%

bboxgrid = 1

setwindow maxof

workgrid = %ingrid% + bboxgrid

setwindow workgrid

kill workgrid

/* naming:

/*  h - the value being smoothed/interpolated

/*  vg - total variance of group of data, or of individual measurement

/*  v_bg - variance between groups

/*  v_wg - variance within groups

/*  wa - weighting for aggregation, based on total variance

/*  vm - variance of the calculated mean

/*  mv - mean of finer scale variances

/*  n - effective number of measurements

/* 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

h0 = %ingrid%

v0 = con(^ isnull(h0), sqr(%sd%))

vg0 = v0

w0 = con(isnull(v0), 0, 1.0 / v0)

wsq0 = sqr(w0)

n0 = con(^ isnull(h0), 1, 0)

&describe v0

&sv bigvar %grd$zmax%

setcell minof

/* aggregate to broader scales

&sv i 1

&sv done .false.

&describe h0

&do &until %done%

  &sv j [calc %i% - 1]

  &type Aggregate from %j% to %i%

  &describe h%j%

  &sv cell3 [calc %grd$dx% * 3]

  &describe h0

  &sv nx0 [round [calc %grd$xmin% / %cell3% - 0.5]]

  &sv ny0 [round [calc %grd$ymin% / %cell3% - 0.5]]

  &sv nx1 [round [calc %grd$xmax% / %cell3% + 0.5]]

  &sv ny1 [round [calc %grd$ymax% / %cell3% + 0.5]]

  &sv x0 [calc ( %nx0% - 0.5 ) * %cell3%]

  &sv y0 [calc ( %ny0% - 0.5 ) * %cell3%]

  &sv x1 [calc ( %nx1% + 0.5 ) * %cell3%]

  &sv y1 [calc ( %ny1% + 0.5 ) * %cell3%]

  setwindow %x0% %y0% %x1% %y1%

  w%i% = aggregate(w%j%, 3, sum)

  wsq%i% = aggregate(wsq%j%, 3, sum)

  n%i% = aggregate(n%j%, 3, sum)

  neff%i% = w%i% * w%i% / wsq%i%

  h%i% = aggregate(w%j% * h%j%, 3, sum) / w%i%

  vbg%i% = aggregate(w%j% * sqr(h%j% - h%i%), 3, sum) / w%i%

  &if %i% eq 1 &then vwg%i% = n%i% - n%i% /* zero, but with window and cell size set for us

  &else vwg%i% = aggregate(w%j% * vg%j%, 3, sum) / w%i%

  vg%i% = vbg%i% + vwg%i%

  vm%i% = 1.0 / w%i%

  mv%i% = n%i% / w%i%

  chisq%i% = 1 + %chisqa% / sqrt(neff%i% - 1) + %chisqb% / (neff%i% - 1)

  v%i% = con(vg%i% / chisq%i% < mv%i%, vm%i%, vg%i%)

  /* remove everything except h%i% and v%i%

  kill w%j%

  kill wsq%j%

  kill n%j%

  kill neff%i%

  kill vbg%i%

  kill vwg%i%

  kill vg%j%

  kill vm%i%

  kill mv%i%

  kill chisq%i%

  &sv done %i% eq 4

  &sv i [calc %i% + 1]

&end

&sv maxstep [calc %i% - 1]

&sv bigvar [calc %bigvar% * 10]

kill w%maxstep%

kill wsq%maxstep%

kill n%maxstep%

kill vg%maxstep%

/* smooth, refine and combine each layer in turn

copy h%maxstep% hs%maxstep%

copy v%maxstep% vs%maxstep%

kill h%maxstep%

kill v%maxstep%

setcell hs%maxstep%

setwindow hs%maxstep%

&do j := %maxstep% &to 1 &by -1

  &sv i [calc %j% - 1]

  &type Refine from %j% to %i%

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

  setcell h%i%

  setwindow maxof

  /* create smoothed higher resolution versions of h and v_bg, hopefully with no nulls!

  hs%j%_%i% = focalmean(hs%j%, circle, 2)

  vs%j%_%i% = focalmean(vs%j%, circle, 2)

  setcell h%i%

  &describe h%i%

  &sv cellsize %grd$dx%

  &describe bboxgrid

  setwindow [calc %grd$xmin% - 4 * %cellsize%] [calc %grd$ymin% - 4 * %cellsize%] [calc %grd$xmax% + 4 * %cellsize%] [calc %grd$ymax% + 4 * %cellsize%] h%i%

  /* create no-null version of finer h and v

  h%i%_c = con(isnull(h%i%), 0, h%i%)

  v%i%_c = con(isnull(v%i%), %bigvar%, v%i%)

  /* combine two values using least variance

  hs%i% = (h%i%_c / v%i%_c + hs%j%_%i% / vs%j%_%i% ) / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)

  vs%i% = 1 / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)

  kill v%i%_c  

  kill h%i%_c

  kill v%i%

  kill h%i%

  kill vs%j%_%i%

  kill hs%j%_%i%

  kill hs%j%

  kill vs%j%

&end

/* result is hs0, with variance vs0

kill bboxgrid