clear
set more off
set mem 1g
# delimit ;

/*****************************************************************************************************************;
**This dofile calculates the 'simulated' grid measures of arsenic contamination and generates hh_simulatedmindistance_2004;
pre-generated datasets this dofile uses:
hhXwelldata_uncontaminated_2004
hhXwelldata_contaminated_2004

*****************************************************************************************************************/;

*****************************************************************************************************************;
*Paths;
*****************************************************************************************************************;

*access files in subdirectory generated by data creating dofiles;
local indatasetU "dtafiles/hhXwelldata_uncontaminated_2004";
local indatasetC "dtafiles/hhXwelldata_contaminated_2004";
local interdataset "dtafiles/hh_simulatedmindistance_2004_temp";
local outdataset "dtafiles/hh_simulatedmindistance_2004";

local path=subinstr(c(pwd),"Documents","Dropbox\BF\data replication",.);
cap cd "`path'";
cap cd "C:\Users\gshastry\Dropbox\BF\data replication";

*****************************************************************************************************************;
*Main
*****************************************************************************************************************;

*install necessary ado files;
ssc inst _gwtmean, replace;

cap log close;
log using SimulateExpectedMinDistance_2004.txt, text replace;

/*clear;
set obs 12760;
*set obs 20;

set seed 12345678;

gen x=uniform()*2-1;
gen y=uniform()*2-1;

gen distance=sqrt(x^2+y^2);
drop if distance>1;
drop distance;

count;
local numrandom=r(N);
gen point=_n;
sort point;
save random, replace;
*/;

use "`indatasetU'", clear;

*mathematically, only need to add wells within X+10.1 miles of the cluster center;
*where X is the distance from the cluster center to the closest uncontaminated well (call it well A);
*this way, point on other side of cluster is closer to the well A (distance X+5) than it is to;
*any potential wells just outside the radius (distance X+10.1-5);
egen mindistance=min(distance), by(dhsid);
drop if distance>max(15.1,mindistance+10.1);

for var lathhid-arsenic: rename X XU;

bysort dhsid: gen count=_n;
sum count;
local maxwellU=r(max);
drop distance mindistance;
reshape wide wellcountU latwellU longwellU arsenicU, i(dhsid) j(count);
order dhsid lathhid longhhid;

count;
local numU=r(N);

gen clusternum=_n;
sort clusternum;

tempfile U;
save `U', replace;

use "`indatasetC'", clear;

*mathematically, only need to add wells within X+10.1 miles of the cluster center;
*where X is the distance from the cluster center to the closest uncontaminated well (call it well A);
*this way, point on other side of cluster is closer to the well A (distance X+5) than it is to;
*any potential wells just outside the radius (distance X+10.1-5);
egen mindistance=min(distance), by(dhsid);
drop if distance>max(15.1,mindistance+10.1);

for var lathhid-arsenic: rename X XC;

bysort dhsid: gen count=_n;
sum count;
local maxwellC=r(max);
drop distance mindistance;
reshape wide wellcountC latwellC longwellC arsenicC, i(dhsid) j(count);
order dhsid lathhid longhhid;

count;
local numC=r(N);

gen clusternum=_n;
sort clusternum;

tempfile C;
save `C', replace;

** code is written assuming numC = numU - need to fix if not true;
assert `numC'==`numU';


*** generate a grid of points;

clear;
local oneside=114;

set obs `oneside';
gen x=-1+2*(_n-1)/(sqrt(_N*_N)-1);
expand `oneside';
sort x;
gen y=-1+2*(mod(_n,sqrt(_N))+sqrt(_N)*(mod(_n,sqrt(_N))==0)-1)/(sqrt(_N)-1);

gen distance=sqrt((x-x[_n-1])^2+(y-y[_n-1])^2) if mod(_n,sqrt(_N))!=1;
gen distance2=sqrt((x-x[_n-sqrt(_N)])^2+(y-y[_n-sqrt(_N)])^2);

sum dis*;
drop dis*;
gen distance=sqrt(x^2+y^2);
drop if distance>1;
drop distance;

count;
local numgrid=r(N);
gen point=_n;
sort point;

*save dtafiles/grid, replace;

expand `numU';

bysort point: gen clusternum=_n;

merge m:1 clusternum using `U';
assert _m==3; drop _m;

gen pointlatU=lathhidU+x*5/69.11;
gen pointlongU=longhhidU+y*5/(69.11*cos(lathhidU*_pi/180));

/*checks;
gen latdis=(lathhid-pointlat);
gen longdis=(longhhid-pointlong);

vincenty pointlat longhhid lathhid longhhid, h(distance_x); 
vincenty lathhid pointlong lathhid longhhid, h(distance_y); 

table dhsid, contents(mean lathhid min latdis max latdis max distance_x);
table dhsid, contents(mean lathhid mean longhhid min longdis max longdis max distance_y);
*/;

vincenty pointlat pointlong lathhid longhhid, h(distance_clusterU); 

forvalues i=1(1)`maxwellU' {;
	di "`i' of `maxwellU'";
	vincenty pointlat pointlong latwellU`i' longwellU`i', h(distance_wellU`i');
};

egen mindistanceU_grid=rmin(distance_wellU*);
gen invdistanceU=1/distance_cluster;

egen mean_mindistanceU=mean(mindistanceU_grid), by(dhsid);
egen wmean_mindistanceU=wtmean(mindistanceU_grid), weight(invdistance) by(dhsid);

gen mindistanceU_1=mindistanceU_grid<=1 if mindistanceU_grid!=.;
egen fraction_mindistanceU_1m=mean(mindistanceU_1), by(dhsid);
egen wfraction_mindistanceU_1m=wtmean(mindistanceU_1), weight(invdistance) by(dhsid);

/*gen mindistanceU_2=mindistanceU_grid<=2 if mindistanceU_grid!=.;
egen fraction_mindistanceU_2m=mean(mindistanceU_2), by(dhsid);
egen wfraction_mindistanceU_2m=wtmean(mindistanceU_2), weight(invdistance) by(dhsid);
*/;
/*temp stuff to debug: restore, not; sort dhsid; br dhsid mindistance_grid-fraction_mindistance_2;glu;*/;

** do the same with contaminated;
merge m:1 clusternum using `C';
assert _m==3; drop _m;

gen pointlatC=lathhidC+x*5/69.11;
gen pointlongC=longhhidC+y*5/(69.11*cos(lathhidC*_pi/180));

/*checks;
gen latdis=(lathhid-pointlat);
gen longdis=(longhhid-pointlong);

vincenty pointlat longhhid lathhid longhhid, h(distance_x); 
vincenty lathhid pointlong lathhid longhhid, h(distance_y); 

table dhsid, contents(mean lathhid min latdis max latdis max distance_x);
table dhsid, contents(mean lathhid mean longhhid min longdis max longdis max distance_y);
*/;

vincenty pointlatC pointlongC lathhidC longhhidC, h(distance_clusterC); 

forvalues i=1(1)`maxwellC' {;
	di "`i' of `maxwellC'";
	vincenty pointlatC pointlongC latwellC`i' longwellC`i', h(distance_wellC`i');
};

egen mindistanceC_grid=rmin(distance_wellC*);

gen mindistanceC_1=mindistanceC_grid<=1 if mindistanceC_grid!=.;
egen fraction_mindistanceC_1m=mean(mindistanceC_1), by(dhsid);
egen wfraction_mindistanceC_1m=wtmean(mindistanceC_1), weight(invdistance) by(dhsid);

/*gen mindistanceC_2=mindistanceC_grid<=2 if mindistanceC_grid!=.;
egen fraction_mindistanceC_2m=mean(mindistanceC_2), by(dhsid);
egen wfraction_mindistanceC_2m=wtmean(mindistanceC_2), weight(invdistance) by(dhsid);

gen mindistance_grid=min(mindistanceC_grid, mindistanceU_grid);
gen closest_contaminated=mindistance_grid==mindistanceC_grid;

gen closest_arsenic1=.;
gen closest_arsenic2=.;
gen closest_arsenic3=.;
gen closest_arsenic4=.;
gen closest_arsenic5=.;
forvalues i=1(1)`maxwellC' {;
	replace closest_arsenic5=arsenicC`i' if round(mindistanceC_grid,.0001)
		==round(distance_wellC`i',.0001) & closest_arsenic5==. & closest_arsenic1!=. & closest_contaminated!=0 & closest_arsenic2!=. & closest_arsenic3!=. & closest_arsenic4!=.;
	replace closest_arsenic4=arsenicC`i' if round(mindistanceC_grid,.0001)
		==round(distance_wellC`i',.0001) & closest_arsenic4==. & closest_arsenic1!=. & closest_contaminated!=0 & closest_arsenic2!=. & closest_arsenic3!=.;
	replace closest_arsenic3=arsenicC`i' if round(mindistanceC_grid,.0001)
		==round(distance_wellC`i',.0001) & closest_arsenic3==. & closest_arsenic1!=. & closest_contaminated!=0 & closest_arsenic2!=.;
	replace closest_arsenic2=arsenicC`i' if round(mindistanceC_grid,.0001)
		==round(distance_wellC`i',.0001) & closest_arsenic2==. & closest_arsenic1!=. & closest_contaminated!=0;
	replace closest_arsenic1=arsenicC`i' if round(mindistanceC_grid,.0001)
		==round(distance_wellC`i',.0001) & closest_arsenic1==. & closest_contaminated!=0;
};

forvalues i=1(1)`maxwellU' {;
	replace closest_arsenic5=arsenicU`i' if round(mindistanceU_grid,.0001)
		==round(distance_wellU`i',.0001) & closest_arsenic5==. & closest_arsenic1!=. & closest_contaminated==0 & closest_arsenic2!=. & closest_arsenic3!=. & closest_arsenic4!=.;
	replace closest_arsenic4=arsenicU`i' if round(mindistanceU_grid,.0001)
		==round(distance_wellU`i',.0001) & closest_arsenic4==. & closest_arsenic1!=. & closest_contaminated==0 & closest_arsenic2!=. & closest_arsenic3!=.;
	replace closest_arsenic3=arsenicU`i' if round(mindistanceU_grid,.0001)
		==round(distance_wellU`i',.0001) & closest_arsenic3==. & closest_arsenic1!=. & closest_contaminated==0 & closest_arsenic2!=.;
	replace closest_arsenic2=arsenicU`i' if round(mindistanceU_grid,.0001)
		==round(distance_wellU`i',.0001) & closest_arsenic2==. & closest_arsenic1!=. & closest_contaminated==0;
	replace closest_arsenic1=arsenicU`i' if round(mindistanceU_grid,.0001)
		==round(distance_wellU`i',.0001) & closest_arsenic1==. & closest_contaminated==0;
};

sum closest_arsenic*;
** if this assert fails, need to add a fourth option;
assert closest_arsenic5==.;
egen closest_arsenic=rmean(closest_arsenic1-closest_arsenic5);

egen fraction_closest_contam=mean(closest_contaminated), by(dhsid);
egen wfraction_closest_contam=wtmean(closest_contaminated), weight(invdistance) by(dhsid);
*/;
** distance between point and cluster shouldn't depend on whether well is C or U;
assert distance_clusterU==distance_clusterC;

/*egen mean_closest_arsenic=mean(closest_arsenic), by(dhsid);
egen wmean_closest_arsenic=wtmean(closest_arsenic), weight(invdistance) by(dhsid);

gen temp=mindistanceU_grid if closest_contaminated==1;
egen mean_mindistanceU_ofCpoints=mean(temp), by(dhsid);
egen wmean_mindistanceU_ofCpoints=wtmean(temp), weight(invdistance) by(dhsid);

gen temp3=mindistanceU_grid if mindistanceC_1==1;
egen mean_mindistanceU_inC1radius=mean(temp3), by(dhsid);
egen wmean_mindistanceU_inC1radius=wtmean(temp3), weight(invdistance) by(dhsid);

gen temp2=mindistanceU_grid if mindistanceC_2==1;
egen mean_mindistanceU_inC2radius=mean(temp2), by(dhsid);
egen wmean_mindistanceU_inC2radius=wtmean(temp2), weight(invdistance) by(dhsid);

*/;
compress;
*save "`interdataset'_grid", replace;

keep dhsid mean_* wmean_* fraction* wfraction*;
label variable mean_mindistanceU "Average Distance to Closest Uncontaminated Well";
label variable wmean_mindistanceU "Weighted Average Distance to Closest Uncontaminated Well";
label variable fraction_mindistanceU_1m "Fraction of Points within 1 m of Uncontaminated Well";
label variable wfraction_mindistanceU_1m "Weighted Fraction of Points within 1 m of Uncontaminated Well";
*label variable fraction_mindistanceU_2m "Fraction of Points within 2 m of Uncontaminated Well";
*label variable wfraction_mindistanceU_2m "Weighted Fraction of Points within 2 m of Uncontaminated Well";

label variable fraction_mindistanceC_1m "Fraction of Points within 1 m of Contaminated Well";
label variable wfraction_mindistanceC_1m "Weighted Fraction of Points within 1 m of Contaminated Well";
*label variable fraction_mindistanceC_2m "Fraction of Points within 2 m of Contaminated Well";
*label variable wfraction_mindistanceC_2m "Weighted Fraction of Points within 2 m of Contaminated Well";

*label variable fraction_closest_contam "Fraction of Points where closest well is Contaminated";
*label variable wfraction_closest_contam "Weighted Fraction of Points where closest well is Contaminated";

/*label variable mean_closest_arsenic "Average Contamination Level of Closest Well";
label variable wmean_closest_arsenic "Weighted Average Contamination Level of Closest Well";

label variable mean_mindistanceU_ofCpoints "Average Distance to Closest Uncontaminated Well for Points where Closest is Contaminated, Missing if No Contaminated Points";
label variable wmean_mindistanceU_ofCpoints "Weighted Average Distance to Closest Uncontaminated Well for Points where Closest is Contaminated, Missing if No Contaminated Points";

label variable mean_mindistanceU_inC1radius "Average Distance to Closest Uncontaminated Well for Points within 1 m of Contaminated Well, Missing if No Contaminated Points";
label variable wmean_mindistanceU_inC1radius "Weighted Average Distance to Closest Uncontaminated Well for Points within 1 m of Contaminated Well, Missing if No Contaminated Points";

label variable mean_mindistanceU_inC2radius "Average Distance to Closest Uncontaminated Well for Points within 2 m of Contaminated Well, Missing if No Contaminated Points";
label variable wmean_mindistanceU_inC2radius "Weighted Average Distance to Closest Uncontaminated Well for Points within 2 m of Contaminated Well, Missing if No Contaminated Points";
*/;
duplicates drop;

sort dhsid;
save "`outdataset'_grid", replace;



glu;

/*
(uniform()*2-1)*5/69.11;

1° of latitude = about 69.11 miles 

1° of longitude = about 69.11 miles along the equator. But all of the longitudes 
come together at the poles, so the farther you are from the equator, the fewer 
miles there are in one degree. 

Number of miles in 1° of longitude = (69.11) x (cosine of the latitude)

Read more: http://wiki.answers.com/Q/How_many_miles_are_in_a_degree_of_longitude_or_latitude#ixzz211i6u75C
*/;

