*! version 1.0.0  17apr2008
program _haplologit_d1_rarex
	version 10.0
	args todo b lnf g
	// dependent variable
	local y    $ML_y1
	// haplotype indices
	local hap1 $MY_hap1
	local hap2 $MY_hap2
	// subject ID
	local by $MY_panel
	tempvar xb L last
	// xb
	mleval `xb' = `b', eq(1)
	qui gen double `L' = .
	qui by `by': gen byte `last' = (_n==_N)
	// haplotype effects
	local neffect = colsof($MY_riskhapmat)
	local nriskhparms = 0
	forvalues i=1/`neffect' {
		tempvar riskhap`i'
		mleval `riskhap`i'' = `b', eq(`=`i'+1')
		local riskhapxb `riskhapxb' `riskhap`i''
		local eq ML_k`=`i'+1'
		local nriskhparms = `nriskhparms'+ $`eq'
		local nrhparms `nrhparms' $`eq'
		local eq ML_x`=`i'+1'
		local xrhvars `xrhvars' $`eq'
	}
	// haplotype frequencies
	tempname numat
	mat `numat' = 0
	local n = `neffect'+2
	forvalues i = $ML_n(-1)`n' {
		local j = `i'-2
		tempname nu`j' 
		mleval `nu`j'' = `b', eq(`i') scalar
		mat `numat' = `nu`j'', `numat'
	}
	tempname scores
	tempvar g1
	qui gen double `g1' = .
	mata: _haplologit_nrx(&$MY_inhermodel(), `"`scores'"')
	
	mlsum `lnf' = `L' if `last'
	if (`todo'==0 | `lnf' >= .) exit
	
	// g1 - gradient for xb
	tempname db
	mlvecsum `lnf' `db' = `g1' if `last'==1, eq(1)
	// matrix scores contains derivatives wrt other parameters
	// computed in _haplologit_nr()
	matrix `g' = (`db',`scores')
	
end