*! joyplot v1.71 (03 Oct 2023)
*! Asjad Naqvi (asjadnaqvi@gmail.com)

* v1.71 (03 Oct 2023): Fixed a bug in single density joyplots for incorrect locals. Fixed a bugs in lines. Added n() option.
* v1.7  (14 Jul 2023): xline(), saving(), peaks, peaksize() options added. ridgeline duplicates joyplot.
* v1.62 (28 May 2023): change over() to by() to align it with other packages. added offset() and laboffset()
* v1.61 (01 Mar 2023): ylabel in densities fixed. normalize in densities fixed.
* v1.6  (05 Nov 2022): bug fixes
* v1.5  (03 Sep 2022): bandwidth fixed. xlabel is now passthru. defaults updated.
* v1.42 (22 Jun 2022): y-axis was bugged in the over plot
* v1.41 (20 Jun 2022): installations fix, numerical over fix.
* v1.4  (26 Apr 2022): axes reverse options added. various optimizations
* v1.3  (24 Apr 2022): stacked densities added. label placement optimized.
* v1.21 (15 Apr 2022): xsize/ysize added. ylabels on right option added.
* v1.2  (13 Apr 2022): xlabel angle, local normalization, lines only option added
* v1.1  (07 Apr 2022): several options added
* v1.0  (13 Dec 2021): first release

**********************************
* Step-by-step guide on Medium   *
**********************************

// if you want to go for even more customization, you can read this guide:

* COVID-19 visualizations with Stata Part 8: Ridgeline plots (Joy plots) (30 Oct, 2020)
* https://medium.com/the-stata-guide/covid-19-visualizations-with-stata-part-8-joy-plots-ridge-line-plots-dbe022e7264d


cap program drop joyplot


program joyplot, sortpreserve rclass

version 15

	syntax varlist(min=1 max=2 numeric) [if] [in], by(varname) [overlap(real 6) BWIDth(real 0.5) palette(string) alpha(real 80)  lines		] ///
		[ LColor(string) LWidth(string) YLABColor(string) YLABSize(real 1.7) YLABPOSition(string) 									] ///
		[ YLine YLColor(string) YLPattern(string) YLWidth(real 0.04) YREVerse XREVerse 												] ///
		[ xtitle(passthru) ytitle(passthru) xlabel(passthru) title(passthru) subtitle(passthru) note(passthru)	     				] ///
		[ scheme(passthru) name(passthru) aspect(passthru) xsize(passthru) ysize(passthru)											] ///
		[ NORMalize(str) rescale droplow  ] ///  // v1.6 options
		[ LABOFFset(real 0) OFFset(real 0) ]  ///  // v1.62 options
		[ xline(passthru) saving(passthru) PEAKs peaksize(real 0.2) ] ///    // v1.7 options
		[ n(real 50) ]  // v1.71 options

	// check dependencies
	capture 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
	}


	capture findfile labmask.ado
	if _rc != 0 {
		qui ssc install labutil, replace
		exit
	}

	if `overlap' < 1 {
		display as error "overlap() should be >= 1"
		exit
	}

	// local options

	if "`lcolor'" 	 == "" local lcolor 	white
	if "`lwidth'" 	 == "" local lwidth 	0.15
	if "`ylabcolor'" == "" local ylabcolor 	black
	if "`ylcolor'"	 == "" local ylcolor  	black
	if "`ylpattern'" == "" local ylpattern  solid
	if "`xreverse'"  != "" local xreverse 	xscale(reverse)
	if "`palette'" == "" {
		local palette CET C1
	}
	else {
		tokenize "`palette'", p(",")
		local palette  `1'
		local poptions `3'
	}


	marksample touse, strok

	local length : word count `varlist'


////////////////////
// two variables  //
////////////////////

if `length' == 2 {

	gettoken yvar xvar : varlist


qui {
	preserve

	keep if `touse'

	gen ones = 1
	bysort `by': egen counts = sum(ones)
	egen tag = tag(`by')
	summ counts, meanonly

	if r(min) < 10 {
		if "`droplow'" == "" {
			count if counts < 10 & tag==1
			di as error "Groups with errors:"
			noi list `by' if counts < 10 & tag==1
			di as error "`r(N)' over group(s) (`by') have fewer than 10 observations. Either clean them manually or use the {it:droplow} option to automatically filter them out."
			exit
		}
		else {
			count
			local obstot = r(N)

			count if counts < 10
			local obsdrop = r(N)

			drop if counts < 10

			local obsdiff = `obstot' - `obsdrop'
			di in yellow "Out of `obstot' observations, `obsdrop' dropped. `obsdiff' observations remaining."

		}
	}


	count
	if r(N) == 0 {
		di as error "All observations dropped. No groups fulfill the criteria for {stata help joyplot:joyplot}."
		exit
	}


	drop ones tag counts

	sort `by' `xvar'
	cap drop _fillin
	fillin `by' `xvar'

	cap drop _fillin


	tempvar myvar // duplicate
	gen `myvar' = `yvar'
	replace `myvar' = . if `myvar' < 0   // TODO: see how to deal with negative values. v1.6: Suggest rescale option


	cap confirm numeric var `by'
		if _rc!=0 {
			tempvar over2
			encode `by', gen(`over2')
			local by `over2'
		}
		else {
			tempvar tempov over2
			egen   `over2' = group(`by')

			if "`: value label `by''" != "" {
				decode `by', gen(`tempov')
				labmask `over2', val(`tempov')
			}
			local by `over2'
		}


	if "`yreverse'" != "" {

		clonevar over2 = `by'

		summ `by', meanonly
		replace `by' = r(max) - `by' + 1

		if "`: value label over2'" != "" {
			tempvar group2
			decode over2, gen(`group2')
			replace `group2' = string(over2) if `group2'==""
			labmask `by', val(`group2')
		}
		else {
			labmask `by', val(over2)
		}
	}


	// rescale v1.6 added

	if "`rescale'" != "" {
		summ `myvar', meanonly
		local gmin = r(min)
		local gmax = r(max)

		return local jmin `"`gmin'"'
		return local jmax `"`gmax'"'

		tempvar myvar2
		gen double `myvar2' = `myvar' - `gmin'
		local myvar `myvar2'

		noi di in yellow "Ridgelines have been rescaled to start from the minimum value. Values range from `gmin' to `gmax' (see {stata return list})."
	}


	// normalization v1.6 fix

	tempvar norm
	gen double `norm' = .

	if "`normalize'" == "global" |  "`normalize'" == ""  {

		levelsof `by', local(lvls)

		foreach x of local lvls {
			 summ `myvar', meanonly
			 replace `norm' = `myvar' / r(max) if `by'==`x'
		}
	}

	else  {

		levelsof `by', local(lvls)

		foreach x of local lvls {
			 summ `myvar' if `by'==`x', meanonly
			 replace `norm' = `myvar' / r(max) if `by'==`x'
		}
	}



	// y labels

	tempvar tag xpoint ypoint y0

	gen `xpoint' = .
	gen `y0' = 0

	egen `tag' = tag(`by')


	if ("`ylabposition'" == "") | ("`ylabposition'" == "left")  {
		qui summ `xvar', meanonly
		replace `xpoint' = r(min) - ((r(max) - r(min)) * 0.10) + `laboffset' if `tag'==1
	}
	if ("`ylabposition'" == "right")  {
		qui summ `xvar', meanonly
		replace `xpoint' = r(max) + ((r(max) - r(min)) * 0.01) + `laboffset' if `tag'==1
	}


	gen `ypoint' = .

	local mygraph
	local yli


	levelsof `by', local(lvls)
	local items = `r(r)'

	foreach x of local lvls {


	summ `by', meanonly
	local newx = r(max) + 1 - `x'

	tempvar y`newx'

    lowess `norm' `xvar' if `by'==`newx', bwid(`bwidth') gen(`y`newx'') nograph

    tempvar ybot`newx' ytop`newx'

	 gen double `ybot`newx'' =  `newx'/ `overlap'
	 gen double `ytop`newx'' = `y`newx'' + `ybot`newx''

	 if "`peaks'" != "" {
	 	summ `ytop`newx'', meanonly
		gen peak`newx' = 1 if `ytop`newx''==r(max)
	 }

	}

	summ `ybot1', meanonly
	local shift = r(mean)

	levelsof `by', local(lvls)

	foreach x of local lvls {
		replace `ybot`x'' = `ybot`x'' - `shift'
		replace `ytop`x'' = `ytop`x'' - `shift'




	}


	levelsof `by', local(lvls)
	local items = `r(r)'

	foreach x of local lvls {


		summ `by', meanonly
		local newx = r(max) + 1 - `x'

		if "`lines'" != "" {
			colorpalette `palette', n(`items') nograph `poptions'
			local mygraph `mygraph' line `ytop`newx'' `xvar', lc("`r(p`newx')'") lw(`lwidth') ||
		}
		else {
			colorpalette `palette', n(`items') nograph `poptions'
			local mygraph `mygraph' rarea  `ytop`newx'' `ybot`newx'' `xvar', fc("`r(p`newx')'%`alpha'") fi(100) lw(none) ||  line `ytop`newx'' `xvar', lc(`lcolor') lw(`lwidth') ||

		}

		replace `ypoint' = (`ybot`newx'' + 0.02) if `tag'==1 & `by'==`newx'

		if "`yline'" != "" {
			local yli `yli' (line `ybot`newx'' `xvar' if `by'==`newx', lp(`ylpattern') lc(`ylcolor') lw(`ylwidth')) ||

		}


		if "`peaks'"!= "" {
			local mypeaks `mypeaks' (scatter `ytop`newx'' `xvar' if peak`newx'==1 , msym(circle) msize(`peaksize') mcolor(black)) ||

		}

	}


	if "`ylabposition'" == "" | "`ylabposition'" == "left" {
		summ `xvar', meanonly
		local x1 = r(min) - ((r(max) - r(min)) * (`offset' / 100)) + `laboffset'
		local x2 = r(max)
	}

	if "`ylabposition'" == "right"  {
		summ `xvar', meanonly
		local x1 = r(min)
		local x2 = r(max) + ((r(max) - r(min)) * (`offset' / 100))  + `laboffset'

	}

	twoway  		///
		`yli'		///
		`mygraph'  	///
		`mypeaks'     ///
		(scatter `ypoint' `xpoint', mcolor(none) mlabcolor(`ylabcolor') mlabel(`by') mlabsize(`ylabsize') ) ///
		, ///
			`xlabel' xscale(range(`x1' `x2')) `xreverse' ///
			ylabel(, nolabels noticks nogrid) yscale(noline) ///
				legend(off) ///
				`title' `subtitle' `xtitle' `ytitle' `note' `xline' ///
				`aspect' `xsize' `ysize' `scheme' `name' `saving'

	restore
	}
}



///////////////////
// one variable  //
///////////////////

if `length' == 1 {


qui {
preserve

	keep if `touse'
	keep `varlist' `by'

	gen ones = 1
	bysort `by': egen counts = sum(ones)
	egen tag = tag(`by')
	summ counts, meanonly

	if r(min) < 10 {
		if "`droplow'" == "" {
			count if counts < 10 & tag==1
			di as error "Groups with errors:"
			noi list `by' if counts < 10 & tag==1
			di as error "`r(N)' over group(s) (`by') have fewer than 10 observations. Either clean them manually or use the {it:droplow} option to automatically filter them out."
			exit
		}
		else {
			count
			local obstot = r(N)

			count if counts < 10
			local obsdrop = r(N)

			drop if counts < 10

			local obsdiff = `obstot' - `obsdrop'
			noi di in yellow "Out of `obstot' observations, `obsdrop' dropped. `obsdiff' observations remaining."

		}
	}

	count
	if r(N) == 0 {
		di as error "All observations dropped. No groups fulfill the criteria for {stata help joyplot:joyplot}."
		exit
	}

	drop ones tag counts


	// finetune the over variable

	cap confirm numeric var `by'

	if _rc!=0 {  // if string
		tempvar over2
		encode `by', gen(`over2')
		local by `over2'
	}
	else {

		tempvar tempov over2
		egen   `over2' = group(`by')   // group the categories


		if "`: value label `by''" != "" {
			decode `by', gen(`tempov')
			labmask `over2', val(`tempov')
		}
		else {
			tempvar bynames
			gen `bynames' = string(`by')
			labmask `over2', val(`bynames')
		}

		local by `over2'
	}


	if "`yreverse'" != "" {

		clonevar over2 = `by'

		summ `by', meanonly
		replace `by' = r(max) - `by' + 1

		if "`: value label over2'" != "" {
			tempvar group2
			decode over2, gen(`group2')
			replace `group2' = string(over2) if `group2'==""
			labmask `by', val(`group2')
		}
		else {
			labmask `by', val(over2)
		}
	}




	// counters
	local dmax  = 0

	summ `varlist', meanonly

	local xrmin = r(min)
	local xrmax = r(max)


	levelsof `by', local(lvls)
	local items = `r(r)'

	foreach x of local lvls {


		summ `by', meanonly
		local newx = r(max) + 1 - `x'   // reverse the sorting

		kdensity `varlist' if `by'==`x', generate(x`newx' y`newx') bwid(`bwidth') n(`n') nograph

		summ y`newx', meanonly
		if r(max) > `dmax' local dmax = r(max)   // global max

		summ x`newx', meanonly
			if r(min) < `xrmin' local xrmin = r(min)
			if r(max) > `xrmax' local xrmax = r(max)
	}

	// pad the minmax to avoid tight axes

	local pad = (`xrmax' -`xrmin') * 0.02

	local xrmin = `xrmin' - `pad'
	local xrmax = `xrmax' + `pad'


	// normalization
	if "`normalize'" == "" | "`normalize'" == "global"  {

		levelsof `by', local(lvls)

		foreach x of local lvls {
			 replace y`x' = y`x' / `dmax'
		}
	}

	if "`normalize'" == "local"  {

		levelsof `by', local(lvls)

		foreach x of local lvls {
			 summ y`x', meanonly
			 replace y`x' = y`x' / r(max)
		}
	}

	// y labels

	tempvar tag xpoint ypoint y0

	gen `xpoint' = .
	gen `y0' = 0

	egen `tag' = tag(`by')


	if ("`ylabposition'" == "") | ("`ylabposition'" == "left") {
		replace `xpoint' = `xrmin' - ((`xrmax' - `xrmin') * 0.10) + `offset' if `tag'==1
	}
	if ("`ylabposition'" == "right")  {
		replace `xpoint' = `xrmax' + ((`xrmax' - `xrmin') * 0.01) + `offset' if `tag'==1
	}


	gen `ypoint' = .

	local mygraph
	local yli


	levelsof `by', local(lvls)
	local items = `r(r)'

	foreach x of local lvls {

		summ `by', meanonly
		local newx = r(max) + 1 - `x'


		tempvar ybot`newx' ytop`newx'

		 gen double `ybot`newx'' =  `newx' / `overlap'     // if  y`newx'!=.
		 gen double `ytop`newx'' = y`newx' + `ybot`newx''  // if  y`newx'!=.

		if "`peaks'"!="" {
			summ `ytop`newx'', meanonly
			gen peak`newx' = 1 if `ytop`newx''== r(max)
		}
	}


	summ `ybot1', meanonly
	local shift = r(mean)

	levelsof `by', local(lvls)

	foreach x of local lvls {
		replace `ybot`x'' = `ybot`x'' - `shift'
		replace `ytop`x'' = `ytop`x'' - `shift'
	}




	levelsof `by', local(lvls)
	local items = `r(r)'


	foreach x of local lvls {


		summ `by', meanonly
		local newx = r(max) + 1 - `x'

		if "`lines'" != "" {
			colorpalette `palette', n(`items') nograph `poptions'
			local mygraph `mygraph' line `ytop`newx'' x`newx', lc("`r(p`newx')'") lw(`lwidth') ||
		}
		else {
			colorpalette `palette', n(`items') nograph `poptions'
			local mygraph `mygraph' rarea  `ytop`newx'' `ybot`newx'' x`newx', fc("`r(p`newx')'%`alpha'") fi(100) lw(none) ||  line `ytop`newx'' x`newx', lc(`lcolor') lw(`lwidth') ||

		}

		qui replace `ypoint' = (`ybot`newx'' + 0.02) if `tag'==1 & `by'==`x'

		if "`yline'" != "" {

			summ `ybot`newx'', meanonly
			local yvalue = r(mean)

			local yli `yli' (pci `yvalue' `xrmin' `yvalue' `xrmax', lp(`ylpattern') lc(`ylcolor') lw(`ylwidth')) ||

		}

		if "`peaks'"!= "" {

			local mypeaks `mypeaks' (scatter `ytop`newx'' x`newx' if peak`newx'==1 , msym(circle) msize(`peaksize') mcolor(black)) ||

		}


	}


	if "`ylabposition'" == "" | "`ylabposition'" == "left" {
		local x1 = `xrmin' - ((`xrmax' - `xrmin') * 0.12) + `offset'
		local x2 = `xrmax'
	}

	if "`ylabposition'" == "right"  {
		local x1 = `xrmin'
		local x2 = `xrmax' + ((`xrmax' - `xrmin') * 0.12) + `offset'
	}


	twoway  ///
		`yli'	    ///
		`mygraph'   ///
		`mypeaks'		///
		(scatter `ypoint' `xpoint', mcolor(none) mlabcolor(`ylabcolor') mlabel(`by') mlabsize(`ylabsize') ) ///
		, ///
			`xlabel' xscale(range(`x1' `x2'))  ///
			ylabel(, nolabels noticks nogrid) yscale(noline) ///
			legend(off) ///
			`title' `subtitle'  `xtitle' `ytitle' `xrev' `note' `xline'  ///
			`aspect' `xsize' `ysize' `scheme' `name' `saving'


	restore
	}


}


end



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