*! version 1.1.0 , feb 98, Guy van Melle  STB-45 gr30
*!
*! 3-D summary of data, interface to 3D-ados <hidlin> or <altitude>
*!
*! syntax: makfun x y 
*!		[, x(# or mat) y(# or mat) z(exp) GRound(#) SMooth ]
*!
*! data preserved but 3 mats created: hlXcut, hlYcut, hlZfun
*! 	for use by function <hlmakfun> in hidlin or altitude call
*!
*
*--------------
prog def makfun
*--------------
	version 5.0
	loc varlist "req exi min(2) max(2)"
	loc options "x(string) y(string) z(string) GRound(real 0) SMooth"

	parse "`*'"
	parse "`varlist'", parse(" ")
	loc xv `1'
	loc yv `2'
	cap memory
	loc matsz= cond(_rc,40,_result(12))-1

*--- check x & y info --- 

	if"`x'"=="" { loc x . }
	verif x `x' hlXcut `matsz'
	loc  mx  $S_1

	if"`y'"=="" { loc y . }
	verif y `y' hlYcut `matsz'
	loc  my  $S_1

*--- check zexp ---

	loc zx `z'

	if ("`z'"=="") | ("`z'"=="freq") { loc zx count `yv' }
	loc wz: word count `zx'
	if `wz'!=2 {
		di in red "z-expression should be <statistic varname>"
		exit
	}
	loc z1: word 1 of `zx'
	loc z2: word 2 of `zx'
	confirm var `z2'
	loc j 1
	while `xv'[`j']==. & `j'<_N { loc j=`j'+1 }
	if `xv'[`j']==. { 
		di in blue "no observations"
		exit 2000
	}
	cap table `xv' in `j', c(`zx')
	if _rc {
		di in red "illegal z-expression z(`z')"
		exit
	}

*--- make hlXcut & hlYcut mat

	preserve
	tempvar xp, xt, yp, yt

	if `mx'==0 {
		qui g `xp'=.
		loc r 0		
		while  `r'< rowsof(hlXcut) {
			loc r=`r'+1
			qui replace `xp'=hlXcut[`r',1] in `r'
		}
	}
	else{
		pctile `xp'=`xv', nq(`mx') 

		loc mx1=`mx'-1
		mkmat  `xp' in 1/`mx1', mat(hlXcut)
	}
	xtile `xt'=`xv', cut(`xp')

	if `my'==0 {
		qui g `yp'=.
		loc r 0		
		while  `r'< rowsof(hlYcut) {
			loc r=`r'+1
			qui replace `yp'=hlYcut[`r',1] in `r'
		}
	}
	else{
		pctile `yp'=`yv', nq(`my') 
		loc my1=`my'-1
	    mkmat  `yp' in 1/`my1', mat(hlYcut)
	}
	xtile `yt'=`yv', cut(`yp')

	collapse (`z1') `z2', by(`xt' `yt') 

	su `z2'

*--- make hlZfun mat
 
	loc rx = rowsof(hlXcut)
	loc rx1=`rx'+1
	loc ry = rowsof(hlYcut)
	loc ry1=`ry'+1

	mat hlZfun=J(`rx1',`ry1',`ground')	
	loc L  0
	while `L'<_N {
		loc L=`L'+1
		loc ix=`xt'[`L']
		loc iy=`yt'[`L']
		mat hlZfun[`ix',`iy']=`z2'[`L']
	}

*--- smooth Z if requested

	if "`smooth'"=="" { exit}
	mat hlZfu1= hlZfun
	mat hlZfun= J(`rx1',`ry1',`ground')
	loc i 0
	while `i'<`rx1' {
		loc i = `i'+1
		loc j 0
		while `j'<`ry1' {
			loc j = `j'+1
			loc s 0
			loc w 0
			loc a=`i'-2
			while `a'<`i'+1 {
			    loc a=`a'+1
			    if `a'>0 & `a'<=`rx1' {
				loc b=`j'-2
				while `b'<`j'+1 {
			 	    loc b=`b'+1
				    if `b'>0 & `b'<=`ry1' {
					if hlZfu1[`a',`b']!=`ground' {
					    loc k= 1 + (`a'==`i' & `b'==`j')
					    loc s=`s' + `k'* hlZfu1[`a',`b']
					    loc w=`w' + `k'
					}
				  }
				}
			  }
			}
			if `w'>0 { mat hlZfun[`i',`j']=`s'/`w' }
		}
	}
end	
*
prog def verif
*
	loc nm `1'
	loc u  `2'
	loc M  `3'
	loc sz `4'
	
	loc m -1
	if "`u'"=="." {	loc m= min(max(int((_N^.5)/3),2),`sz')}
	else {
		cap confirm number `u'
		if _rc==0 {	loc m= max(int(abs(`u')),2)}
		else {
			cap mat li `u'
			if _rc==0 { 
				if rowsof(`u') > 1 & colsof(`u') > 1 {
				    di in red "matrix `u' must be 1xk or kx1"
				    error 198
				}
				mat `M'=`u'
				if rowsof(`M')==1 { mat `M' = `M'' }
				loc m= 0 
			}
		}
	}
	if `m'==-1 {
		di " `nm'(`u') ! " _cont
		di in red "'`u'' not a number, not a matrix... " _cont
		error 198
	}
	glo S_1 `m'
end
*