// Copyright Brendan Halpin 2016 brendan.halpin@ul.ie
// Distribution is permitted under the terms of the GNU General Public Licence
// Apr  1 2016
// Calculate and plot silhouette width for a cluster solution and its pairwise distance matrix


program define silhouette, sortpreserve
version 10.0
syntax varlist(min=1 max=1) , DISTmat(string) IDvar(varname) [SILH(string) GEN(string) NAME(string) *]
tempname gv groupsize gsize gm so distmat2

if (("`silh'" != "") & ("`gen'" != "")) {
  di as error "Options SILH and GEN are synonyms; use one only"
  exit 499
}

if (("`silh'" != "") & ("`gen'" == "")) {
  local gen "`silh'"
}

if ("`name'"!="") {
  local name "name(`name')"
}

gen `so' = _n

// Check dist-mat exists
qui matlist `distmat'[1,1]

// Insist on correct sort order, and that it is unique
qui des, varlist
local so `r(sortlist)'
local mainsort : word 1 of `so'
if ("`mainsort'" != "`idvar'") {
  di in red "Error: data must be sorted by same ID variable as used for defining distances"
  error 5
}
isid `idvar'

//marksample touse
qui tab `varlist' /* if `touse' */, matcell(`groupsize')
bysort `varlist': gen `gsize' = _N
sort `idvar'

/* mkmat `touse', matrix(`touse') */
/* mata: `touse'=st_matrix("`touse'") */

//qui putmata `gv'=(`varlist'), replace

mata: st_matrix("`gv'", st_data(.,"`varlist'"))
mata: `gv'=st_matrix("`gv'")
mata: `distmat2' = st_matrix("`distmat'")
mata: silhgroup(`distmat2',`gv',st_matrix("`groupsize'"),"`gm'")
mata: mata drop `gv'

qui su `varlist' /* if `touse' */
local ngroups `=r(max)'
forvalues x = 1/`ngroups' {
  tempvar `gm'`x'
}

svmat `gm', names(`gm')

tempname ai bi sil idx zero

forvalues x = 1/`ngroups' {
  capture gen `ai' = .
  capture gen `bi' = . // intentionally using missing as bignum
  qui replace `ai' = `gm'`x'*(`gsize'/(`gsize'-1)) if `x'==`varlist'
  qui replace `bi' = `gm'`x' if `gm'`x' < `bi' & `x'!=`varlist'
}

qui gen `sil' = (`bi' - `ai')/max(`ai', `bi')
label var `sil' "Silhouette"

if "`gen'" != "" {
  gen `gen' = `sil'
}

table `varlist' /* if `touse' */, c(n `sil' min `sil' mean `sil' max `sil') format(%5.2f)

gsort `varlist' -`sil'

gen `idx' = _n //+ 10*(`varlist'-1)
gen `zero' = 0
label var `sil' "Silhouette"

//local gcmd "pcspike `zero' `idx' `sil' `idx' if `varlist'==1 & `touse'"
local gcmd "pcspike `zero' `idx' `sil' `idx' if `varlist'==1, `options'"

forvalues x = 2/`ngroups' {
//  local gcmd "`gcmd' || pcspike `zero' `idx' `sil' `idx' if `varlist'==`x' & `touse'"
  local gcmd "`gcmd' || pcspike `zero' `idx' `sil' `idx' if `varlist'==`x',  `options'"
}

twoway `gcmd' legend(off) ytitle("Silhouette width") xtitle("Cases") `name'

end

mata:
  real matrix silhgroup(real matrix dist, real vector groupvar, real vector groupsize, string gm) {
    real scalar ngroups, i, N
    real matrix selector, groupmeans

    ngroups  = rows(groupsize)
    N=rows(dist)

    groupmeans = J(rows(dist), 0, .)

    // Order the distance matrix by the grouping variable, columns only

    for (i=1; i<=ngroups; i++) {
      selector = select(range(1,N,1), (groupvar :== i))
      groupmeans = groupmeans, transposeonly(mean(transposeonly(dist[.,selector])))
    }

    st_matrix(gm, groupmeans)
  }
end