/* Copyright 2007-2017 Brendan Halpin brendan.halpin@ul.ie
   Distribution is permitted under the terms of the GNU General Public Licence
*/
#delimit ;

program omamatv3, plugin;

   mata:;
   real matrix function expandpwdist(real matrix raw, real matrix seqid, real matrix nd)
     {
     output = mm_expand(raw,nd,nd,1)
     return(output[invorder(seqid), invorder(seqid)])
     }
     end;
   
program define sdhollister;
version 9;
   syntax varlist [if] [in] [using/] ,
     SUBSmat(string)
     PWDist(string)
     LENgth(string)
     TIMEcost(real)  /* x in formula */
     LOCalcost(real) /* y in formula */
     [WORkspace STAndard(string) DUps];

  local norm 0;
   if ("`standard'"=="longer") {;
      local norm 1;
      };
   else if (inlist("`standard'","none")) {;
      local norm 0;
      };
   if (`norm'==1) {;
      di "Normalising distances with respect to length";
      };
   else {;
      di "Not normalising distances with respect to length";
      };
   
   
   marksample touse, novarlist; // novarlist mean keep cases with missing vars

   tempname twtype;
   local twtype 4;    
   
   tempvar idvar;
   tempvar lengthvar;
   gen `lengthvar' = `length';
   tempname indelcost;
   scalar `indelcost' = 0.0;

   /* tempvar hol_x; */
   /* scalar `hol_x' = `timecost'; */
   /* tempvar hol_y; */
   /* scalar `hol_y' = `localcost'; */
   local hol_x = `timecost';
   local hol_y = `localcost';
   

   
   local printworkspace 0;
   if "`workspace'" ~= "" {;
      local printworkspace 1;
      };

   
   local adjdur 1;
   local facexp 1;

   tempname ndups;
   tempname first;

   preserve;

   gen `idvar'=_n;

   sort `varlist';
   //                   mkmat `idvar';
   mata: st_matrix("`idvar'", st_data(.,"`idvar'"));
   by `varlist': gen `ndups' = _N;
   by `varlist': gen `first' = _n==1;

   qui count if `first';
   di "`r(N)' unique observations";
   qui keep if `first';

   // mkmat `ndups';
   mata: st_matrix("`ndups'", st_data(.,"`ndups'"));
   

   // matrix `pwdist' = J(_N,_N,0);
   mata: st_matrix("`pwdist'",   J(`=_N',`=_N',0));
   
     /* Arguments hardcoded into omamatv3:
     0: substitution matrix name
     1: indel cost
     2: output matrix
     3: dimensions of subsmatrix
     4: adjust for duration?
     5: show workspace?
     6: exponent
     */
   /* Checks?
   1: is subsmat a matrix, a square matrix, with dimension >= n-states
   2: is indel an integer? relate to max subscost?
   3: let pwdist be a name only
*/

   scalar subsrows = rowsof(`subsmat');
   scalar subscols = colsof(`subsmat');
   if subsrows!=subscols {;
      di "Error: non square substitution matrix";
      exit;
      };

   plugin call omamatv3 `idvar' `lengthvar' `varlist',
     `subsmat' `indelcost' subsrows `pwdist' `adjdur' `printworkspace' `facexp' `twtype' 0 0 0
     `hol_x' `hol_y' `norm'
     ;
   capture mata mata which mm_expand();
   if _rc {;                                                                              
      di as error "mm_expand() from -moremata- is required; type -ssc install moremata- to obtain it";
      exit 499;           
      };
   mata: `pwdist'= expandpwdist(st_matrix("`pwdist'"),st_matrix("`idvar'"),st_matrix("`ndups'"));
   mata: st_matrix("`pwdist'",`pwdist');
   restore;
   
   
end;