/* Copyright 2007 Brendan Halpin brendan.halpin@ul.ie Distribution is permitted under the terms of the GNU General Public Licence Jul 19 2008 08:23:29 copy of timewarp.ado $Id: twed.ado,v 1.14 2017/02/06 10:41:20 brendan Exp $ $Log: twed.ado,v $ Revision 1.14 2017/02/06 10:41:20 brendan Summary: refstr to ref Revision 1.13 2015/02/14 11:25:23 brendan Summary: Bypass Stata matrix to lift 11k limit, add DUPS option Revision 1.12 2014/03/30 21:43:09 brendan Summary: Regularised duplicates message Revision 1.11 2012/06/28 23:10:27 brendan Made version string optional, fixed numeric/string issue with normalisation parameter Revision 1.10 2012/06/21 20:28:19 brendan Standardising message shouldn't be a warning Revision 1.9 2012/06/20 14:58:48 brendan Included mm_expand check/hint Revision 1.8 2012/06/16 22:49:16 brendan Correcting "standard not recognised" message Revision 1.7 2012/06/15 23:43:00 brendan Normalisation Revision 1.6 2012/06/15 21:42:10 brendan Added normalisation and proper handling of `touse' Revision 1.5 2011/11/01 21:13:04 brendan Removed unneeded save of "hold" matrix, copy of pwdist Revision 1.4 2011/10/03 22:08:10 brendan Check \$id\$ */ #delimit ; capture program drop twed; capture program drop omamatv3; 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 twed; version 9; syntax varlist [if] [in] [using/] , SUBSmat(string) LAMbda(real) NU(real) PWDist(string) LENgth(string) [REF(integer 0) WORkspace DUps STAndard(string) VERsion] ; if ("`version'" != "") {; di "TWED version: \$Id: twed.ado,v 1.14 2017/02/06 10:41:20 brendan Exp $"; }; 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"; }; if ("`dups'"=="") {; local dups 0; }; else {; local dups 1; }; if (`ref'!=0) {; local dups 1; }; marksample touse, novarlist; // novarlist mean keep cases with missing vars preserve; keep if `touse'; tempname twtype; local twtype 3; tempvar idvar; tempvar lengthvar; gen `lengthvar' = `length'; tempname indelcost; scalar `indelcost' = 0; local printworkspace 0; if "`workspace'" ~= "" {; local printworkspace 1; }; local adjdur 1; local facexp 1; tempname ndups; tempname first; gen `idvar'=_n; if (`dups'==0) {; 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'")); }; //this setting needs to be here since _N has changed local ncols `ref'; if `ncols' == 0 {; local ncols = _N; }; //matrix `pwdist' = J(_N,`ncols',0); mata: st_matrix("`pwdist'", J(`=_N',`ncols',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' `lambda' `nu' `ref' 1 1 `norm'; if (`dups'==0) {; capture mata mata which mm_expand(); if _rc {; di as error "mm_expand() from -moremata- is required; type -ssc install moremata- to obtain it"; di as error "Alternatively, use the {cmd:dups} option to treat duplicate sequences"; exit 499; }; mata: `pwdist'= expandpwdist(st_matrix("`pwdist'"),st_matrix("`idvar'"),st_matrix("`ndups'")); mata: st_matrix("`pwdist'",`pwdist'); mata: mata drop `pwdist'; // Drop the mata copy of the PWdist matrix, no longer needed & potentially large }; restore; end;