*! version 1.0.0  22apr2008
program _haplo_em_phased, eclass
	version 10.0
	if replay() {
		if "`e(cmd)'" != "_haplo_em_phased" {
			error 301
		}
		Display `0'
		exit
	}
	cap noi Estimate `0'
	eret local cmdline `"_haplo_em_phased `0'"'
	eret scalar rc = _rc
	exit _rc
end

program Estimate, eclass
	version 10.0
	syntax [namelist(min=3 max=4)] [using/] ,   NSNP(integer)	   ///
						  [ EMITERate(int 500)	   ///
						    EMTOLerance(real 1e-6) ///
						    SAVing(string)	   ///
						    CUToff(real 1e-5)	   ///
						    EMLOG EMDOTs noHEADer  ///
						    noHAPTable BASE10	   ///
						    ONLYCUTOFF 		   ///
						    eminit(string)	   ///
						    SORT ]
	if `"`eminit'"' != "" {
		cap confirm matrix `eminit'
		local k1 = 2^`nsnp'-1
		if _rc {
			di as err "eminit(): 1x`=string(`k1',"%9.0g")'" ///
			   " matrix with numbers between 0 and 1 is required"
			exit 198
		}
		if colsof(`eminit')!=`k1' | rowsof(`eminit')!=1 {
			di as err "eminit(): 1x`=string(`k1',"%9.0g")'" ///
			   " matrix with numbers between 0 and 1 is required"
			exit 198
		}
	}
	if `"`using'"' != "" & `"`namelist'"' == "" {
		di as err "you must specify names for variables containing"
		di as err "subject identifiers, and haplotype pairs"
		exit 198
	}
	if `"`saving'"' != "" {
		_prefix_saving `saving'
		local replace   `"`s(replace)'"'
		if `"`replace'"' == "" {
			confirm new file `"`s(filename)'"'
		}
	}
	if `"`emdots'"' != "" {
		if `"`emlog'"' != "" {
			di as err "emlog and emdots are not allowed together"
			exit 198
		}
		local emdots .
	}
	local phasehaplo `base10'

	preserve
	qui use `namelist' using `"`using'"', clear
	local varlist `namelist'
	syntax [varlist] [using] [, *]
	marksample touse
	tokenize `varlist'
	if `"`4'"' != "" {
		qui replace `touse' = 0 if `4'==0
	}
	qui keep if `touse'
	if `c(N)' == 0 {
		error 2000
	}
	
	tempname freq
	mata: em_haplo( `nsnp', `"`freq'"', `emtolerance', `emiterate', ///
			"`emlog'", "`emdots'", "`base10'")
	ereturn matrix freq `freq'
	eret local cmd _haplo_em_phased
	eret local snpvars `2' `3'
	Display, `header' cutoff(`cutoff') `base10' ///
		  saving(`saving') `replace' `haptable' `sort'
	restore
end

program Display
	version 10.0
	syntax [, SAVing(string) REPLACE CUToff(real 1e-5)   ///
		  noHEADer noHAPTable BASE10 SORT ]

	di _n "{txt}Haplotype frequency EM estimation"
	local eqalign = 24
	if `"`header'"' == "" {
di _n "{txt}Number of obs{ralign `=`eqalign'-13': = }{res}" %10.0g e(N)
di   "{txt}Number of subjects{ralign `=`eqalign'-18': = }{res}" %10.0g e(N_s)
di   "{txt}Obs per subject:{ralign `=`eqalign'-16':min = }{res}" %10.0g e(N_min)
di   "{txt}{ralign `eqalign':max = }" as res %10.0g e(N_max)
	}
	di as txt _n "Number of iterations  = " as res %10.0g e(iter)
	di as txt    "Sample log-likelihood = " as res %10.0g e(ll)

	tempvar touse haplo hapfreq
	tempname hap freq
	qui gen byte `touse' = 1
	if `"`base10'"' == "" {
		qui gen str`=e(N_snp)' `haplo' = ""
		local ifmiss if `haplo' != ""
	}
	else {
		qui gen long `haplo' = .
		local ifmiss if `haplo' != .
	}
	mat `freq' = e(freq)'
	qui svmat double `freq', n(`hapfreq')
	rename `hapfreq'1 `hapfreq'
	mata: _di_haplotypes("`haplo'", "`base10'")
	if `"`haptable'"' == "" { 
		char `haplo'[varname] haplotype
		if `cutoff' == 0 {
			local star
			local starnote
		}
		else {
			local star *
			local starnote di as txt ///
			  "  * frequencies > " as res string(`cutoff',"%9.0g")
		}
		char `hapfreq'[varname] frequency`star'
		format `hapfreq' %8.0g
		if `"`sort'"' != "" {
			gsort -`hapfreq' `haplo'
		}
		`haptable' list `haplo' `hapfreq' ///
				`ifmiss' & `hapfreq'>=`cutoff', ///
			subvarname abbrev(10) noobs sep(0) div string(20)
		`starnote'
	}
	if (`e(converged)' == 0) {
		di as txt "Warning: " as err "convergence not achieved"
	}
	qui if `"`saving'"' != "" {
		_prefix_saving `saving'
		local saving    `"`s(filename)'"'
		local replace   `"`s(replace)'"'
		if `"`replace'"' == "" {
			confirm new file `"`s(filename)'"'
		}
		cap preserve
		keep `haplo' `hapfreq'
		keep `ifmiss' & `hapfreq' >= `cutoff'
		rename `haplo' haplotype
		rename `hapfreq' frequency
		save `"`saving'"', `replace'
		cap restore
	}
	qui keep if `touse'
end