*! Alluvial v1.21 (19 Oct 2023).
*! Asjad Naqvi

* v1.21 (19 Oct 2023): Fixed the showmiss bug (reported by Matthias Schonlau)
* v1.2  (04 Apr 2023): Minor fixes. If/in added back in.
* v1.1  (15 Jan 2023): fix label pass through. Weights added. offset added. valcond is just numeric. missing now has a color.
* v1.0  (10 Dec 2022): Beta release.

cap program drop alluvial


program alluvial, sortpreserve

version 15

	syntax varlist [if] [in] [aw fw pw iw/],  ///
		[ palette(string) colorby(string) smooth(numlist >=1 <=8) gap(real 2) RECENter(string) SHAREs showmiss alpha(real 75) ]  ///
		[ CATGap(string) CATSize(string)    ]  ///
		[ LABAngle(string) LABSize(string) LABPOSition(string) LABGap(string) SHOWTOTal  ] ///
		[ VALSize(string)   VALFormat(string) VALGap(string) NOVALues  ]  ///
		[ LWidth(string) LColor(string)  ]  ///
		[ VALCONDition(real 0) offset(real 0) ]  ///  // v1.1
		[ BOXWidth(string) ] /// // v1.2
		[ title(passthru) subtitle(passthru) note(passthru) scheme(passthru) name(passthru) xsize(passthru) ysize(passthru)		]


	// check dependencies
	cap findfile colorpalette.ado
	if _rc != 0 {
		display as error "The palettes package is missing. Install the {stata ssc install palettes, replace:palettes} and {stata ssc install colrspace, replace:colrspace} packages."
		exit
	}


	marksample touse, strok

qui {
preserve

	*keep if `touse'   // this is dropping missing values. DO NOT ENABLE
	keep `varlist' `exp'

	foreach x of local varlist {

		cap confirm numeric var `x'  // convert from numeric to string.
		if _rc!=0 {
			decode `x', gen(`x'2)
			drop `x'
			ren `x'2 `x'
		}
		else {
			qui levelsof `x'

			if `r(r)' > 20 {
				di as error "`x' has more than 20 categories. The variable might be continuous."
				di as error "Please simplify or drop the variable."
				exit
			}
		}
	}


	// store labels for later passthru
	foreach x of local varlist {
		if "`: var label `x''" != "" {
			local mylab`x' : var label `x'
		}
		else {
			local mylab`x' `x'
		}
	}


	gen temp = 1
	local items : word count `varlist'
	local items2 = `items' - 1


	tokenize `varlist'

	local obs = _N


	forval i = 1/`items2' {

		local j = `i' + 1
		gen f`i' = ``i''
		gen t`i' = ``j''

		gen catf`i' = "``i''"
		gen catt`i' = "``j''"


		if "`: value label ``i'''" != "" {
			decode ``i'', gen(labf`i')
		}
		else {
			gen labf`i' = string(``i'')
		}

		if "`: value label ``j'''" != "" {
			decode ``j'', gen(labt`i')
		}
		else {
			gen labt`i' = string(``j'')
		}

	}

	drop `varlist'

	drop temp
	gen id = _n
	order id

	reshape long f t labf labt catf catt, i(id) j(layer)

	drop id
	sort layer f
	gen value = 1

	if "`showmiss'" != "" {
		replace labf = "_missing" if f==.
		replace labt = "_missing" if t==.

		replace f = 99999 if f==.
		replace t = 99999 if t==.

	}
	else {
		drop if f==.
		drop if t==.
	}

	if "`weight'" != "" {
		local myweight  [`weight' = `exp']
	}


	collapse (sum) value `myweight', by(f t layer labt labf catf catt)
	sort layer f t


	if "`shares'" != "" {
		if "`weight'" == "" {
			replace value =  (value / `obs')
			format value %5.2f
		}
		else {
			bysort layer: egen mysum = sum(value)
			replace value =  (value / mysum)
		}
	}


	**** SANKEY ROUTINE BELOW

	ren layer x1
	summ x1, meanonly
	replace x1 = x1 - r(min) // rebase to 0

	gen x2 = x1 + 1

	ren value val1
	gen val2 = val1

	ren catf cat1
	ren catt cat2

	ren labf lab1
	ren labt lab2

	sort x1 f t  // this affects the draw order

	gen id = _n
	order id

	egen grp1 = group(x1 f)  // draw order
	egen grp2 = group(x1 t)

	sort x1 grp1 grp2
	by x1: gen y1 = sum(val1)

	sort x2 grp2 grp1
	by x2: gen y2 = sum(val2)

	sort x1 f t

	gen layer = x1 + 1

	ren f flo1
	ren t flo2


	reshape long x flo val grp y cat lab, i(id layer) j(tt)
	drop tt


	order grp id lab x
	sort id x lab


	sort layer x y

	by layer x: gen y1 = y[_n-1]
	recode y1 (.=0)
	gen y2 = y
	drop y

	order layer grp id lab x y1 y2 val

	// mark the highest value and the layer

	levelsof x, local(lvls)

	local hival = 0 // track the value

	foreach x of local lvls {

		summ y2 if x==`x', meanonly
		if r(max) > `hival' {
			local hilayer = `x'
			local hival = r(max)
		}
	}


	*** add gaps

	local propgap = `hival' * `gap' / 100

	sort x layer flo id
	egen tag = tag(x flo)
	bysort x: replace tag = sum(tag)
	cap drop offset
	gen offset = (tag - 1) * `propgap'
	replace y1 = y1 + offset
	replace y2 = y2 + offset

	drop tag offset
	encode lab, gen(order)


	*** transform the groups to be at the mid points

	sort x id y1 y2
	gen y1t = .
	gen y2t = .


	levelsof layer
	local tlayers = r(r) - 1


	forval i = 1/`tlayers' {

	local left  = `i'
	local right = `i' + 1


	levelsof flo if layer== `left', local(lleft)  // y:   to in the  first cut
	levelsof flo if layer==`right', local(lright) // x: from in the second cut


	foreach y of local lleft {  // left
		foreach x of local lright {      // right


			if `x' == `y' & (`x'!=99999 & `y'!=99999) {  // check if the groups are equal

				// in layer range
				summ y1 if flo==`x' & layer==`left' & x==`left', meanonly
				if r(N) > 0 {
					local y1max `r(max)'
					local y1min `r(min)'
				}
				else {
					local y1max 0
					local y1min 0
				}

				summ y2 if flo==`x' & layer==`left' & x==`left', meanonly
				if r(N) > 0 {
					local y2max `r(max)'
					local y2min `r(min)'
				}
				else {
					local y2max 0
					local y2min 0
				}

				local l1max = max(`y1max',`y2max')
				local l1min = min(`y1min',`y2min')

				// out layer range
				summ y1 if flo==`x' & layer==`right' & x==`left', meanonly
				if r(N) > 0 {
					local y1max `r(max)'
					local y1min `r(min)'
				}
				else {
					local y1max 0
					local y1min 0
				}

				summ y2 if flo==`x' & layer==`right' & x==`left', meanonly
				if r(N) > 0 {
					local y2max `r(max)'
					local y2min `r(min)'
				}
				else {
					local y2max 0
					local y2min 0
				}

				local l2max = max(`y1max',`y2max')
				local l2min = min(`y1min',`y2min')


				// calculate the displacement
				local displace = ((`l1max' - `l1min') - (`l2max' - `l2min')) / 2


				// displace the top and bottom parts
				replace y1t = y1 + `displace' + `l1min' - `l2min' if flo==`x' & layer==`right' & x==`left'
				replace y2t = y2 + `displace' + `l1min' - `l2min' if flo==`x' & layer==`right' & x==`left'

				}
			}
		}
	}


	replace y1t = y1 if y1t==.
	replace y2t = y2 if y2t==.

	drop y1 y2

	ren y1t y1
	ren y2t y2


	*** recenter

	levelsof x, local(lvls)

	foreach x of local lvls {

		qui summ y1 if x==`x', meanonly
		local ymin = r(min)
		qui summ y2 if x==`x', meanonly
		local ymax = r(max)


		if "`recenter'" == "bot" {
			local displace = 0
		}


		if "`recenter'" == "" | "`recenter'" == "mid" {
			local displace = (`hival' - `ymax') / 2
		}

		if "`recenter'" == "top" {
			local displace = `hival' - `ymax'
		}

		replace y1 = y1 + `displace' if x==`x'
		replace y2 = y2 + `displace' if x==`x'

	}

	cap drop sums
	bysort layer x lab: egen sums = sum(val)


	*** generate the curves
	local newobs = 30
	expand `newobs'
	sort id x
	cap drop xtemp
	bysort id: gen xtemp =  (_n / (`newobs' * 2))

	if "`smooth'" == "" local smooth = 4

	gen ytemp =  (1 / (1 + (xtemp / (1 - xtemp))^-`smooth'))

	gen y1temp = .
	gen y2temp = .


	levelsof layer	, local(cuts)
	levelsof id		, local(lvls)

	foreach x of local lvls {   // each id is looped over

		foreach y of local cuts {

			summ ytemp if id==`x' & layer==`y', meanonly

			if r(N) > 0 {
				local ymin = r(min)
				local ymax = r(max)
			}
			else {
				local ymin = 0
				local ymax = 0
			}

			sum x if layer==`y', meanonly
				local x0 = r(min)
				local x1 = r(max)


			summ y1 if id==`x' & x==`x0' & layer==`y', meanonly
			if r(N) > 0 {
				local y1min = r(min)
			}
			else {
				local y1min = 0
			}

			summ y1 if id==`x' & x==`x1' & layer==`y', meanonly
			if r(N) > 0 {
				local y1max = r(max)
			}
			else {
				local y1max = 0
			}

			replace y1temp = (`y1max' - `y1min') * (ytemp - `ymin') / (`ymax' - `ymin') + `y1min' if id==`x' & layer==`y'

			summ y2 if id==`x' & x==`x0' & layer==`y', meanonly
			if r(N) > 0 {
				local y2min = r(min)
			}
			else {
				local y2min = 0
			}

			summ y2 if id==`x' & x==`x1' & layer==`y', meanonly
			if r(N) > 0 {
				local y2max = r(max)
			}
			else {
				local y2max = 0.0000001
			}

			replace y2temp = (`y2max' - `y2min') * (ytemp - `ymin') / (`ymax' - `ymin') + `y2min' if id==`x' & layer==`y'
		}
	}



	replace xtemp = xtemp + layer - 1


	***** mid points for wedges


	egen tag = tag(x flo)

	cap gen midy = .

	levelsof x, local(lvls)
	foreach x of local lvls {

	levelsof flo if x ==`x', local(odrs)

		foreach y of local odrs {

		summ y1 if x==`x' & flo==`y', meanonly
		local min = r(min)

		summ y2 if x==`x' & flo==`y', meanonly
		local max = r(max)

		replace midy = (`min' + `max') / 2 if tag==1 & x==`x' & flo==`y'

		}
	}


	***** mid points for sankey labels

	egen tagp = tag(id)
	cap gen midp   = .  // outgoing
	cap gen midpin = .  // incoming
	cap gen xin = .


	levelsof id, local(lvls)
	foreach x of local lvls {

	// outbound values
		summ x if id==`x', meanonly
		local xval = r(min)

		summ y1 if id==`x' & x==`xval', meanonly
		local min = r(min)

		summ y2 if id==`x'  & x==`xval', meanonly
		local max = r(max)

		replace midp = (`min' + `max') / 2 if id==`x' & tagp==1

	// inbound values
		summ x if id==`x', meanonly
		local xval = r(max)

		summ y1 if id==`x' & x==`xval', meanonly
		local min = r(min)

		summ y2 if id==`x'  & x==`xval', meanonly
		local max = r(max)

		replace midpin = (`min' + `max') / 2 if id==`x' & tagp==1
		replace xin = x + 1 if id==`x' & tagp==1
	}


	*** fix boxes

	sort layer grp x y1 y2
	bysort x order: egen ymin = min(y1)
	bysort x order: egen ymax = max(y2)

	egen wedge = group(x flo)
	egen tagw = tag(wedge)

	replace lab = "{it:missing}" if flo== 99999

	*** tag the category labels

	summ ymin, meanonly
		local catmin = r(min)
	summ ymax, meanonly
		local catmax = r(max)

	egen tagc = tag(cat)
	if "`catgap'"		== "" local catgap 2

	gen grpy = `catmin' - ((`catmax' - `catmin') * `catgap' / 100) if tagc==1

	gen grpnm = ""

	foreach x of local varlist {
		replace grpnm = "`mylab`x''" if cat=="`x'" & tagc==1

	}


	sort x flo

	if "`boxwidth'"    == "" local boxwidth 3.2


	********************
	*** final plot   ***
	********************

		if "`palette'" == "" {
			local palette tableau
		}
		else {
			tokenize "`palette'", p(",")
			local palette  `1'
			local poptions `3'
		}


		if "`colorby'" == "layer" | "`colorby'" == "level" {
			local switch 1
		}
		else {
			local switch 0
		}


	sort layer x order xtemp y1temp y2temp


	// boxes

	local boxes

	levelsof wedge, local(lvls)
	local items = r(r)


	foreach x of local lvls {


		if `switch' == 1 { // by layer
			summ x if wedge==`x', meanonly
			local clr = r(mean) + 1
		}
		else {  			// by category
			summ order if wedge==`x', meanonly
			local clr = r(mean)
		}


		colorpalette `palette' , nograph `poptions'
		local boxes `boxes' (rspike ymin ymax x if wedge==`x' & tagw==1, lcolor("`r(p`clr')'%100") lw(`boxwidth')) ||

	}

	// arcs

	if "`lcolor'"  == "" local lcolor white
	if "`lwidth'"  == "" local lwidth none

	levelsof wedge
	local groups = r(r)

	local shapes


	levelsof id, local(lvls)

	foreach x of local lvls {

		if `switch' == 1 { 	// by layer
			qui sum layer if id==`x'
		}
		else {  			// by category
			qui sum x if id==`x'
			qui sum order if id==`x' & x == r(min)
		}



		if r(N) > 0 {
			local clr = r(mean)
			colorpalette `palette' , nograph `poptions'

			local shapes `shapes' (rarea y1temp y2temp xtemp if id==`x', lc(`lcolor') lw(`lwidth') fi(100) fcolor("`r(p`clr')'%`alpha'"))  ||
		}
	}


	**** PLOT EVERYTHING ***

	if "`catsize'"		== "" local catsize 2.3

	if "`labangle'" 	== "" local labangle 90
	if "`labsize'" 		== "" local labsize 2
	if "`labposition'"  == "" local labposition 0
	if "`labgap'" 		== "" local labgap 0


	if "`valsize'"  == "" local valsize 1.5

	local labcon "if val >= `valcondition'"


	if "`valformat'" == "" {
		if "`shares'"=="" {
			local valformat "%12.0fc"
		}
		else {
			local valformat "%4.3f"
		}
	}



	format val `valformat'

	if "`valgap'" 	 == "" local valgap 2

	summ ymax, meanonly
	local yrange = r(max)

	if "`showtotal'" != "" {
		gen lab2 = lab + " (" + string(sums, "`valformat'") + ")"
		local lab lab2
	}
	else {
		local lab lab
	}

	if "`novalues'" == "" {
		local values `values' (scatter midp   x   `labcon', msymbol(none) mlabel(val) mlabsize(`valsize') mlabpos(3)             mlabgap(`valgap')                       mlabcolor(black)) ///

		local values `values' (scatter midpin xin `labcon', msymbol(none) mlabel(val) mlabsize(`valsize') mlabpos(9)             mlabgap(`valgap')                       mlabcolor(black)) ///

	}


	// offset

	summ x, meanonly
	local xrmin = r(min)
	local xrmax = r(max) + ((r(max) - r(min)) * `offset' / 100)


	// FINAL PLOT //


	twoway ///
		`shapes' ///
			`boxes' ///
			(scatter midy x if  tag==1, msymbol(none) mlabel(`lab') mlabsize(`labsize') mlabgap(`labgap') mlabpos(`labposition') mlabangle(`labangle') mlabcolor(black)) ///
			`values' ///
			(scatter grpy x if tagc==1, msymbol(none) mlabel(grpnm) mlabgap(`catgap') mlabsize(`catsize') mlabpos(0) mlabcolor(black)) ///
			, ///
				legend(off) ///
					xlabel(, nogrid) ylabel(0 `yrange', nogrid)     ///
					xscale(off range(`xrmin' `xrmax')) yscale(off)	 ///
					`title' `subtitle' `note' `scheme' `name'	///
					`xsize' `ysize'
*/
restore
}


end




*********************************
******** END OF PROGRAM *********
*********************************