capture program drop dudahart
program define dudahart, rclass
syntax, DISTmat(name) IDvar(varname) [NGroups(integer 15) GRaph(string) NAME(name) *]

version 10

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

// 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'


mata: pwd = st_matrix("`distmat'")

tempname resmat

tempvar d1
gen `d1' = 1

tempvar index dh dht
qui {
  gen `index' = .
  gen `dh' = .
  gen `dht' = .
}

forvalues x = 2/`=1+`ngroups'' {
  tempvar d`x'
  cluster gen `d`x'' = groups(`x'), ties(fewer) `name'
}


mata: result = J(0,2,.)
forvalues soln=1/`ngroups' {
  local next = `soln'+1
  local whichsplit 0
  local new1 0
  local new2 0
  qui su `d`soln''
  forvalues x = 1/`r(max)' {
    qui su `d`next'' if `d`soln'' == `x'
    if (`r(min)'!=`r(max)') {
      local whichsplit `x'
      local new1 `r(min)'
      local new2 `r(max)'
    }
  }
  //di "`soln'/`next': `whichsplit' `new1' `new2'"

  tempvar j1 j2
  qui gen `j1' = `d`soln''==`whichsplit'
  qui gen `j2' = `d`next''==`new1'
  qui replace `j2' = 2 if `d`next''==`new2'

  qui putmata j1=`j1' j2=`j2', replace

  mata: j11 = j2 :==1
  mata: j12 = j2 :==2
  mata: SSt = sum(select(select(pwd, j1), j1'))/sum(j1)
  mata: SS1 = sum(select(select(pwd, j11), j11'))/sum(j11)
  mata: SS2 = sum(select(select(pwd, j12), j12'))/sum(j12)
  mata: dh = (SS1+SS2)/SSt
  mata: t = (1/dh - 1)*(sum(j1) - 2)
  mata: result = result\(dh, t)

}
mata: st_matrix("`resmat'",result)

local maxdigs 3
local maxwid 9
local maxwid2 9

di
di as txt "{c TLC}{hline 13}{c TT}{hline 27}{c TRC}"
di as txt "{c |}             {c |}   Duda/Hart on distances  {c |}"
di as txt "{c |}  Number of  {c |}             {c |}  pseudo     {c |}"
di as txt "{c |}  clusters   {c |} Je(2)/Je(1) {c |}  T-squared  {c |}"
di as txt "{c LT}{hline 13}{c +}{hline 13}{c +}{hline 13}{c RT}"
forvalues i = 1/`ngroups' {
  di as txt "{c |} " _c
  di as res "{center 11:{ralign `maxdigs':`i'}}" _c
  di as txt _col(15) "{c |} " _c
  local tmp : di %9.4f `resmat'[`i',1]
  local tmp `tmp'
  di as res "{center 11:{ralign `maxwid':`tmp'}}" _c
  di as txt _col(29) "{c |} " _c
  local tmp : di %9.2f `resmat'[`i',2]
  local tmp `tmp'
  di as res "{center 11:{ralign `maxwid2':`tmp'}}" _c
  di as txt _col(43) "{c |}"
  return scalar duda_`i' = `resmat'[`i',1]
  return scalar dudat2_`i' = `resmat'[`i',2]
  qui {
    replace `index' = `i' in `i'
    replace `dh' = `resmat'[`i',1] in `i'
    replace `dht' = `resmat'[`i',2] in `i'
  }

}
di as txt "{c BLC}{hline 13}{c BT}{hline 13}{c BT}{hline 13}{c BRC}"

if ("`graph'"!="") {
  if ("`graph'"=="dh") {
    line `dh' `index', title("Duda-Hart index") ytitle("DH index") xtitle("N-Clusters")||scatter `dh' `index', legend(off) `options'
  }
  else if ("`graph'"=="dht") {
    line `dht' `index', title("Duda-Hart T-squared") ytitle("DH Tsq") xtitle("N-Clusters")||scatter `dht' `index', legend(off) `options'
  }
  else {
    line `dh' `dht' `index', title("Duda-Hart index") ytitle("DH index") xtitle("N-Clusters")||scatter `dh' `dht' `index', legend(off) `options'
  }
}

end