program clfuzz, rclass
syntax newvarlist(min=1 max=1), DISTmat(string) IDvar(string) K(real) [ FUZZiness(real 3.0) ]
tempname medvar uij medoids V minimand

mata: result = main("`distmat'", `k', `fuzziness', "`medvar'")
mata: st_matrix("`uij'", result.U)
mata: st_matrix("`V'", result.V)
mata: st_numscalar("`minimand'", result.minimand)
mata: shipvar(result.V, "`medvar'",`=_N')


qui getgroup `medvar', dist(`distmat') id(`idvar') gen(`varlist')

return matrix uij = `uij'
return matrix V = `V'
return scalar minimand = `minimand'

end



mata

struct mainval {
  real matrix V
  real matrix U
  real scalar minimand
}

// From pam.ado
real matrix function random_medoids(real Nk, real matrix distmat) {
  N = rows(distmat)
  // Initialise with random medoids (Nk)
  medoids = uniqrows(floor(runiform(Nk,1):*N) :+ 1)
  // Replace initial medoids if there are duplicates (zero on off-diag)
  while ( (length(medoids)<Nk) | (min(distmat[medoids,medoids] + I(length(medoids))) == 0)) {
    medoids = uniqrows(floor(runiform(Nk,1):*N) :+ 1)
  }
  return(medoids[order(medoids,1)])
}

real matrix makeU(real matrix distmat, real matrix medoids, real scalar m) {
  // The membership weight NxNk matrix
  Nk = length(medoids)
  N = rows(distmat)
  dic = distmat[,medoids]
  dicf = (1 :/ dic) :^ (1/(m-1))
  dicfs = rowsum(dicf)
  U = dicf :/ dicfs
  /* U = editmissing(U,0) */

  for (i=1; i<=Nk; i++) {
      // "1 0 0 0" etc for the medoid cases
      // Should deal with duplicates too, probably by looking at distmat==0
      U[medoids[i],] = J(1,Nk,0) //1/(Nk-1))
      U[medoids[i],i]= 1
    }

  //  U,dic
  return(U)
}

real matrix reassignV(real matrix distmat, real matrix medoids, real matrix U, real scalar m) {
  N = rows(distmat)
  Nk = length(medoids)
  newmedoids = medoids

  for (k = 1; k<=Nk; k++) {
    ddk = transposeonly(distmat :* U[,k]:^m)
    ddksum =rowsum(ddk)
    mindist = min(ddksum)
    for (i = 1; i<=N; i++) {
      if ((ddksum[i]==mindist)) {
        newmedoids[k]=i
        break
      }
    }
  }
  newmedoids = newmedoids[,order(transposeonly(newmedoids),1)]
  return(newmedoids[order(newmedoids,1)])
}
end

mata
void shipvar(real matrix medoids, string medstr, real scalar N) {
  Nk = rows(medoids)
  varnames = (medstr)
  medvar = J(N,1,0)
  for (i=1; i<=Nk; i++) {
    medvar[medoids[i]]=1
  }
  View = .
  idx = st_addvar("double",medstr)
  st_view(View,.,medstr)
  View[1..N,.] = medvar
}
end


mata
struct mainval main(string distmat, real scalar Nk, real scalar m, string medstr) {
  struct mainval scalar retval
  dd = st_matrix(distmat)
  N = rows(dd)
  V = random_medoids(Nk, dd)

  Vhist = transposeonly(V)
  for(i=1;i<=100;i++) {
    U = makeU(dd, V, m)
    Vold = V
    V = reassignV(dd, V, U, m)
    Vhist = Vhist \ transposeonly(V)
    sprintf("Iter %1.0f: %10.7f", i, sum(dd[,V] :* U:^m))
    if (Vold == V) break
  }
  "Medoids history"
  Vhist
  if (Vold != V) {
    "Warning: end of iterations reached without convergence"
  }

  // Minimand: sum of U^m times distance to medoids
  minimand = sum(dd[,V] :* U:^m)

  retval.V = V
  retval.U = U
  retval.minimand = minimand
  return(retval)
}

end