// May  1 2017 20:52:22
// Make dendrohmap.ado self contained


capture program drop dhm
program define dhm, sortpreserve
syntax, DISTmat(string) IDvar(varname) [BY(varname)  ORDer(varname) EXPonent(real 1.0) LEVels(real 100) PALette(real 1) BOX RGBperm(real 1) NAME(string) *]

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

if ("`name'"!="") {
  local name name(`name')
}
tempname cumu ncumu dord rawx
if ("`by'"!="") {
  qui tab `by', matcell(`rawx')
  mata: cumu = cumulate("`rawx'")
  mata: st_matrix("`cumu'",cumu)

  local xcumulist "0.5"
  local lablist ""
  qui su `by'
  forvalues x = 1/`=r(max)' {
    local xcumulist `xcumulist' `=`cumu'[`x',1]+0.5'
    local lablist `lablist' `=`cumu'[`x',1]+0.5-`rawx'[`x',1]/2' `"`x'"'
  }
  local ncumulist `"ytick(`xcumulist')"'
  local xcumulist `"xtick(`xcumulist')"'
  local xlablist `"xlabel(`lablist', noticks)"'
  local ylablist `"ylabel(`lablist', noticks)"'
}

tempvar so g999
if ("`order'"=="") {
  cluster generate `g999'=groups(`=_N'), ties(fewer) `name'
}
else {
  gen `so' = _n
  sort `by' `order'
  generate `g999'= _n
  sort `so'
}

mata: st_matrix("`dord'", dendrohmap("`distmat'", "`g999'"))

//pmc, mat(`dord') legend(off) xsize(4) ysize(4) `xcumulist' `ncumulist' `addboxlist' `options'

dhmhmapmat, mat(`dord') levels(`levels') exponent(`exponent') pal(`palette') rgb(`rgbperm') `box' `xcumulist' `ncumulist' `xlablist' `ylablist' `options'

end

program dhmpalette, rclass
// Take a number in 0-1, return 3 0-255 RGB values
args zlev palette permute
local rgb1 = 0
local rgb2 = 0
local rgb3 = 0

if (`palette' == 1) {
  local rgb1 = min(1, 3*(`zlev'))
  local rgb2 = min(1, 3*max(0, (`zlev'-1/3)))
  local rgb3 = min(1, 3*max(0, (`zlev'-2/3)))
}
else if (`palette' == 2) {
  if inrange(`zlev',0,1/6) {
    local rgb1 = `zlev'*6
    local rgb2 = 0
    local rgb3 = 0
  }
  if inrange(`zlev',1/6,2/6) {
    local rgb1 = 1
    local rgb2 = 6*(`zlev'-1/6)
    local rgb3 = 0
  }
  if inrange(`zlev',2/6,3/6) {
    local rgb1 = 6*(3/6 - `zlev')
    local rgb2 = 1
    local rgb3 = 6*(`zlev'-2/6)
  }
  if inrange(`zlev',3/6,4/6) {
    local rgb1 = 0
    local rgb2 = 6*(4/6 - `zlev')
    local rgb3 = 1
  }
  if inrange(`zlev',4/6,5/6) {
    local rgb1 = 3*(`zlev'-4/6)
    local rgb2 = 0
    local rgb3 = 6*(5/6 - `zlev')
  }
  if inrange(`zlev',5/6,6/6) {
    local rgb1 = 3*(`zlev'-4/6)
    local rgb2 = 6*(`zlev' - 5/6)
    local rgb3 = 6*(`zlev' - 5/6)
  }
}
else if (`palette' == 3) {
  if inrange(`zlev',0,1/2) {
    local rgb1 = `zlev'
    local rgb2 = 1 - 2*`zlev'
    local rgb3 = 2*`zlev'
  }
  if inrange(`zlev',1/2,1) {
    local rgb1 = `zlev'
    local rgb2 = 2*`zlev' - 1
    local rgb3 = 1
  }
}
else if (`palette' == 4) {
local rgb1 = `zlev'
if inrange(`zlev',0,1/3) {
  local rgb2 = 3*`zlev'
  local rgb3 = 0
}
if inrange(`zlev',1/3,3/3) {
  local rgb2 = 1
  local rgb3 = 1.5*(`zlev'-1/3)
}
}
else if (`palette' == 5) {
local rgb1 = `zlev'
local rgb2 = 0
local rgb3 = 0
if (`zlev'>2/3) {
  local rgb2 = (`zlev'-1/3)*3/2
}
if (`zlev'>2/3) {
  local rgb3 = (`zlev'-2/3)*3
}
}

if (`permute' == 1) {
  return scalar R = int(`=255*`rgb1'')
  return scalar G = int(`=255*`rgb2'')
  return scalar B = int(`=255*`rgb3'')
}
else if (`permute' == 2) {
  return scalar R = int(`=255*`rgb1'')
  return scalar B = int(`=255*`rgb2'')
  return scalar G = int(`=255*`rgb3'')
}
else if (`permute' == 3) {
  return scalar G = int(`=255*`rgb1'')
  return scalar R = int(`=255*`rgb2'')
  return scalar B = int(`=255*`rgb3'')
}
else if (`permute' == 4) {
  return scalar G = int(`=255*`rgb1'')
  return scalar B = int(`=255*`rgb2'')
  return scalar R = int(`=255*`rgb3'')
}
else if (`permute' == 5) {
  return scalar B = int(`=255*`rgb1'')
  return scalar R = int(`=255*`rgb2'')
  return scalar G = int(`=255*`rgb3'')
}
else if (`permute' == 6) {
  return scalar G = int(`=255*`rgb1'')
  return scalar B = int(`=255*`rgb2'')
  return scalar R = int(`=255*`rgb3'')
}

end

program define dhmhmapmat
syntax, MATrix(string) PALette(real) RGBperm(real) [LEVels(real 100) BOX EXPonent(real 1.0) *]

local asprat = rowsof(`matrix')/colsof(`matrix')

preserve
tempname v
clear
mata: d = hmapmat("`matrix'")

getmata (v*)=d

qui su v3
local scalemax = r(max)
gen v4 = int((`levels'*(v3/`r(max)')^`exponent'))

local addplot ""
if (("`box'"!="") & ("`by'"!="")) { // ignore box if by not given
  mata: boxcoords(cumu)
  local addplot "addplot(line x1 x2, lcolor(black) || line x3 x4, lcolor(black) legend(off))"
}


tempname dhmmat dhmscale
dhmmyhmap v2 v1 v4, nosc xlabel(none) ylabel(none) aspect(`asprat') reverse name(`dhmmat') `addplot' pal(`palette') rgb(`rgbperm') `options'
dhmdrawscale, bottom(0) top(`scalemax') fysize(50) fxsize(10) vertical name(`dhmscale') exponent(`exponent') pal(`palette') rgb(`rgbperm')
graph combine `dhmmat' `dhmscale'
graph drop `dhmmat' `dhmscale'

restore

end


program define dhmmyhmap
version 8.2
*! based on hmap version 1.0 16mar2010 by Austin Nichols <austinnichols@gmail.com>
syntax varlist(min=3 max=3) [if] [in] ,  PALette(real) RGBperm(real) [ XLabel(passthru) YLabel(passthru) ADDPlot(string) REVerse MONochrome noSCatter * ]
marksample touse
tokenize `varlist'

xunits `1' `touse'
local xdist=r(units)
xunits `2' `touse'
local ydist=r(units)
if (`"`xlabel'"'=="") {
  su `1', mean
  local xlabel `"xlabel(`r(min)'(`xdist')`r(max)', val angle(90) nogrid)"'
}

if (`"`ylabel'"'=="") {
  su `2', mean
  local ylabel `"ylabel(`r(min)'(`ydist')`r(max)', val angle(0) nogrid)"'
}
if (`"`addplot'"'!="") {
  local addplot `"||`addplot'"'
}

preserve

tempvar y0 y1
gen `y1'=`2'+`ydist'/2
label value `y1' `:val lab `2''
gen `y0'=`2'-`ydist'/2
qui levelsof `3' if `touse', loc(lev)
su `3' if `touse', mean
local minlev=r(min)
local maxlev=r(max)
if `maxlev'==`minlev' {
  di as err "variable `3' does not vary; graph would be a solid color"
  err 198
}

local gcom
foreach x of local lev {
  local zlev=(`x'-`minlev')/(`maxlev'-`minlev')
  if ("`reverse'"!="") {
    local zlev=1-`zlev'
  }
  if ("`monochrome'"!="") {
    local R=int(255*`zlev')
    local G=int(255*`zlev')
    local B=int(255*`zlev')
  }
  else {
    dhmpalette `zlev' `palette' `rgbperm'
    local R=r(R)
    local G=r(G)
    local B=r(B)
  }
  local gcom `"`gcom'||rbar `y1' `y0' `1' if (abs(`3'-`x')<0.000000001), barw(`xdist') col("`R' `G' `B'") fi(inten100) lw(none)"'
}
if "`scatter'"=="" local sc `"sc `2' `1', mlab(`3') mlabpos(0) msy(i) `scopt'"'
local gcom `"`gcom' legend(off)  yscale(reverse) `ylabel' `xlabel'  xtitle("") ytitle("")"'
// di `"G: `gcom'"'
// di `"O: `options'"'
timer clear
timer on 1
twoway `gcom' plotr(fc(gray) m(zero)) `options' `addplot' ||`sc'
timer off 1
timer list
restore
end


program xunits, rclass
args v touse
qui summ `v' if `touse', mean
local p = 1
capture assert float(`v') == float(round(`v',1)) if `touse'
if _rc == 0 {
  while _rc == 0 {
    local p=`p'*10
    capture assert float(`v') == float(round(`v',`p')) if `touse'
  }
  local p=`p'/10
}
else {
  while _rc {
    local p=`p'/10
    capture assert float(`v') == float(round(`v',`p')) if `touse'
  }
}
return scalar units = `p'
end

// ----


program dhmdrawscale
syntax , [ BOTtom(real 0.0) TOP(real 1.0) VERtical EXPonent(real 1.0) NBLocks(integer 11) *] PALette(real) RGBperm(real)

local min = `bottom'
local max = `top'

preserve
qui {
  clear
  set obs `nblocks'
  gen n = _n
  gen scalevar =   `min' +  (`max' - `min')*(_n-1)/`=`nblocks'-1'
  gen z = round(scalevar^`exponent'/`max', 0.000001)
  gen x = 1

  if ("`sqrt'" != "") {
    replace z=sqrt(abs(z))
  }

  gen lab = string(scalevar, "%5.1f")
  replace lab = `"""' + lab + `"""'

  local lab ""
  forvalues x = 1/`nblocks' {
    local lab = `"`lab' `x' `=lab[`x']'"'
  }
  if ("`vertical'" != "") {
    noi dhmmyhmap x n z, reverse nosc aspect(`nblocks') ylabel(`lab', labsize(medium) angle( 0)) xlabel(none) title("Scale", size(medium)) pal(`palette') rgb(`rgbperm') `options'
  }
  else {
    noi dhmmyhmap n x z, reverse nosc aspect(`=1/`nblocks'') xlabel(`lab', labsize(medium) angle(90)) ylabel(none) title("Scale", size(medium)) pal(`palette') rgb(`rgbperm') `options'
  }
}

restore
end


mata:
  real matrix dendrohmap(string pwd, string dendorder) {
    romat = st_matrix(pwd)
    dord = st_data(.,dendorder)
    romat = romat[order(dord,1), order(dord,1)]
    return(romat)
  }

  real matrix cumulate(string freqs) {
    c = st_matrix(freqs)
    for (i=2; i<=rows(c); i++) {
      c[i] = c[i-1]+c[i]
    }
    return(c)
  }

  real matrix hmapmat(string mat) {
    d = st_matrix(mat)
    v = vec(d)

    i = mm_expand(1::rows(d), cols(d))
    j = floor((((1::rows(v)) :- 1) :/ rows(d)) :+ 1)

    v = i, j, v
    return(v)
  }

real function boxcoords (real matrix cumul) {
  n = 3*rows(cumul)
  x = J(n,4,0)

  epsilon = 0.5
  prior = 0

  for (i=1; i<=rows(cumul); i++) {
    x[(i-1)*3+1,] = J(1,4,prior)
    x[(i-1)*3+2,] = (prior, cumul[i], cumul[i], prior)
    x[(i-1)*3+3,] = (cumul[i], cumul[i], cumul[i], cumul[i])
    prior = cumul[i]
  }

  x = x :+ epsilon

  idx=st_addvar("double",("x1", "x2", "x3", "x4"))
  st_view(V=.,.,.)
  V[1..n,idx]= x
}




end