*! Version : 1.30  (TSJ-1: st0003)
*! Date    : 19 Oct 2001
*! Author  : Adrian Mander
*! Email   : adrian.mander@mrc-bsu.cam.ac.uk

/* 2/3/01 bug fix in haplo expansion for haplotypes
 *********************************************************************
 * 8/11/99 Now I want to implement some log-linear modelling within the
 * hapfreq3.ado algorithm. Basically the expansion of the dataset
 * is still carried out but instead of a simple calculating
 * algorithm I shall use the log-linear approach hence any model
 * can be specified
 *************************************************************************************
 * 15/11/99 About to add the LDIM as the varlist to calculate the likelihood over.
 * This is important otherwise the algorithm will use the minimal continguency
 * table and when comparing tables you need the same dimension tables. The varlist
 * defines the table. Not sure about the missing data and cells :(
 *************************************************************************************/

program define hapipf, rclass
version 7.0
syntax [varlist] [using/] [if] ,[ LDIM(varlist) IPF(string) START DISplay EXPect KNOWN PHASE(varname) ACC(real 0.001) IPFACC(real 0.000001) DEBUG(integer 0) NOLOG MODEL(integer 0) LRTEST(numlist) RARE(real 0) CONVARS(string) CONFILE(string) QUIET NOISE MV MVDEL MENU]
                                                                 
if "`menu'"~="" {
  cap _qmenu
  if _rc==3000 {
    di in red "Command will not be run"
    exit
  }
  if _rc==3001 {
    di in blue "Running $Qcommand..."
    $Qcommand
    exit
  }
}
                                                                 
tokenize `varlist'

cap which ipf
if _rc~=0 {
  di in red "YOU MUST INSTALL ipf.ado which was first introduced in STB55"
  di "This function performs the log-linear modelling"
  exit
}

tempfile origin
qui save "`origin'",replace

if "`if'"~="" { keep `if'}

global all_name = "`varlist'"
global nloc 0
global nsub= _N
                                                                
/***************************************************
 * Delete all lines with missing marker information
 * for each marker in the varlist
 ***************************************************/
                                                              
  if "`mv'"~="" & "`known'"~="" { di as error "If you want to impute missing data then phase is not known for those people with missing typings" }
  
  if "`mvdel'"~="" {
    di
    local len = length("You have selected the deletion of missing lines")
    di as text "You have selected the deletion of missing lines"
    di in smcl "{dup `len':{c -}}"
    
    di as result "Variable [variable label] (#lines deleted)"
    while "`1'"~="" {
      qui count if `1'==.
      local lab : variable label `1'
      if "`lab'"~="" { di  _continue "`1' [`lab'] (" %4.0f `r(N)' ")" }
      else { di  _continue "`1'  (" %4.0f `r(N)' ")" }
      qui drop if `1'==.
      mac shift 1
      qui count if `1'==.
      local lab : variable label `1'
      if "`lab'"~="" { di  _col(29) "`1' [`lab'] ("  %4.0f `r(N)' ")" }
      else { di  _col(29) "`1'  ("  %4.0f `r(N)' ")" }
      qui drop if `1'==.
      mac shift 1
    }
  }
tokenize "`varlist'"
                                                            
/*********************************************
 * Pair the data from the varlist
 * and check that they are numeric
 *********************************************/

di
di as text "Marker information"
di in smcl "{dup 18:{c -}}"

local i 1
while "`1'"~="" {
   if "`2'"=="" { di in red "There must be paired data" }
   global nloc =$nloc+1
   local wc1=2*$nloc-1
   local wc2=`wc1'+1
   local root1: word `wc1' of $all_name
   local root2: word `wc2' of $all_name
   cap confirm string variable `root1'
   if _rc==0 {
     di "`root1' cannot be a string variable"
     exit(7)
   }
   cap confirm string variable `root2'
   if _rc==0 {
     di "`root2' cannot be a string variable"
     exit(7)
   }
   if "`quiet'"=="" {
     local lab1 : variable label `root1'
     local lab2 : variable label `root2'
     if "`lab1'"~="" & "`lab2'"~="" { di as res "Alleles for l$nloc are (`root1' , `root2') [`lab1' , `lab2']" }
     else { di as res "Alleles for l$nloc are (`root1' , `root2') " }
   }
   return local loci`i' = "`root1' `root2'"
   local i = `i'+1
mac shift 2
}
                                                                 
di

local nloc = $nloc
if "`nloc'"~="1" {
  local haptext "Haplotype"
  local haptvar "Haplo"
}
else {
  local haptext "Allele"
  local haptvar "Allele"
}

/**********************************
 * Check out the constraint file
 **********************************/

tempfile temp

qui save "`temp'",replace

if "`confile'"~="" {
   if "`convars'"=="" { di as error "You must specify the convars() option when using constraint files" }
   use "`confile'",replace
   tokenize "`convars'"
   cap confirm variable Efreqold
   if _rc~=0 {
     cap confirm variable Ifreq
     if _rc~=0 {
       di as error "Constraints file doesnt contain Efreqold or Ifreq!"
       exit(198)
     }
     else { rename Ifreq Efreqold }
   }
   sort `convars'
   qui save "`confile'",replace
   
   while "`1'"~="" {
     mac shift 1
   }
}

/*********************************************
 * Make sure about including the ipf() option
 *********************************************/

  if "`ipf'"=="" {
    di in red "You must specify the loglinear model in the ipf() option "
    exit(198)
  }

use "`temp'",replace

/**********************************************************************
 * haplotyp expands the dataset into the phases
 * Takes the varlist and constructs the haplotypes
 * locus1 and locus2
 * For phase unknown locus1 and locus2 contain all possible haplotypes
 **********************************************************************
 * Additionally the user can specify that there are missing values
 * in the DNA typings with the MV option
 */
 
if "`known'"~="" & "`phase'"=="" { haplotyp, known `mv'}
if "`known'"=="" & "`phase'"==""{ haplotyp, `mv'}
if "`phase'"~="" { haplotyp,known phase(`phase') `mv'}

/* What is the next line? Probably gets rid of the chance of duplicates */

qui {
drop `varlist'
compress
}

/***********************************************************
 * Start frequencies f1 - future thing to do is to 
 * implement a weighting system on the start points
 *
 * Note that if a subject has been expanded by phase then
 * the start frequencies are spread evenly by default
 ***********************************************************/

tempvar f1 tef1
if "`start'"=="" {
   sort subject
   qui by subject: gen double `f1'=1/_N 
}
else {
   sort subject
   qui gen `f1'=uniform()
   qui by subject: gen double `tef1'= sum(`f1')
   qui by subject: replace `f1'= `f1'/`tef1'[_N]
   drop `tef1'
}

/*********************************************
 *   EM algorithm
 *********************************************/

local cont 1
local it 0
local saveimp 0 

      cap confirm new variable locus
      if _rc~=0 {
        di in red "Rename the variable locus"
        exit(111)
      }

while(`cont'==1) { /*Start of while loop */
                       
   cap drop `pr' `prs'
   tempvar pr prs stupid plocus

   cap drop llh
   local it = `it'+1

/************************************************************
 *       # New estimate of haplotype relative frequencies
 ************************************************************/

   cap drop freq

      gen long `stupid' = _n /* this will preserve line numbers */
      qui _stack locus1 locus2
      qui rename locus1 locus
      qui drop locus2

     /* generate individual loci variables */

      local vlist " "
      qui gen str40 `plocus' = locus
      forval i=1/$nloc {
         tempvar len len2
         qui gen `len' = index(`plocus',".")-1
         qui gen l`i'=real(substr(`plocus',1,`len'))
         if `i'==$nloc { qui replace l`i'=real(`plocus') }
         else { qui gen `len2' = index(`plocus',".")+1
         qui replace `plocus'=substr(`plocus',`len2',.) }
         drop `len'
         cap drop `len2'
         local vlist "`vlist' l`i'"
      }
                   
   if `debug'>0 {
     tab locus
     list l* if locus==""
   }

   drop `plocus'

/************************************************************
 * Use the iterative proportional fitting algorithm and save
 * the expected frequencies in fit.dta
 * The variables present in the model are output into the
 * r(vlist) as this may contain more variables than loci
 ************************************************************/

tempfile fit
local file = substr("`fit'",1,index("`fit'",".")-1)


 if "`ldim'"=="" {
      if "`convars'"=="" {
         if "`noise'"=="" {  qui  ipf [fw=`f1'],exp fit(`ipf') save("`file'") acc(`ipfacc') }
         else{ ipf [fw=`f1'],exp fit(`ipf') save("`file'") acc(`ipfacc')  }
      }
      else {
        qui ipf [fw=`f1'],exp convars(`convars') confile("`confile'") fit(`ipf') save(`file') acc(`ipfacc')
      }
       
      local df=r(df)
      local nparms=r(parms)
      local ncells=r(ncells)
      sort `r(vlist)'
      global ipflist = "`r(vlist)'"
      merge `r(vlist)' using `file'   /* merge in the new frequencies */
 }
 else {
    qui ipf `ldim' [fw=`f1'], fit(`ipf') save(`file') acc(`ipfacc')
    local df=r(df)
    sort `ldim'
    global ipflist = "`ldim'"
    merge `ldim' using `file'   /* merge in the new frequencies */
 }

/******************************************************
 * Sometimes fit.dta has an unobserved line of 
 ******************************************************/

local rvlist "`r(vlist)'"
qui count if locus==""

if `debug'>1 {
  l locus caco Efreq prob if locus==""
  di "`r(N)' loci missing"
}

if `r(N)' > 0 {
  tempvar temprep
  gen `temprep' = locus==""
  forval i=1/$nloc {
    qui replace locus = cond( locus=="", locus + string(l`i'), locus +"."+ string(l`i')) if `temprep'==1
  }
}


if `debug'>2 {
  l locus caco Efreq prob 
}

/****************************************
 * Saving the expected Frequencies.
 ****************************************/

  if `saveimp'==1 {
    tempfile imputef now
    qui save "`now'"
    sort `rvlist'
    qui by `rvlist': keep if _n==1
    qui save "`imputef'"
    use "`now'",clear
  }

      drop `vlist' Efreq Ofreq _merge
 

      sort _stack `stupid'
      qui _ustack locus prob, by(_stack) val(1 2)
      drop _stack merg _merge

      drop `stupid' 
      cap drop npar
	
      qui gen npar = 1
      local npar = npar[_N]

/*************************************************************
 *     # New probabilities (per phase and per subject)
 ************************************************************/

   qui gen double `pr' = prob1*prob2
   drop prob1 prob2 

/******************************************************************************
 *     # Calculate log likelihood
 * 1/11/99 I think in the next part within subject `by' shouldnt change
 * hence could ignore `by' UNLESS in each by strata individuals ids are
 * determined
 ******************************************************************************/
 
   cap confirm new variable llh llhh
   if _rc~= 0 { di in red "rename the variable below:"
      confirm new variable llh llhh
   }
 
      sort subject  
      qui by subject : gen double `prs'=sum(`pr')
      qui by subject : gen double llh = log(`prs'[_N]) 
      qui by subject : replace `pr' = `pr'/`prs'[_N]

      qui by subject :gen double llhh= cond(_n==_N,llh[_n],0)

      qui replace llhh = sum(llhh)
      local llh = llhh[_N]
      drop llhh


      if "`nolog'"=="" { di as res "Iteration `it' loglhd = `llh'" }

/**************************************
 * Use the new weights for the ipf
 **************************************/
  
      qui replace `f1' = `pr'

/**************************************
 *       # Convergence test
 **************************************/

   if (`it'>1) {
     local cont = (`llh' - `lastllh')>`acc' 
     if (`llh' < `lastllh') {
       local tmpdiff = (`lastllh'-`llh')
       di in smcl "{err}Likelihood not increased in hapipf. Decrease= `tmpdiff'"
     }
   }

if `saveimp'==0 {
  if `cont'==0 {
    local saveimp 1
    local cont 1
  }
}

local lastllh = `llh'

}  /*end of while loop */


_stack locus1 locus2 `pr' `pr'
rename locus1 locus

/********************************************************************************
 * In order to display the results I need the variables included
 * in the ipf model so that the expected frequencies are split
 * in groups of this as well
 ********************************************************************************/

     local vlist "$ipflist"
     forval i=1/$nloc {
         local length= length("l`i'")
         local ind = index("`vlist'","l`i'")
         if `ind'>0 { local vlist = substr("`vlist'",1,`ind'-1)+substr("`vlist'",`ind'+`length',.) }
      }

sort locus `vlist'
qui by locus `vlist' : gen double freq = sum(`pr')
qui _unique subject
qui gen double eprob=freq/(2*r(unique))
qui by locus `vlist' : keep if _n==_N

sort `vlist' locus

/****************************************************************
 * A routine to remove rare haplotypes but not zero haplotypes
 * and also remove the degrees of freedom
 ****************************************************************/

if `rare'~=0 {
  di "Removing rare haplotypes...."

/*
  gen rare= eprob<`rare'
  count if rare==1
  local nrare=r(N)
                      
      local i 1
      local slist ""
      qui gen str40 plocus = locus
      while `i'<=$nloc {
         qui gen len = index(plocus,".")-1
         qui gen l`i'=real(substr(plocus,1,len))
         if `i'==$nloc { qui replace l`i'=real(plocus) }
         else { qui gen len2 = index(plocus,".")+1
         qui replace plocus=substr(plocus,len2,.) }
         drop len
         cap drop len2
   local slist "`slist' l`i'"
         local i=`i'+1
      }
     drop plocus
  gen Efreqold = cond(rare==1, 0 ,.)
  sort `slist'
  save constrain,replace
*/

  gen rare= eprob<`rare' & eprob>0
  count if rare==1
  local nrare=r(N)

  drop if rare==1

  list locus `vlist' freq eprob if rare==1, noobs


}

/**********************************
 * Display the testing expressions
 **********************************/

if "`quiet'"=="" {
  di as text ""
  local len = length("`haptext' Frequency Estimation by EM algorithm")
  di "`haptext' Frequency Estimation by EM algorithm"
  di in smcl "{dup `len':{c -}}"
  
  
  di as res " No. loci ", _col(20) "= $nloc"
  di " Log-Likelihood ", _col(20) "= `llh'"
  *di " Log-Likelihood under null ", _col(20) "= `llh0'"
  *di " 2*LogLikelihoodRatio  ", _col(20) "=", 2*(`llh'-`llh0')
  di " Df ", _col(20) "=" , `df'
  di " No. parameters ", _col(20) "=" , `nparms'
  di " No. cells ", _col(20) "=" , `ncells'
  /*di " No of rare parameters dropped ", _col(20) "=", `nrare'
  *di " p-value ", _col(20) "=", chiprob(`npar'-`npar0',2*(`llh'-`llh0')) */
  if "`confile'"~="" {
    di in smcl as text "{c TLC}{dup 45:{c -}}{c TRC}"
    di in smcl as text "{c |}" as res "WARNING: df wrong when using constraints file" as text "{c |}"
    di in smcl as text "{c BLC}{dup 45:{c -}}{c BRC}"
  }
  di
}

/************************************************
 * Should display the expected frequencies and
 * the imputed frequencies.
 ************************************************/

if "`display'"~="" {
   sort locus `vlist'
   di in gr "Imputed Frequencies"
   qui gen str80 `haptvar'=locus
   qui compress `haptvar'
   list `haptvar' `vlist' freq eprob, noobs
   qui drop `haptvar'
   di in gr "Expected Frequencies"
}

 qui use "`imputef'",clear
 sort locus `vlist'
 rename Efreq freq
 tempvar total
 gen double `total'=sum(freq)
 gen double eprob=freq/`total'[_N]
 sort locus `vlist'
         
if "`display'"~="" {
   qui gen str80 `haptvar'=locus
   qui compress `haptvar'
   list `haptvar' `vlist' freq eprob, noobs
   qui drop `haptvar'
   * LOOK at the total frequency
   qui gen double totf = sum(freq)
   di 
   di "TOTAL FREQ is ", totf[_N]
}

global loglik=`llh'
global df =`df'

keep locus `vlist' freq eprob
if "`using'"~="" {

/* Must create the l1 l2 l3...  variables from locus for the profile likelihood */

  qui gen str40 plocus = locus
  forval i=1/$nloc {
      qui gen len = index(plocus,".")-1
      qui gen l`i'=real(substr(plocus,1,len))
      if `i'==$nloc { qui replace l`i'=real(plocus) }
      else { qui gen len2 = index(plocus,".")+1
      qui replace plocus=substr(plocus,len2,.) }
      drop len
      cap drop len2
  }
  drop plocus

  cap save "`using'", replace
}

/****************************************************************
 * Saving the model ipf string and the loglikelihood and
 * degrees of freedom
 ****************************************************************/

global hapmod`model'="`ipf'"
global hapdf`model'="`df'"
global hapllhd`model'="`llh'"

if "`lrtest'"~="" {
   tokenize "`lrtest'"
   local m1 "hapmod`1'"
   local m2 "hapmod`2'"
   local l1 "hapllhd`1'"
   local l2 "hapllhd`2'"
   local d1 "hapdf`1'"
   local d2 "hapdf`2'"

   /******************************************
    * Check whether these MACROS exist or not
    ******************************************/

   local error 0
   if "$`l2'" == "" {
      di in red "GLOBAL: `l2' does not exist is there a model `2'?"
      local error 1
   }
   if "$`l1'" == "" { 
      di in red "GLOBAL: `l1' does not exist is there a model `1'?"
      local error 1
    }
   if "$`m2'" == "" { 
      di in red "GLOBAL: `m2' does not exist is there a model `2'?"
      local error 1
    }
   if "$`m1'" == "" { 
      di in red "GLOBAL: `m1' does not exist is there a model `1'?" 
      local error 1
   }
   if "$`d2'" == "" { 
      di in red "GLOBAL: `d2' does not exist is there a model `2'?" 
      local error 1
   }
   if "$`d1'" == "" { 
      di in red "GLOBAL: `d1' does not exist is there a model `1'?" 
      local error 1
   }
   if `error'~=1 {
     di
     local len =length("Likelihood Ratio Test Comparing Model $`m2' to $`m1'")
     di as text "Likelihood Ratio Test Comparing Model $`m2' to $`m1'"
     di in smcl "{dup `len':{c -}}"
     di as res " llhd2 (df2)           =", $`l2', $`d2'
     di " llhd1 (df1)           =", $`l1', $`d1'
     di
     di "-2*(llhd2-llhd1)  =",-2*($`l2'-$`l1')
     di "Change in df      = ", $`d2'-$`d1'
     local lrt = -2*($`l2'-$`l1')
     local chdf = $`d2'-$`d1'
     local pv=chiprob(`chdf',`lrt')
     di "p-value           = ", chiprob(`chdf',`lrt')
     if `chdf'<0 {
       di as error "WARNING: negative chi-squared statistic"
       di as error "       The order of models in lrtest() is wrong"
     }                   
     else {
       if `pv'<0.05 { di as text "Accept Model $`m1' at 5% significance level" }
       else {  di as text "Accept Model $`m2' at 5% significance level" }
     }
     return scalar lrtpv = `pv'
     return scalar lrtdf = `chdf'
     return scalar lrtchi = `lrt'
   }
}


use "`origin'",clear
         
end

/*********************************************************
 * Tabulate for both locus the expected frequency
 * The 3rd and fourth variables contain the counting 
 * vectors
 *********************************************************/

program define tabduo
version 7.0
syntax	varlist(min=2 max=100), [BY(string)]
tokenize `varlist'

if "`by'"=="" {
	stack `1' `3' `2' `4', into(temp one) clear
	sort temp
	qui by temp: gen freq=sum(one)
	qui by temp: keep if _n==_N
}
if "`by'"~="" {
	stack `1' `3' `by' `2' `4' `by', into(temp one `by') clear
	sort temp `by'
	qui by temp `by': gen freq=sum(one)
	qui by temp `by': keep if _n==_N
}
end

/*********************************************************
 * Take the 2 haplotypes from each subject 
 * and do some basic tabulating saving the resulting 
 * tabulate in merge and merge2, these two represent 
 * the different sorted variables
 **********************************************************/

program define tabhap
version 7.0
syntax varlist(min=2 max=100)
tokenize `varlist'

cap confirm new file merge.dta
if _rc~=0 { di in red "merge.dta will be deleted within the program" }
cap confirm new file merge2.dta
if _rc~=0 { di in red "merge2.dta will be deleted within the program" }

stack `1' `2', into(haplo) clear
qui egen hapgrp = group(haplo)
sort haplo
qui by haplo: keep if _n==_N
keep haplo hapgrp
cap save merge,replace
sort hapgrp
cap save merge2,replace
end

/*************************************************************************
 * Calculate all the 2*($nloc-1) possible haplotypes 
 * given the 
 * observed data for each phase and person.
 * have locus1 and locus2 contain the string versions of the haplotypes
 * and the global unknown having the number of subjects with phase unknown
 **************************************************************************
 * 8/11/99 - The strings are out and separate variables are in
 **************************************************************************/

program define haplotyp
version 7.0
syntax [varlist] [,NOISE KNOWN PHASE(string) MV]
tokenize `varlist'

cap confirm new variable locus1 locus2 subject
if _rc~=0 {
  di in red "You must rename the following variable"
  confirm new variable locus1 locus2 subject
}
tempvar subj expand

* NOTE that `subj' contains the actual subject ids and at the end it gets renamed
* as subject

qui gen long subject = _n
qui gen long `subj' = _n

if "`mv'"~="" {
  di "EXPANDING MISSING DATA....."

  local li 1
  while `li'<=$nloc {
    local wc1=2*`li'-1
    local wc2=`wc1'+1
    local root1: word `wc1' of $all_name
    local root2: word `wc2' of $all_name

    qui count if `root1'==.
    local c1 = `r(N)'
    qui count if `root2'==.
    local cnt = `c1'+`r(N)'

    if `cnt'>0 {
      di "There are `cnt' missing values at locus `li'"
      qui tab `root1' , matrow(row)
      qui tab `root2' , matrow(col)

      /* NEED to make one matrix with all the values at locus i called unique */

      mat values = row \ col  
      local j 1
      while `j'<=rowsof(values) {
        if `j'==1 {
          mat unique = values[1,1] 
          local j =`j'+1
        }
        local not 0
        local jj 1
        while `jj'<`j' {
          if values[`j',1]==values[`jj',1] {
            local not 1
            local jj=`j'
          }
          local jj=`jj'+1
        }
        if `not'==0 { mat unique = unique \ values[`j',1] }

       local j = `j'+1
      }
      
    /* The matrix unique now contains just one copy of the alleles
     * Replace missing with one from each of unique
     * FOR TWO loci this would be the number of unique*(unique+1)/2
     */
      
      tempvar temp
      sort subject
      gen long `temp'=_n
      local missex = rowsof(unique)
      local missex = `missex'*(`missex'+1)/2
      qui expand  `missex' if `root1'==. | `root2'==.
      sort `temp'
      
      /*Create the new phenotypes */

      local lineno 1

      local rowsofunique = rowsof(unique)
      forval i=1/`rowsofunique' {
        forval j=1/`rowsofunique' {
           qui by `temp': replace `root1' = cond(`root1'==. & _n==`lineno', unique[`i',1],`root1')
           qui by `temp': replace `root2' = cond(`root2'==. & _n==`lineno', unique[`j',1],`root2')
           local lineno = `lineno'+1
         }
      }

 
    } /* end of if dealing with missing */


    local li =`li'+1  /*loop of loci */
  }

}
  di


if "`known'"~="" & "`phase'"=="" {
	qui gen str40 locus1=""
	qui gen str40 locus2=""
	/*******************************************
	 * Construct locus strings from alleles
	 *******************************************/
	sort subject
	global unknown = 0
	forval i=1/$nloc {
		local wc1=2*`i'-1
		local wc2=`wc1'+1
		local root1: word `wc1' of $all_name
		local root2: word `wc2' of $all_name

		if `i'==1 {
			qui replace locus1 =string(`root1')
			qui replace locus2 =string(`root2')
		}
		else {
			qui replace locus1 = locus1+"."+string(`root1')
			qui replace locus2 = locus2+"."+string(`root2')
		}

	}
}
if "`known'"=="" & "`phase'"=="" {

	local root1: word 1 of $all_name
	local root2: word 2 of $all_name
	qui gen str40 locus1=string(`root1')
	qui gen str40 locus2=string(`root2')
	qui gen `expand'=.

	/*******************************************
	 * Construct locus strings from alleles and 
	 * also expand to all possibilities
	 *******************************************/

        if "`mv'"=="" {	local i 2 }
	if "`mv'"~="" { local i 1 }

	while `i'<=$nloc {
          sort `subj'
          /* Gets the two loci names for locus i */
          local wc1=2*`i'-1
          local wc2=`wc1'+1
          local root1: word `wc1' of $all_name
          local root2: word `wc2' of $all_name
		

        /* Expand out heterozygotes */
                qui replace `expand'= 2*(`root1'~=`root2')

          /* old subject no is _n */
                  
		qui replace subject=_n
		qui expand `expand'
		sort subject

/* OLD WAY OF DOING MISSING DATA
                qui by subject: replace locus1= locus1+"."+string(cond(_n==2,`root2',`root1'))
               qui by subject: replace locus2= locus2+"."+string(cond(_n==2,`root1',`root2'))
*/

                if `i'>1 {
                  qui by subject: replace locus1= locus1+"."+string(cond(_n==2,`root2',`root1'))
		  qui by subject: replace locus2= locus2+"."+string(cond(_n==2,`root1',`root2'))
                }
                if `i'==1 {
                  qui by subject: replace locus1= string(cond(_n==2,`root2',`root1'))
		  qui by subject: replace locus2= string(cond(_n==2,`root1',`root2'))
                }

		
	local i=`i'+1
	}
	qui drop subject `expand'
	qui rename `subj' subject
}

if "`phase'"~="" {
	local root1: word 1 of $all_name
	local root2: word 2 of $all_name
	qui gen str40 locus1=string(`root1')
	qui gen str40 locus2=string(`root2')
	qui gen `expand'=.
	cap confirm variable `phase'
	if _rc~=0 {
		di in red "variable `phase' does not exist!"
		exit 111
	}

	/*******************************************
	 * Construct locus strings from alleles and 
	 * also expand to all possibilities
	 *******************************************/

        if "`mv'"=="" {	local i 2 }
	if "`mv'"~="" { local i 1 }

	while `i'<=$nloc {
		sort `subj'
		local wc1=2*`i'-1
		local wc2=`wc1'+1
		local root1: word `wc1' of $all_name
		local root2: word `wc2' of $all_name
		
                /* THIS BIT IS FOR NOT EXPANDING FOR WHOLE DATASET*/

	 qui replace `expand'= cond(`phase'==0,2*(`root1'~=`root2'),1)

		qui replace subject=_n
		qui expand `expand'
		sort subject

                if `i'>1 {
                  qui by subject: replace locus1= locus1+"."+string(cond(_n==2,`root2',`root1'))
		  qui by subject: replace locus2= locus2+"."+string(cond(_n==2,`root1',`root2'))
                }
                if `i'==1 {
                  qui by subject: replace locus1= string(cond(_n==2,`root2',`root1'))
		  qui by subject: replace locus2= string(cond(_n==2,`root1',`root2'))
                }
		
	local i=`i'+1
	}
	qui drop subject `expand'
	qui rename `subj' subject

}
qui compress

end

program define _stack
version 7.0
syntax varlist(min=1) [using/] [if],[NOISE START KNOWN PHASE(string) BY(string) ACC(real 0.001) DEBUG]
tokenize `varlist'

gen _stack=1
tempfile stack 
qui save "`stack'"

local i 1
while "``i''"~="" {
	local p=`i'
	local i = `i'+1
	if "``p''"~="``i''" {
		drop ``p'' 
		rename ``i'' ``p''
	}
	local i=`i'+1
}

qui replace _stack=2

qui append using `stack'


end

program define _ustack
version 7.0
syntax varlist(min=1) [using/] [if],BY(string) VALues(string) [NOISE START KNOWN PHASE(string) ]
tokenize `varlist'

tokenize "`values'", parse(" ")
tempfile first merge

qui by `by': gen merg=_n
save "`first'"

keep if `by'==`1'
keep `varlist' merg
foreach var of varlist `varlist' {
	rename `var' `var'1
}
sort merg
save "`merge'"

use "`first'"
tokenize `values', parse(" ")
keep if `by'==`2'
foreach var of varlist `varlist' {
	rename `var' `var'2
}

sort merg
merge merg using `merge'

end

/*******************************************
 * This programs unique function
 *******************************************/

program define _unique, rclass
   version 7.0
   syntax varlist(min=1 max=1)
   tokenize "`varlist'"
   local var `1'

preserve
sort `var'
qui by `var': keep if _n==1
return scalar N=_N
global S_1 =_N
qui drop if `var'==.
global S_2 =_N
return scalar unique=_N

restore

end

/******************************************
 * BUILDING the menu system for easy use
 ******************************************/

prog def _qmenu
version 7.0
syntax [varlist]

win c clear
local statnloc 10

/* Window width and height textheight */

local winw  210
local winh 240
local texth 7
local leftmargin 5
local rightmargin =`winw'-5

/* Initialising globals */

forval i=1/`statnloc' {
  global statloc`i' ""
}
global Qadd ""
global Qipf ""
global Qvar ""
global Qtvar ""
global infoqt ""
global infoqt1 ""
global Qcommand ""

/* Varlist options */
  
global Qsel "Select variables for:"
global Qmany "Loci"
global Qqt "Disease Outcome"
global Qvarlist "`varlist'"
window control static Qsel 5 10 110 `texth' center
window control static Qqt 60 20 50 `texth' center
window control ssimple Qvarlist 60 30 50 100 Qtvar
window control static Qmany 5 20 50 10 center
window control msimple Qvarlist 5 30 50 100 Qvar

/* BUTTONS */

local tey = `winh'-11-15
global Qvchk1 "hapupdate 0"
global Qvchk2 "hapupdate 1"
window control button "Apply" 5 `tey' 30 11 Qvchk1
window control button "Put in review window" 40 `tey' 70 11 Qvchk2

/* Get rid of window */

global Qrun "exit 3001"
global Qexit "exit 3000"
window control button "Run command" 115 `tey' 55 11 Qrun
window control button "Exit" 175 `tey' 30 11 Qexit

/* INFORMATION
 * Display info qdix and qdiy control the top left corner */

local qdix=`winw'-95
local qdiy 30
local qdisx = `qdix'+35
local qdiw = `winw'-5-`qdix'

global Qdi "Information"
window control static Qdi `qdix' 30 `qdiw' 100 blackframe
window control static Qdi `qdisx' 27 30 7 center

/* Put in pairs of vars */
  
tokenize $Qvar
forval i=1/`statnloc' {
  local tey = `qdiy'+10
  local tex = `qdix'+10
  local tew = `qdiw'-20
  local staty = `tey' + 7*(`i'-1)
  win c static statloc`i' `tex' `staty' `tew' `texth' left
}
local tey=`staty'+10
win c static infoqt1 `tex' `tey' `tew' `texth'
local tey=`tey'+`texth'
win c static infoqt `tex' `tey' `tew' `texth'

/* Display three models with check box and ipf() syntax box */

global Qllm "Loglinear Model"
local tey=`winh'-11-15-75
local teyy = `tey'-3
local tey2 = `tey'+9
local tew = `rightmargin'-`leftmargin'
local tew2= `tew'-20
local tex = `leftmargin'+10
local temid = int((`rightmargin'-`leftmargin')/2)-30
window control static Qllm `leftmargin' `tey' `tew' 20 blackframe
window control static Qllm `temid' `teyy' 60 8 center
window control static Qipf `tex' `tey2' `tew2' `texth' center

/* Display three models with check box and ipf() syntax box */

global Qmod1 ""
global Qmod2 ""
global Qmod3 ""
global Qmod4 ""

global Qqm "4 common models"
local tey=`winh'-11-15-50
local teyy = `tey'-3
local tew = `rightmargin'-`leftmargin'
local temid = int((`rightmargin'-`leftmargin')/2)-30
local tew2= `tew'-20
local tex = `leftmargin'+10
local tey2 = `tey'+9
window control static Qqm `leftmargin' `tey' `tew' 45 blackframe
window control static Qqm `temid' `teyy' 64 `texth' center
local tey3 = `tey'+7
local tex2 = `tex'+80
win c radbegin "Haplotype frequencies" `tex' `tey3' 75 `texth' Qrad
win c static Qmod1  `tex2' `tey3' 90 `texth'
local tey3 = `tey3'+9
win c radio "Linkage equilibrium" `tex' `tey3' 90 8 Qrad
win c static Qmod2  `tex2' `tey3' 90 `texth'
local tey3 = `tey3'+9
win c radio "Disease association" `tex' `tey3' 90 `texth' Qrad
win c static Qmod3  `tex2' `tey3' 90 `texth'
local tey3 = `tey3'+9
win c radend "No association"  `tex' `tey3' 90 `texth' Qrad
win c static Qmod4  `tex2' `tey3' 90 `texth'

/* Display the syntax */

global Qdis ""
local tey=`winh'-1-`texth'
local tex=`leftmargin'+31
window control static Qdis `leftmargin' `tey' 30 `texth' left
window control static Qcommand `tex' `tey' 200 `texth' left

/* ENd of windowing */

window dialog "A window to help with the syntax" . . `winw' `winh'

end