program prtab, rclass
    version 11

	syntax varlist (numeric min=2 max=2) [if] [in] , [prec_at(real .2) COMpare(varname) INTERpolate rank THREShold fscore NOGraph NOREFline *]

	_get_gropts, graphopts(`options')
	local twowayopts `"`s(graphopts)'"'

	gettoken y_name classifier_name : varlist

 /* Call preserve because we will need to create new variables */
	preserve

	tempvar y classifier
	qui gen `y' = `y_name' `if' `in'
	qui gen `classifier' = `classifier_name' `if' `in'
	qui drop if missing(`y') | missing(`classifier')
	gsort -`classifier'

 /* Checks on inputs */
	cap assert `y'==0 | `y'==1
	if  _rc~=0 {
	  noi di in red "true status of variable `y_name' must be 0 or 1"
	  exit 198
	}
	if "`rank'"!="" & "`threshold'"!="" {
	  noi di in red "you cannot specify both rank and threshold"
	  exit 198
	}
	if "`compare'"!="" & "`threshold'"!="" {
	  noi di in red "you cannot use compare with threshold"
	  exit 198
	}
    if "`rank'"!="" & "`fscore'"!="" {
	  noi di in red "you cannot specify both rank and fscore"
	  exit 198
	}
	if "`rank'"!="" & "`interpolate'"!="" {
	  noi di in red "you cannot specify both rank and interpolate"
	  exit 198
	}
    if "`prec_at'" > "1" & "`rank'"=="" {
	  noi di in red "you have provided prec_at with a value of recall that is greater than 1"
	  exit 198
	}
 /* Finished checks on inputs */



  if "`compare'" == "" {
    qui sum `y'
	local rnd_rate = r(mean)
	local obs = r(N)
	local total_pos = `rnd_rate'*`obs'

    if "`rank'"=="" {
	  tempvar running_total prec recall ranking fscore_val
	  gen `running_total' = sum(`y')
	  gen `ranking' = _n
	  gen `prec' = `running_total'/`ranking'
      gen `recall' = `running_total'/`total_pos'
      gen `fscore_val' = 2*`prec'*`recall'/(`prec'+`recall')

	  tempvar unique
	  gsort `classifier' -`recall'
	  by `classifier': gen `unique'= _n==1
	  qui sum `unique'
	  local unique_obs = r(N)*r(mean)
	  gsort -`classifier'
      qui drop if `unique'==0

      tempvar obs_number
	  gen `obs_number' = _n

 /* Display table */

      di as txt _n "{col 12}Number of observations       =  " as result    `obs'
      di as txt "{col 12}Unique values of classifier  =  "  as result   `unique_obs'
      di as txt "{col 12}Number of positive cases     =  " as result    `total_pos'
      di as txt  "{col 12}Portion of positive cases    =" as result  %8.04f `rnd_rate'

	  if `unique_obs' >= 3 {
       tempvar firsttime1 firsttime2 firsttime3
       egen `firsttime1' = min(cond(`recall' >= `prec_at'-.1, `obs_number', .))
       egen `firsttime2' = min(cond(`recall' >= `prec_at', `obs_number', .))
       if `firsttime1'==`firsttime2'{
	     qui replace `firsttime2'=`firsttime1' + 1
	   }
	   egen `firsttime3' = min(cond(`recall' >= `prec_at'+.1, `obs_number', .))
       if `firsttime3'<=`firsttime2'{
	     qui replace `firsttime3'=`firsttime2' + 1
	   }
	   di as txt _n "{hline 50}
       di as txt "{col 5} Recall =" %8.04f `recall'[`firsttime1']  _col(25)  %8.04f `recall'[`firsttime2'] _col(37)  %8.04f `recall'[`firsttime3']
       di as txt "{hline 50}
       di as res "{col 2}Precision" as result  _col(14) %8.04f `prec'[`firsttime1'] _col(25) %8.04f `prec'[`firsttime2'] _col(37) %8.04f `prec'[`firsttime3']
       di as txt "{hline 50}
	  }

 /* Interpolate */
	  if "`interpolate'" !=""{
	    local new = _N + 1
		qui set obs `new'
		qui replace `recall'=0 if missing(`recall')
		qui replace `prec'=0 if missing(`prec')
	    gsort -`recall'
      	tempvar inter_prec
        qui gen `inter_prec'=`prec'[1] in 1
		qui replace `inter_prec'=max(`prec', `inter_prec'[_n-1]) if missing(`inter_prec')
		qui replace `prec'=`inter_prec'
		sort `recall'
	  }
 /* Display graph */
      if "`nograph'"=="" {
	    if "`threshold'"=="" & "`fscore'"==""{
          if "`norefline'"==""{
            graph twoway line `prec' `recall' , yline(`rnd_rate') ytitle("Precision") xtitle("Recall") `twowayopts'
	      }
	      else{
	        graph twoway line `prec' `recall', ytitle("Precision") xtitle("Recall") `twowayopts'
	      }
        }
	    else if "`fscore'"==""{
          if "`norefline'"==""{
		    graph twoway line `prec' `classifier', xsc(reverse) yline(`rnd_rate') ytitle("Precision") xtitle("Threshold") `twowayopts'
		  }
		  else{
		    graph twoway line `prec' `classifier', xsc(reverse) ytitle("Precision") xtitle("Threshold") `twowayopts'
		  }
	    }
		else{
		  if "`threshold'"=="" {
            graph twoway line `fscore_val' `recall', ytitle("F-score") xtitle("Recall") `twowayopts'
		  }
		  else{
		    graph twoway line `fscore_val' `classifier', xsc(reverse) ytitle("F-score") xtitle("Threshold") `twowayopts'
		  }
		}
	  }
 /* Calculate and display AUC */
	qui integ `prec' `recall'
	local auc = r(integral)
	di as txt _n"{col 2}Area under precision-recall curve:"%8.04f `auc'

 /* Return values */
	return scalar AUC = `auc'
	return scalar N = `obs'
	return scalar unique_val = `unique_obs'
	return scalar pos_cases = `total_pos'
    return scalar prct_pos = `rnd_rate'
	}
	else {
	  tempvar running_total prec ranking avg_pos
	  gen `ranking' = _n
      bysort `classifier': egen `avg_pos'= mean(`y')
	  gsort -`classifier'
	  gen `running_total' = sum(`avg_pos')
	  gen `prec' = `running_total'/`ranking'

	  tempvar unique
	  bysort `classifier': gen `unique'= _n==1
	  qui sum `unique'
	  local unique_obs = r(N)*r(mean)
	  gsort -`classifier'

      tempvar obs_number
	  gen `obs_number' = _n

 /* Display table */
      di as txt _n "{col 12}Number of observations       =  " as result    `obs'
      di as txt "{col 12}Unique values of classifier  =  "  as result   `unique_obs'
      di as txt "{col 12}Number of positive cases     =  " as result    `total_pos'
      di as txt  "{col 12}Portion of positive cases    =" as result  %8.04f `rnd_rate'
       local top_1 = round(.01*`obs',1)
       local top_5 = round(.05*`obs',1)
	   local top_10 = round(.10*`obs',1)

	   if "`prec_at'"==".2" {
	     di as txt _n"{hline 50}
         di as txt  "{col 15} Top 1% {col 28} Top 5% {col 40} Top 10%"
         di as txt "{hline 50}
         di as res "{col 2}Precision" _col(14) as result  %8.04f `prec'[`top_1'] _col(27) %8.04f `prec'[`top_5'] _col(39) %8.04f `prec'[`top_10']
         di as txt "{hline 50}
       }
	   else if floor(`prec_at') <= `obs'{
	   	 di as txt _n"{hline 50}
         di as txt  "{col 23} At rank " floor(`prec_at')
         di as txt "{hline 50}
         di as res "{col 12}Precision" _col(24)as result   %8.04f `prec'[floor(`prec_at')]
         di as txt "{hline 50}
	   }
	   else{
	   	 di as txt _n"{hline 50}
         di as txt  "{col 23} At rank " floor(`obs')
         di as txt "{hline 50}
         di as res "{col 12}Precision" _col(24)as result   %8.04f `prec'[floor(`obs')]
         di as txt "{hline 50}
		  di as res _n"Note: The provided value of prec_at is greater than the"
	      di as res "number of observations."
	   }

 /* Display graph */
      if "`nograph'"=="" {
        if "`norefline'"==""{
          graph twoway line `prec' `ranking', yline(`rnd_rate') ytitle("Precision") xtitle("Rank") `twowayopts'
	    }
	    else{
	      graph twoway line `prec' `ranking', ytitle("Precision") xtitle("Rank") `twowayopts'
	    }
	    if `unique_obs'<`obs' {
	      di as res _n"Note: There are multiple observations with the same classifer"
	      di as res "value. Average precision over possible rankings was used."
	    }

      }

 /* Return values */
	return scalar N = `obs'
	return scalar unique_val = `unique_obs'
	return scalar pos_cases = `total_pos'
    return scalar prct_pos = `rnd_rate'
	}

  }
 /* Now handle comparisons of two classifiers */
  else {
     tempvar classifier2
	 qui gen `classifier2' = `compare' `if' `in'
	 qui drop if missing(`classifier2')
	 qui sum `y'
	 local rnd_rate = r(mean)
	 local obs = r(N)
	 local total_pos = `rnd_rate'*`obs'

	 if "`rank'"=="" {

	 /*  Main classifier */
	  tempvar running_total prec recall ranking fscore_val
	  qui gen `running_total' = sum(`y')
	  qui gen `ranking' = _n
	  qui gen `prec' = `running_total'/`ranking'
      qui gen `recall' = `running_total'/`total_pos'
      qui gen `fscore_val' = 2*`prec'*`recall'/(`prec'+`recall')

	  tempvar unique
	  gsort `classifier' -`recall'
	  by `classifier': gen `unique'= _n==1

  /*  Comparison classifier */
	  tempvar running_total2 prec2 recall2 ranking2 fscore_val2
      gsort -`classifier2'
	  qui gen `running_total2' = sum(`y')
	  qui gen `ranking2' = _n
	  qui gen `prec2' = `running_total2'/`ranking2'
      qui gen `recall2' = `running_total2'/`total_pos'
      qui gen `fscore_val2' = 2*`prec2'*`recall2'/(`prec2'+`recall2')

	  tempvar unique2
	  gsort `classifier2' -`recall2'
	  by `classifier2': gen `unique2'= _n==1

 /* Interpolate */
	  if "`interpolate'" !=""{
	    tempvar temp_prec new_prec
	    qui gen `temp_prec'=`prec' if `unique'
	    qui bysort `classifier': egen `new_prec'=max(`temp_prec')

	    local new = _N + 1
		qui set obs `new'
		qui replace `recall'=0 if missing(`recall')
		qui replace `new_prec'=0 if missing(`prec')
		qui replace `unique'=1 if missing(`unique')
		qui sum `classifier'
		qui replace `classifier'=r(max)+1 if missing(`classifier')
		gsort `classifier' -`recall'
      	tempvar inter_prec
        qui gen `inter_prec'=`new_prec'[1] in 1
		qui replace `inter_prec'=max(`new_prec', `inter_prec'[_n-1]) if  missing(`inter_prec')
		qui replace `prec'=`inter_prec'

	    tempvar temp_prec2 new_prec2
	    qui gen `temp_prec2'=`prec2' if `unique2'
	    qui bysort `classifier2': egen `new_prec2'=max(`temp_prec2')

		qui replace `recall2'=0 if missing(`recall2')
		qui replace `prec2'=0 if missing(`prec2')
		qui sum `classifier2'
		qui replace `classifier2'=r(max)+1 if missing(`classifier2')
        gsort -`recall2' `classifier2'
      	tempvar inter_prec2
        qui gen `inter_prec2'=`new_prec2'[1] in 1
		qui replace `inter_prec2'=max(`new_prec2', `inter_prec2'[_n-1]) if  missing(`inter_prec2')
		qui replace `prec2'=`inter_prec2'
	  }
 /* Display graph */
      if "`nograph'"=="" {
	  	gsort -`classifier'
	    tempvar neg_c2
	    gen `neg_c2' = -`classifier2'

	    if "`fscore'"==""{
          if "`norefline'"==""{
		    graph twoway line `prec' `recall' if `unique', yline(`rnd_rate') ytitle("Precision") xtitle("Recall") `twowayopts' legend(order(1 "`classifier_name'" 2 "`compare'"))|| line `prec2' `recall2' if `unique2', sort(`neg_c2')
		  }
	      else{
	        graph twoway line `prec' `recall' if `unique', ytitle("Precision") xtitle("Recall") `twowayopts' legend(order(1 "`classifier_name'" 2 "`compare'"))|| line `prec2' `recall2' if `unique2', sort(`neg_c2')
	      }
        }
	    else {
		  graph twoway line `fscore_val' `recall' if `unique', sort xsc(reverse) ytitle("F-score") xtitle("Recall") `twowayopts' legend(order(1 "`classifier_name'" 2 "`compare'")) || line `fscore_val2' `recall2' if `unique2', xsc(reverse) sort(`neg_c2')
		}
	  }
 /* Calculate AUC */
	qui integ `prec' `recall' if `unique'
	local auc = r(integral)
	qui integ `prec2' `recall2' if `unique2'
	local auc2 = r(integral)

 /* Display table */
    local shortname1 = substr("`classifier_name'", 1, 17)
	local shortname2 = substr("`compare'", 1, 17)
    di as txt _n "{col 12}Number of observations       =  " as result   `obs'
    di as txt "{col 12}Number of positive cases     =  "   as result  `total_pos'
    di as txt  "{col 12}Portion of positive cases    ="  as result %8.04f `rnd_rate'
	di as txt _n"{hline 50}
	di as txt "{col 14}`shortname1'{col 34}`shortname2'"
	di as txt "{hline 50}
	di as txt "{col 2}PR AUC{col 12}"as result  %8.04f `auc' "{col 32}" %8.04f `auc2'
	di as txt "{hline 50}

 /* Return values */
	return scalar AUC = `auc'
	return scalar AUC2 = `auc2'
	return scalar N = `obs'
	return scalar pos_cases = `total_pos'
    return scalar prct_pos = `rnd_rate'
	 }
	 else{
	  tempvar running_total prec ranking avg_pos
	  gen `ranking' = _n
      bysort `classifier': egen `avg_pos'= mean(`y')
	  gsort -`classifier'
	  gen `running_total' = sum(`avg_pos')
	  gen `prec' = `running_total'/`ranking'

	  tempvar unique
	  bysort `classifier': gen `unique'= _n==1
	  qui sum `unique'
	  local unique_obs = r(N)*r(mean)

	  tempvar running_total2 prec2 ranking2 avg_pos2
	  gsort -`classifier2'
	  gen `ranking2' = _n
      bysort `classifier2': egen `avg_pos2'= mean(`y')
	  gsort -`classifier2'
	  gen `running_total2' = sum(`avg_pos2')
	  gen `prec2' = `running_total2'/`ranking2'

	  tempvar unique2
	  bysort `classifier2': gen `unique2'= _n==1
	  qui sum `unique2'
	  local unique_obs2 = r(N)*r(mean)

      tempvar obs_number
	  gen `obs_number' = _n

 /* Start table */
	  di as txt _n "{col 12}Number of observations       =  " as result   `obs'
      di as txt "{col 12}Number of positive cases     =  "   as result  `total_pos'
      di as txt  "{col 12}Portion of positive cases    ="  as result %8.04f `rnd_rate'

 /* Display graph */
      if "`nograph'"=="" {
	  	gsort -`classifier'
        if "`norefline'"==""{
          graph twoway line `prec' `ranking', sort yline(`rnd_rate') ytitle("Precision") xtitle("Rank") legend(order(1 "`classifier_name'" 2 "`compare'")) `twowayopts' || line `prec2' `ranking2', sort(`ranking2')
	    }
	    else{
	      graph twoway line `prec' `ranking', sort ytitle("Precision") xtitle("Rank") legend(order(1 "`classifier_name'" 2 "`compare'")) `twowayopts' || line `prec2' `ranking2', sort(`ranking2')
	    }
      }

 /* Calculate AUC to return and display */
    drop `prec' `prec2' `ranking' `ranking2' `running_total' `running_total2'

	/*  Main classifier */
	  tempvar running_total prec recall ranking
	  qui gen `running_total' = sum(`y')
	  qui gen `ranking' = _n
	  qui gen `prec' = `running_total'/`ranking'
      qui gen `recall' = `running_total'/`total_pos'

	  tempvar unique
	  gsort `classifier' -`recall'
	  by `classifier': gen `unique'= _n==1

     /*  Comparison classifier */
	  tempvar running_total2 prec2 recall2 ranking2
      gsort -`classifier2'
	  qui gen `running_total2' = sum(`y')
	  qui gen `ranking2' = _n
	  qui gen `prec2' = `running_total2'/`ranking2'
      qui gen `recall2' = `running_total2'/`total_pos'

     /* Interpolate */
	  if "`interpolate'" !=""{
  	    sort `classifier'
      	tempvar inter_prec
        qui gen `inter_prec'=`prec'[1] in 1
		qui replace `inter_prec'=max(`prec', `inter_prec'[_n-1]) if missing(`inter_prec')
		qui replace `prec'=`inter_prec'

	    sort `classifier2'
      	tempvar inter_prec2
        qui gen `inter_prec2'=`prec2'[1] in 1
		qui replace `inter_prec2'=max(`prec2', `inter_prec2'[_n-1]) if missing(`inter_prec2')
		qui replace `prec2'=`inter_prec2'
	  }
     /* Calculate AUC */
	qui integ `prec' `recall' if `unique'
	local auc = r(integral)
	qui integ `prec2' `recall2' if `unique2'
	local auc2 = r(integral)

     /* Display table */
    local shortname1 = substr("`classifier_name'", 1, 17)
	local shortname2 = substr("`compare'", 1, 17)
    di as txt _n"{hline 50}
	di as txt "{col 14}`shortname1'{col 34}`shortname2'"
	di as txt "{hline 50}
	di as txt "{col 2}PR AUC{col 12}"as result  %8.04f `auc' "{col 32}" %8.04f `auc2'
	di as txt "{hline 50}

	 if "`nograph'"=="" & (`unique_obs2'<`obs' | `unique_obs'<`obs') {
	      di as res _n"Note: There are multiple observations with the same classifer"
	      di as res "value. Average precision over possible rankings was used."
	    }



 /* Return values */
	return scalar AUC = `auc'
	return scalar AUC2 = `auc2'
	return scalar N = `obs'
	return scalar pos_cases = `total_pos'
    return scalar prct_pos = `rnd_rate'
	}




 }



 /* Call restore to remove created variables */
	restore

	end