#!/usr/bin/env perl

# USER INPUT #######################################################################################
$inputframe=$ARGV[0];
$inputfile=$ARGV[1];
if ($inputfile eq "" or $inputframe eq ""){die "wrong arguments";}

#####################################################################################################

#process sampling frame
open (DAT, "$inputframe") or die "missing sampling frame"; 
$header = (<DAT>); @sample_frame=<DAT>; close (DAT); foreach (@sample_frame) {chomp;}
$maparea=0; $mappixels=0; $nstrata=@sample_frame;
for ($i=0; $i<$nstrata; $i++){
($strata,$area,$count)=split('\t',@sample_frame[$i]);
if ($strata==0 or $strata eq "" or $area==0 or $area eq "" or $count==0 or $count eq ""){die "incorrect sampling frame";}
$matrix[$strata][0]=$strata; $matrix[$strata][1]=$area; $matrix[$strata][2]=$count; 
$maparea=$maparea+$area; $mappixels=$mappixels+$count;
}
print "N of strata: $nstrata\n";

#read interpretation table
open (DAT, "$inputfile") or die "missing sampling table"; $header = (<DAT>); @samplesi=<DAT>; close (DAT); foreach (@samplesi) {chomp;}
($name,$ext)=split('\.',$inputfile);
open (OUT, ">$name\_report.txt");

#check sampling frame
&sampframe;

$outype="TREE_COVER_2000"; print "$outype\n";
@samples=();
for $sample (@samplesi){
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
$proportion=$TC2000; $type=$FT2000;
push(@samples,"$id,$strata,$proportion,$type");
}
&process;

$outype="TREE_COVER_2016"; print "$outype\n";
@samples=();
for $sample (@samplesi){
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
$proportion=$TC2016; $type=$FT2016;
push(@samples,"$id,$strata,$proportion,$type");
}
&process;

$outype="GROSS_LOSS"; print "$outype\n";
@samples=();
for $sample (@samplesi){
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
if (($MAX-$TC2016)>($TC2000-$MIN)){$proportion=$MAX-$TC2016;} else {$proportion=$TC2000-$MIN;}
$type=$TLOSS;
push(@samples,"$id,$strata,$proportion,$type");
}
&process;

$outype="GROSS_GAIN"; print "$outype\n";
@samples=();
for $sample (@samplesi){
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
if (($MAX-$TC2016)>($TC2000-$MIN)){$proportion=$MAX-$TC2016;} else {$proportion=$TC2000-$MIN;}
$proportion=$TC2016-$TC2000+$proportion; $type=$FT2016;
push(@samples,"$id,$strata,$proportion,$type");
}
&process;

$outype="NET_CHANGE"; print "$outype\n";
@samples=();
for $sample (@samplesi){
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
$proportion=$TC2016-$TC2000; $type=$FT2000;
push(@samples,"$id,$strata,$proportion,$type");
}
&process;

close (OUT);

sub process {
#fill the matrix
@matrix1=(); @matrix2=(); @matrix3=(); @matrix4=(); @matrix5=(); @matrix6=(); 
for ($i=1; $i<=$nstrata; $i++){$matrix1[$strata][0]=0;}
%types=(); $lastkey=0;
for $sample (@samples){
($id,$strata,$proportion,$type)=split(',',$sample);
#check consistency
if ($strata<1 or $strata>$nstrata) {die "$sample";}
#check type 
$check=0; while (($key, $value) = each(%types)){ if ($value eq $type){$check=1; $typekey=$key;}} 
if ($check==0){$typekey=$lastkey; $types{$lastkey}=$type; $lastkey++;}
#fill values
$matrix1[$strata][0]=$matrix1[$strata][0]+1;
$matrix1[$strata][1]=$matrix1[$strata][1]+($proportion/100);
$matrix1[$strata][2]=$matrix1[$strata][2]+($proportion/100)*($proportion/100);
$matrix2[$strata][$typekey]=$matrix2[$strata][$typekey]+($proportion/100);
$matrix3[$strata][$typekey]=$matrix3[$strata][$typekey]+($proportion/100)*($proportion/100);
if ($proportion>=50){$matrix4[$strata][0]=$matrix4[$strata][0]+1;}
}

#area total
$prop=0; for ($i=1; $i<=$nstrata; $i++){
if ($matrix1[$i][0]>0){$prop=$prop+($matrix[$i][2]/$mappixels)*$matrix1[$i][1]/$matrix1[$i][0];}}
$totarea=$prop*$maparea;

#area by type
@typearea=();
for ($k=0; $k<$lastkey; $k++){
$prop=0; for ($i=1; $i<=$nstrata; $i++){
if ($matrix1[$i][0]>0){$prop=$prop+($matrix[$i][2]/$mappixels)*$matrix2[$i][$k]/$matrix1[$i][0];}}
$typearea[$k]=$prop*$maparea;
}

#error total
for ($i=1; $i<=$nstrata; $i++){if ($matrix1[$i][0]>0){$matrix1[$i][3]=($matrix1[$i][2]-($matrix1[$i][1]*$matrix1[$i][1])/$matrix1[$i][0])/($matrix1[$i][0]-1);}}
for ($i=1; $i<=$nstrata; $i++){if ($matrix1[$i][0]>0){$matrix1[$i][4]=$matrix[$i][2]*$matrix[$i][2]*(1-$matrix1[$i][0]/$matrix[$i][2])*($matrix1[$i][3]/$matrix1[$i][0]);}}
$t1=0; for ($i=1; $i<=$nstrata; $i++){$t1=$t1+$matrix1[$i][4];}
$setot=$maparea*sqrt($t1/($mappixels*$mappixels));

#error by type
@setype=();
for ($k=0; $k<$lastkey; $k++){
for ($i=1; $i<=$nstrata; $i++){if ($matrix1[$i][0]>0){$matrix5[$i][0]=($matrix3[$i][$k]-($matrix2[$i][$k]*$matrix2[$i][$k])/$matrix1[$i][0])/($matrix1[$i][0]-1);}}
for ($i=1; $i<=$nstrata; $i++){if ($matrix1[$i][0]>0){$matrix6[$i][0]=$matrix[$i][2]*$matrix[$i][2]*(1-$matrix1[$i][0]/$matrix[$i][2])*($matrix5[$i][0]/$matrix1[$i][0]);}}
$t1=0; for ($i=1; $i<=$nstrata; $i++){$t1=$t1+$matrix6[$i][0];}
$setype[$k]=$maparea*sqrt($t1/($mappixels*$mappixels));
}

#print report
print OUT "\n$outype\n";
print OUT "Type\tarea\tSE\tConf95%\n";
$UI=$setot*1.96; 
print OUT "Total\t$totarea\t$setot\t$UI\n";
for ($k=0; $k<$lastkey; $k++){
if ($types{$k} ne "0"){
$UI=$setype[$k]*1.96;
print OUT "$types{$k}\t$typearea[$k]\t$setype[$k]\t$UI\n";
}}
}

sub sampframe {
#fill the matrix
@matrix1=();
for $sample (@samplesi){
(@t)=split('\t',$sample); $num=@t; if ($num != 29){die "wrong sample\n$sample\n";}
($id,$xcoord,$ycoord,$strata,$country,$stype,$maploss,$mapyear,$mapgain,$maptc2000,$maptc2016,$P2000,$TC2000,$FT2000,$FST2000,$P2016,$TC2016,$FT2016,$FST2016,$PMIN,$MIN,$PMAX,$MAX,$YLOSS,$LOSS,$TLOSS,$VHR1,$VHR2,$CERT,$NOTES)=split('\t',$sample);
#check consistency
if ($strata<1 or $strata>$nstrata) {die "$sample";}
#fill values
$matrix1[$strata][0]=$matrix1[$strata][0]+1; 
$matrix1[$strata][1]=$matrix1[$strata][1]+($proportion/100);
$matrix1[$strata][2]=$matrix1[$strata][2]+($proportion/100)*($proportion/100);
}

#print report
print OUT "Strata table\n";
print OUT "code\tarea\tcount\tN_samples\n";
for ($i=1; $i<=$nstrata; $i++){
print OUT "$matrix[$i][0]\t$matrix[$i][1]\t$matrix[$i][2]\t$matrix1[$i][0]\n";
}
}