*! version 1.0.0  01aug2002  Matthias Schonlau   (SJ2-4: st0028)
program define clustergram
	version 7
	syntax varlist [if] [in], CLuster(varlist) [ fill FRaction(real 0.2) /*
*/ XLAbel(numlist) YLAbel(numlist sort) COnnect(str) Symbol(str)  Gap(int 5) * ]

	* trap -connect()-, -symbol()-
	foreach opt in connect symbol {
		if "``opt''" != "" {
			di as err "`opt'() not supported: please see help"
			exit 198
		}
	}

	* observations to use
	marksample touse
	qui count if `touse'
	if r(N) == 0 { error 2000 }

	* check whether jth variable defines j clusters
	local nc = 0
	foreach v of var `cluster' {
		local nc = `nc' + 1
		tab `v' if `touse', nofreq
		if `nc' != `r(r)' {
			di as txt "Warning: `v' does not define `nc' clusters"
		}
	}

	preserve
	qui keep if `touse'

	tempvar id clustery clmean clmeanlag width clus_nlag
	gen long `id'  = _n
	local N = _N

	* stack the y variables into one variable clustery.
	tempname clustery
	local k : word count `varlist'
	tokenize `varlist'
	forval i = 1 / `k' {
		rename ``i'' `clustery'`i'
	}
	qui reshape long `clustery', i(`id')

	local max : word count `cluster'

	* for each obs replace the cluster number with the cluster mean
	quietly foreach v of var `cluster' {
		tempvar clmean
		bysort `v' : egen `clmean' = mean(`clustery')
		local clmeans "`clmeans' `clmean'"
	}

	* data now contain (#id's * #y's) observations
	* the cluster mean is the same across all y's, can collapse
	* (mean) clustery is not used
	collapse (mean) `clustery', by(`id' `clmeans')

	* data now contain (#id's) observations and variables:
	* id, the clmeans, clustery

	* clus_n is the index of clmean
	foreach v of varlist `clmeans' {
		local args "`args' `v' `id'"
	}
	stack `args', into(`clmean' `id') clear
	local clus_n "_stack"
	qui {
		bysort `id' (`clus_n') : gen `clmeanlag' = `clmean'[_n+1]
		drop if `clmeanlag' == .
	}

	* clus_n is important in case there are the same means at
	* different cluster sizes
	contract `clmean' `clmeanlag' `clus_n'

	* preparation of graph
	label var `clus_n' "Number of clusters"
	label var `clmean' "Cluster mean"

	* create -ylabel()- -xlabel()- if none provided
	* set `range' and use it to standardize `width'
	if "`ylabel'" == "" {
		summarize `clmeanlag', meanonly
		local round = 1
		local miny = round(r(min),`round')
		local maxy = round(r(max),`round' )
		local midy = round(r(mean),`round')
		local ylabel "`miny' `midy' `maxy'"
		local range = r(max) - r(min)
	}
	else {
	* if `ylabel' supplied, reset `range' to depend on the range of `ylabel'
	* this causes different clustering algorithms with the same `ylabel'
	* to be standardized in the same way w.r.t width
		local ylamin : word 1 of `ylabel'
		local ylan : word count `ylabel'
		local ylamax : word `ylan' of `ylabel'
		local range = `ylamax' -`ylamin'
	}

	gen `width' =_freq / `N' * `range' * `fraction'

	if "`xlabel'" == "" { local xlabel "1 `max'" }

	graph `clmean' `clus_n' , /*
	*/ s(.) ylab(`ylabel') xlab(`xlabel')  gap(`gap') `options'

	* user typed -fill- or not => fill is 0 or 1
	local fill = "`fill'" != ""
	* note the following line requires that the clusters increase one
	* at a time
	gen `clus_nlag' = `clus_n' + 1
	Graphbox `clus_n' `clus_nlag' `clmean' `clmeanlag' `width' `fill' `range'
end

program define Graphbox
* draw a parallelogram (graph must exist already)
* xstart xend ystart yend 	(vectors)
*	start and end of the line segments for the center of the parallelogram
* width 	width of the parallelogram
* fill 	(scalar) if fill>0 then the parallelogram is filled with lines
* to make it appear solid
	args xstart xend ystart yend width fill range
	tempvar x1 x2 y1 y2 y3 y4
	gph open
	graph
	local ay = r(ay)
	local by = r(by)
	local ax = r(ax)
	local bx = r(bx)
	gen `x1' = `ax' * `xstart' + `bx'
	gen `x2' = `ax' * `xend' + `bx'

	* graph clockwise , starting at top left
	gen `y1' = `ay' * (`ystart' + `width' / 2) + `by'
	gen `y2' = `ay' * (`yend'   + `width' / 2) + `by'
	gen `y3' = `ay' * (`yend'   - `width' / 2) + `by'
	gen `y4' = `ay' * (`ystart' - `width' / 2) + `by'
	gph vpoly `y1' `x1' `y2' `x2' `y3' `x2' `y4' `x1' `y1' `x1'

	* summarize `yend', meanonly
	* local range = r(max) - r(min)
	if `fill' {
		* increment determines how densely each box is plotted
		local increment = .001 * `range'
		local N = _N

		forval i = 1 / `N' {
	   		local add = 0
   			while (`add' < `width'[`i']) {
				local y5 = /*
*/ `ay' * (`ystart'[`i'] - `width'[`i'] / 2 + `add') + `by'
				local y6 = /*
*/ `ay' * (`yend'[`i'] - `width'[`i'] / 2 + `add') + `by'
				local startx1 = `x1'[`i']
				local endx1 = `x2'[`i']
				gph line `y5' `startx1' `y6' `endx1'
				local add = `add' + `increment'
   			}
		}
	}
	gph close
end