program define simplex
        version 4.0


        local mldepn : word 1 of `*'
        global S_mldepn = "`mldepn'"
        pause Before Parse.
        local varlist "req ex"
	local if "opt"
	local in "opt"
	local options "FROM(string) SPAN(real .1) TOLerance(real .00001) NOCONStant ITERate(real 16000) LF(string) n_obs(real 0) STEP(real .1) TRACE"

        macro shift 1

        pause * is now `*'.

        parse "`*'"

        pause Parsed arguments.  Varlist is `varlist', mldepn is `mldepn', if is `if', in is `in'.

	tempname missing sampwgt vmat

	qui egen `missing' = rmiss(`varlist')
	qui generate `sampwgt' = 0
	qui replace `sampwgt' = 1 `if' `in'
	qui replace `sampwgt' = 0 if `missing' > 0
	global S_mlwgt "`sampwgt'"


        if "`noconst'" == "" { local varlist "`varlist' _cons" }
        local numvs : word count `varlist'
        local vecsize = colsof(`from') + 1

        pause Duh.

	if `numvs' != `vecsize' - 1 {
                display "Variable list does not match FROM list in dimension."
                exit(1)

                }

        pause Vecsize checks out as `vecsize', numvs as `numvs'.

	/* Note: Starting matrix variables are named by variable list */



	tempname normval curtl binit tomax
	local iter_no = 1
	scalar `normval' = 0

	matrix `binit' = `from'

	matrix colnames `binit' = `varlist'


	/* Take the vector in `from' , add `span' to one component, */
	/* and call it a new vector.   Do this for each component,  */
	/* and the resulting vectors plus the original will define  */
	/* the simplex   */


	local simp "`binit'"
	local i = 1
	while `i' <= `numvs' {
		tempname simp`i'
		matrix `simp`i'' = `binit'
		local ith_elt = el(`from',1,`i')
		local ith_elt = `ith_elt' + `span'
		matrix `simp`i''[1,`i'] = `ith_elt'
		local simp "`simp' `simp`i''"
		local i = `i' + 1
		}

	



	pause from is `from', trace is `trace', lf is `lf', simp is `simp'.

	norm `normval' , simplex(`simp') `trace' lf(`lf')


	pause After norm.  Normval is `normval'.

	display _n

	/* Begin the optimization loop */

	while (`normval' > `toleran') & `iter_no' <= `iterate' {

		pause Beginning of loop.

		local i = 1
		local min_pos = 1
		local min_lf = 2000000
		local max_pos = 1
		local max_lf = -2000000

		while `i' <= `vecsize' {

			local ith_vec : word `i' of `simp'
			`lf' `ith_vec' `curtl'
			if `curtl' < `min_lf' {
				local min_pos = `i'
				local min_lf = `curtl'
				}
			if `curtl' > `max_lf' {
				local max_pos = `i'
				local max_lf = `curtl'
				}

			local i = `i' + 1
			}

			pause min_pos is `min_pos'.
			display "Iteration `iter_no': Minimum likelihood `min_lf', maximum norm is " `normval'


		local lambda = `step'
		tempname concomb maxcomb combml ithtemp
		local maxml = -2000000000


		while `lambda' <= 1 {

			if "`trace'" != "" { display _n "Lambda: " `lambda' }

			matrix `concomb' = J(1,`numvs',0)
			local factor1 = (1 - `lambda') / (`vecsize' - 1)
			local factor2 = (1 + `lambda') / (`vecsize' - 1)

			local i = 1


			while `i' <= `vecsize' {

				local ith_vec : word `i' of `simp'

				if `i' != `min_pos' {
					matrix `ithtemp' = `factor1'* `ith_vec'
					}

				else {matrix `ithtemp' = `lambda' * `ith_vec' }

				matrix `concomb' = `concomb' + `ithtemp'

				local i = `i' + 1

				}


			`lf' `concomb' `combml'

			if "`trace'" != "" {
				matrix list `concomb'
				display _n "Likelihood: " `combml' _n
				}

			if `combml' > `maxml' {
				matrix `maxcomb' = `concomb'
				local maxml = `combml'
				}

			matrix `concomb' = J(1,`numvs',0)
			local i = 1

			while `i' <= `vecsize' {

				local ith_vec : word `i' of `simp'

				if `i' != `min_pos' {
					matrix `ithtemp' = `factor2'* `ith_vec'
					}
				else {matrix `ithtemp' = -`lambda' * `ith_vec'}

				matrix `concomb' = `concomb' + `ithtemp'

				local i = `i' + 1

				}


			`lf' `concomb' `combml'

			if "`trace'" != "" {
				matrix list `concomb'
				display _n "Likelihood: " `combml' _n
				}

			if `combml' > `maxml' {
				matrix `maxcomb' = `concomb'
				local maxml = `combml'
				}


			local lambda = `lambda' + `step'
			}

		if "`trace'" != "" {
			display "New matrix for step " `iter_no' ":"
			matrix list `maxcomb'
			}

		local min_vec : word `min_pos' of `simp'
		local max_vec : word `max_pos' of `simp'


		if `maxml' > `min_lf' {
			matrix `min_vec' = `maxcomb'
			}

		else {

			display "Nonconcavity encountered.  Flipping and shrinking simplex."


			local i = 1
			while `i' <= `vecsize' {
				local ith_vec : word `i' of `simp'
				matrix `tomax' = `max_vec' - `ith_vec'
				matrix `tomax' = .5*`tomax'
				matrix `ith_vec' = `max_vec' + `tomax'
				local i = `i' + 1
				}



			}



		norm `normval' , simplex(`simp') `trace' lf(`lf')

		pause After norm.  normval is `normval' , tolerance is `toleran' , iter_no is `iter_no', iterate is `iterate', simp is `simp'


		local iter_no = `iter_no' + 1

		}

	if `iter_no' > `iterate' {
		display "stopped at `iterate'th iteration." _n
		}

	display _n "Estimate using simplex, tolerance `toleran':"

	qui sum $S_mldepn if $S_mlwgt
	local n_obs = _result(1)

	matrix `vmat' = J(`numvs',`numvs',0)
	matrix rownames `vmat' = `varlist'
	matrix colnames `vmat' = `varlist'

	matrix post `max_vec' `vmat' , dep($S_mldepn) obs(`n_obs')

	display "Likelihood value  `max_lf'" 
	display "No. Observations  `n_obs'" _n
	matrix mlout

end