*! Description : Proportional Venn Diagrams
*! Version 1.6 22May2022 by Mead Over
*! Enhanced version of pvenn.ado by Wenfeng Gong; Jan Osterman (Duke Center for Health Policy)

program define pvenn2
version 11.0
preserve
set more off

syntax [varlist(default=none max=3)] [if] [, plabel(namelist min=1 max=3) ///
	pA(real 0) pB(real 0) pC(real 0) pAB(real 0) pAC(real 0) pBC(real 0) ///
	ptotal(real 0) ANGle(real 0) pos(integer 1) DRAWTotal(integer 1) DRAWLabel(integer 1) ///
	lc(string) lp(string) lw(string) n(integer 2000)   /// Following options added by M. Over
	pABC(real 0) nolabel text(passthru) legend(passthru) saving(passthru) name(passthru)  ///
	note(passthru) title(passthru) subtitle(passthru) scheme(passthru) plotregion(passthru) NOTEXoffset(passthru) debug UNDERscore *]

*	Graph options not spcified in syntax
local gopt "`options'"

if "`debug'"~="" {
	di as txt "options: " as res `"`title' `subtitle' `note' `legend' `text'  `name' `saving' `plotregion' `scheme' `debug' `options' "'
}

*	-legdefault- does not accept spaces in the names of the sets
foreach i of numlist 1/3 {
	local ABorC = word("A B C",`i')
	local name`ABorC' = cond(word("`plabel'",`i')~="",word("`plabel'",`i'),"Population_`ABorC'")
}

if "`debug'"~="" di as txt "nameA is " as res "`nameA'"

if "`varlist'"=="" {
	if `pA' + `pB' + `pC' + `pAB' + `pAC' + `pBC' + `pABC' == 0 {
		di as err "User must choose to execute -pvenn2- with"  ///
			_n "EITHER a -varlist- OR by specifying the sizes" ///
			_n "of the sets and their intersections."  ///
			_n "Here the varlist and population sizes are all blank."
		exit 198
	}
	local popA=`pA'
	local popB=`pB'
	local popC=`pC'
	local popAB=`pAB'
	local popAC=`pAC'
	local popBC=`pBC'
	local popABC=`pABC'  //  Added version 1.1
	local popt=`ptotal'
}
else {  //  When -varlist- is present
	if `pA' + `pB' + `pC' + `pAB' + `pAC' + `pBC' + `pABC' ~= 0 {
		di as err "User must choose to execute -pvenn2- with"  ///
			_n "EITHER a -varlist- OR by specifying the sizes" ///
			_n "of the sets and their intersections."  ///
			_n "Here the varlist and population sizes are both specified."
		exit 198
	}
	// Initialize pop sizes
	local popA= 0
	local popB= 0
	local popC= 0
	local popAB=0
	local popAC=0
	local popBC=0
	local popABC=0  //  Added version 1.1
	local popt=`ptotal'
	marksample touse
	local i 0
	foreach v of local varlist {
		local vtab`++i'= "`v'"
		qui check01 `v' if `touse'

		if "`debug'"~="" di as txt "vtab`i'= " as res "`vtab`i''"
	}
	if `i'>=1 {
		tab `vtab1' if `touse', matcell(cellvals)
		if "`debug'"~="" matrix list cellvals
		local popA=cellvals[2,1]
		local popt=`r(N)'
		local asid1 as txt " as identified by the variable " as res "`vtab1'"
	}
	if `i'>=2 {
		tab `vtab1' `vtab2' if `touse', matcell(cellvals)
		local popB=cellvals[1,2]+cellvals[2,2]
		local popAB=cellvals[2,2]
		local asid2 as txt " as identified by the variable " as res "`vtab2'"
	}
	if `i'==3 {
		tab `vtab1' `vtab3' if `touse', matcell(cellvals)
		local popC=cellvals[1,2]+cellvals[2,2]
		local popAC=cellvals[2,2]
		local asid3 as txt " as identified by the variable " as res "`vtab3'"
		tab `vtab2' `vtab3' if `touse', matcell(cellvals)
		local popBC=cellvals[2,2]
		qui count if `touse' & ///
			(`vtab1' >0 & `vtab1' <.) &  ///
			(`vtab2' >0 & `vtab2' <.) &  ///
			(`vtab3' >0 & `vtab3' <.)
		local popABC = `r(N)'
	}
}

di _n as txt "Sizes of subsets:"
di as txt "Set A: " as res "`nameA'" `asid1' as txt ": " as res "`popA'"
di as txt "Set B: " as res "`nameB'" `asid2' as txt ": " as res "`popB'"
di as txt "Set C: " as res "`nameC'" `asid3' as txt ": " as res "`popC'"
di as txt "AB Overlap: " as res "`popAB'"
di as txt "AC Overlap: " as res "`popAC'"
di as txt "BC Overlap: " as res "`popBC'"
di as txt "ABC Overlap: " as res "`popABC'"
di as txt "Total population: " as res "`popt'"
di as txt "Draw total population rectangle (0,not; 1, yes): " as res "`drawtotal'"
di as txt "Rotation angle of circle B: " as res "`angle'"
di as txt "Relative position of circle C (1,clockwise; -1,counterclockwise): " as res "`pos'"

local poptot = max(`ptotal', `popt')

if `drawtotal'==0 {
	local popt=0
}

if "`plotregion'"=="" {
	if `drawtotal'==0 {
		local plotregion plotregion(margin(vlarge))
	}
	else {
		local plotregion plotregion(margin(medlarge))
	}
}

* Set locals which determine graph colors according to active or selected scheme
tw (function y = x, range(0 10)), nodraw `scheme'

local p1 `.__SCHEME.color.p1'
local p2 `.__SCHEME.color.p2'
local p3 `.__SCHEME.color.p3'
*local p4 `.__SCHEME.color.p4'
local p4 `.__SCHEME.color.axisline'

/* Set up line color, pattern, width */
*if "`lc'"=="" local lc="cranberry emerald brown lavender"
if "`lc'"=="" local lc="`p1' `p2' `p3' `p4'"

if "`lw'"=="" local lw="0.8 0.8 0.8 0.8"
if "`lp'"=="" local lp="solid solid solid solid"

forvalues z=1/4  {
	local xopt`z'=" n(`n') `gopt' `plotregion' "
	local conlc`z': word `z' of `lc'
	if "`debug'"=="debug" di as txt "DEBUG: conlc`z' " as res "`conlc`z''"
	if "`conlc`z''"~="" local xopt`z' "`xopt`z'' lc(`conlc`z'')"
		else local xopt`z' "`xopt`z'' lc(`conlc1')"
	local conlp`z': word `z' of `lp'
	if "`conlp`z''"~="" local xopt`z' "`xopt`z'' lp(`conlp`z'')"
		else local xopt`z' "`xopt`z'' lp(`conlp1')"
	local conlw`z': word `z' of `lw'
	if "`conlw`z''"~="" local xopt`z' "`xopt`z'' lw(`conlw`z'')"
		else local xopt`z' "`xopt`z'' lw(`conlw1')"
	if "`debug'"=="debug" di as txt "DEBUG: xopt`z'= " as res "`xopt`z''"
}


/* Adjust scale of parameters */
local r1= 1                    // radius of circle 1 is 1
local r2= ((`popB'/`popA')^(0.5))* `r1'
local r3= ((`popC'/`popA')^(0.5))* `r1'
local p1= _pi * `r1' * `r1'  // area of a circle 1 is pi * r1^2
local p2= (`popB'/`popA') * `p1'
local p3= (`popC'/`popA') * `p1'
local p12= (`popAB'/`popA') * `p1'
local p13= (`popAC'/`popA') * `p1'
local p23= (`popBC'/`popA') * `p1'
local popt= (`popt'/`popA') * `p1'

if "`debug'"=="debug" {
	di as err "Line 152"
	foreach c1 of numlist 1/3 {
		di as txt "DEBUG: radius of circle `c1': r`c1'= " as res "`r`c1''"
		di as txt "DEBUG: Multiple of circle `c1''s to circle 1's area: p`c1'= " as res "`p`c1''"
		local fromc1 = `c1'
		foreach c2 of numlist `fromc1'/3 {
			di as txt "DEBUG: p`c1'`c2' = " as res  `p`c1'`c2''
		}
		local fromc1 = `c1' + 1
	}
}

if `popA'<=0 {
  local r1= 0
  local r2= 1
  local r3= ((`popC'/`popB')^(0.5))* `r2'
  local p1= 0
  local p2= _pi * `r2' * `r2'
  local p3= (`popC'/`popB') * `p2'
  local p12= (`popAB'/`popB') * `p2'
  local p13= (`popAC'/`popB') * `p2'
  local p23= (`popBC'/`popB') * `p2'
  local popt= (`popt'/`popB') * `p2'
  }
if `popA'<=0 & `popB'<=0 {
  local r1=0
  local r2=0
  local r3=1
  local p1= 0
  local p2= 0
  local p3= _pi * `r3' * `r3'
  local p12= (`popAB'/`popC') * `p3'
  local p13= (`popAC'/`popC') * `p3'
  local p23= (`popBC'/`popC') * `p3'
  local popt= (`popt'/`popC') * `p3'
}

if "`debug'"=="debug" {
	di as err "Line 190"
	foreach c1 of numlist 1/3 {
		di as txt "DEBUG: radius of circle `c1': r`c1'= " as res "`r`c1''"
		di as txt "DEBUG: Multiple of circle `c1''s to circle 1's area: p`c1'= " as res "`p`c1''"
		local fromc1 = `c1'
		foreach c2 of numlist `fromc1'/3 {
			di as txt "DEBUG: p`c1'`c2' = " as res  `p`c1'`c2''
		}
		local fromc1 = `c1' + 1
	}
}

if "`debug'"=="debug"  {
	di as txt "DEBUG: r1= " as res "`r1'"
	di as txt "DEBUG: r2= " as res "`r2'"
	di as txt "DEBUG: r3= " as res "`r3'"
	di as txt "DEBUG: p1= " as res "`p1'"
	di as txt "DEBUG: p2= " as res "`p2'"
	di as txt "DEBUG: p3= " as res "`p3'"
	di as txt "DEBUG: p12= " as res "`p12'"
	di as txt "DEBUG: p13= " as res "`p13'"
	di as txt "DEBUG: p23= " as res "`p23'"
}

/* Estimate a1 a2*/
local a1=`r1'
local a2=`r2'
local tp12=0
if `p12'>0 {
	if `r1'>= `r2' {
		while (abs(`tp12'-`p12')/`p12')> 0.01 {
			local a2= `a2'-`r2'/1000
			local a1= (`r1'^2-`r2'^2+`a2'^2)^(0.5)
			if `a2'< ((-1)*`r2') {
				di as err "cannot find appropriate location for Circle B"
				exit
			}
			else {
					local tp12 = (`r1'^2)*(acos(`a1'/`r1'))- `a1'*((`r1'^2-`a1'^2)^(0.5)) + (`r2'^2)*(acos(`a2'/`r2'))- `a2' * ((`r2'^2-`a2'^2)^(0.5))
			}
		}
	}
	else{
		while (abs(`tp12'-`p12')/`p12')> 0.01 {
			local a1= `a1'-`r1'/1000
			local a2= (`r2'^2-`r1'^2+`a1'^2)^(0.5)
			if `a1'< ((-1)*`r1') {
				di as err "cannot find appropriate location for Circle A"
				exit
			}
			else {
					local tp12 = (`r1'^2)*(acos(`a1'/`r1'))- `a1'*((`r1'^2-`a1'^2)^(0.5)) + (`r2'^2)*(acos(`a2'/`r2'))- `a2' * ((`r2'^2-`a2'^2)^(0.5))
			}
		}
	}
}

if "`debug'"=="debug"  {
	di as txt "a1= "  as res "`a1'"
	di as txt "a2= "  as res "`a2'"
	di as txt "tp12= "  as res "`tp12'"
}

if `r3'>0 & `r3'!=. {
	/* Estimate a13 a31*/
	local a13=`r1'
	local a31=`r3'
	local tp13=0
	if `p13'>0 {
		if `r1'>= `r3' {
			while (abs(`tp13'-`p13')/`p13')> 0.01 {
				local a31= `a31'-`r3'/1000
				local a13= (`r1'^2-`r3'^2+`a31'^2)^(0.5)
				if `a31'< ((-1)*`r3') {
					di as err "cannot find appropriate location for Circle C"
					exit
				}
				else {
						local tp13 = (`r1'^2)*(acos(`a13'/`r1'))- `a13'*((`r1'^2-`a13'^2)^(0.5)) + (`r3'^2)*(acos(`a31'/`r3'))- `a31' * ((`r3'^2-`a31'^2)^(0.5))
				}
			}
		}
		else{
			while (abs(`tp13'-`p13')/`p13')> 0.01 {
				local a13= `a13'-`r1'/1000
				local a31= (`r3'^2-`r1'^2+`a13'^2)^(0.5)
				if `a13'< ((-1)*`r1') {
					di as err "cannot find appropriate location for Circle C"
					exit
				}
				else {
						local tp13 = (`r1'^2)*(acos(`a13'/`r1'))- `a13'*((`r1'^2-`a13'^2)^(0.5)) + (`r3'^2)*(acos(`a31'/`r3'))- `a31' * ((`r3'^2-`a31'^2)^(0.5))
				}
			}
		}
	}

if "`debug'"=="debug"  {
	di as txt "a13= " as res "`a13'"
	di as txt "a31= " as res "`a31'"
	di as txt "tp13= " as res "`tp13'"
}

	/* Estimate a23 a32*/
	local a23=`r2'
	local a32=`r3'
	local tp23=0
	if `p23'>0 {
		if `r2'>= `r3' {
			while (abs(`tp23'-`p23')/`p23')> 0.01 {
				local a32= `a32'-`r3'/1000
				local a23= (`r2'^2-`r3'^2+`a32'^2)^(0.5)
				if `a32'< ((-1)*`r3') {
					di as err "cannot find appropriate location for Circle C"
					exit
				}
				else {
						local tp23 = (`r2'^2)*(acos(`a23'/`r2'))- `a23'*((`r2'^2-`a23'^2)^(0.5)) + (`r3'^2)*(acos(`a32'/`r3'))- `a32' * ((`r3'^2-`a32'^2)^(0.5))
				}
			}
		}
		else{
			while (abs(`tp23'-`p23')/`p23')> 0.01 {
				local a23= `a23'-`r2'/1000
				local a32= (`r3'^2-`r2'^2+`a23'^2)^(0.5)
				if `a23'< ((-1)*`r2') {
					di as err "cannot find appropriate location for Circle C"
					exit
				}
				else {
						local tp23 = (`r2'^2)*(acos(`a23'/`r2'))- `a23'*((`r2'^2-`a23'^2)^(0.5)) + (`r3'^2)*(acos(`a32'/`r3'))- `a32' * ((`r3'^2-`a32'^2)^(0.5))
				}
			}
		}
	}
	if "`debug'"=="debug"  {
		di as txt "a23= " as res "`a23'"
		di as txt "a32= " as res "`a32'"
		di as txt "tp23= " as res "`tp23'"
	}
}

/* locate the center of circles */


local cx1=0
local cy1=0
local cx2=0
local cy2=0
local cx3=0
local cy3=0
local l12=0
local l13=0
local l23=0
local angle213=0
local l12= `a1'+`a2'
local cx2= `cx1'+(`l12') * cos(_pi*`angle'/180)
local cy2= `cy1'+(`l12') * sin(_pi*`angle'/180)
if `r3'>0 & `r3'!=.{
	local l13= `a13'+`a31'
	local l23= `a23'+`a32'
	if `l12'*`l13'!=0 local angle213= acos((`l12'^2+`l13'^2-`l23'^2)/(2*`l12'*`l13'))
	local angle2=_pi*`angle'/180-`pos'*`angle213'
	local cx3=`cx1'+(`l13') * cos(`angle2')
	local cy3=`cy1'+(`l13') * sin(`angle2')
}

if "`debug'"~="" {
	di as txt "`l12' `l13' `l23'"
	di as txt "angle213= " as res "`angle213'"
	di as txt "angle2= " as res "`angle2'"
	di as txt "cx3= " as res "`cx3'"
	di as txt "cy3= " as res "`cy3'"
}

/* locate the place for label */
local lnameA=""
local lnameB=""
local lnameC=""
local xlabel1=`cx1'
local ylabel1=`cy1'
local xlabel2=`cx2'
local ylabel2=`cy2'
local xlabel3=`cx3'
local ylabel3=`cy3'
if `drawlabel'==1 {
*	if `r1'>0 local lnameA="`nameA'"
*	if `r2'>0 local lnameB="`nameB'"
*	if `r3'>0 local lnameC="`nameC'"
	local ll12=`l12'-`r2'
	local ll13=`l13'-`r3'
	local ll21=`l12'-`r1'
	local ll23=`l23'-`r3'
	local ll31=`l13'-`r1'
	local ll32=`l23'-`r2'

	local angles1=(_pi*`angle'/180-`pos'*`angle213'/2)
	local temp= min(`ll12',`ll13',0)*cos(`angles1')
	if `temp'!=0 local xlabel1=`cx1'+ `temp' +0.3 *`temp'/abs(`temp')
	local temp= min(`ll12',`ll13',0)*sin(`angles1')
	if `temp'!=0 local ylabel1=`cy1'+ `temp' +0.3 *`temp'/abs(`temp')

	local angle123=0
	if `l12'*`l23'!=0 local angle123= acos((`l12'^2+`l23'^2-`l13'^2)/(2*`l12'*`l23'))
	local angles2=(_pi*`angle'/180+`pos'*`angle123'/2)
	local temp= min(`ll21',`ll23',0)*cos(`angles2')
	if `temp'!=0 local xlabel2=`cx2'- `temp' -0.3 *`temp'/abs(`temp')
	local temp= min(`ll21',`ll23',0)*sin(`angles2')
	if `temp'!=0 local ylabel2=`cy2'- `temp' -0.3 *`temp'/abs(`temp')

	local angle132=0
	if `l13'*`l23'!=0 local angle132= acos((`l13'^2+`l23'^2-`l12'^2)/(2*`l13'*`l23'))
	local angles3=_pi*`angle'/180+`pos'*(`angle213'+ `angle132'/2)
	local temp= min(`ll31',`ll32',0)*cos(`angles3')
	if `temp'!=0 local xlabel3=`cx3'+ `temp' +0.3 *`temp'/abs(`temp')
	local temp= min(`ll31',`ll32',0)*sin(`angles3')
	if `temp'!=0 local ylabel3=`cy3'+ `temp' +0.3 *`temp'/abs(`temp')
	*di as txt "`ll12',`ll13'"
	*di as txt "angles1:" as res "`angles1'"
	*di as txt "xlabel1:" as res "`xlabel1'"
	*di as txt "`ll21',`ll23'"
	qui local angles2=180* `angles2'/_pi
	*di as txt "angles2:" as res "`angles2'"
	*di as txt "xlabel2:" as res "`xlabel2'"
	*di as txt "xlabel3:" as res "`xlabel3'"

}

/*draw the circles*/
foreach circles in 1 2 3 {
	local rangel`circles'= (`cx`circles''-`r`circles'')
	local rangeh`circles'= (`cx`circles''+`r`circles'')
	*local erangel`circles'=floor(2*`rangel`circles'') /2
	*local erangeh`circles'=ceil(2*`rangeh`circles'') /2
}
local minX=min(`rangel1',`rangel2', `rangel3')
local maxX=max(`rangeh1',`rangeh2', `rangeh3')
*di as txt "minX= " as res "`minX'"
*di as txt "maxX= " as res "`maxX'"
local Xrange=`maxX'-`minX'
foreach circles in 1 2 3 {
	local rangeb`circles'= (`cy`circles''-`r`circles'')
	local ranget`circles'= (`cy`circles''+`r`circles'')
	*local erangeb`circles'=floor(2*`rangeb`circles'')/2
	*local eranget`circles'=ceil(2*`ranget`circles'') /2
}
local minY=min(`rangeb1',`rangeb2', `rangeb3')
local maxY=max(`ranget1',`ranget2', `ranget3')
*di as txt "minY= " as res "`minY'"
*di as txt "maxY= " as res "`maxY'"
local Yrange=`maxY'-`minY'
local aspectratio=`Yrange'/`Xrange'
local rangearea=`Xrange'*`Yrange'
*di as txt "rangearea: `rangearea'"
if `rangearea'<=`popt' & `popt'>0 {
	local totalX=sqrt(`popt'/`rangearea')*`Xrange'
	local totalY=sqrt(`popt'/`rangearea')*`Yrange'
}
else {
	if `popt'>0 di as res "Margin (i.e. subjects outside the union of the sets A, B, and C) is not large enough;"  ///
		_n "total population cannot be represented by a rectangle outside the circles representing A, B & C."
	local totalX=`Xrange'
	local totalY=`Yrange'
}
local totalminX=`minX'-(`totalX'-`Xrange')/2
local totalmaxX=`maxX'+(`totalX'-`Xrange')/2
local totalminY=`minY'-(`totalY'-`Yrange')/2
local totalmaxY=`maxY'+(`totalY'-`Yrange')/2
local minX=floor(2*`totalminX')/2
local maxX=ceil(2*`totalmaxX')/2
local minY=floor(2*`totalminY')/2
local maxY=ceil(2*`totalmaxY')/2
local ylabel4=`totalmaxY'-0.1
local xlabel4=`totalminX'+0.1

if "`debug'"~= ""  {
	di as txt "aspectratio= " as res "`aspectratio'"
	if "`debug'"=="debug" {
		foreach r of numlist 1/3 {
			foreach lh in l h {
				di as txt "DEBUG: The value of range`lh'`r' is: " as res  `range`lh'`r''
			}
		}
	}
}

/*	Add variable labels instead of variable names

if "`nolabel'"=="nolabel" {

	di as txt "label is `label'"
	di as txt "lnameA is: `lnameA'"

	local lnameA : variable label `vtab1'
	local lnameB : variable label `vtab2'
	local lnameC : variable label `vtab3'

}

*	Add counts instead of variable names
local lnameA : display %6.0f `popA'
local lnameB : display %6.0f `popB'
local lnameC : display %6.0f `popC'

*/

if "`label'"=="nolabel" {
	local lnameA
	local lnameB
	local lnameC
}
else {
	if "`underscore'"=="" {
		local lnameA = subinstr("`nameA'","_"," ",.)
		local lnameB = subinstr("`nameB'","_"," ",.)
		local lnameC = subinstr("`nameC'","_"," ",.)
	}
	else {
		local lnameA "`nameA'"
		local lnameB "`nameB'"
		local lnameC "`nameC'"
	}
}

*	Next step, get rid of the variable labels in the middle of the circles...

if `rangearea'<=`popt' {
	if `r3'>0 & `r3'!=.{
		if `"`legend'"'== ""  {

			//  New in version 1.5 25May2021
			legdefault `nameA' `nameB' `nameC' ,   ///
				popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
				popA(`popA') popB(`popB') popC(`popC')  ///
				popAB(`popAB') popAC(`popAC') popBC(`popBC')  popABC(`popABC') `underscore'
			local legend `r(legdefault)'

		}
		twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off)  ///
				xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio')  ///
			|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ///
			|| function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' `"`lnameB'"') `xopt2'  ///
			|| function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2'  ///
			|| function y = sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') text(`ylabel3' `xlabel3' `"`lnameC'"') `xopt3'  ///
			|| function y = -sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') `xopt3'  ///
			|| function y = `totalminX',range(`totalminY' `totalmaxY') horizontal xtitle("") ytitle("")  `xopt4'  ///
			|| function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal  `xopt4'  ///
			|| function y = `totalminY',range(`totalminX' `totalmaxX')  `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) ///
			|| function y = `totalmaxY',range(`totalminX' `totalmaxX')  `xopt4' ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'
	}
	else {
		if `r2'>0 & `r2'!=.{
			if `"`legend'"'== ""  {

				//  New in version 1.5 25May2021
				legdefault `nameA' `nameB' ,   ///
					popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
					popA(`popA') popB(`popB')  ///
					popAB(`popAB') `underscore'
				local legend `r(legdefault)'

			}
			twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off)  ///
					xtitle("") ytitle("") `xopt1'  aspectratio(`aspectratio') ///
				|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ///
				|| function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' `"`lnameB'"') `xopt2'   ///
				|| function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2'   ///
				|| function y = `totalminX',range(`totalminY' `totalmaxY') xtitle("") ytitle("") horizontal `xopt4'  ///
				|| function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal `xopt4'  ///
				|| function y = `totalminY',range(`totalminX' `totalmaxX')  `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) ///
				|| function y = `totalmaxY',range(`totalminX' `totalmaxX')  `xopt4' ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'
		}
		else {
			if `"`legend'"'== ""  {

				//  New in version 1.5 25May2021
				legdefault `nameA' ,   ///
					popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
					popA(`popA') `underscore'
				local legend `r(legdefault)'

			}
			twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off) ///
					xtitle("") ytitle("") `xopt1'  aspectratio(`aspectratio') ///
				|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ///
				|| function y = `totalminX',range(`totalminY' `totalmaxY') xtitle("") ytitle("") horizontal `xopt4'  ///
				|| function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal `xopt4'  ///
				|| function y = `totalminY',range(`totalminX' `totalmaxX')  `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) ///
				|| function y = `totalmaxY',range(`totalminX' `totalmaxX')  `xopt4' ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'
			}
	}

	exit
}

if `r3'>0 & `r3'!=.{
	if "`debug'"=="debug" {
		foreach r of numlist 1/3 {
			foreach lh in l h {
				di as txt "DEBUG: The value of range`lh'`r' is: " as res  `range`lh'`r''
			}
		}
	}
	if `"`legend'"'== ""  {

		//  New in version 1.5 25May2021
		legdefault `nameA' `nameB' `nameC' ,   ///
			popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
			popA(`popA') popB(`popB') popC(`popC')  ///
			popAB(`popAB') popAC(`popAC') popBC(`popBC') popABC(`popABC') `underscore'
		local legend `r(legdefault)'

	}
	twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off)  ///
			xtitle("") ytitle("")  `xopt1' aspectratio(`aspectratio') `legend_was_here'  ///
		|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ///
		|| function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' `"`lnameB'"') `xopt2'   ///
		|| function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2'  ///
		|| function y = sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') text(`ylabel3' `xlabel3' `"`lnameC'"') `xopt3'  ///
		|| function y = -sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') `xopt3' ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'

}
else {
	if `r2'>0 & `r2'!=.{
		if `"`legend'"'== ""  {
*			local legend legend(order(1 "`nameA'" 3 "`nameB'"))
*			legdefault `vtab1' `vtab2' ,   ///

			//  New in version 1.5 25May2021
			legdefault `nameA' `nameB' ,   ///
				popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
				popA(`popA') popB(`popB')  ///
				popAB(`popAB')  `underscore'
			local legend `r(legdefault)'

		}
		twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off)  ///
				xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio')  ///
			|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ///
			|| function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' `"`lnameB'"') `xopt2'  ///
			|| function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2' ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'
	}
	else {
		if `"`legend'"'== ""  {

			//  New in version 1.5 25May2021
			legdefault `nameA' ,   ///
				popt(`popt') poptot(`poptot') r2(`r2') r3(`r3') `notexoffset'   ///
				popA(`popA') `underscore'
			local legend `r(legdefault)'

		}
		twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' `"`lnameA'"') xlabel(none) ylabel(none) xscale(off) yscale(off)  ///
				xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio')  ///
			|| function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1')  `xopt1'  ,  `text' `legend' `saving' `name' `note' `title' `subtitle' `schemeopt'
	}
}

end

prog define legdefault , rclass
	syntax [namelist(max=3)], popt(real) poptot(real) r2(real) r3(real)  ///
		[popA(integer -9) popB(integer -9) popC(integer -9)  ///
		popAB(integer -9) popAC(integer -9) popBC(integer -9) popABC(integer -9) notexoffset(real 0) UNDERscore]

*	These legends could be executed from the command line syntax such as:
*	local order order(1 "\`vtab1' = `popA')" 3 "\`vtab2' = `popB')" 5 "\`vtab3' = `popC')" )

	if "`namelist'"~=""  tokenize `namelist'

	if `poptot'> 0 {
	    local totnote `""Total = \`poptot'""'
	}

	//  Convert set sizes to Stata's fc format
	foreach p in popA popB popC popAB popAC popBC popABC poptot {
	    if ``p'' > 0 {
		    local `p' : display %25.0fc ``p''
			local `p' = trim("``p''")
		}
		else {
		    local `p' 0
		}
	}

	if `notexoffset'>0 {
		local xoffset xoffset(`notexoffset')
	}

	if `r3'>0 & `r3'!=. {
		if "`namelist'"~=""  {
			local legtitle title("Set A:  `1'" "Set B:  `2'" "Set C:  `3'", size(*.8))
			if "`underscore'"=="" local legtitle = subinstr(`"`legtitle'"',"_"," ",.)
		}
		local order order(1 "A = `popA'" 3 "B = `popB'" 5 "C = `popC'" )
		local legnote note(`totnote' "A {&intersect} B = `popAB'" "A {&intersect} C = `popAC'" "B {&intersect} C = `popBC'" "A {&intersect} B {&intersect} C = \`popABC'", size(*1.3) justification(right) `xoffset')
		local legend legend(`legtitle' `order' `legnote' pos(3) col(1) region(lstyle(none)) )
    }

	else {
		if `r2'>0 & `r2'!=. {
			if "`namelist'"~=""  {
				local legtitle title("Set A:  `1'" "Set B:  `2'", size(*.8))
				if "`underscore'"=="" local legtitle =  subinstr(`"`legtitle'"',"_"," ",.)
			}
			local order order(1 "A = `popA'" 3 "B = `popB'" )
			local legnote note(`totnote' "A {&intersect} B = `popAB'", size(*1.3) justification(right) `xoffset')
			local legend legend(`legtitle' `order' `legnote' pos(3) col(1) region(lstyle(none)) )
		}
		else {
			if "`namelist'"~=""  {
				local legtitle title("Set A:  `1'", size(*.8))
				if "`underscore'"=="" local legtitle =  subinstr(`"`legtitle'"',"_"," ",.)
			}
			local order order(1 "A = `popA'" )
			local legnote note(`totnote', size(*1.3) justification(right) )
			local legend legend(`legtitle' `order' `legnote' pos(3) col(1) region(lstyle(none)) `xoffset')
		}
	}

	return local legdefault `legend'

end /* End of subroutine legdefault */

prog define check01

	syntax varname if/

	confirm numeric variable `varlist'

	count if `if' & `varlist' < 0
		local lt0 `r(N)'
	count if `if' & `varlist' == 0
		local eq0 `r(N)'
	count if `if' & `varlist' > 0 & `varlist' < 1
		local btwn01 `r(N)'
	count if `if' & `varlist' == 1
		local eq1 `r(N)'
	count if `if' & `varlist' > 1 & `varlist' <.
		local gt1 `r(N)'
	count if `if' & `varlist' >=.
		local msng `r(N)'

	local badvalues = `lt0' + `btwn01' + `gt1'
	if `badvalues' > 0 {
		di as err "When generating a Venn diagram from a dataset, each of the variables in the -varlist- must be binary."  ///
			_n as err "All values must be either 0, 1 or missing.  Here the variable -`varlist'- has `badvalues' which are neither 0, 1 nor missing."
		exit 198
	}
	if `msng' > 0 {
		di as txt "The variable " as res "`varlist'" as txt " has " as res "`msng'" as txt " missing values"
	}


end  /* End of subroutine check01 */


/*

Future ideas by W. Gong
1. treat problem when a population is involved in the overlap part of
   the other two populations. eg. pvenn pA(45) pB(25) pC(50) pAB(18) pAC(24) pBC(3) ang(120)
2. fill color in circles.

Future ideas by Mead Over
3.	Look for inspiration to the article:
	eulerAPE: Drawing Area-Proportional 3-Venn Diagrams Using Ellipses
	by Luana Micallef,  Peter Rodgers
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0101717

	and from:
Alexandersson, A. 1998. gr32: Confidence ellipses. Stata Technical Bulletin 46: 10-13. In Stata Technical Bulletin Reprints, vol. 8,
        54-57. College Station, TX:  Stata Press.
------. 2004. Graphing confidence ellipses: An update of ellip for Stata 8. Stata Journal 4(3): 242-256.

4.  Fix broken example from help file:
		pvenn2, pA(3) pB(2) pC(3) pAB(1) pAC(2) pBC(2) ptotal(5)
		Is this the same as Gong's issue #1 above?
5.	Add the number/percentage of observations in each intersection (as in Lauritsen's -venndiag-)
6.	Add the [if] and [in] options
7.	Figure out how to close the gaps in the circle when using 2000 rather than 20,000 points
		.i.e. Perhaps allocate a constant number of points per degree around the circles?
*/

*	Version 1.0 26 Oct 2011 by Wenfeng Gong; Jan Osterman (Duke Center for Health Policy)
*				Email: gongwenf@gmail.com
*	Version 1.1 21/6/2015 by Mead Over, Center for Global Development
*	Version 1.2 4Dec2017  by Mead Over, Center for Global Development
*		Cosmetic changes to output to the results window
*		Fix variable label placement on the graph
*		When user specifies the -legend(order)- option, s/he must use
*		order(1 "hkhk" 3 "hkhk" 5 "hkhk" 7 "Total population"),
*		not order(1 "hkhk" 2 "hkhk" 3 "hkhk" 4 "Total population")
*		The -addtext- option did not work.  Replace with -text- option.
*		For text option, point (0,0) is center of rectangle,
*		(1,1) is northeast corner and (-1,-1) is southwest corner.
*	Version 1.21 18Dec2018 Change default number of points from 20,000 to 2,000
*	Version 1.4 23Feb2021
*		Some debug added. But the broken example in the original help file for
*		-pvenn.ado- remains broken.  Try command:  pvenn2, pA(3) pB(2) pC(3) pAB(1) pAC(2) pBC(2) ptotal(5)
*		Better algoritm which uses ellipses and a hill-climbing solution algorithm is that in:
*		Micallef L, Rodgers P (2014) eulerAPE: Drawing Area-Proportional 3-Venn Diagrams Using Ellipses. PLoS ONE 9(7): e101717. https://doi.org/10.1371/journal.pone.0101717
*	Version 1.5 25May2021 Add the subroutine -legdefault- to provide a default legend
*		Suppress the display of useless x- and y-axes. Allow default colors to be set by the scheme.
*		Add undocumented option -notexoffset(#)-. It is only relevant with the default legend.
*		Its effect is to shift the legend's "note" to the right or left,
*		so that it is better aligned with the keys above it.
*		Retain separate -title_options- so as to preserve their ability
*		to incorporate quotes, and thus multiple lines.
*	Version 1.6 22May2022 Allow the [if] syntax.
*		Set scheme colors correctly on first execution of -pvenn2-.
*		Allow option -plabel()- to provide names for sets A, B and C
*		The option NOTEXTOffset() remains undocumented.
*		Exit with error message if the optional -varlist- is empty and the optional
*		"immediate" options (-pA()-, -pB()-, etc.) for specifying the population sizes are also all empty.
*		Add undocumented option -Underscore- which retains underscores rather than replacing them with spaces in legend