*! version 0.2 2023-01-23 Niels Henrik Bruun
*! Support: Niels Henrik Bruun, niels.henrik.bruun@gmail.com
*! 2023-02-21 v0.2  Adding graph options are now working (Thank you to Eric Melse)
*! 2023-02-21 v0.2  functionality for weights added (Thank you to Eric Melse)
* 2023-02-10 v0.1   created


* 2002 Brandon - Two Versions of the Contrasting-Groups Standard-Setting Method: A Review
/*
1 /(sqrt(2*pi)*s1)*exp(-0.5*(x-m1)^2 / s1^2) = 1 /(sqrt(2*pi)*s2)*exp(-0.5*(x-m2)^2 / s2^2)

1 /s1*exp(-0.5*(x-m1)^2 / s1^2) = 1 /s2*exp(-0.5*(x-m2)^2 / s2^2)

-log(s1) + -0.5*(x-m1)^2 / s1^2 = -log(s2) + -0.5*(x-m2)^2 / s2^2

0 = log(s1) + 0.5*(x-m1)^2 / s1^2 - log(s2) - 0.5*(x-m2)^2 / s2^2

0 = 2*log(s1/s2) + (x-m1)^2 / s1^2 - (x-m2)^2 / s2^2

0 = 2*log(s1/s2)*s1^2*s2^2 + s2^2*(x-m1)^2 - s1^2*(x-m2)^2

0 = (s2^2 - s1^2)*x^2 - 2*(m1*s2^2-m2*s1^2)*x + 2*log(s1/s2)*s1^2*s2^2 + m1^2*s2^2-m2^2*s1^2


same sigmas:
0 = s^2*(x-m1)^2 - s^2*(x-m2)^2

0 = (2*s^2*(m1-m2))*x - (m1^2-m2^2)*s^2

x = (m1^2-m2^2) / (m1-m2) / 2 = (m1+m2) / 2
*/
capture program drop cgssmi
program define cgssmi, rclass

  version 12
  
  syntax , MSDrowmatrix(string) /*
    */[ /*
      */Refname(string) /*
      */Compname(string) /*
      */Header(string) /*
      */CPFormat(string) /*
      */PFormat(string) /*
      */Graph /*
      */ * /*
    */]

  _get_gropts , graphopts(`options')

  if "`cpformat'" != "" {
    capture confirm format `cpformat'
    if _rc mata: _error("Format for Cut-point is not a format")
  }
  else local cpformat %6.0f
  if "`pformat'" != "" {
    capture confirm format `pformat'
    if _rc mata: _error("Format for p-values is not a format")
  }
  else local pformat %5.1f
  
  capture matrix _msd = `msdrowmatrix' //First row is for references/experts
  if !_rc {
    local m1 = `=_msd[1,1]'
    local sd1 = `=_msd[1,2]'
    local m2 = `=_msd[2,1]' 
    local sd2 = `=_msd[2,2]'
  }
  else mata _error("Option msdrowmatrix is not a matrix")
  
  local lb = min(`m1'-3*`sd1', `m2'-3*`sd2')
  local ub = max(`m1'+3*`sd1', `m2'+3*`sd2')    
  
  if "`refname'" == "" local refname "Row 1 (ref.)"
  if "`compname'" == "" local compname "Row 2 (comp.)"
  if "`header'" != "" local header `"- "`header'""'

  *mata: ncut(`m1', `sd1', `m2', `sd2')
  mata: st_local("ncut", strofreal(ncut(`m1', `sd1', `m2', `sd2')))
  if "`ncut'" != "." {
    scalar fn = normal(`=(`ncut' - `m1') / `sd1'')
    scalar fp = 1 - normal(`=(`ncut' - `m2') / `sd2'')
    if `m1' < `m2' {
      scalar fp = normal(`=(`ncut' - `m2') / `sd2'')
      scalar fn = 1 - normal(`=(`ncut' - `m1') / `sd1'')
    }
    local xline xline(`ncut')
    local pfcut `""Pass-Fail cut-off = `=string(`ncut', "`cpformat'")'""'
    local fp `""Theoretical false positive = `=string(100*fp, "`pformat'")'%""'
    local fn `""Theoretical false negative  = `=string(100*fn, "`pformat'")'%""'
  }
  else {
      scalar fp = .
      scalar fn = .    
  }

  if "`graph'`options'" != "" | "`graph'" != "" {
    local gr_txt twoway (function y=normalden(x,`m1',`sd1'), range(`lb' `ub')) ///
      (function y=normalden(x,`m2',`sd2'), range(`lb' `ub')), `xline' ///
      subtitle(`pfcut' `fp' `fn', ring(0) position(11) size(small)) ///
      legend(order(`header' 1 "`refname'" 2 "`compname'")) ///
      yscale(off) `s(graphopts)'
    `gr_txt'
    return local graph_text = `"`gr_txt'"'
  }
  return matrix meansd = _msd
  return scalar cutoff = `ncut'
  return scalar fn = fn
  return scalar fp = fp
end
  
capture mata mata drop ncut()
mata:
  function ncut(m1, sd1, m2, sd2){
    real scalar a, b, c, d
    
    if ( min((sd1, sd2)) <= 0 ) _error("SDs must be positive")
    if ( sd1 == sd2 ) return((m1 + m2) / 2)
    else {
      a = sd2 == sd1 ? 0.000001 : sd2^2 - sd1^2
      b = -2 * (m1 * sd2^2 - m2 * sd1^2)
      c = 2 * log(sd1/sd2) * sd1^2 * sd2^2 + m1^2 * sd2^2 - m2^2 * sd1^2
      d = b^2 - 4 * a * c
      if ( d < 0 ) _error("No solution")
      if (m1 < m2) return((-b + sqrt(d)) / 2 / a)
      else return((-b - sqrt(d)) / 2 / a)
    }
  }
end