program define survsim_user
	version 14.2

	syntax newvarname(min=2 max=2), 								///
										MAXTime(string)				///
																	///
									[								///
										LTruncated(string)			///
																	///
										LOGHazard(string)			///
										Hazard(string)				///
										LOGCHazard(string)			///
										CHazard(string)				///
										NODES(int 30)				///
																	///
										MIXTURE						///
										Distribution(string)		///
										Lambdas(numlist)			///
										Gammas(numlist)				///
										PMIX(real 0.5)				///
																	///
										COVariates(string)			///
										TDE(string)					///
										TDEFUNCtion(string)			///
																	///
										MARGinal					///
									]								//

		local stime : word 1 of `varlist'
		local died  : word 2 of `varlist'

		cap which lmoremata.mlib
		if _rc {
			display in yellow "You need to install the moremata package. This can be installed using,"
			display in yellow ". {stata ssc install moremata}"
			exit 198
		}

		//====================================================================================================================//
		//mixture - build hazard/chazard and continue as user

			if "`mixture'"!="" {

				local Nls  : word count `lambdas'
				local Ngs  : word count `gammas'

				if `Nls'!=2 | `Ngs'!=2  {
					di as error "Number of lambdas() and gammas() must be 2 under a mixture model"
					exit 198
				}

				local l = length("`distribution'")
				if `l'==0 {
					di as error "distribution() required"
					exit 198
				}
				if 		substr("exponential",1,max(1,`l'))=="`distribution'" 	local dist "exp"
				else if substr("gompertz",1,max(3,`l'))=="`distribution'" 		local dist "gompertz"
				else if substr("weibull",1,max(1,`l'))=="`distribution'" 		local dist "weibull"
				else {
					di as error "Unknown distribution"
					exit 198
				}

				local l1 : word 1 of `lambdas'
				local l2 : word 2 of `lambdas'
				if "`dist'"=="weibull" | "`dist'"=="gompertz" {
					local g1 : word 1 of `gammas'
					local g2 : word 2 of `gammas'
				}

				if "`tde'"=="" {

					if "`dist'"=="exp" {
						local chazard "-log(`pmix':*exp(-`l1':*{t}) :+ (1:-`pmix'):*exp(-`l2':*{t}))"
					}
					else if "`dist'"=="weibull" {
						local chazard "-log(`pmix':*exp(-`l1':*{t}:^(`g1')) :+ (1:-`pmix'):*exp(-`l2':*{t}:^(`g2')))"
					}
					else {
						local chazard "-log(`pmix':*exp((`l1':/`g1'):*(1:-exp(`g1':*{t}))) :+  (1:-`pmix'):*exp((`l2':/`g2'):*(1:-exp(`g2':*{t}))))"
					}

				}
				else {

					if 		"`dist'"=="exp" {
						local base_surv "(`pmix':*exp(-`l1':*{t}) :+ (1:-`pmix'):*exp(-`l2':*{t}))"
						local numer 	"(`l1':*`pmix':*exp(-`l1':*{t}) :+ `l2':*(1:-`pmix'):*exp(-`l2':*{t}))"
						local hazard 	"`numer' :/ `base_surv'"
					}
					else if "`dist'"=="weibull" {
						local base_surv "(`pmix':*exp(-`l1':*{t}:^(`g1')) :+ (1:-`pmix'):*exp(-`l2':*{t}:^(`g2')))"
						local numer 	"(`l1':*`g1':*`pmix':*{t}:^(`g1':-1):*exp(-`l1':*{t}:^(`g1')) :+ `l2':*`g2':*(1:-`pmix'):*{t}:^(`g2':-1):*exp(-`l2':*{t}:^(`g1')))"
						local hazard 	"`numer' :/ `base_surv'"
					}
					else {
						local base_surv "(`pmix':*exp((`l1':/`g1'):*(1:-exp(`g1':*{t}))) :+  (1:-`pmix'):*exp((`l2':/`g2'):*(1:-exp(`g2':*{t}))))"
						local numer 	"`pmix':*exp((`l1':/`g1'):*(1:-exp(`g1':*{t}))) :* (-`l1':*exp(`g1':*{t})) + (1:-`pmix'):*exp((`l2':/`g2'):*(1:-exp(`g2':*{t}))) :* (-`l2':*exp(`g2':*{t}))"
						local hazard 	"(`numer') :/ `base_surv'"
					}

				}

				//now continues as user function
			}


		//====================================================================================================================//
		//baseline covariates

			tempvar expxb
			if "`covariates'"!="" {

				tokenize `covariates'
				local ncovlist : word count `covariates'
				local ncovvars = `ncovlist'/2
				cap confirm integer number `ncovvars'
				if _rc>0 {
					di as error "Variable/number missing in covariates"
					exit 198
				}
				local ind = 1
				local error = 0
				forvalues i=1/`ncovvars' {
					cap confirm numeric var ``ind'', exact
					if _rc {
						local errortxt "invalid covariates(... ``ind'' ``=`ind'+1'' ...)"
						local error = 1
					}
					cap confirm num ``=`ind'+1''
					if _rc {
						local errortxt "invalid covariates(... ``ind'' ``=`ind'+1'' ...)"
						local error = 1
					}
					tempvar vareffect`i'
					qui gen double `vareffect`i'' = ``ind''*``=`ind'+1''
					local ind = `ind' + 2
				}
				if `error' {
					di as error "`errortxt'"
					exit 198
				}
				local cov_linpred "`vareffect1'"
				if `ncovvars'>1 {
					forvalues k=2/`ncovvars' {
						local cov_linpred "`cov_linpred' + `vareffect`k''"
					}
				}

				qui gen double `expxb' = exp(`cov_linpred')
			}
			else qui gen byte `expxb' = 1

		//====================================================================================================================//
		// Time-dependent effects

			tempvar tdexb
			if "`tde'"!="" {
				tokenize `tde'
				local ntde : word count `tde'
				local ntdevars = `ntde'/2
				cap confirm integer number `ntdevars'
				if _rc>0 {
					di as error "Variable/number missing in tde"
					exit 198
				}

				local ind = 1
				local error = 0
				forvalues i=1/`ntdevars' {
					cap confirm var ``ind'', exact
					if _rc {
						local errortxt "invalid tde(... ``ind'' ``=`ind'+1'' ...)"
						local error = 1
					}
					cap confirm num ``=`ind'+1''
					if _rc {
						local errortxt "invalid tde(... ``ind'' ``=`ind'+1'' ...)"
						local error = 1
					}
					tempvar tdeeffect`i'
					qui gen double `tdeeffect`i'' = ``ind''*``=`ind'+1''

					local ind = `ind' + 2
				}
				if `error' {
					di as error "`errortxt'"
					exit 198
				}
				local tde_linpred "`tdeeffect1'"
				if `ntdevars'>1 {
					forvalues k=2/`ntdevars' {
						local tde_linpred "`tde_linpred' + `tdeeffect`k''"
					}
				}

				qui gen double `tdexb' = `tde_linpred'
			}
			else qui gen double `tdexb' = 0

		//====================================================================================================================//
		//setup stuff

			quietly {

				gen double `stime' = .

				tempvar logu
				gen double `logu' = log(runiform())

				cap drop _survsim_rc
				gen _survsim_rc = 0

			}

		//====================================================================================================================//
		//core

			if "`loghazard'"!="" | "`hazard'"!="" {

				gaussquad_ss, n(`nodes')	//Gauss-Legendre nodes and weights

				if "`loghazard'"!="" {
					local hazard exp(`loghazard')
				}

				mata: st_local("matahazard",subinstr("`hazard'","{t}","tnodes"))

				local matahazard "(`matahazard') :* expxb"

				if "`tde'"!="" {
					if "`tdefunction'"=="" {
						local tdefunction tnodes
					}
					else {
						mata: st_local("tdefunction",subinstr("`tdefunction'","{t}","tnodes"))
					}

					local matahazard "`matahazard' :* exp(tdexb :* (`tdefunction'))"
				}

				//push variables found in loghazard()/hazard() into Mata
                local Nhvars = 0
				if "`mixture'"=="" {

					macro drop overallsyntax1 overallsyntax2 mmrootsyntax1 mmrootsyntax2

					gettoken first rest : matahazard, parse("[ ,\^\*\(\)-\+/:<>=]")
					while "`rest'"!="" {
						if trim("`first'")!="," {
							cap confirm var `first', exact
							if !_rc {
								local test1 = 0
								local hvarindex = 0
								foreach var in `covlist' {
									local hvarindex = `hvarindex' + 1
									if "`first'"=="`var'" {
										local test1 = 1
										local refhvarindex = `hvarindex'
									}
								}
								if `test1'==0 {
									local covlist `covlist' `first'		//contains a list of all varnames specified, be they time-indep or time-dependent
									local Nhvars = `Nhvars' + 1
									mata: st_local("matahazard",subinstr("`matahazard'","`first'","hvars[,`Nhvars']"))
								}
								else {
									mata: st_local("matahazard",subinstr("`matahazard'","`first'","hvars[,`refhvarindex']"))
								}
							}
						}
						gettoken first rest : rest, parse("[ ,\^\*\(\)-\+/:<>=]")
					}
					if `Nhvars' {
						mata: st_view(hvars=.,.,tokens("`covlist'"))
						global overallsyntax1 real matrix hvars
						global overallsyntax2 hvars
						if "`marginal'"=="" global mmrootsyntax1 , hvars[i,]
						else				global mmrootsyntax1 , hvars
						global mmrootsyntax2 , hvars
					}
				}

				//test the hazard function
				if `Nhvars'==0 {
					mata: tnodes = expxb = tdexb = 0.1
					cap mata: test1 = `matahazard'
					if _rc {
						di as error "Error in loghazard()/hazard()"
						exit 198
					}
					mata mata drop tnodes expxb tdexb test1
				}

				//push function into global
				global chaz `matahazard'
				cap pr drop survsim_user_core

				survsim_user_core, 	stime(`stime') 				///
									maxtime(`maxtime') 			///
									logu(`logu') 				///
									expxb(`expxb')				///
									tdexb(`tdexb')				///
									ltruncated(`ltruncated') 	///
									`marginal'					//

				//tidy up
				cap macro drop cumhaz
				if `Nhvars' {
					cap macro drop overallsyntax1 overallsyntax2 mmrootsyntax1 mmrootsyntax2
				}
			}

			//cumulative hazard
			else {

				if "`logchazard'"!="" {
					local chazard exp(`logchazard')
				}
				mata: st_local("matachazard",subinstr("`chazard'","{t}","t"))
				local matachazard "(`matachazard') :* expxb"
				if "`ltruncated'"!="" {
					mata: st_local("matachazard0",subinstr("`chazard'","{t}","lt"))
					local matachazard0 "(`matachazard0') :* expxb"
				}

				//tde's
				if "`tde'"!="" {
					if "`tdefunction'"=="" {
						local tdefunction1 t
					}
					else {
						mata: st_local("tdefunction1",subinstr("`tdefunction'","{t}","t"))
					}
					local matachazard "`matachazard' :* exp(tdexb :* (`tdefunction1'))"
					if "`ltruncated'"!="" {
						if "`tdefunction'"=="" {
							local tdefunction0 lt
						}
						else {
							mata: st_local("tdefunction0",subinstr("`tdefunction'","{t}","lt"))
						}
						local matachazard0 "`matachazard0' :* exp(tdexb :* (`tdefunction0'))"
					}
				}

				//any variables in ch() or logch()
				local nhazvars = 0
				macro drop overallsyntax1 overallsyntax2 mmrootsyntax1 mmrootsyntax2
				gettoken first rest : matachazard, parse("[ ,\^\*\(\)-\+/:<>=]")
				while "`rest'"!="" {
					if trim("`first'")!="," {
						cap confirm var `first', exact
						if !_rc {
							local test1 = 0
							foreach var in `covlist' {
								if "`first'"=="`var'" local test1 = 1
							}
							if `test1'==0 {
								local covlist `covlist' `first'		//contains a list of all varnames specified, be they time-indep or time-dependent
								local nhazvars = `nhazvars' + 1
								mata: st_local("matachazard",subinstr(st_local("matachazard"),"`first'","hvars[1,`nhazvars']"))
								if "`ltruncated'"!="" {
									mata: st_local("matachazard0",subinstr(st_local("matachazard0"),"`first'","hvars[1,`nhazvars']"))
								}
							}
						}
					}
					gettoken first rest : rest, parse("[ ,\^\*\(\)-\+/:<>=]")
				}
				if `nhazvars' {
					mata: st_view(hvars=.,.,tokens("`covlist'"))
					global overallsyntax1 real matrix hvars
					global overallsyntax2 hvars
					if "`marginal'"=="" global mmrootsyntax1 , hvars[i,]
					else				global mmrootsyntax1 , hvars
					global mmrootsyntax2 , hvars
				}

				//test
				if `nhazvars'==0 {
					mata: t = expxb = tdexb = 0.1
					cap mata: test1 = `matachazard'
					if _rc {
						di as error "Error in logchazard()/chazard()"
						exit 198
					}
					mata mata drop t expxb tdexb test1
				}

				//push function into global
				global chaz `matachazard'
				if "`ltruncated'"!="" {
					global chaz ${chaz} :- `matachazard0'
				}
				cap pr drop survsim_user_core

				survsim_user_core, 	stime(`stime') 				///
									maxtime(`maxtime') 			///
									logu(`logu') 				///
									expxb(`expxb')				///
									tdexb(`tdexb')				///
									ltruncated(`ltruncated') 	///
									chazard						///
									`marginal'					//

				//tidy up
				cap macro drop chaz
				if `nhazvars' {
					cap macro drop overallsyntax1 overallsyntax2 mmrootsyntax1 mmrootsyntax2
				}

			}

	//done

end

program define gaussquad_ss, rclass
	syntax [, N(int 30)]
	tempname weights nodes
	mata ss_gq("`weights'","`nodes'")
	return matrix weights = `weights'
	return matrix nodes = `nodes'
end

mata:
void ss_gq(string scalar weightsname, string scalar nodesname)
{
	n 	=  strtoreal(st_local("n"))
	i 	= range(1,n,1)'
	i1 	= range(1,n-1,1)'

	muzero = 2
	a = J(1,n,0)
	b = i1:/sqrt(4 :* i1:^2 :- 1)

	A= diag(a)
	for(j=1;j<=n-1;j++){
		A[j,j+1] = b[j]
		A[j+1,j] = b[j]
	}
	symeigensystem(A,vec,nodes)
	weights = (vec[1,]:^2:*muzero)'
	weights = weights[order(nodes',1)]
	nodes = nodes'[order(nodes',1)']
	st_matrix(weightsname,weights)
	st_matrix(nodesname,nodes)
}

end