*! ML program that acts as the middle man between the stata ado lcmixlogit and the mata command lcmixl_ll.
program define mixmixlogit_d1
	version 12
	args todo b lnf g

	mata: mixmix_ll("`b'")
	scalar `lnf' = r(ll)

	if (`todo'==0 | `lnf'>=.) exit
	matrix `g' = r(gradient)
end

*! The mata command lcmixl_ll has the main function of calculating the log likelihood and gradient vector at each iteration.
version 12
mata:
void mixmix_ll(string scalar B_s)
{
	/* Import Variables from the Stata .ado and assign them internal names */
	external data_X
	external data_Y
	external data_ccov
	external data_CSID
	external macro_nrep
	external macro_np
	external macro_k
	external macro_burn
	external macro_class
	external macro_kclass
	external lognormal
	external macro_corr
	external macro_T
	external macro_ccovnum
	external macro_cutnum
	external gradpick

	nrep = macro_nrep
	np = macro_np
	k = macro_k
	burn = macro_burn
	ccovnum = macro_ccovnum
	cutnum = macro_cutnum
	classes = macro_class
	kclass = macro_kclass
	classesp = classes - 1
	B = st_matrix(B_s)'
	corr = macro_corr

	/* Split up the passed B matrix into the mean coefficients, standard errors, and class probabilities */

	MRND = B[|1,1\kclass,1|]

	if (corr == 1) {
		external macro_cho
		cho = macro_cho
		SRND = J(kclass,kclass,0)
		for(j=1; j<=classes; j++) {
			SRND[|(((j-1)*k)+1),(((j-1)*k)+1)\(j*k),(j*k)|] = invvech(B[|(kclass+(((j-1)*(cho/classes))+1)),1\(kclass+(j*(cho/classes))),1|]) :* lowertriangle(J(k,k,1))
		}
		if (ccovnum > 0) gammas = B[|(kclass+cho+1),1\(kclass+cho+ccovnum),1|]
		cuts = B[|(kclass+cho+ccovnum+1),1\(kclass+cho+ccovnum+cutnum),1|]
	}
	else {
		SRND = diag(B[|kclass+1,1\(kclass*2),1|])
		if (ccovnum > 0) gammas = B[|(kclass*2)+1,1\(kclass*2)+ccovnum,1|]
		cuts =   B[|(kclass*2)+ccovnum+1,1\(kclass*2)+ccovnum+cutnum,1|]
	}

	/* Set up empty matrices for the LL and gradient vector */
	p = J(np,classes,0)
	P = J(np,1,0)

	if (corr == 1) {
		G = J(np,(kclass+cho),0)
	}
	else {
		G = J(np,(kclass*2),0)
	}

	Gp = J(np,ccovnum+cutnum,0)
	Gpre = J(1,classes,0)
	COMP = J(k,nrep,0)

	/* Initial Halton sequence for later shuffling */
	rseed(1234567)
	firstrow = invnormal(halton(nrep,1,(1+burn)))
	ERR = firstrow'
	for (z=2; z<=kclass; z++) {
		if (MRND[z,1] == MRND[z-1,1] & SRND[z,1] == SRND[z-1,1]) {
			addz = ERR[|z-1,1\z-1,nrep|]
			ERR = ERR \ addz
		}
		else {
			addz = jumble(firstrow)
			ERR = ERR \ addz'
		}
	}

	i = 1
	/* Calculate choice probabilities and gradients for each individual */
	for (n=1; n<=np; n++) {

		/* Shuffle the Halton sequence */
		ERR[|1,1\1,nrep|] = jumble(ERR[|1,1\1,nrep|]')'
		for (z=2; z<=kclass; z++) {
			if (MRND[z,1] == MRND[z-1,1] & SRND[z,1] == SRND[z-1,1]) {
				ERR[|z,1\z,nrep|] = ERR[|z-1,1\z-1,nrep|]
			}
			else {
				ERR[|z,1\z,nrep|] = jumble(ERR[|z,1\z,nrep|]')'
			}
		}

		/* Calculate the simulated beta vector using the halton series, the means, and the standard errors of the coefficients. */
		BETA = MRND :+ (SRND*ERR)

		/* If some beta distributions are specified by the user to be positive or negative log-normal, transform them here. */
		for (j=1; j<=kclass; j++) {
			if (lognormal[1,j] != 0) {
				BETA[j,.] = lognormal[1,j]:*exp(BETA[j,.])
			}
		}

		/* Calculate some more stuff */
		if (corr == 1) {
			M = J((kclass+cho),nrep,0)
		}
		else {
			M = J((kclass*2),nrep,0)
		}
		R = J(1,nrep,0)
		nc = macro_T[i,1]
		istart = i

		/* Loop for each class */
		for(j=1; j<=classes; j++) {
			RJ = J(1,nrep,1)
			i = istart

			for (t=1; t<=nc; t++) {
				/* Grab the data */
				YMAT = data_Y[|i,1\(i+data_CSID[i,1]-1),cols(data_Y)|]
				XMAT = data_X[|i,1\(i+data_CSID[i,1]-1),cols(data_X)|]

				/* Choice probability for that class, time period, and individual */
				EV = exp(XMAT*BETA[|(((j-1)*k)+1),1\(j*k),nrep|])
				EV = (EV :/ colsum(EV))

				/* Multiply with other time periods for that specific class and individual */
				RJ = RJ :* colsum(YMAT :* EV)

				/* Calc M for gradient vector */
				PMAT = YMAT :- EV
				for (s=1; s<=k; s++) {
					COMP[s,.] = colsum(PMAT :* XMAT[.,s])
				}

				for(s=1; s<=k; s++) {
					if (lognormal[1,(((j-1)*k)+s)] != 0) {
						M[(((j-1)*k)+s),.] = M[(((j-1)*k)+s),.] :+ (COMP[s,.] :* BETA[(((j-1)*k)+s),.])
					}
					else {
						M[(((j-1)*k)+s),.] = M[(((j-1)*k)+s),.] :+ COMP[s,.]
					}
				}
				if (corr == 1) {
					num = 1
					for (l=1; l<=k; l++) {
						for (s=l; s<=k; s++) {
							if (lognormal[1,(((j-1)*k)+s)] != 0) {
								M[(((j-1)*(cho/classes))+num+kclass),.] = M[(((j-1)*(cho/classes))+num+kclass),.] :+ (COMP[s,.] :* BETA[(((j-1)*k)+s),.] :* ERR[(((j-1)*k)+l),.])
							}
							else {
								M[(((j-1)*(cho/classes))+num+kclass),.] = M[(((j-1)*(cho/classes))+num+kclass),.] :+ (COMP[s,.] :* ERR[(((j-1)*k)+l),.])
							}
							num = num + 1
						}
					}
				}
				else {
					for(s=1; s<=k; s++) {
						if (lognormal[1,(((j-1)*k)+s)] != 0) {
							M[(((j-1)*k)+s+kclass),.] = M[(((j-1)*k)+s+kclass),.] :+ (COMP[s,.] :* BETA[(((j-1)*k)+s),.] :* ERR[(((j-1)*k)+s),.])
						}
						else {
							M[(((j-1)*k)+s+kclass),.] = M[(((j-1)*k)+s+kclass),.] :+ (COMP[s,.] :* ERR[(((j-1)*k)+s),.])
						}
					}
				}

				i = i + data_CSID[i,1]
			}

			/* Calculate class probability denominator for this specific individual */
			if (ccovnum > 0) {
				CCOVMAT = mean(data_ccov[|istart,1\(i-1),cols(data_ccov)|],1)
				if (j == 1) {
					COMP4 = exp(cuts[1,1] - CCOVMAT*gammas)
					p[n,j] = COMP4 / (1 + COMP4)
				}
				else if (j == classes) {
					COMP5 = exp(cuts[classesp,1] - CCOVMAT*gammas)
					p[n,j] = 1 - COMP5 / (1 +COMP5)
				}
				else {
					COMP4 = exp(cuts[j-1,1] - CCOVMAT*gammas)
					COMP5 = exp(cuts[j,1] - CCOVMAT*gammas)
					p[n,j] = COMP5 / (1 + COMP5) - COMP4 / (1 + COMP4)
				}

			}
			else {
				if (j == 1) {
					COMP4 = exp(cuts[1,1])
					p[n,j] = COMP4 / (1 + COMP4)
				}
				else if (j == classes) {
					COMP5 = exp(cuts[classesp,1])
					p[n,j] = 1 - COMP5 / (1 +COMP5)
				}
				else {
					COMP4 = exp(cuts[j-1,1])
					COMP5 = exp(cuts[j,1])
					p[n,j] = COMP5 / (1 + COMP5) - COMP4 / (1 + COMP4)
				}
			}

			/* Aggregate over classes to obtain the choice probability for that individual (average over simulation) */
			P[n,1] = P[n,1] :+ (p[n,j] :* mean(RJ',1))

			/* Start calc for gradient vector for this specific class */
			if (corr == 1) {
				num = 1
				for (l=1; l<=k; l++) {
					G[n,(((j-1)*k)+l)] = p[n,j] :* mean((RJ :* M[(((j-1)*k)+l),.])',1)
					for (s=l; s<=k; s++) {
						G[n,(((j-1)*(cho/classes))+num+kclass)] = p[n,j] :* mean((RJ :* M[(((j-1)*(cho/classes))+num+kclass),.])',1)
						num = num + 1
					}
				}
			}
			else {
				for(s=1; s<=k; s++) {
					G[n,(((j-1)*k)+s)] = p[n,j] :* mean((RJ :* M[(((j-1)*k)+s),.])',1)
					G[n,(((j-1)*k)+s+kclass)] = p[n,j] :* mean((RJ :* M[(((j-1)*k)+s+kclass),.])',1)
				}

			}
			Gpre[1,j] = mean(colsum(RJ)',1)
		}

		/* Finish G calc now that statistics for all classes have been obtained */
		COMP6 = (1 :/ P[n,1])
		G[n,.] = COMP6 :* G[n,.]
		for (c=1; c<= (ccovnum+cutnum); c++) {
			if (c <= ccovnum) {
				for (j=1; j<=classesp; j++) {
					if (ccovnum > 0) COMPZ = exp(cuts[j,1] - CCOVMAT*gammas)
					else COMPZ = exp(cuts[j,1])
					Gp[n,c] = Gp[n,c] + ((CCOVMAT[1,c] * COMPZ) / (1 + COMPZ)^2) :* COMP6 :* Gpre[1,j+1] - ((CCOVMAT[1,c] * COMPZ) / (1 + COMPZ)^2) :* COMP6 :* Gpre[1,j]
				}
			}
			else {
				if (ccovnum > 0) COMPZ = exp(cuts[c-ccovnum,1] - CCOVMAT*gammas)
				else COMPZ = exp(cuts[c-ccovnum,1])
				Gp[n,c] = ((COMPZ)/(1 + COMPZ)^2) :* COMP6 :* Gpre[1,c-ccovnum] - ((COMPZ)/(1 + COMPZ)^2) :* COMP6 :* Gpre[1,c-ccovnum+1]
			}
		}
	}

	G = G , Gp

	for(j=1; j<=classes; j++) {
		if (colsum(p[.,j]) < 0.015) {
			ll_penalty = -99999999
		}
	}

	/* Output the log-likelihood and gradient vector. */
	st_numscalar("r(ll)", colsum(ln(P)))
	st_matrix("r(gradient)", colsum(G))
}

end