#!/usr/bin/perl 
use English;
# r.downscale.pl - GRASS GIS script that sets up r.mapcalc to downscale 
# resolution for raster variables such as average temperature or degree-days,
# that tend to be highly correlated to elevation. 
#
# example use: you have low resolution (4 km) avg temperature and elev. layers,
#              and a high res. elev. (100m) layer. On the first pass, build the 
#              regression function for the low res. temp. (y) vs. elev. (x) using wt.d 
#              linear regression (5 x 5 cell local areas). On the second
#              pass, after changing to the higher resolution, apply the model
#              to the higher res elev data, using a weighted average
#              over 5 x 5 cell local areas to provide a smoothing. Wt. values are
#              selectable for both passes, allowing one to emulate 1/distance,
#              1/distance2 and other relationships.
#
# This perl program selects and sends an equation to r.mapcalc, by:
#	1. remove comments: s/#.*$// (or s/#.*?$//mg, if it's in one string)
#       2. drop blank lines; r.mapcalc treats them as EOFs
#	3. substitute input parameters: print in double-quotes
#       4. feed to r.mapcalc's stdin
#
# Developed and released to the public under a GNU public license by 
# Leonard Coop, Oregon State University, coopl@bcc.orst.edu
# Perl to r.mapcalc script orig. by Dan Upper
#
$allparams=<<END;
  ######################################################################
  ##                    CHANGE THE DEFAULTS HERE                      ##
  ##    This is not perl code -- it's key=value pairs, parsed below   ##
  ##                    Use 'NULL' in place of ""                     ##
  ######################################################################
  #
  #parameter =default    # description (all filenames are GRASS raster maplayers) 
  #--------------------------------------------------------------------------------
  calc       =BUILD      # req: function performed (BUILD  or EST)
  input      =H_50       # req: BUILD low res filename to be downscaled
  elev_lo    =us_25m_dem # req: BUILD low res elev filename
  elev_hi    =DEM30m     # req: EST high res elev filename
  output     =HIDD5d     # req: EST output filename with downscaled values
  b1         =BONE       # opt: BUILD & EST slope (x1000) filename
  b0         =BZERO      # opt: BUILD & EST intercept (x100) filename
  rval       =RVAL       # opt: BUILD corr. coeff. (r) (x100) filename 
  resid      =RESID      # opt: BUILD residual (Yi - YHATi) (x10) filename
  ssxy       =SSXY       # opt: BUILD intermediate filename SSXY (sum of squares xy)
  ssxx       =SSXX       # opt: BUILD intermediate filename SSXX (sum of squares xx)
  ssyy       =SSYY       # opt: BUILD intermediate filename SSYY (sum of squares yy)
  d0         =1.0        # opt: BUILD wt value for 5 x 5 center cell      
  d1         =1.0        # opt: BUILD wt value at dist. = 1 cell from center
  d2         =1.0        # opt: BUILD wt value at dist. = 2 cells from center
  e0         =8.0        # opt: EST wt value for 5 x 5 center cell      
  e1         =2.0        # opt: EST wt value at dist. = 1 cell from center
  e2         =1.0        # opt: EST wt value at dist. = 2 cells from center
  #--------------------------------------------------------------------------------
  ######################################################################
END

# now, stuff all the above key=value pairs into a hash
%dscaleparams = ();
@allparams = split '\n', $allparams;
while(@allparams) {
    $_ = shift @allparams;
    s/#.*$//;  # delete comments
    if (($key, $value) = /(\S+)\s*=\s*(\S+)\s*$/) {
        if ($value eq "NULL") { $value = "" }
        $dscaleparams{$key} = $value;
    }
}

if ($#ARGV < 0 or $#ARGV == 0)  {   # try to process options
    if (($#ARGV < 0) or ($ARGV[0] =~ /(-h)|(--help)/))  {
        print<<END;
$0: Perl script to use r.mapcalc to downscale from low to high 
    resolution using 5 x 5 local regressions between variables
    like rainfall or temperature (dep. variable) and elevation 
    (indep. variable). 
    The program expects to be run twice: 
    Once at the lower resolution to build models using calc=BUILD,
    then at a higher resolution using calc=EST to calculate new 
    high resolution values. Example usage:

    GRASS >g.region res=0:02:30.00                             # set lower res
    GRASS >r.downscale.pl input=H_50 elev_lo=lodem calc=BUILD  # build model
    GRASS >g.region res=0:00:09.00                             # set higher res
    GRASS >r.downscale.pl elev_hi=DEM30m output=HI_50 calc=EST # estimate 

    For help:    
    --help or -h        print this message and exit.
    --keys or -k        print the list of keys and exit.
    --defaults or -d    print the list of keys and defaults and exit.
   Other arguments are interpreted as parameter=value pairs

   Notes: For GRASS 4.x, all stats are scaled up: B1 (x1000), B0 (x100), 
          r (x100), residuals (x100).  Keep this in mind when checking
          results
END
#print $#ARGV, "\n";
exit;        
    } elsif ($ARGV[0] =~ /(-k)|(--keys)/) {
        $count=0;
        for $key (sort keys %dscaleparams) {
            print "$key ";
            if (++$count > 5) {
                print "\n";
                $count=0;
            }
        }
        print "\n";
        exit;
        
    } elsif ($ARGV[0] =~ /(-)|(--defaults)/) {
        foreach $_ (split '\n', $allparams) {
            next if /##/;
            print "$_\n";
        }
        exit;
        
    }
}         # end of option processing


# now loop over command-line arguments
while($_ = shift @ARGV) {
    # parse them to key=value,
    ($key, $value) = /(.*)=(.*)/;

    # is the key defined?
    if (defined($dscaleparams{$key})) {
        if ($value eq "NULL") { $value = "" }     # "" can get dropped, so...
        $dscaleparams{$key}=$value;                # replace default with value.
    } else {
        # If not, it's probably a typo.
        die "$0: bad parameter '$key'\n";
    }
}

# DEBUG: print the hash
#for $key (sort keys %dscaleparams) {
#    print "$key\t=$dscaleparams{$key}\n";
#}
#exit;

# input parameters
$x= $dscaleparams{'elev_lo'};
$y= $dscaleparams{'input'};
$hidem= $dscaleparams{'elev_hi'};
$hidd= $dscaleparams{'output'};
$calc= $dscaleparams{'calc'};
$ssxy = $dscaleparams{'ssxy'};
$ssxx = $dscaleparams{'ssxx'};
$ssyy = $dscaleparams{'ssyy'};
$rval = $dscaleparams{'rval'};
$resid = $dscaleparams{'resid'};
$b1 = $dscaleparams{'b1'}; 
$b0 = $dscaleparams{'b0'};
# weights to use for wtd regression - 0 for middle cell 1 for dist.=1 cell, etc
# so for $w0=8, $w1=2 $w2=1, the mid cell gets 8X the weight of a cell 
#     2 cells distant, which would work for 1/D2, and $w0=3, $w1=1.5, 
#     and $w2=1 for 1/D
$w0 = $dscaleparams{'d0'};
$w1 = $dscaleparams{'d1'};
$w2 = $dscaleparams{'d2'};
# and set of wts. for estimation using the model, to provide smoothing 
#  of results, these wts should possibly be different than model building,
# so for $e0=8, $e1=2 $e2=1, the mid cell gets 8X the weight of a 
#   cell 2 cells distant
$e0 = $dscaleparams{'e0'};
$e1 = $dscaleparams{'e1'};
$e2 = $dscaleparams{'e2'};

# constants - save a little typing
$m12 = "[-2,-2]";
$m11 = "[-2,-1]";
$m10 = "[-2,0]";
$m9 = "[-2,1]";
$m8 = "[-2,2]";
$m7 = "[-1,-2]";
$m6 = "[-1,-1]";
$m5 = "[-1,0]";
$m4 = "[-1,1]";
$m3 = "[-1,2]";
$m2 = "[0,-2]";
$m1 = "[0,-1]";
$m0 = "[0,0]";
$p1 = "[0,1]";
$p2 = "[0,2]";
$p3 = "[1,-2]";
$p4 = "[1,-1]";
$p5 = "[1,0]";
$p6 = "[1,1]";
$p7 = "[1,2]";
$p8 = "[2,-2]";
$p9 = "[2,-1]";
$p10 = "[2,0]";
$p11 = "[2,1]";
$p12 = "[2,2]";

# use hash to specify different stored r.mapcalc procedures
%methods=(); 

#  test simple avg
$methods{R}='$prefix ( $x [-1,-1] + $x [0,0] + $x [1,1])/3';
# ratio method DDhires = ELEVhires x DDlores / ELEVlores
$methods{R2}=<<'EOL';
RATIO = $y / $x
EOL

# model building using weighted regression method on 5 x 5 local neighborhood
$methods{BUILD}=<<'EOL';
$ssxy = eval(  eoln
xmean = float(     \ eoln
       if($x $m12<=0,$x $m0,$x $m12) * $w2 +  \ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 +  \ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 +  \ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 +  \ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 +  \ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 +  \ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 +  \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +  \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +  \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +  \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 +  \ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 +  \ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 +  \ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 +  \ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 +  \ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 +  \ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2 \ eoln
       )/float($w2 * 16 + $w1 * 8 + $w0), \ eoln
sm12 = if($x $m12 <=0,$x $m0 ,$x $m12 ) - xmean, \ eoln
sm11 = if($x $m11 <=0,$x $m0 ,$x $m11 ) - xmean, \ eoln
sm10 = if($x $m10 <=0,$x $m0 ,$x $m10 ) - xmean, \ eoln
sm9 = if($x $m9 <=0,$x $m0 ,$x $m9 ) - xmean, \ eoln
sm8 = if($x $m8 <=0,$x $m0 ,$x $m8 ) - xmean, \ eoln
sm7 = if($x $m7 <=0,$x $m0 ,$x $m7 ) - xmean, \ eoln
sm6 = if($x $m6 <=0,$x $m0 ,$x $m6 ) - xmean, \ eoln
sm5 = if($x $m5 <=0,$x $m0 ,$x $m5 ) - xmean, \ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4 ) - xmean, \ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3 ) - xmean, \ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean, \ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean, \ eoln
s0 = $x [0,0] - xmean, \ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1 ) - xmean, \ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2 ) - xmean, \ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3 ) - xmean, \ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4 ) - xmean, \ eoln
s5 = if($x $p5 <=0,$x $m0 ,$x $p5 ) - xmean, \ eoln
s6 = if($x $p6 <=0,$x $m0 ,$x $p6 ) - xmean, \ eoln
s7 = if($x $p7 <=0,$x $m0 ,$x $p7 ) - xmean, \ eoln
s8 = if($x $p8 <=0,$x $m0 ,$x $p8 ) - xmean, \ eoln
s9 = if($x $p9 <=0,$x $m0 ,$x $p9 ) - xmean, \ eoln
s10 = if($x $p10 <=0,$x $m0 ,$x $p10 ) - xmean, \ eoln
s11 = if($x $p11 <=0,$x $m0 ,$x $p11 ) - xmean, \ eoln
s12 = if($x $p12 <=0,$x $m0 ,$x $p12 ) - xmean, \ eoln
       (sm12 * if($y $m12 <=0,$y $m0 ,$y $m12 ) * $w2 + \ eoln
       sm11 * if($y $m11 <=0,$y $m0 ,$y $m11 ) * $w2 + \ eoln
       sm10 * if($y $m10 <=0,$y $m0 ,$y $m10 ) * $w2 + \ eoln
       sm9 * if($y $m9 <=0,$y $m0 ,$y $m9 ) * $w2 + \ eoln
       sm8 * if($y $m8 <=0,$y $m0 ,$y $m8 ) * $w2 + \ eoln
       sm7 * if($y $m7 <=0,$y $m0 ,$y $m7 ) * $w2 + \ eoln
       sm6 * if($y $m6 <=0,$y $m0 ,$y $m6 ) * $w1 + \ eoln
       sm5 * if($y $m5 <=0,$y $m0 ,$y $m5 ) * $w1 + \ eoln
       sm4 * if($y $m4 <=0,$y $m0 ,$y $m4 ) * $w1 + \ eoln
       sm3 * if($y $m3 <=0,$y $m0 ,$y $m3 ) * $w2 + \ eoln
       sm2 * if($y $m2 <=0,$y $m0 ,$y $m2 ) * $w2 + \ eoln
       sm1 * if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 + \ eoln
       s0 * ($y $m0 ) * $w0 + \ eoln
       s1 * if($y $p1 <=0,$y $m0 ,$y $p1 ) * $w1 + \ eoln
       s2 * if($y $p2 <=0,$y $m0 ,$y $p2 ) * $w2 + \ eoln
       s3 * if($y $p3 <=0,$y $m0 ,$y $p3 ) * $w2 + \ eoln
       s4 * if($y $p4 <=0,$y $m0 ,$y $p4 ) * $w1 + \ eoln
       s5 * if($y $p5 <=0,$y $m0 ,$y $p5 ) * $w1 + \ eoln
       s6 * if($y $p6 <=0,$y $m0 ,$y $p6 ) * $w1 + \ eoln
       s7 * if($y $p7 <=0,$y $m0 ,$y $p7 ) * $w2 + \ eoln
       s8 * if($y $p8 <=0,$y $m0 ,$y $p8 ) * $w2 + \ eoln
       s9 * if($y $p9 <=0,$y $m0 ,$y $p9 ) * $w2 + \ eoln
       s10 * if($y $p10 <=0,$y $m0 ,$y $p10 ) * $w2 + \ eoln
       s11 * if($y $p11 <=0,$y $m0 ,$y $p11 ) * $w2 + \ eoln
       s12 * if($y $p12 <=0,$y $m0 ,$y $p12 ) * $w2) \ eoln
)
CR
$ssxx = eval( eoln
xmean = float(     \ eoln
       if($x $m12<=0,$x $m0,$x $m12) * $w2 +  \ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 +  \ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 +  \ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 +  \ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 +  \ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 +  \ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 +  \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +  \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +  \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +  \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 +  \ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 +  \ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 +  \ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 +  \ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 +  \ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 +  \ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2 \ eoln
       )/float($w2 * 16 + $w1 * 8 + $w0), \ eoln
sm12 = if($x $m12 <=0,$x $m0 ,$x $m12 ) - xmean, \ eoln
sm11 = if($x $m11 <=0,$x $m0 ,$x $m11 ) - xmean, \ eoln
sm10 = if($x $m10 <=0,$x $m0 ,$x $m10 ) - xmean, \ eoln
sm9 = if($x $m9 <=0,$x $m0 ,$x $m9 ) - xmean, \ eoln
sm8 = if($x $m8 <=0,$x $m0 ,$x $m8 ) - xmean, \ eoln
sm7 = if($x $m7 <=0,$x $m0 ,$x $m7 ) - xmean, \ eoln
sm6 = if($x $m6 <=0,$x $m0 ,$x $m6 ) - xmean, \ eoln
sm5 = if($x $m5 <=0,$x $m0 ,$x $m5 ) - xmean, \ eoln
sm4 = if($x $m4 <=0,$x $m0 ,$x $m4 ) - xmean, \ eoln
sm3 = if($x $m3 <=0,$x $m0 ,$x $m3 ) - xmean, \ eoln
sm2 = if($x $m2 <=0,$x $m0 ,$x $m2) - xmean, \ eoln
sm1 = if($x $m1 <=0,$x $m0 ,$x $m1) - xmean, \ eoln
s0 = $x [0,0] - xmean, \ eoln
s1 = if($x $p1 <=0,$x $m0 ,$x $p1 ) - xmean, \ eoln
s2 = if($x $p2 <=0,$x $m0 ,$x $p2 ) - xmean, \ eoln
s3 = if($x $p3 <=0,$x $m0 ,$x $p3 ) - xmean, \ eoln
s4 = if($x $p4 <=0,$x $m0 ,$x $p4 ) - xmean, \ eoln
s5 = if($x $p5 <=0,$x $m0 ,$x $p5 ) - xmean, \ eoln
s6 = if($x $p6 <=0,$x $m0 ,$x $p6 ) - xmean, \ eoln
s7 = if($x $p7 <=0,$x $m0 ,$x $p7 ) - xmean, \ eoln
s8 = if($x $p8 <=0,$x $m0 ,$x $p8 ) - xmean, \ eoln
s9 = if($x $p9 <=0,$x $m0 ,$x $p9 ) - xmean, \ eoln
s10 = if($x $p10 <=0,$x $m0 ,$x $p10 ) - xmean, \ eoln
s11 = if($x $p11 <=0,$x $m0 ,$x $p11 ) - xmean, \ eoln
s12 = if($x $p12 <=0,$x $m0 ,$x $p12 ) - xmean, \ eoln
      (sm12 * sm12 * $w2 + sm11 * sm11 * $w2 + sm10 * sm10 * $w2 + \ eoln
      sm9 * sm9 * $w2 + sm8 * sm8 * $w2 + sm7 * sm7 * $w2 + \ eoln
      sm6 * sm6 * $w1 + sm5 * sm5 * $w1 + sm4 * sm4 * $w1 + \ eoln
      sm3 * sm3 * $w2 + sm2 * sm2 * $w2 + sm1 * sm1 * $w1 + \ eoln
      s0 * s0 * $w0 + \
      s1 * s1 * $w1 + s2 * s2 * $w2 + s3 * s3 * $w2 + s4 * s4 * $w1 + \ eoln
      s5 * s5 * $w1 + s6 * s6 * $w1 + s7 * s7 * $w2 + s8 * s8 * $w2 + \ eoln
      s9 * s9 * $w2 + s10 * s10 * $w2 + s11 * s11 * $w2 + s12 * s12 * $w2) \ eoln
)
CR 
$b1 = 1000 * (float($ssxy) / float($ssxx)) 
CR
$b0 = eval( eoln
xmean = float(     \ eoln
       if($x $m12<=0,$x $m0,$x $m12) * $w2 +  \ eoln
       if($x $m11<=0,$x $m0,$x $m11) * $w2 +  \ eoln
       if($x $m10<=0,$x $m0,$x $m10) * $w2 +  \ eoln
       if($x $m9<=0,$x $m0,$x $m9) * $w2 +  \ eoln
       if($x $m8<=0,$x $m0,$x $m8) * $w2 +  \ eoln
       if($x $m7<=0,$x $m0,$x $m7) * $w2 +  \ eoln
       if($x $m6<=0,$x $m0,$x $m6) * $w1 +  \ eoln
       if($x $m5<=0,$x $m0,$x $m5) * $w1 +  \ eoln
       if($x $m4<=0,$x $m0,$x $m4) * $w1 +  \ eoln
       if($x $m3<=0,$x $m0,$x $m3) * $w2 +  \ eoln
       if($x $m2<=0,$x $m0,$x $m2) * $w2 +  \ eoln
       if($x $m1<=0,$x $m0,$x $m1) * $w1 +  \ eoln
       $x $m0 * $w0 +  \ eoln
       if($x $p1<=0,$x $m0,$x $p1) * $w1 +  \ eoln
       if($x $p2<=0,$x $m0,$x $p2) * $w2 +  \ eoln
       if($x $p3<=0,$x $m0,$x $p3) * $w2 +  \ eoln
       if($x $p4<=0,$x $m0,$x $p4) * $w1 +  \ eoln
       if($x $p5<=0,$x $m0,$x $p5) * $w1 +  \ eoln
       if($x $p6<=0,$x $m0,$x $p6) * $w1 +  \ eoln
       if($x $p7<=0,$x $m0,$x $p7) * $w2 +  \ eoln
       if($x $p8<=0,$x $m0,$x $p8) * $w2 +  \ eoln
       if($x $p9<=0,$x $m0,$x $p9) * $w2 +  \ eoln
       if($x $p10<=0,$x $m0,$x $p10) * $w2 +  \ eoln
       if($x $p11<=0,$x $m0,$x $p11) * $w2 +  \ eoln
       if($x $p12<=0,$x $m0,$x $p12) * $w2 \ eoln
       )/float($w2 * 16 + $w1 * 8 + $w0), \ eoln
ymean = float(     \ eoln
       if($y $m12 <=0,$y $m0 ,$y $m12 ) * $w2 +  \ eoln
       if($y $m11 <=0,$y $m0 ,$y $m11 ) * $w2 +  \ eoln
       if($y $m10 <=0,$y $m0 ,$y $m10 ) * $w2 +  \ eoln
       if($y $m9 <=0,$y $m0 ,$y $m9 ) * $w2 +  \ eoln
       if($y $m8 <=0,$y $m0 ,$y $m8 ) * $w2 +  \ eoln
       if($y $m7 <=0,$y $m0 ,$y $m7 ) * $w2 +  \ eoln
       if($y $m6 <=0,$y $m0 ,$y $m6 ) * $w1 +  \ eoln
       if($y $m5 <=0,$y $m0 ,$y $m5 ) * $w1 +  \ eoln
       if($y $m4 <=0,$y $m0 ,$y $m4 ) * $w1 +  \ eoln
       if($y $m3 <=0,$y $m0 ,$y $m3 ) * $w2 +  \ eoln
       if($y $m2 <=0,$y $m0 ,$y $m2 ) * $w2 +  \ eoln
       if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 +  \ eoln
       ($y $m0 ) * $w0 +  \ eoln
       if($y $p1 <=0,$y $m0 ,$y $p1 ) * $w1 +  \ eoln
       if($y $p2 <=0,$y $m0 ,$y $p2 ) * $w2 +  \ eoln
       if($y $p3 <=0,$y $m0 ,$y $p3 ) * $w2 +  \ eoln
       if($y $p4 <=0,$y $m0 ,$y $p4 ) * $w1 +  \ eoln
       if($y $p5 <=0,$y $m0 ,$y $p5 ) * $w1 +  \ eoln
       if($y $p6 <=0,$y $m0 ,$y $p6 ) * $w1 +  \ eoln
       if($y $p7 <=0,$y $m0 ,$y $p7 ) * $w2 +  \ eoln
       if($y $p8 <=0,$y $m0 ,$y $p8 ) * $w2 +  \ eoln
       if($y $p9 <=0,$y $m0 ,$y $p9 ) * $w2 +  \ eoln
       if($y $p10 <=0,$y $m0 ,$y $p10 ) * $w2 +  \ eoln
       if($y $p11 <=0,$y $m0 ,$y $p11 ) * $w2 +  \ eoln
       if($y $p12 <=0,$y $m0 ,$y $p12 ) * $w2 \ eoln
       )/float($w2 * 16 + $w1 * 8 + $w0), \ eoln
((ymean - ((float($b1) /1000) * xmean))) * 100  \ eoln
)
CR
$ssyy = eval( eoln
ymean = float(     \ eoln
       if($y $m12 <=0,$y $m0 ,$y $m12 ) * $w2 +  \ eoln
       if($y $m11 <=0,$y $m0 ,$y $m11 ) * $w2 +  \ eoln
       if($y $m10 <=0,$y $m0 ,$y $m10 ) * $w2 +  \ eoln
       if($y $m9 <=0,$y $m0 ,$y $m9 ) * $w2 +  \ eoln
       if($y $m8 <=0,$y $m0 ,$y $m8 ) * $w2 +  \ eoln
       if($y $m7 <=0,$y $m0 ,$y $m7 ) * $w2 +  \ eoln
       if($y $m6 <=0,$y $m0 ,$y $m6 ) * $w1 +  \ eoln
       if($y $m5 <=0,$y $m0 ,$y $m5 ) * $w1 +  \ eoln
       if($y $m4 <=0,$y $m0 ,$y $m4 ) * $w1 +  \ eoln
       if($y $m3 <=0,$y $m0 ,$y $m3 ) * $w2 +  \ eoln
       if($y $m2 <=0,$y $m0 ,$y $m2 ) * $w2 +  \ eoln
       if($y $m1 <=0,$y $m0 ,$y $m1 ) * $w1 +  \ eoln
       ($y $m0 ) * $w0 +  \ eoln
       if($y $p1 <=0,$y $m0 ,$y $p1 ) * $w1 +  \ eoln
       if($y $p2 <=0,$y $m0 ,$y $p2 ) * $w2 +  \ eoln
       if($y $p3 <=0,$y $m0 ,$y $p3 ) * $w2 +  \ eoln
       if($y $p4 <=0,$y $m0 ,$y $p4 ) * $w1 +  \ eoln
       if($y $p5 <=0,$y $m0 ,$y $p5 ) * $w1 +  \ eoln
       if($y $p6 <=0,$y $m0 ,$y $p6 ) * $w1 +  \ eoln
       if($y $p7 <=0,$y $m0 ,$y $p7 ) * $w2 +  \ eoln
       if($y $p8 <=0,$y $m0 ,$y $p8 ) * $w2 +  \ eoln
       if($y $p9 <=0,$y $m0 ,$y $p9 ) * $w2 +  \ eoln
       if($y $p10 <=0,$y $m0 ,$y $p10 ) * $w2 +  \ eoln
       if($y $p11 <=0,$y $m0 ,$y $p11 ) * $w2 +  \ eoln
       if($y $p12 <=0,$y $m0 ,$y $p12 ) * $w2 \ eoln
       )/float($w2 * 16 + $w1 * 8 + $w0), \ eoln
sm12 = if($y $m12 <=0,$y $m0 ,$y $m12 ) - ymean, \ eoln
sm11 = if($y $m11 <=0,$y $m0 ,$y $m11 ) - ymean, \ eoln
sm10 = if($y $m10 <=0,$y $m0 ,$y $m10 ) - ymean, \ eoln
sm9 = if($y $m9 <=0,$y $m0 ,$y $m9 ) - ymean, \ eoln
sm8 = if($y $m8 <=0,$y $m0 ,$y $m8 ) - ymean, \ eoln
sm7 = if($y $m7 <=0,$y $m0 ,$y $m7 ) - ymean, \ eoln
sm6 = if($y $m6 <=0,$y $m0 ,$y $m6 ) - ymean, \ eoln
sm5 = if($y $m5 <=0,$y $m0 ,$y $m5 ) - ymean, \ eoln
sm4 = if($y $m4 <=0,$y $m0 ,$y $m4 ) - ymean, \ eoln
sm3 = if($y $m3 <=0,$y $m0 ,$y $m3 ) - ymean, \ eoln
sm2 = if($y $m2 <=0,$y $m0 ,$y $m2 ) - ymean, \ eoln
sm1 = if($y $m1 <=0,$y $m0 ,$y $m1 ) - ymean, \ eoln
s0 = if($y $m0 <=0,$y $m0 ,$y $m0 ) - ymean, \ eoln
s1 = if($y $p1 <=0,$y $m0 ,$y $p1 ) - ymean, \ eoln
s2 = if($y $p2 <=0,$y $m0 ,$y $p2 ) - ymean, \ eoln
s3 = if($y $p3 <=0,$y $m0 ,$y $p3 ) - ymean, \ eoln
s4 = if($y $p4 <=0,$y $m0 ,$y $p4 ) - ymean, \ eoln
s5 = if($y $p5 <=0,$y $m0 ,$y $p5 ) - ymean, \ eoln
s6 = if($y $p6 <=0,$y $m0 ,$y $p6 ) - ymean, \ eoln
s7 = if($y $p7 <=0,$y $m0 ,$y $p7 ) - ymean, \ eoln
s8 = if($y $p8 <=0,$y $m0 ,$y $p8 ) - ymean, \ eoln
s9 = if($y $p9 <=0,$y $m0 ,$y $p9 ) - ymean, \ eoln
s10 = if($y $p10 <=0,$y $m0 ,$y $p10 ) - ymean, \ eoln
s11 = if($y $p11 <=0,$y $m0 ,$y $p11 ) - ymean, \ eoln
s12 = if($y $p12 <=0,$y $m0 ,$y $p12 ) - ymean, \ eoln
      (sm12 * sm12 * $w2 + sm11 * sm11 * $w2 + sm10 * sm10 * $w2 + \ eoln
        sm9 * sm9 * $w2 + sm8 * sm8 * $w2 + sm7 * sm7 * $w2 + \ eoln
        sm6 * sm6 * $w1 + sm5 * sm5 * $w1 + sm4 * sm4 * $w1 + \ eoln
        sm3 * sm3 * $w2 + sm2 * sm2 * $w2 + sm1 * sm1 * $w1 + \ eoln
      s0 * s0 * $w0 + \ eoln
      s1 * s1 * $w1 + s2 * s2 * $w2 + s3 * s3 * $w2 + s4 * s4 * $w1 + \ eoln
      s5 * s5 * $w1 + s6 * s6 * $w1 + s7 * s7 * $w2 + s8 * s8 * $w2 + \ eoln
      s9 * s9 * $w2 + s10 * s10 * $w2 + s11 * s11 * $w2 + s12 * s12 * $w2) \ eoln
) 
CR
$rval = ((float($ssxy) / 1000)/ sqrt((float($ssxx)/1000) * (float($ssyy)/1000))) * 100 \ eoln
CR
$resid = ($y - (float($b0) /100  + $x * (float($b1) / 1000))) \ eoln
EOL

#calculate estimated DDs from regression params - simple
$methods{EST3}=<<'EOL';
$hidd = (float($b0) / 100) + ($hidem * (float($b1) / 1000)) \ eoln
EOL

#calculate estimated DDs from regression params - weighted avg neighbor 5x5 only
# weighting approximate discrete version of 1/distance-squared 8-2-1 for distance
# 1-2-3
$methods{EST}=<<'EOL';
$hidd = float(if($b0 $m12 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m12) / 100) + ($hidem * (float($b1 $m12) / 1000))) * $e2 + \ eoln
        if($b0 $m11 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m11) / 100) + ($hidem * (float($b1 $m11) / 1000))) * $e2 + \ eoln
        if($b0 $m10 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m10) / 100) + ($hidem * (float($b1 $m10) / 1000))) * $e2 + \ eoln
        if($b0 $m9 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m9) / 100) + ($hidem * (float($b1 $m9) / 1000))) * $e2 + \ eoln
        if($b0 $m8 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m8) / 100) + ($hidem * (float($b1 $m8) / 1000))) * $e2 + \ eoln
        if($b0 $m7 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m7) / 100) + ($hidem * (float($b1 $m7) / 1000))) * $e2 + \ eoln
        if($b0 $m6 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m6) / 100) + ($hidem * (float($b1 $m6) / 1000))) * $e1 + \ eoln
        if($b0 $m5 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m5) / 100) + ($hidem * (float($b1 $m5) / 1000))) * $e1 + \ eoln
        if($b0 $m4 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m4) / 100) + ($hidem * (float($b1 $m4) / 1000))) * $e1 + \ eoln
        if($b0 $m3 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m3) / 100) + ($hidem * (float($b1 $m3) / 1000))) * $e2 + \ eoln
        if($b0 $m2 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m2) / 100) + ($hidem * (float($b1 $m2) / 1000))) * $e2 + \ eoln
        if($b0 $m1 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $m1) / 100) + ($hidem * (float($b1 $m1) / 1000))) * $e1 + \ eoln
              ((float($b0 $m0) / 100) + ($hidem * (float($b1 $m0) / 1000))) * $e0 + \ eoln
        if($b0 $p1 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p1) / 100) + ($hidem * (float($b1 $p1) / 1000))) * $e1 + \ eoln
        if($b0 $p2 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p2) / 100) + ($hidem * (float($b1 $p2) / 1000))) * $e2 + \ eoln
        if($b0 $p3 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p3) / 100) + ($hidem * (float($b1 $p3) / 1000))) * $e2 + \ eoln
        if($b0 $p4 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p4) / 100) + ($hidem * (float($b1 $p4) / 1000))) * $e1 + \ eoln
        if($b0 $p5 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p5) / 100) + ($hidem * (float($b1 $p5) / 1000))) * $e1 + \ eoln
        if($b0 $p6 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p6) / 100) + ($hidem * (float($b1 $p6) / 1000))) * $e1 + \ eoln
        if($b0 $p7 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p7) / 100) + ($hidem * (float($b1 $p7) / 1000))) * $e2 + \ eoln
        if($b0 $p8 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p8) / 100) + ($hidem * (float($b1 $p8) / 1000))) * $e2 + \ eoln
        if($b0 $p9 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p9) / 100) + ($hidem * (float($b1 $p9) / 1000))) * $e2 + \ eoln
        if($b0 $p10 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p10) / 100) + ($hidem * (float($b1 $p10) / 1000))) * $e2 + \ eoln
        if($b0 $p11<=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p11) / 100) + ($hidem * (float($b1 $p11) / 1000))) * $e2 + \ eoln
        if($b0 $p12 <=0,($b0 $m0 / 100) + ($hidem * ($b1 $m0 / 1000)), \ eoln
              (float($b0 $p12) / 100) + ($hidem * (float($b1 $p12) / 1000))) * $e2) \ eoln
        /float($e2 * 16 + $e1 * 8 + $e0) \ eoln
EOL

# clean up r.mapcalc input template: $methods{$calc}
$methods{$calc} =~ s/#.*?$//mg;
$methods{$calc} =~ s/\n//mg;
$methods{$calc} =~ s/\s\s//mg;
$methods{$calc} =~ s/CR/\n/mg;
$methods{$calc} =~ s/eoln/\\ eoln\n/mg;

# ----- set up r.mapcalc pipe -----
# sed required because perl cannot allow final slash marks expected by r.mapcalc
open RMAPCALC, "| sed -e 's|eoln|\\\\|g' | r.mapcalc "
  or die "Can't launch r.mapcalc";

 # trickery to interpolate the parameters into the command string
    # i.e. to replace "$var" in with the current value of $var in the string
   $calcstring =  eval '"' . "$methods{$calc}" . '\n"';

# use for debug:
# print "$calcstring\n";

# ----- actual call to r.mapcalc passing in the calc method -----
print RMAPCALC $calcstring;
close RMAPCALC;
# end of r.downscale.pl