program define gllam_ll
version 6.0
args todo bo lnf alv probs res
*matrix list `bo'
if $HG_dots {
	noi disp in gr "." _c
}
/* ----------------------------------------------------------------------------- */
/* set up variables and macros needed */

local toplev = $HG_tplv
local topi = $HG_tpi
tempname pkpl b mzlc
local clus $HG_clus
sort `clus'
* reset the the clock and reset znow
local i = 1
matrix M_ip[1,1] = 1 /* in case topi=0 */
while (`i' <= `topi'){
	matrix M_ip[1,`i'] = 1
	local i = `i' + 1
}

/* -------------------------------------------------------------------------------- */
/* set up lprod1 ... lint1 */
local i=1
tempvar extra
quietly gen double `extra'=0
while `i' <= `toplev'{               
	tempvar  lint`i'
	gen double `lint`i'' = 0.0   /* used to integrate, therefore must be zero */
	if (`i'>1){
		tempvar  lprod`i'
		quietly gen double `lprod`i'' = .
		tempvar lfac`i'
		quietly gen double `lfac`i'' = 0  
	}
	local i = `i' + 1
}
/* set up names for HG_xb`i' */
local i = 1
while (`i' <= $HG_tpff){
	tempname junk
	global HG_xb`i' "`junk'"
	local i = `i' + 1
}
/* set up names for HG_s`i' */
local i = 1
while (`i'<=$HG_tprf){
	tempname junk
	global HG_s`i' "`junk'"
	local i = `i'+1
}

qui remcor6 "`bo'"

if $HG_error==1{
	disp in re "error in remcor6"
	scalar `lnf' = .
	exit
}

* disp "HG_xb1 after remcor: $HG_xb1[$which]:" $HG_xb1[$which]
*disp "HG_s2 after remcor: $HG_s2[$which]:" $HG_s2[$which]

local i = 2
while `i' <= `toplev'{
	quietly zip `i'
	local i = `i' + 1
}
*disp "STARTING LOOP"
/* --------------------------------------------------------------------------------- */
/* recursive loop through r.effs (levels)    */
/* topi nested loops: irf from 1 ro nip(rf) */
/* ip is "clock": (irf) stages of loops      */

*local lc = ln(10000) /* term to add to each lpyz */
*local lc = 0
local levno = `toplev'
local rf = `topi'
while (`rf' <= `topi') { /* for each r.eff */

/* ----------------------------------------------------------------------------------*/
/* reset ip to 1st point for all lower r.effs... */
/* update znow                                   */  
	*disp "reset ip up to random effect " `rf'
	while (`rf' > 1) {
		local rf = `rf' - 1
		*disp `rf'
		matrix M_ip[1,`rf'] = 1
	}
	while (`levno' > 1){ 
	/* update znow for all new ips    */
		quietly zip `levno'
		local levno = `levno' - 1
	}
/* --------------------------------------------------------------------------------- */
/* cycling through lev1: set lint1 to lpyz for new znow  */
		
	local rf = 1
	local levno = 1
	local sortlst `clus' /* cluster variables to aggregate over */
	*disp "rf = " `rf'
	while (M_ip[1,`rf'] <= M_nip[1,`rf']) { /* for each int. point for curr. r.eff */
		if (M_nip[1,`rf'] > 1) {quietly zip `levno'} /* no need to call zip for level 1 */
		/* replace lint1 by current cond. log likelihood for lev1 unit*/
		*matrix list M_ip

		qui lpyz `lint1'

		* noi disp " after lpyz lint1 = " `lint1'[$which]
		quietly count if `lint1'==.& $HG_ind==1
		if r(N) > 0{
			* overflow problem
			noi disp "overflow at level 1 ( " r(N) " missing values)"
			*list $HG_clus if `lint1'==.
			*matrix list `bo'
			*lpyz `lint1'
			if "`res'" == ""{
				scalar `lnf' = .
				exit
			}	
		}
		quietly replace `lint1' = `lint1' * $HG_wt1
		matrix M_ip[1,`rf'] = M_ip[1,`rf'] + 1
		*next ip
	}
/* --------------------------------------------------------------------------------- */
/* update lint for all completed levels up to */
/* highest completed level, reset lower lints */
/* to zero (for models including a random effect) */

	* now ip too high for current `rf'
	matrix M_ip[1, `rf'] = M_nip[1, `rf']    /* reset ip at rf to highest value */
	while (M_ip[1,`rf'] == M_nip[1,`rf'] & `rf' <= `topi'){ 
	* digit equals its max => increment next digit	
		if (`rf' == M_nrfc[1,`levno'] & `levno' < `toplev'){  
		* done last r.eff of current level
			* disp "********** level " `levno' " complete ************"
			* next level
			local levno = `levno'+1 
 
			/* change sortlst */         
			local l = `toplev' - `levno' + 2
			tokenize "`sortlst'"
			* take away var. of level to sum over          	
			local `l' " "                 
			local sortlst "`*'"

/*------------------------------------------------------------------------------------ */
/* update lint and lprod */
		quietly{
			*!! disp "!!! update level " `levno'

			/* get pkpl: product of r.effs at level */            
			* update znow for levno
			zip `levno'               
			zprob `levno'
			scalar `pkpl' = $S_1

			/* set previous lint to ln(lint) */
			local lprev = `levno' - 1
			if(`levno' > 2){
				*!! disp " replace lint" `lprev' " by ln(lint" `lprev' ")"
				quietly count if `lint`lprev'' < 1e-308
				if r(N) > 0{
					/* overflow problem */
					noi disp "overflow at level " `lprev'
					if "`res'" == ""{
						scalar `lnf' = .
						exit
					}	
				}
				quietly replace `lint`lprev'' = ln(`lint`lprev'')
				quietly replace `lint`lprev'' = (`lint`lprev''-`lfac`lprev'')*${HG_wt`lprev'}
			}

			/* sum previous lprod within cluster at current level */
			*!! disp " "
			*!! disp "by `sortlst': replace lprod" `levno' "=cond(_n==N, sum(lint" `lprev' "))"
			quietly by `sortlst': replace `lprod`levno'' = /*
				*/ cond(_n==_N,sum(`lint`lprev''),.)
			*!! noi disp "lprod" `levno' " = " `lprod`levno''[$which]
			/* avoid underflow: extra added to lfac, lint multiplied by exp(extra) */
			
			*local nxtlev = `levno' + 1
			* noi matrix list M_ip
			local lstrf = M_nrfc[1,`levno']
			local frstrf = M_nrfc[1,`lprev'] + 1
			* noi disp "lstrf = `lstrf' and firstrf = `frstrf'"
			*!! disp "checking that M_ip[1,`frstrf' to `lstrf'] = 1"
			local frst = 1
			while(`frstrf' <= `lstrf'){
				if M_ip[1,`frstrf'] ~= 1{
					local frst = 0
				}
				local frstrf = `frstrf' + 1
			}	
			if `frst' == 1{ /* first term in integral */
				*!! noi disp "first term for level `levno', random effect `rf'"
				quietly replace `extra' = 0
				quietly replace `lfac`levno'' = -`lprod`levno''
				* noi disp "lfac`levno' = " `lfac`levno''[$which]
				if "`res'" ~= ""{
					local lst
					local sofar = 0
					local comp = M_nrfc[1,2]
					local found = 0
					tempvar rest resl
					tempname mrot
					gen double `rest' = `lint`levno''
					gen double `resl' = 0
					local ll = M_nrfc[2,1]
					while `ll'<M_nrfc[2,2]{
						gen double `res'm`ll' = 0
						* disp "          `res'm`ll'"
						local ll =  `ll' + 1
					}
				} 
			}
			else{   
				quietly replace `extra' = /*
				*/cond((`lprod`levno''+`lfac`levno''-500)>0, /*
				*/ -(`lprod`levno''+`lfac`levno''-500),0)
				* noi disp "extra = " `extra'[$which]
				quietly replace `lfac`levno''=`lfac`levno''+`extra'
				*!! noi disp "lfac`levno' = " `lfac`levno''[$which]
			}
			/* increment lint at current level using lprod at previous level */
			*!! disp " "
			*!! disp "increase lint" `levno' " by exp(lprod" `levno' ")*pkpl"
			quietly replace `lint`levno''=/*
				*/ exp(`extra')*`lint`levno''+exp(`lprod`levno''+`lfac`levno'')*`pkpl'
			*!! noi disp "increase by " exp(`lprod`levno''[$which]+`lfac`levno''[$which])*`pkpl' " to "  `lint`levno''[$which]
			if "`res'" ~= ""{
			   noi{ /* for current ip value: */
				if `levno' == 2{
					*!! disp "sortlst = `sortlst'"

					local sofar = `sofar' + 1
					* disp "updating `rest' the `sofar'th time"
					qui replace `rest' = `lint`levno''
					qui replace `resl' = exp(`lprod`levno''+`lfac`levno'')*`pkpl'

					* disp "resl = " `resl'[$which]
					local first = M_nrfc[1,1] + 1 /* loop */
					local ll = `first'
					
					if $HG_free==1{
						matrix `mrot' = M_znow
					}
					else{
						*matrix `mrot' = M_chol*M_znow'
						matrix `mrot' = CHmat*M_znow'
						matrix `mrot' = `mrot''
						* matrix list `mrot'
					}
					local kk = M_nrfc[2,1] /* r. eff */
					while `ll' <= M_nrfc[1,`levno']{
					while `kk' < M_nrfc[2,`levno']{
						local kkp = `kk'+1
						* disp "loop " `ll'
						* disp "random effect " `kk'
						local junk = M_ip[1,`ll']
						local w = M_nip[2,`kkp']
						*scalar `mzlc' = M_zlc`w'[1,`junk']
						scalar `mzlc' = `mrot'[1,`kk']
						* disp "mzlc = " `mzlc'
						qui replace `res'm`kk' = `res'm`kk'*exp(`extra')+/*
							*/ `mzlc'*`resl' /* *${HG_s`kkp'} */
						* disp "after adding M_znow[1," `kk' "] = " M_znow[1,`kk'] ", `res'm`kk' = " `res'm`kk'[$which]
						if M_ip[1,`comp'] == M_nip[1,`comp']{
							* disp "loop " `comp' " is complete"
							if `comp' == `first'{
							  *matrix list M_ip
							  * disp "last update: divide `res'm`kk' by `rest'"
							  qui replace `res'm`kk' = `res'm`kk'/`rest'
							}
							else if `found' == 0{
								* disp "reduce comp by 1"
								local comp = `comp' -1
								local found = 1
							}
						}
						local kk = `kk' + 1
						}
						local ll = `ll' + 1
					}
					if `probs'>0 {
					/* create and update posterior probabilities */
					local lab
					local ll = M_nrfc[1,1]+1
					while `ll' <= M_nrfc[1,2]{ /* each loop */
						local junk = M_ip[1,`ll'] 
						*!! disp "ip = " `junk'
						local lab "`lab'`junk'"
						local ll = `ll'+1
					}
					*!! disp "lab is `lab'"
					qui gen double `res'`lab' = `resl'

					if M_ip[1,`comp'] == M_nip[1,`comp']&`comp'==`first'{
						*!! disp "last update: divide `res'`lab' by `rest'"
						qui replace `res'`lab'=`res'`lab'/`rest'
					}

					*!! disp "`res'`lab' = " `res'`lab'[$which]
					tokenize "`lst'"
					while "`1'" ~= ""{ /* for all posterior probabilities */
						* disp "multiplying `1' by "  exp(`extra'[$which])
						qui replace `1' = `1'*exp(`extra')
						if M_ip[1,`comp'] == M_nip[1,`comp']&`comp'==`first'{
							*!! disp "dividing `1' by pt" 
							qui replace `1' = `1'/`rest'
						}
						mac shift
					}
					local lst "`lst' `res'`lab'"
					} /* endif */
				}
			   } /* noi */
			}

			/* reset previous lint to zero */
			if `levno'>2{
				*!! disp "setting lint" `lprev' " to zero"
				quietly replace `lint`lprev'' = 0
				quietly replace `lfac`lprev'' = 0
			}
		 } /* qui */
/* ------------------------------------------------------------------------------------------- */	
		}
		* next digit
		local rf = `rf' + 1
	}
	* rf is first r.eff that is not complete
	* increase clock in lowest incomplete digit
	* disp "update rf = " `rf'
	matrix M_ip[1,`rf'] = M_ip[1,`rf'] + 1
	*matrix list M_ip
}/*}*/
quietly{ 
	*now rf too high
	*!! disp "********** level " `toplev' " complete ************"
	*!! noi disp "lint" `toplev' " = " `lint`toplev''[$which]
	if(`toplev'>1){

		disp "taking log of lint" `toplev' " = " `lint`toplev''[$which]
		disp "subtracting " `lfac`toplev''[$which]
		quietly replace `lint`toplev'' = (ln(`lint`toplev'')-`lfac`toplev'')* ${HG_wt`toplev'}
	}
	* noi display "lnf[" $which "] = " `lint`toplev''[$which]
	qui by `sortlst': replace  `extra' = cond(_n==_N,1,0)
	*mlsum `lnf' = `lint`toplev'' if `extra' == 1 /* can only use this when program called by ML */
	count if `extra' == 1
	local n = r(N)
	summarize `lint`toplev'', meanonly
	if `n' > r(N) {
		noi disp "there are " r(N) " values of likelihood, should be " `n'
		noi list `sortlst' if `extra' == 1& `lint`toplev''==.
		noi disp "lnf equal to missing in last step"
		scalar `lnf' = .
		exit
	}
	scalar `lnf' = r(sum)
	* noi display "total lnf = " `lnf'
	* capture drop lint`toplev'
	* gen double lint`toplev' = `lint`toplev''
} /* qui */
end

program define zip
	version 6.0
	* updates znow 
	* matrix list M_ip
	args levno
/* -----------------------------------------------------------------------------*/
/* do we need to update all r.effs at current level?   */

	disp "in zip, levno is " `levno'
	local i = M_nrfc[2,`levno'-1] + 1 

	*!! disp "update" 
	local k = M_nrfc[1,`levno']
	local k = M_ip[1,`k']
	local last = M_nrfc[2,`levno']
	while `i' <= `last'{
		if $HG_free{
			* same class for all random effects
			local which = `k'
		}
		else{
			local which = M_ip[1,`i']
		}
		local npt = M_nip[2,`i']
		local im = `i' - 1
	 	disp "     "`im' "th z to " `which' "th location"
		disp " using M_zlc`npt' "
		matrix M_znow[1,`im'] = M_zlc`npt'[1,`which']
		*!! disp M_znow[1,`im'] 
		local i = `i' + 1
	}
end

program define zprob
	version 6.0
	* returns product of pk needed for integration at level lev for current ip
	args levno
	tempname pkpl mzps
	disp "in zprob, levno is " `levno'
	scalar `pkpl' = 1.0

	local i=M_nrfc[1,`levno'-1] + 1 

	*!! disp "-----------pkpl: product of" 
	local last = M_nrfc[1,`levno'] 
	local k = M_ip[1,`last']
	while `i' <= `last'{
		if $HG_free{
			local which = `k'
		}
		else{
			local which = M_ip[1,`i']
		}
		local npt = M_nip[2,`i']
		disp "     prob for " `i' "th r.eff: " `which' "th weight"
		disp " using M_zps`npt' "
		* product of probabilities
		scalar `mzps' = M_zps`npt'[1,`which']
		scalar `pkpl' = `pkpl'*`mzps'
		local i=`i'+1
	}
	global S_1 = `pkpl'
	*!! disp "pkpl = " `pkpl'
end


program define lpyz
	version 6.0
* returns log of prob of obs. given znow
	args lpyz

	disp "-----------------called lpyz"

	tempname mznow
	tempvar zu xb mu /* linear predictor and zu: r.eff*design matrix for r.eff */

/* ----------------------------------------------------------------------------- */
*quietly{
/* level 1 was done in remcor*/

/* levels 2 to toplev*/
	

	if $HG_tprf>1{
		matrix score double `zu' = M_znow
	}
	else{
		qui gen double `zu' = 0
	}

	matrix list M_znow
	* disp "ML_y1: $ML_y1 " $ML_y1[$which]

	*matrix list M_ip

	disp " xb1 = " $HG_xb1[$which]
	disp " zu = " `zu'[$which]

	*replace `zu' = 0

	if "$HG_link"=="mlogit"|"$HG_link"=="smlogit"{
		tempvar s
		if "$HG_link"=="smlogit"{
			gen double `s' = $HG_s1
		}
		else{
			gen double `s' = 1
		}
		if $HG_exp==1&$HG_expf==0{
			tempvar mu
			*qui gen double `mu' = 1 if $ML_y1==M_resp[1,1]
			qui gen double `mu' = exp(`zu'/`s') if $ML_y1==M_resp[1,1]
			local n=rowsof(M_resp)
			local i=2
			while `i'<=`n'{
				local prev = `i' - 1 
				* disp "xb`prev':" ${HG_xb`prev'}[$which]
				qui replace `mu' = exp((${HG_xb`prev'} + `zu')/`s') if $ML_y1==M_resp[`i',1]
				local i = `i' + 1
			}

			sort $HG_clus $HG_ind
			qui by $HG_clus: replace `lpyz'=cond(_n==_N,sum(`mu'),.)
			qui replace `lpyz' = ln(`mu'/`lpyz') 
		}
		else if $HG_exp==1&$HG_expf==1{
			tempvar mu
			qui gen double `mu' = exp(($HG_xb1 + `zu')/`s')
			list $HG_xb1 `mu' a ind in 1/4
			sort $HG_clus $HG_ind
			qui by $HG_clus: replace `lpyz'=cond(_n==_N,sum(`mu'),.)
			disp "denom = " `lpyz'[$which]
			qui replace `lpyz' = ln(`mu'/`lpyz') 
		}
		else{
			tempvar den mu tmp
			local n=rowsof(M_resp)
			local i = 2
			qui gen double `mu' = 1 if $ML_y1==M_resp[1,1]
			qui gen double `den' = 1
			qui gen double `tmp' = 0
			while `i'<= `n'{
				local prev = `i' - 1 
				qui replace `tmp' = exp((${HG_xb`prev'} + `zu')/`s')
				qui replace `mu' =  `tmp' if $ML_y1==M_resp[`i',1]
				replace `den' = `den' + `tmp'
				local i = `i' + 1
			}
			replace `lpyz' = ln(`mu'/`den')
		}
	}

	else if $HG_olog==0|"HG_lv"~=""{
		local and 
		if "$HG_lv"~=""&$HG_olog>0{
			local and & $HG_fv ~= $HG_olog
		}
		* disp in re "and=`and'"
		quietly gen double `mu' = 0
		link "$HG_link" `mu' $HG_xb1 `zu' $HG_s1
		disp " mu = " `mu'[$which]
		family "$HG_famil" `lpyz' `mu' $HG_s1 "`and'"
	}

	if $HG_olog~=0{
		local lnk: word $HG_olog of $HG_link
		if "`lnk'"=="ologit"{
			local func logitl
		}
		else if "`lnk'"=="oprobit"{
			local func probitl
		}
		else if "`lnk'"=="ocll"{
			local func cll
		}
		tempvar mu p1 p2
		local and 
		if "$HG_lv"~=""&$HG_olog>0{
			local and & $HG_fv == $HG_olog
		}
		local n=rowsof(M_resp)
		qui gen double `p1' = 0
		qui gen double `mu' = $HG_xb1+`zu'-$HG_xb2
		`func' `mu' `p1'
		qui replace `lpyz' = ln(1-`p1') /*
			*/ if $ML_y1==M_resp[1,1] `and'
		qui gen double `p2' = `p1'
		local i = 2
		while `i' < `n'{
			local nxt = `i' + 1
			qui replace `mu' = $HG_xb1+`zu'-${HG_xb`nxt'}
			`func' `mu' `p2'
			qui replace `lpyz' = ln(`p1' -`p2') /*
				*/ if $ML_y1==M_resp[`i',1] `and'
			qui replace `p1' = `p2'
			local i = `i' + 1
		}
		qui replace `lpyz' = ln(`p2') /*
			*/ if $ML_y1==M_resp[`n',1] `and'
		qui replace `lpyz' = -100 if `lpyz'==.
	}
	noi disp "lnz = " `lpyz'[$which]
*} /* qui */
end

program define logitl
	version 6.0
	args mu p
	qui replace `p' = 1/(1+exp(-`mu'))
end

program define cll
	version 6.0
	args mu p
	qui replace `p' = 1-exp(-exp(`mu'))
end

program define probitl
	version 6.0
	args mu p
	qui replace `p' = normprob(`mu')
end


program define link
	version 6.0
* returns mu for requested link
	args which mu xb zu s1
	*disp " in link, which is `which' "

	tokenize "`which'"
	local i=1
	while "`1'"~=""{
		if "$HG_lv" == ""{
			local ifs 
		}
		else{
			local ifs if $HG_lv==`i'
		}
		
		if ("`1'" == "logit"){
		noi disp "quietly replace `mu' = 1/(1+exp(-`xb'-`zu')) `ifs'"
			quietly replace `mu' = 1/(1+exp(-`xb'-`zu')) `ifs'
		}
		else if ("`1'" == "probit"){
			*disp "doing probit "
			quietly replace `mu' = normprob((`xb'+`zu')) `ifs'
		}
		else if ("`1'" == "sprobit"){
			quietly replace `mu' = normprob((`xb'+`zu')/`s1') `ifs'
		}
		else if ("`1'" == "log"){
			*disp "doing log "
			quietly replace `mu' = exp(`xb'+`zu') `ifs'
		}
		else if ("`1'" == "recip"){
			*disp "doing recip "
			quietly replace `mu' = 1/(`xb'+`zu') `ifs'
		}
		else if ("`1'" == "cll"){
			*disp "doing cll "
			quietly replace `mu' = 1 - exp(-exp(`xb'+`zu')) `ifs'
		}
		else{
		noi disp "quietly replace `mu' = `xb'+`zu' `ifs'"
			quietly replace `mu' = `xb'+`zu' `ifs'
		}
		local i = `i' + 1
		mac shift
	}

end

program define family
	version 6.0
	args which lpyz mu s1 and

	tokenize "`which'"
	local i=1
	local und `and'
	while "`1'"~=""{
		if "$HG_fv" == ""{
			local ifs 
			local und
		}
		else{
			local ifs if $HG_fv == `i'
		}
		if ("`1'" == "binom"){
			famb `lpyz' `mu' "`ifs'" "`und'"
		}
		else if ("`1'" == "poiss"){
			famp `lpyz' `mu' "`ifs'" "`und'"
		}
		else if ("`1'" == "gauss") {
			famg `lpyz' `mu' `s1' "`ifs'" "`und'"  /* get log of conditional prob. */
		}
		else if ("`1'" == "gamma"){
			famga `lpyz' `mu' `s1' "`ifs'" "`und'"
		}
		else{
			disp in re "unknown family in gllam_ll"
			exit 198
		}
		local i = `i' + 1
		mac shift
	}
end
	
program define famg
	version 6.0
* returns log of normal density conditional on r.effs
	args lpyz mu s1 if and
	disp "running famg `if' `and'"
	disp "s1 = " `s1'[$which] ", mu = " `mu'[$which] " and Y = " $ML_y1[$which]
      	quietly replace `lpyz' = /*
		*/ -(ln(2*_pi*`s1'^2) + (($ML_y1-`mu')/`s1')^2)/2 `if' `and'
end

program define famb
	version 6.0
* returns log of binomial density conditional on r.effs
* $HG_denom is denominator
	args lpyz mu if and
	disp "running famb `if' `and'"
	disp "mu = " `mu'[$which] " and Y = " $ML_y1[$which]
	quietly replace `lpyz' = /*
		*/ $ML_y1*ln(`mu')+($HG_denom-$ML_y1)*cond(ln(1-`mu')~=.,ln(1-`mu'),-100) `if' `and'
	disp "done famb"
end

program define famp
	version 6.0
* returns log of poisson density conditional on r.effs
	args lpyz mu if and
	*!! disp "running famp `if'"
	*disp in re "if and: `if' `and'"
	quietly replace `lpyz' = /*
		*/ $ML_y1*(ln(`mu'))-`mu'-lngamma($ML_y1+1) `if' `and'
	* disp "done famp"
end

program define famga
	version 6.0
* returns log of gamma density conditional on r.effs
	args lpyz mu s1 if and
	*!! disp "running famg `if'"
	*!! disp "mu = " `mu'[$which]
	*!! disp "s1 = " `s1'[$which]
	qui replace `mu' = 0.0001 if `mu' <= 0
	tempvar nu
	qui gen double `nu' = `s1'^(-2)
      	quietly replace `lpyz' = /*
		*/ `nu'*(ln(`nu')-ln(`mu')) - lngamma(`nu')/*
		*/ + (`nu'-1)*ln($ML_y1) - `nu'*$ML_y1/`mu' `if' `and'
end