* 1.1		RAR 21 January 2004, with improved coding
* 1.0.1	RAR 12 December 2003, with suggestions from NJC
* 1.0		RAR 22 November 2003 - England won RU World Cup
program cvxhull, sortpreserve
	version 8.0
	syntax varlist(min=2 max=2) [if] [in] ///
	[ , Hulls(numlist int min=0 max=2 >0) MDPlot(int 8) ///
	 GROup(varname) SELect(numlist int min=0 >0 sort) ///
	 MEAns PREfix(string) noREPort noRETain ///
	 SCATopt(string) noGRAph SAVing(string)]

	loc reporting = ("`report'" != "noreport")
	if "`retain'" == "noretain" & "`graph'" == "nograph" {
		di as err "That combination of options is silly: quitting now"
		exit 498
	}

* Nick Cox thinks this parsing is still flakey 
	if `"`saving'"' != "" { //  Check saving filename   
		tokenize `"`saving'"', parse(",")
		args saving comma replace
		if index(`"`saving'"', ".") == 0 loc saving `"`saving'.gph"'
		capture confirm file `saving'
		if _rc == 0 { /* existing file */
			if "`replace'" != "replace" {
				di as err "Invalid save file - did you mean to replace?"
				exit 198
			}
			loc saving `"`saving',replace"'
		}
		else confirm new file `saving' // invalid file - unless new
		loc saving `"saving("`saving'")"'
	}

	if "`hulls'" == "" {
		loc hulls = 1 
		loc hullgap = 1 // Set default values
	}
	else {
		tokenize "`hulls'"
		loc hulls = `1'
		if   "`2'" == ""  loc hullgap = 1
		else              loc hullgap = `2'
	}

	marksample touse // Deal with `if', `in' & missing
	qui count if `touse'
	if r(N) == 0 error 2000

	tempvar sample grp scratch count pts hullno
	qui gen `sample'=`touse' // & remember all points considered
	if "`group'" == "" {     // group variable may be absent, string or numeric
		qui gen `grp' = 1
		loc select "1"
		loc maxgrp = 1
	}
	else {
		qui egen `grp' = group(`group') if `touse', label        
		if "`select'" != "" {
			loc j = 0
			foreach i in `select' {
				loc ++j
			}
			loc maxgrp = `j'  // no of items
		}
		else {
			su `grp', meanonly
			loc maxgrp = r(max)
			loc select "1 / `maxgrp'"
		}
		qui {
			egen `scratch' = eqany(`grp'), values(`select')
			replace `touse' = `touse' * `scratch'
		}
		if `reporting' di as txt "Codes for groups based on " as res "`group'"
		label list `grp'
	}

* Double check, in case group selection has dropped all obs
	qui count if `touse'
	if r(N) == 0 error 2000

	loc retd = 1 + int((`hulls'-1)/`hullgap')  // number of hulls
	loc hull = plural(`retd', "hull")
	loc gs = plural(`maxgrp', "group")
	if `reporting'  di as txt "Up to `retd' `hull' to be saved for `maxgrp' `gs'"

	if "`prefix'" == "" loc prefix "_cvxh"
	capture drop `prefix'*l
	capture drop `prefix'*r
	capture drop `prefix'grp
	capture drop `prefix'hull
	capture drop `prefix'pts
	capture drop `prefix'cnt
	capture drop `prefix'*mindex
	capture drop `prefix'*maxdex
	capture drop `scratch'

	tokenize `varlist'
	args y x 
	sort `grp' `x' `y' // within group, by ascending x then y

	qui {
		gen `count' = 0
		gen `pts' = 0
		gen `hullno' = 0
	}
	tempname gap
	loc gap = `hullgap' -1      // force retain hull 1
	loc maxhull = 0

***** Main loop ******************************************
	tempvar leftpath rightpath onhull  // re-use scratch for points on segment
	qui {
		gen `leftpath' = .    // set up outside loop to use replace
		gen `rightpath' = .   // within loop & save RAM allocation time
		gen `onhull' = 0
		gen `scratch' = 0
	}
	forvalues h = 1 / `hulls' { // b1 : brackets numbered for sanity! 
		qui count if `touse'    // stop when all points peeled
		if r(N) > 0  { // b2
		 	qui {
				replace `leftpath' = .
				replace `rightpath' = .
				replace `onhull' = 0
			}
			loc sp = 1          // Starting point
			loc notstarted = 1
			
			while `sp' <= _N  { // b3 scan all observations
				loc curgrp = `grp'[`sp']
				if !`touse'[`sp']  {
					loc ++sp    // ignore point & increment pointer
				}
				else {          // b4 - point processing ...
					while `grp'[`sp'] == `curgrp'  { // b5 within group
						if `notstarted' { // b6 mark first point of current group
							qui {
								replace `leftpath' = `y' in `sp'
								replace `rightpath' = `y' in `sp'
								replace `onhull' = (_n==`sp')
								replace `scratch' = 0
							} 
							loc leftcur = `sp'
							loc rightcur = `sp'
							loc sp1 = `sp'
							loc curgrp = `grp'[`sp']
							loc mindex = 1   
							loc notstarted = 0
							loc ++sp
						} 
						else  { // 6
							loc j = `leftcur' + 1
							loc maxL = -1
							loc maxD = 0
							loc found = 0 
							while `j'<=_N & `grp'[`j']==`curgrp' { // b7 find next left */
								if `touse'[`j']==1 {                 // b8
									loc d = sqrt( (`y'[`j'] - `y'[`leftcur'])^2 + (`x'[`j'] - `x'[`leftcur'])^2 )
									if float(`d') > 0 {  // b9
										loc cosa = (`y'[`j'] - `y'[`leftcur'])/ `d'
										if float(`cosa') > float(`maxL') { // b10 new leftmost direction
											loc maxL = `cosa'
											loc next = `j'
											qui replace `scratch' = (_n==`j')
											qui replace  `scratch' = 1 in `sp1'
											loc found = 1
										} 
										else {
											if float(`cosa') == float(`maxL') { // collinear
												if float(`d') > float(`maxD')  loc next = `j'
												qui replace `scratch' = 1 in `j'
												loc found = 1 
											}
										}
									} // 9
									else { // 9
										qui replace `onhull' = 1 in `j' // coincident with current point
										loc next = `j'
									} //9
								} // 8
								loc ++j
							} // 7
							qui replace `onhull' = (`onhull' | `scratch') 
							if `found' {
								qui {
									replace `leftpath' = `y' in `next'
									loc ++mindex
									loc leftcur = `next'
								} 
							} 
							
							loc j = `rightcur' + 1
							loc maxR = 1
							loc maxD = 0
							loc found = 0 
							while `j'<=_N & `grp'[`j']==`curgrp' { // b7 find next right
								if `touse'[`j']==1  { //8
									loc d = sqrt( (`y'[`j'] - `y'[`rightcur'])^2 + (`x'[`j'] - `x'[`rightcur'])^2 )
									if float(`d') > 0 { //9
										loc cosa = (`y'[`j'] - `y'[`rightcur']) / `d'
										if float(`cosa') < float(`maxR') { //10
											loc maxR = `cosa'
											loc next = `j'
											qui replace `scratch' = (_n==`j')
											qui replace  `scratch' = 1 in `sp1'
											loc found = 1 
										} 
										else {
											if float(`cosa') == float(`maxR')  { // collinear
												if float(`d') > float(`maxD')  loc next = `j'
												qui replace `scratch' = 1 in `j'
												loc found = 1
											} 
										}
									} //9
									else { //9
										qui replace `onhull' = 1 in `j'  // coincident
										loc next = `j'
									} //9
								} //8
								loc ++j
							} // 7
							qui replace `onhull' = (`onhull' | `scratch') 
							if `found' {
								qui {
									replace `rightpath' = `y' in `next' 
									loc ++mindex
									loc rightcur = `next'
								} 
							} 	
						
* Check if found a hull, or continue						
							loc sp = `leftcur'
							if `rightcur' < `leftcur' loc sp = `rightcur' 
							if `leftcur' == `rightcur'   { // hull closed
								qui {
									replace `hullno' = -`h' in `sp'
									replace `hullno' = `h' in `sp1'
									replace `count' = `mindex' - 1 in `sp1'  // last point double counted
									qui count if `onhull'
									replace `count' = r(N) in `sp'
									replace `pts' = `h' if `onhull'
									replace `touse' = 0 if `onhull'
								} 
								loc notstarted = 1
								while `grp'[`sp'] == `curgrp' {
									loc ++sp
								}
							} 
						} // e6
					} // e5 found end of group

					if `mindex' == 1  { // special case - single point hull
						qui {
							qui count if `onhull'
							replace `count' = r(N) in `sp1'
							replace `hullno' = `h' in `sp1'
                    		loc notstarted = 1
						} 
					} 
				} // e4 end point processing 
			} // e3 end loop over data
				
* Save hull if requested
			loc ++gap
			if `gap' == `hullgap' {
				loc ++maxhull
				qui {
				gen `prefix'`h'l = `leftpath'
				gen `prefix'`h'r = `rightpath'
				} 
				if `maxhull' <= `mdplot' loc hulllist "`prefix'1l-`prefix'`h'r"
				loc gap = 0                // restart count
				if `reporting' di as txt "  hull level `h' calculated and saved"
			} 
		} // e2 do nothing if no data points
	} // e1 exit loop when all hulls marked

***** Optional plot 
	if `"`graph'"' != "nograph" {
		if `reporting' di as txt "Graph will be plotted presently"
* Build and execute graph command ... as a very long macro text!
		loc colours "black black blue blue dkorange dkorange magenta magenta emerald emerald khaki khaki cyan cyan red red"
* Set up point plot and add line plots for each selected group
* RAR's preference for default horizontal y labels
		loc gr `"scatter `y' `x' if `sample',yti("`y'")yla(,angle(0))`saving' `scatopt'"'
		if "`means'" == "means" {
			tempvar ymean xmean
			egen `ymean' = mean(`y') if `select', by(`grp')
			egen `xmean' = mean(`x') if `select', by(`grp')
			loc gr `"`gr'||scatter `ymean' `xmean',ms(T)xti("`x'")"'
		} 
		foreach i of numlist `select' {
			loc gr `"`gr'||line `hulllist' `x' if `grp'==`i',clc(`colours')legend(off)"'
		} 
		`gr' // & execute macro command 
	} 
	
	if "`retain'" != "noretain" {
		qui {
			if "`group'" != "" gen `prefix'grp  = `grp'
			gen `prefix'hull = `hullno'
			gen `prefix'cnt  = `count'
			gen `prefix'pts = `pts'
		} 
	} 
	else {
		capture drop `prefix'*l `prefix'*r 
	}

	if `reporting' di as txt "cvxhull run"

end // of cvxhull