mata:
mata clear

void MYNLExp(real scalar todo, real vector b, 	///
	real vector y, real matrix X,		///
	val, grad, hess)
{
	real vector  r, f
	real matrix  df

	f   = exp(X*b')
	r   = y - f 
	val = -(r:^2)
	df  = f:*X

	if (todo>=1) {
		grad = r:*df
	}
	if (todo==2) {
		hess = -1*quadcross(df, df)
	}
}

void MYNLProbit(real scalar todo, real vector b,  ///
	real vector y, real matrix X,	 	  ///
	val, grad, hess)
{
	real vector  r, f, xb
	real matrix  df

	xb  = X*b'
	f   = normal(xb)
	r   = y - f 
	val = -(r:^2)
	df  = normalden(xb):*X

	if (todo>=1) {
		grad = r:*df
	}
	if (todo==2) {
		hess = -1*quadcross(df, df)
	}
}
 
void MYNLWork( string scalar depvar,  string scalar indeps,
             string scalar touse,   string scalar constant,
             string scalar bname,   string scalar Vname,
             string scalar nname,   string scalar rname,
	     string scalar model)
{

    real vector       y, b
    real matrix       X, V
    real scalar       n, p, ssr
    string scalar     emsg
    pointer(function) f
    transmorphic      S

    if (model=="expm") {
        f = &MYNLExp()
    }
    else if (model=="probit") {
        f = &MYNLProbit()
    }
    else {
        emsg = "{red}model " + model + " invalid\n"
        printf(emsg)
	exit(error(498))
    }
    y = st_data(., depvar, touse)
    n = rows(y)
    X = st_data(., indeps, touse)
    if (constant == "") {
        X = X,J(n, 1, 1)
    }
    p = cols(X)

    S = optimize_init()
    optimize_init_argument(S, 1, y)
    optimize_init_argument(S, 2, X)
    optimize_init_evaluator(S, f)
    optimize_init_params(S, J(1, p, .01))
    optimize_init_evaluatortype(S, "gf2")
    optimize_init_conv_vtol(S, 1e-10)
    b   = optimize(S)
    V   = invsym(-1*optimize_result_Hessian(S))
    ssr = (-1/(n-p))*optimize_result_value(S)
    V   = ssr*V

    st_matrix(bname, b)
    st_matrix(Vname, V)
    st_numscalar(nname, n)
    st_numscalar(rname, p)
}

mata mlib create lmynllib, replace
mata mlib add    lmynllib MYNLWork() MYNLProbit() MYNLExp()

end