// Given a grouping variable and a pairwise distance matrix, generate a
// binary variable marking the group medoids (cases with the minimum
// sum of distances to other variables

program getmedoids
syntax varlist(min=1 max=1), DISTmat(string) IDvar(varname) GENerate(string)

// Check dist-mat exists
qui matlist `distmat'[1,1]
// Insist on sort oder
qui isid `idvar', sort
// Generate medoid variable (initially missing)
qui gen `generate'=.

di "Translating cluster membership variable `varlist' into medoids index variable `generate'"
mata: getmedoids(st_matrix("`distmat'"),st_data(.,"`varlist'"),"`generate'")
end


mata:

// Changes newvar to have 1 for medoid cases, else 0 (by side-effect)
real matrix function getmedoids (real matrix distmat, real matrix group, string newvar) {
  Nk = max(group)
  newmedoids=J(Nk,1,.)
  N = rows(distmat)
  V = .
  st_view(V,.,newvar)
  for (i=1;i<=Nk;i++) {
    // Select the cluster
    selector = select(range(1,N,1), (group :== i))
    clusterdistmat = distmat[selector, selector]
    clustersize = rows(clusterdistmat)
    centroiddist = transposeonly(mean(clusterdistmat)) // For syntax highlighting, "'" causes issues
    mindist = min(centroiddist)
    /* countmedoids = 0 */
    for (row = 1; row <= clustersize; row++) {
      if (abs(mindist - centroiddist[row]) < 0.00000000000001) // "==" should work, since one is a copy of the other
      {
        /* countmedoids = countmedoids + 1 */
        newmedoids[i] = selector[row] // reverse select
      }
    }
    /* if (countmedoids>1) { */
    /*   sprintf("More than one medoid in cluster %3.0f: %3.0f", i, countmedoids) */
    /*   } */
  }
  newmedoids=uniqrows(newmedoids)
  V[.,1] = J(N,1,0)
  V[newmedoids,1] = J(rows(newmedoids),1,1)
}
end