options ls=140 ps=60 pageno=1;

/** Model Weight Updating.SAS
**    Clint Moore, USGS
**    30 September 2011
**    28 August 2013 (fix for 64-bit OS)
**
** This program reads in FWS monitoring data and conducts a bootstrap and MCMC
**  sampling of vegetation transects.  These samples are used to build a
**  probability distribution for annual state transition and model prediction
**  likelihoods for model weight updating.
**
**  Although this program has been used by the U.S. Geological Survey (USGS),
**    no warranty, expressed or implied, is made by the USGS or the
**    U.S. Government as to the accuracy and functioning of the program and
**    related program material nor shall the fact of distribution constitute
**    any such warranty, and no responsibility is assumed by the USGS in
**    connection therewith.
**/

****** CONFIRM CURRENT MANAGEMENT YEAR, GRASS TYPE, AND FILE LOCATIONS! ********
********************************************************************************;
** Set value of current management year **;
  %let myear = 2013;

** Grass type (either m/M/mixed/MIXED or t/T/tall/TALL)  **;
  %let gselect = MIXED;
**  %let gselect = t;

** 64-bit settings for SAS PC Files Server **;
  %let server = CMoore-T1600;
  %let port = 9621;

** Settings for file locations (do not place an ending backslash) **;

  ** Directory containing Access database of monitoring and management data **;
  %let dbpath = C:\Northern Prairie Project\Database\Coordinator Database\NPAM_v8(final)_2013Aug27;

  ** Directory containing transition model parameters **;
  %let modelpath = C:\Northern Prairie Project\Models\Parameterizing\2013\Transition Model;

  ** Directory for storage of bootstrap and MCMC results (SAS database format) **;
  %let workpath = C:\Northern Prairie Project\Models\Updating\2013;

********************************************************************************
********************************************************************************;

  ** Do not alter the next 5 lines! -- they determine strings for file and table selection **;
      %let year1 = %eval(&myear-1);
      %let g0=mixed;
      %let g1=tall;
      %let suffix = %eval(%upcase(%sysfunc(substr(&gselect,1,1)))=T);
      %let gtype = %unquote(%nrstr(&g)&suffix);
  ** --------------------------------------------------------------------------------------**;

** Settings of model file names and Access table names **
  ** These names should not have to be changed if the table naming convention within Access remains unchanged  **;

  ** Name of Excel file containing transition model parameters (do not include .xls extension) **;
  %let modelsheet = trans_probs_&gtype;

  ** Name of Access database file **;
  %let dbfile = NPAM_v8.mdb;

  ** Name of Access table for current-year management data **;
  %let mgmttable = mgmt&myear._&gtype;

  ** Name of Access table for current-year monitoring data **;
  %let mon2table = monitor&myear._&gtype;

  ** Name of Access table for prior-year monitoring data **;
  %let mon1table = monitor&year1._&gtype;

  ** Name of Access table for current-year defoliation index data **;
  %let deftable = index&myear._&gtype;

  ** Name of Access table for current-year (input) model weights **;
  ** fields: ModelYear, MixedWt1,..., MixedWt4, TallWt1,..., TallWt6 **;
  %let mwttable = tblModelWeights;


  %let mwtout = &mwttable;

  ** Name of Access table for OUTPUT model weights **;
  **   Default action is to overwrite the INPUT table in Access;
  **   To create a NEW table in Access, un-comment the following line and provide a valid table name ;
  %let mwtout = tblModelWts2013;


** Parameters for bootstrap and MCMC samplers **
  ** Values are currently set for high performance and do not need to be changed **
  ** For these settings, samplers can take several minutes per unit -- be patient **;
  ** Random number seed **;
  %let rseed = 1234567;

  ** Bootstrap settings **;
    %let runboot = 1;         ** Turns on (1) or off (0) bootstrap sampler **;
    %let bootn = 100000;      ** Number of bootstrap samples -- may crash if bootn=1000K and number of transects>40;

  ** MCMC sampler parameters **;
    %let runmcmc = 1;         ** Turns on (1) or off (0) MCMC sampler **;
    %let fuzz = 1e-200;       ** Small threshold value for Dirichlet generator **;
    %let mcmcreps = 1000000;  ** Number of MCMC samples **;
    %let tunereps = 5000;     ** Number of burn-in samples **;



%let years=&year1&myear;

libname results "&workpath";

** Read data from Access tables **;
proc import out=work.modwts datatable="&mwttable" dbms=accesscs replace;
  database="&dbpath\&dbfile";  server="&server";  port=&port;
  run;
proc import out=work.mon2 datatable="&mon2table" dbms=accesscs replace;
  database="&dbpath\&dbfile";  server="&server";  port=&port;
  run;
proc import out=work.mon1 datatable="&mon1table" dbms=accesscs replace;
  database="&dbpath\&dbfile";  server="&server";  port=&port;
  run;
proc import out=work.indx datatable="&deftable" dbms=accesscs replace;
  database="&dbpath\&dbfile";  server="&server";  port=&port;
  run;
proc import out=work.mgmt datatable="&mgmttable" dbms=accesscs replace;
  database="&dbpath\&dbfile";  server="&server";  port=&port;
  run;

** Create formats for input and output **;
proc format;
  value invfmt 1='SB' 2='Co' 3='KB' 4='RM';
  value npfmt 1='0-30' 2='30-45' 3='45-60' 4='60-100';
  value defolfmt 1='<2' 2='2-4' 3='>4';
  value trtfmt 1="Rest" 2="Graze" 3="Burn" 4="B/G";
  value statefmt 1="{NP 0-30, SB}"
                 2="{NP 0-30, Co}"
                 3="{NP 0-30, KB}"
                 4="{NP 0-30, RM}"
                 5="{NP 30-45, SB}"
                 6="{NP 30-45, Co}"
                 7="{NP 30-45, KB}"
                 8="{NP 30-45, RM}"
                 9="{NP 45-60, SB}"
                10="{NP 45-60, Co}"
                11="{NP 45-60, KB}"
                12="{NP 45-60, RM}"
                13="{NP 60-100, SB}"
                14="{NP 60-100, Co}"
                15="{NP 60-100, KB}"
                16="{NP 60-100, RM}";
  invalue in_invfmt 'SB'=1 'Co'=2 'KB'=3 'RM'=4;
  invalue in_npfmt '0-30'=1 '30-45'=2 '45-60'=3 '60-100'=4;
  invalue in_defolfmt '<2'=1 '2-4'=2 '>4'=3;
  invalue inmtrtfmt "Rest"=1 "Graze"=2 "Burn"=3 "Burn/Graze"=4 "Unclassifiable"=.;
  invalue inttrtfmt "Rest"=1 "Grazew/inwindow"=2 "Burnw/inwindow"=3 "Defoliate"=4 "Unclassifiable"=.;
  invalue inmtrtfmtx "Rest"=1 "Graze"=2 "Burn"=3 "B/G"=4;
  invalue inttrtfmtx "Rest"=1 "GrazeW"=2 "BurnW"=3 "Defol"=4;
  run;

%macro readdata(y1,y2);

  ** Reads
  **   1 file of current-year (t) model weights for focal grassland type
  **     1 record: grassland, w_model1, w_model2,..., w_model6
  **   1 file of prior-year (t-1) monitoring data on units within grassland type:
  **     year, UID, transect, counts of SB, KB, NP, RM
  **   1 file of current-year (t) defoliation history:
  **     year, UID, defoliation history index (0-127) for years t-6, t-5,..., t
  **   1 file of current-year (t) monitoring data:
  **     year, UID, transect, counts of SB, KB, NP, RM
  **   1 file of current-year (t) implemented actions:
  **     year, UID, action taken
  **;

  %global priors;

  ** Read model prior weights and store in macro variable PRIORS **;
  data _null_;
    set modwts (where=(modelyear=&year1));
    w1 = &gtype.wt1;
    w2 = &gtype.wt2;
    w3 = &gtype.wt3;
    w4 = &gtype.wt4;
    %if &gtype=tall %then %do;
        w5 = tallwt5;
        w6 = tallwt6;
    %end;
    %else %do;
        w5 = 0;
        w6 = 0;
    %end;
    length priors $ 106;
    priors = put(w1,16.14) || ", " || put(w2,16.14) || ", " || put(w3,16.14) || ", " ||
             put(w4,16.14) || ", " || put(w5,16.14) || ", " || put(w6,16.14);
    call symput('priors',priors);
    run;

  ** Sets the appropriate informat for reading in models, based on grass type **;
  %if &gtype=mixed %then %let fmt = inmtrtfmt.;
  %if &gtype=tall  %then %let fmt = inttrtfmt.;

  ** Read prior-year monitoring data **;
  data mon1 (keep=year uid transect sb kb np rm);
    length uid $ 40 transect $ 8;
    set mon1 (rename=(year=charyear sb=ch_sb kb=ch_kb np=ch_np rm=ch_rm transect=tr) keep=year transect sb kb np rm complex org unit);
    year = input(charyear,4.);
    if vtype(year)="C" then year = input(charyear,4.);  else year = charyear;
    if vtype(ch_sb)="C" then sb = input(ch_sb,4.);  else sb = ch_sb;
    if vtype(ch_kb)="C" then kb = input(ch_kb,4.);  else kb = ch_kb;
    if vtype(ch_np)="C" then np = input(ch_np,4.);  else np = ch_np;
    if vtype(ch_rm)="C" then rm = input(ch_rm,4.);  else rm = ch_rm;
    transect = compress(tr);
    uid = substr(compress(complex),1,4) || " | " || substr(compress(org),1,4) || " | " || compress(unit);
    run;
  proc sort data=mon1;  by uid transect;  run;

  ** Read current-year defoliation histories **;
  data indx (keep=uid defindex ys);
    length uid $ 40;
    set indx (rename=(index=defindex) keep=index ys complex org unit);
    uid = substr(compress(complex),1,4) || " | " || substr(compress(org),1,4) || " | " || compress(unit);
    run;
  proc sort data=indx;  by uid;  run;

  ** Read current-year monitoring data **;
  data mon2 (keep=year uid transect sb kb np rm);
    length uid $ 40 transect $ 8;
    set mon2 (rename=(year=charyear sb=ch_sb kb=ch_kb np=ch_np rm=ch_rm transect=tr) keep=year transect sb kb np rm complex org unit);
    if vtype(year)="C" then year = input(charyear,4.);  else year = charyear;
    if vtype(ch_sb)="C" then sb = input(ch_sb,4.);  else sb = ch_sb;
    if vtype(ch_kb)="C" then kb = input(ch_kb,4.);  else kb = ch_kb;
    if vtype(ch_np)="C" then np = input(ch_np,4.);  else np = ch_np;
    if vtype(ch_rm)="C" then rm = input(ch_rm,4.);  else rm = ch_rm;
    transect = compress(tr);
    uid = substr(compress(complex),1,4) || " | " || substr(compress(org),1,4) || " | " || compress(unit);
    run;
  proc sort data=mon2;  by uid transect;  run;

  ** Read current-year implemented management actions **;
  data mgmt (keep=year uid treatment);
    length uid $ 40 treatment $ 16;
    set mgmt (rename=(year=charyear final_treatment_classification=trt) keep=year final_treatment_classification complex org unit);
    year = input(charyear,4.);
    treatment = compress(trt);
    uid = substr(compress(complex),1,4) || " | " || substr(compress(org),1,4) || " | " || compress(unit);
    run;
  proc sort data=mgmt;  by uid;  run;

  ** Find which units have monitoring data in both years.  List those that do not **;
  data mon1units;
    set mon1;
    by uid;
    if first.uid;
    run;

  data mon2units;
    set mon2;
    by uid;
    if first.uid;
    run;

  data monmatch (keep=uid) monnomatch (keep=uid missing);
    merge mon1units (in=one) mon2units (in=two);
    by uid;
    if one & two then output monmatch;
    if one & ^two then
      do;
        missing = 'year2';
        output monnomatch;
      end;
    if ^one & two then
      do;
        missing = 'year1';
        output monnomatch;
      end;
    run;

  proc print data=monnomatch;
    title1 "Units lacking monitoring data in a year";
    var uid missing;
    run;
  title1;

  ** Restrict prior-year monitoring data to those units with matches in current year **;
  data mon1;
    merge mon1 (in=one) monmatch (in=two);
    by uid;
    if one & two;
    run;

  ** Restrict current-year monitoring data to those units with matches in prior year **;
  data mon2;
    merge mon2 (in=one) monmatch (in=two);
    by uid;
    if one & two;
    run;

  ** Concatenate both years of monitoring data **;
  data transects (rename=(uid=unit));
    set mon1 mon2;
    run;

  ** Exclude management data to those units with monitoring data in both years.
  **   List those units that lack a management action and those units that do not
  **   have monitoring data **;
  data mgtmatch (keep=year uid treatment trt)
       nomondata (keep=year uid treatment)
       nomgtdata (keep=uid)
       notvalid (keep=year uid treatment);
    merge mgmt (in=one) monmatch (in=two);
    by uid;
    ** Treatment indicator **;
    trt = input(treatment,?? &fmt);
    if trt>. & one & two then output mgtmatch;
    if trt=. & one & two then output notvalid;
    if one & ^two then output nomondata;
    if ^one & two then output nomgtdata;
    run;

  proc print data=nomondata;
    title1 "Units with management data in current year but lacking before-after monitoring data";
    var year uid treatment;
    run;

  proc print data=nomgtdata;
    title1 "Units with before-after monitoring data, but no current-year management action";
    var uid;
    run;

  proc print data=notvalid;
    title1 "Units with a non-usable treatment type for updating";
    var uid treatment;
    run;
  title1;

  ** Merge current-year defoliation history with current-year management data.
  **   List those units that lack a defoliation history **;
  data treatments (keep=year uid treatment defindex trt def ys
                        rename=(uid=unit treatment=trtmnt))
       nodmatch
       notvalid;
    merge mgtmatch (in=one) indx (in=two);
    by uid;
    ** Defoliation class assignment **;
    def = 1*(defindex<2) + 2*(defindex>=2 & defindex<=4) + 3*(defindex>4);
    if one & two then output treatments;
    if one & ^two then output nodmatch;
    run;

  proc print data=nodmatch;
    title1 "Units with treatment and monitoring data but no defoliation history";
    var uid treatment;
    run;
  title1;

%mend;

%readdata(&year1,&myear)


%macro readmodels;

  %let mnames=;
  %let models=;
  %if &gtype=mixed %then
    %do;
      %let mnames=mixed_m1 mixed_m2 mixed_m3 mixed_m4;
      %let models = 4;
      %let fmt = inmtrtfmtx.;
    %end;
  %if &gtype=tall %then
    %do;
      %let mnames=tall_m1 tall_m2 tall_m3 tall_m4 tall_m5 tall_m6;
      %let models = 6;
      %let fmt = inttrtfmtx.;
    %end;

  proc datasets library=work nolist nowarn;
    delete model;
    run;

  %do i=1 %to &models;

    %let modelworksheet = %scan(&mnames,&i);
    ** Read competing models from Excel file **;
    proc import out=work.a
                datafile="&modelpath\&modelsheet..xls"
                dbms=excelcs replace;
      server="&server";  port=&port;
      sheet="&modelworksheet";
      scantext=yes;
      usedate=yes;
      scantime=yes;
      run;

    data a (drop=trtmnt dominv npclass defolx);
      set a;
      model = &i;
      trt = input(trtmnt,&fmt);
      inv0 = input(dominv,in_invfmt.);
      np0 = input(npclass,in_npfmt.);
      defidx = input(defolx,in_defolfmt.);
      v0 = 4*(np0-1)+inv0;
      run;

    proc datasets library=work nolist nowarn;
      append base=model data=a;
      run;

  %end;

  proc sort data=model;
    by model defidx trt np0 inv0;
    run;

%mend;

%readmodels


%macro prepdata;

  proc sort data=transects;
    by unit year transect;
    run;

  proc sort data=treatments;
    by unit year;
    run;

  ** Count up number of unit pairs and transects per unit, compute component proportions **;
  data transects (keep=unit transect year unityear sb kb np rm k tran p1-p4);
    set transects;
    by unit year;
    ** unityear is count of unit-year combinations, tran is count of transects per unit **;
    retain unityear 0 tran;
    if first.year then
      do;
        tran = 0;
        unityear+1;
        ** At end of dataset, the accumulated value of unityear is saved as a macro variable **;
        call symput('unityear',put(unityear,4.));
      end;
    tran+1;
    ** If any count is missing, convert it to 0 **;
    if sb=. then sb = 0;
    if kb=. then kb = 0;
    if np=. then np = 0;
    if rm=. then rm = 0;
    ** Add components to get total number of stops this year **;
    **   (transects designed with 50 stops each, but some years may have less) **;
    k = sb + kb + np + rm;
    ** Calculates component proportions p1 - p4 **;
    p1 = sb/k;
    p2 = kb/k;
    p3 = np/k;
    p4 = rm/k;
    run;

  ** Obtain simple stats on component proportions **;
  proc means data=transects mean noprint;
    by unit year;
    var p1-p4;
    output out=averages (keep=unit year pa1-pa4 pase1-pase4) mean=pa1-pa4 stderr=pase1-pase4;
    run;

  ** Merge vegetation data, treatment actions, and simple stats by unit-year **;
  data transects;
    merge transects treatments averages;
    by unit year;
    run;

  ** Need to re-sort by descending transect number for next step **;
  proc sort data=transects out=transect2;
    by unit year descending tran;
    run;

  ** Make a dataset of back-to-back pairs of units **;
  **   Keep values of years (year1-year2), number of transects (tran1-tran2),
  **   treatment applied in year 2 (trtmnt), defoliation index in year 2 (index),
  **   ID of pairing (pairs), simple stats for both years (avg1-avg8, se1-se8),
  **   unique label for the pairing (name), years since last defoliation (ys)  **;
  data pairs (keep=unit year1 year2 tran1 tran2 trtmnt name pairs
                   avg1-avg8 se1-se8 defindex trt def ys);
    length name $ 60 sameunit $ 60;
    set transect2;
    by unit year;
    retain sameunit " " year1 year2 tran1 tran2 avg1-avg8 se1-se8 pairs 0;
    array avg[2,4] avg1-avg8;
    array se[2,4] se1-se8;
    array pmean[4] pa1-pa4;
    array pse[4] pase1-pase4;
    ** If passing a unit boundary, initialize first year set of variables **;
    if first.unit then
      do;
        sameunit = unit;
        year1 = year;
        tran1 = tran;
        do i=1 to 4;
          avg[1,i] = pmean[i];
          se[1,i] = pse[i];
        end;
      end;
    ** If passing a year boundary within same unit and treatment is valid for use, initialize 2nd year set **;
    else if first.year & sameunit=unit & trt>. then
      do;
        year2 = year;
        tran2 = tran;
        do i=1 to 4;
          avg[2,i] = pmean[i];
          se[2,i] = pse[i];
        end;
        ** Make a number of pairs and put value in macro variable **;
        pairs+1;
        call symput('pairs',put(pairs,3.));
        ** Create a name for the pair for use in file labeling purposes **;
        name = compress(unit || "." || put(year1,4.) || "." || put(year2,4.));
        ** Output the information on this pair **;
        output;
        ** Initialize the current values as the next first year set of variables **;
        **  (only needed if processing more than 2 years of data at a time)  **;
        year1 = year;
        tran1 = tran;
        do i=1 to 4;
          avg[1,i] = pmean[i];
          se[1,i] = pse[i];
        end;
      end;
    run;

%mend;

%prepdata


** This macro performs bootstrap sampling of all unit pairs **;
%macro bootandmcmc;

  ** Clear out any previous results database to receive these results **;
  proc datasets library=work nolist nowarn;
    delete boot mcmc;
    run;

  %* Macro loop over the pairings **;
  %do i=1 %to &pairs;

    data _null_;
      set pairs;
      ** Use loop index to pull current pairing record **;
      if _n_=&i then
        do;
          ** Initialize macro variables from information in current record **;
          call symput('unit',trim(unit));
          call symput('y1',put(year1,4.));
          call symput('y2',put(year2,4.));
          call symput('t1',tran1);
          call symput('t2',tran2);
          call symput('trt',trtmnt);
        end;
      run;

    data pair (keep=unit trtmnt years tran year trial sb kb np rm k);
      ** Subset portion of vegetation file on basis of macro variables just initialized **;
      set transects (where=(unit="&unit" & (year=&y1 | year=&y2)))
                    end=lastone;
      by year;
      if year=&y1 then
        do;
          trial = 1;
          tran = &t1;
        end;
      else
        do;
          trial = 2;
          tran = &t2;
        end;
      trtmnt = "&trt";
      years = compress("&y1-&y2");
      run;

    proc iml;

      *** Procedure definitions ***;

      start stateassign(x);
        ** This procedure returns an assignment of vegetation state (1-16)
        **   from composition proportions of SB, KB, NP, and RM.
        ** Argument
        **  x: n x 4 vector of vegetation composition proportions
        ** Return
        **  y: n x 1 vector -- assignment of state **;
        n = nrow(x);
        r = j(n,4,0);
        np = 1 + (x[,3]>=0.3) + (x[,3]>=0.45) + (x[,3]>=0.6);
        if x[,3]<1.0 then r[,3] = x[,4] / (1 - x[,3]);
        else r[,3] = 0;
        if x[,1]+x[,2]>0 then
          do;
            r[,1] = x[,1] / (x[,1] + x[,2]);
            r[,2] = 1 - r[,1];
          end;
        else
          do;
            r[,1] = 0;
            r[,2] = 0;
          end;
        dom = 4 * (r[,3] >= 2/3) +
              1 * (2/3 >= r[,3]) # (r[,1] >= 2/3) +
              3 * (2/3 >= r[,3]) # (r[,2] >= 2/3) +
              2 * (2/3 >= r[,3]) # (2/3 >= r[,1]) # (2/3 >= r[,2]);
        state = 4*(np-1) + dom;
        return(state);
      finish;

      start bootstrap(x,bootn,rseed,phase,y,samples);
        ** This procedure performs bootstrap sampling of monitoring transects
        **   to compute average state assignments.
        ** Arguments
        **  x : n x 4 matrix of vegetation type counts
        **  bootn: number of bootstrap samples
        **  rseed: random number seed
        **  phase: indicator of which pair this procedure is processing
        ** Returns
        **  y : 1 x 24 vector of output
        **            pos 1-4: vegetation component bootstrap means
        **            pos 5-8: vegetation component bootstrap SDs
        **            pos 8-24: means of state assignments
        **  samples : bootn x 1 vector of state samples  **;
        t = x / (j(1,4,1) @ x[,+]);
        n = nrow(t);
        draws = int(ranuni(j(n,bootn,rseed))*n)+1;
        pavgsum = j(1,4,0);
        pavgsum2 = pavgsum;
        pavgsd = pavgsum;
        h = datetime();
        h0 = h;
        do i=1 to bootn;
          ** Progress monitor (updates log every 15 seconds) **;
          if datetime()-h0>15 then
            do;
              h0 = datetime();
              progress = i/bootn*100;
              remain = ((h0-h)/progress*100 - (h0-h)) / 60;
              put "Bootstrap step " phase 1.0 " of 2, " progress 3.0
                  "% complete, about " remain 4.1
                  " minutes remain for this phase (&unit, &y1-&y2)";
            end;
          ** s is a bootstrap sample of transects drawn with replacement **;
          s = t[draws[,i],];
          pavg = s[:,];
          pavgsum = pavgsum + pavg;
          pavgsum2 = pavgsum2 + (pavg # pavg);
          ** Determine state assignment from average of sample **;
          state = stateassign(pavg);
          samples[i,] = state;
        end;
        ** Summarize over all bootstrap samples **;
        pavgmean = pavgsum/bootn;
        if n>1 then pavgsd = sqrt( ( pavgsum2 - (pavgsum # pavgsum)/bootn ) / bootn );
        m = ((samples @ j(1,16,1)) = (j(bootn,1,1) @ (1:16)))[+,] / bootn;
        ** Return bootstrap samples of state assignment **;
        y = shape(pavgmean,1,4) || shape(pavgsd,1,4) || m;
        put "Bootstrap step " phase 1.0 " of 2 finished (&unit, &y1-&y2)" /;
      finish;

      start GenerateDirich(alpha);
        ** Generates a matrix of Dirichlet random variates of size K **;
        ** Argument
        **   alpha: J x K matrix of Dirichlet parameters
        ** Return
        **   g: J x K matrix of J Dirichlet draws, each of dimension K (each row sums to 1) ;
        g = j(nrow(alpha),ncol(alpha),1);
        **  call randgen(g,'gamma',alpha);  ** This form of call does not work **;
        do r=1 to nrow(alpha);
          do c=1 to ncol(alpha);
            call randgen(z,'gamma',alpha[r,c]);
            g[r,c] = max(z,&fuzz);
          end;
        end;
        g = g / (j(1,ncol(alpha),1) @ g[,+]);
        return(g);
      finish;

      start logFullCondAlphaK(p,alpha,k);
        ** Computes conditional (on kth alpha) log likelihood of Dirichlet
        ** Arguments
        **   p: J x 1 vector of kth-component proportions
        **   alpha: K x 1 vector of Dirichlet parameters
        **   k: scalar referencing focal Dirichlet parameter
        ** Return
        **   y: scalar  ;
        y = ( (alpha[k]-1) * (log(p))[+] ) - log(alpha[k]) +
            ( nrow(p) * (lgamma(alpha[+]) - lgamma(alpha[k])) );
        return(y);
      finish;

      start mcmc(x,mcmcreps,tunetime,rseed,phase,y,samples);
        ** This procedure performs MCMC sampling to compute posterior distribution
        **   of parameters of unit-level Dirichlet and posterior distribution of
        **   vegetation state assignment.
        **   (Based on Gauss code provided by Bill Link, PWRC)
        ** Arguments
        **  x : n x 4 matrix of vegetation type counts
        **  mcmcreps: number of MCMC samples
        **  tunereps: number of tuning samples
        **  rseed: random number seed
        **  phase: indicator of which pair this procedure is processing
        ** Returns
        **  y : 1 x 24 vector of output
        **            pos 1-4: vegetation component simulation means
        **            pos 5-8: vegetation component simulation SDs
        **            pos 8-24: means of state assignments
        **  samples : mcmcreps x 1 vector of state samples  **;
        Nmultis = nrow(x);
        Ncells = ncol(x);
        call randseed(&rseed);
        call randgen(rseed,'uniform');
        rseed = int(rseed*10e8);
        allalpha = j(mcmcreps+tunetime,Ncells,0);
        allalpha[1,] = j(1,Ncells,1)/2;
        maxstepsize = allalpha[1,]/2;
        h = datetime();
        h0 = h;
        i = 1;
        do while(i<mcmcreps + tunetime);
          ** Progress monitor (updates log every 15 seconds) **;
          if datetime()-h0>15 then
            do;
              h0 = datetime();
              progress = i/(mcmcreps+tunetime)*100;
              remain = ((h0-h)/progress*100 - (h0-h)) / 60;
              put "MCMC step " phase 1.0 " of 2, " progress 3.0
                  "% complete, about " remain 4.1
                  " minutes remain for this phase (&unit, &y1-&y2)";
            end;
          i = i+1;
          alphacurr = allalpha[i-1,];
          ** Gibbs sampling of p vectors **;
          pcurr = GenerateDirich( (alphacurr @ j(Nmultis,1,1)) + x );
          ** Gibbs sampling of alphas **;
          k = 0;
          do while(k<Ncells);
            k = k+1;
            alphacand = alphacurr;
            alphacand[k] = alphacand[k] + (2*uniform(rseed)-1) * maxstepsize[k];
            if alphacand[k]<=0 then
              do;
                if i<tunetime then maxstepsize[k] = maxstepsize[k] * 0.98;
              end;
            else
              do;
                logr = logFullCondAlphaK(pcurr[,k],alphacand,k) -
                       logFullCondAlphaK(pcurr[,k],alphacurr,k);
                if log(uniform(rseed))<logr then
                  do;
                    alphacurr = alphacand;
                    if i<tunetime then maxstepsize[k] = maxstepsize[k] * 1.025;
                  end;
                else if i<tunetime then maxstepsize[k] = maxstepsize[k] * 0.98;
              end;
          end;
          allalpha[i,] = alphacurr;
        end;
        ** Remove burn-in samples **;
        allalpha = allalpha[tunetime+1:mcmcreps+tunetime,];
        ** Posterior samples of component proportions **;
        p = allalpha / (j(1,Ncells,1) @ allalpha[,+]);
        ** Posterior samples of state assignments **;
        samples = stateassign(p);
        pavgsum = p[+,];
        pavgsum2 = p[##,];
        pavgmean = pavgsum/mcmcreps;
        pavgsd = sqrt( ( pavgsum2 - (pavgsum # pavgsum)/mcmcreps ) / mcmcreps );
        m = ((samples @ j(1,16,1)) = (j(mcmcreps,1,1) @ (1:16)))[+,] / mcmcreps;
        ** Return MCMC posterior samples of state assignment **;
        y = shape(pavgmean,1,4) || shape(pavgsd,1,4) || m;
        put "MCMC step " phase 1.0 " of 2 finished (&unit, &y1-&y2)" /;
      finish;


      ** Read data and create input matrices **;
      use pair;
      read point 1 var{unit trtmnt years} into unitinfo;
      read all var{sb kb np rm k} where(trial=1) into trial1;
      read all var{sb kb np rm k} where(trial=2) into trial2;

      ** Bootstrap execution **;
      %if &runboot=1 %then
        %do;
          file log;
          m1 = j(1,24,0);
          s1 = j(&bootn,1,0);
          m2 = m1;
          s2 = s1;
          ** Call bootstrap procedure for each year **;
          call bootstrap(trial1[,1:4],&bootn,&rseed,1,m1,s1);
          call bootstrap(trial2[,1:4],&bootn,&rseed,2,m2,s2);
          ** Convert pair of state values (1-16) to 1-256 combined index vector of state pairs **;
          mcomb = (s1-1)*16 + s2;
          free s1 s2;

          ** Compute state pair assignment probabilities from sample frequencies **;
          summ = j(1,256,0);
          test = j(&bootn,1,1);
          do i=1 to 256;
            summ[i] = (mcomb = (test*i))[+] / &bootn;
          end;
          free mcomb test;
          x = (&i) || m1[1:4]` || m2[1:4]` || m1[5:8]` || m2[5:8]` || summ;

          ** Output results to dataset **;
          names = 'pairs' || ('bootm1':'bootm8') || ('bootsd1':'bootsd8') || ('bootx1':'bootx256');
          create bootadd var names;
          append from x;
        %end;

      ** MCMC execution **;
      %if &runmcmc=1 %then
        %do;
          file log;
          m1 = j(1,24,0);
          s1 = j(&mcmcreps,1,0);
          m2 = m1;
          s2 = s1;
          ** Call MCMC procedure for each year **;
          call mcmc(trial1[,1:4],&mcmcreps,&tunereps,&rseed,1,m1,s1);
          call mcmc(trial2[,1:4],&mcmcreps,&tunereps,&rseed,2,m2,s2);
          ** Convert pair of state values (1-16) to 1-256 combined index vector of state pairs **;
          mcomb = (s1-1)*16 + s2;
          free s1 s2;

          ** Compute state pair assignment probabilities from sample frequencies **;
          summ = j(1,256,0);
          test = j(&mcmcreps,1,1);
          do i=1 to 256;
            summ[i] = (mcomb = (test*i))[+] / &mcmcreps;
          end;
          free mcomb test;
          x = (&i) || m1[1:4]` || m2[1:4]` || m1[5:8]` || m2[5:8]` || summ;

          ** Output results to dataset **;
          names = 'pairs' || ('mcmcm1':'mcmcm8') || ('mcmcsd1':'mcmcsd8') || ('mcmcx1':'mcmcx256');
          create mcmcadd var names;
          append from x;
        %end;

      quit;

    ** Append the results for current pair to results for previous pairs **;
    proc datasets library=work nolist nowarn;
      %if &runboot=1 %then %str(append base=boot data=bootadd;);
      %if &runmcmc=1 %then %str(append base=mcmc data=mcmcadd;);
      run;

  %end;

  ** Write temporary datasets to permanent SAS databases **;
  %if &runboot=1 %then %str(data results.&gtype.&years.boot;  set boot;  run;);
  %if &runmcmc=1 %then %str(data results.&gtype.&years.mcmc;  set mcmc;  run;);

%mend;

%bootandmcmc

** This macro reads MCMC and bootstrap output and processes the results **;
%macro read;

  %if &gtype=mixed %then %do;
    %let models = 4;
    %let keepstring = MixedWt1 MixedWt2 MixedWt3 MixedWt4;
  %end;
  %if &gtype=tall %then %do;
    %let models = 6;
    %let keepstring = TallWt1 TallWt2 TallWt3 TallWt4 TallWt5 TallWt6;
  %end;

  ** Only run these statements if both bootstrap and MCMC are requested **;
  %if &runboot=1 and &runmcmc=1 %then %do;

    ** Merge BOOT and MCMC output with pair names and simple statistics **;
    data results (drop=i sbkb other npstate domstate r_rm r_sb r_kb);
      merge pairs results.&gtype.&years.boot results.&gtype.&years.mcmc;
      by pairs;
      array p[2,4] avg1-avg8;
      array s[2] state1-state2;
      ** Compute deterministic state assignment from simple averages **;
      do i=1 to 2;
        sbkb = p[i,1] + p[i,2];  ** Sum of SB & KB only **;
        other = sbkb + p[i,4];   ** Sum of all non-NP components **;
        npstate = 1 + (p[i,3]>=0.3) + (p[i,3]>=0.45) + (p[i,3]>=0.6);
        if other>0 then r_rm = p[i,4] / other;
        else r_rm = 0;
        if sbkb>0 then
          do;
            r_sb = p[i,1] / sbkb;
            r_kb = 1 - r_sb;
          end;
        else
          do;
            r_sb = 0;
            r_kb = 0;
          end;
        domstate = 4 * (r_rm>=2/3) +
                   1 * (r_rm<2/3) * (r_sb>=2/3) +
                   3 * (r_rm<2/3) * (r_kb>=2/3) +
                   2 * (r_rm<2/3) * (r_sb<2/3) * (r_kb<2/3);
        s[i] = 4 * (npstate-1) + domstate;
      end;
      run;

    proc sort data=results;  by pairs name trtmnt def ys tran1 state1 tran2 state2;  run;

    ** Compute likelihoods and posterior probabilities for simple, bootstrap, and Bayesian approaches **;
    **  (Note use of _priors_ macro variable in array below) **;
    data results (drop=i def0);
      set results;
      array z[3,6] lik1_1-lik1_6 lik2_1-lik2_6 lik3_1-lik3_6;
      ** These priors stay fixed over successive units **;
      array pre_fixed[3,6] _temporary_ (&priors, &priors, &priors);
      array post_fixed[3,6] postf1_1-postf1_6 postf2_1-postf2_6 postf3_1-postf3_6;
      ** These priors change over successive units **;
      array pre_dynamic[3,6] _temporary_ (&priors, &priors, &priors);
      array post_dynamic[3,6] postd1_1-postd1_6 postd2_1-postd2_6 postd3_1-postd3_6;
      ** These variables will be read in from model file **;
      array v[16] np1sb np1co np1kb np1rm np2sb np2co np2kb np2rm
                  np3sb np3co np3kb np3rm np4sb np4co np4kb np4rm;
      ** These are posterior distributions on state pair membership **;
      array b[16,16] bootx1-bootx256;
      array m[16,16] mcmcx1-mcmcx256;
      sum1f = 0;  sum1d = 0;
      sum2f = 0;  sum2d = 0;
      sum3f = 0;  sum3d = 0;
      ** Loop over models **;
      do i=1 to &models;
        ** For model 1 (cyclic action model), read from defoliation state 2 if ys<5 **;
        def0 = def;
        if i=1 & ys<5 then def0 = 2;
        ** Compute look-up row in model set based on model, defoliation, and treatment indices
        **  and deterministic assignment of state1 **;
        lookup = (3*4*16)*(i-1) + (4*16)*(def0-1) + 16*(trt-1) + state1;
        ** Direct-read model file to obtain likelihood **;
        set model point=lookup;
        ** Likelihood for deterministic state assignment **;
        z[1,i] = v[state2];
        ** Compute expected likelihoods for state assignments from bootstrap and MCMC analyses **;
        z[2,i] = 0;
        z[3,i] = 0;
        do j=1 to 16;
          ** This pulls out a row (starting state);
          lookup = (3*4*16)*(i-1) + (4*16)*(def0-1) + 16*(trt-1) + j;
          set model point=lookup;
          do k=1 to 16;
            ** This loops over the columns (ending state);
            z[2,i] = z[2,i] + v[k]*b[j,k];
            z[3,i] = z[3,i] + v[k]*m[j,k];
          end;
        end;
        ** Increment accumulators **;
        sum1f = sum1f + z[1,i]*pre_fixed[1,i];
        sum2f = sum2f + z[2,i]*pre_fixed[2,i];
        sum3f = sum3f + z[3,i]*pre_fixed[3,i];
        sum1d = sum1d + z[1,i]*pre_dynamic[1,i];
        sum2d = sum2d + z[2,i]*pre_dynamic[2,i];
        sum3d = sum3d + z[3,i]*pre_dynamic[3,i];
      end;
      ** Compute posterior model weights **;
      do i=1 to &models;
        if sum1f>0 then post_fixed[1,i] = z[1,i]*pre_fixed[1,i]/sum1f;
        if sum2f>0 then post_fixed[2,i] = z[2,i]*pre_fixed[2,i]/sum2f;
        if sum3f>0 then post_fixed[3,i] = z[3,i]*pre_fixed[3,i]/sum3f;
        if sum1d>0 then post_dynamic[1,i] = z[1,i]*pre_dynamic[1,i]/sum1d;
        if sum2d>0 then post_dynamic[2,i] = z[2,i]*pre_dynamic[2,i]/sum2d;
        if sum3d>0 then post_dynamic[3,i] = z[3,i]*pre_dynamic[3,i]/sum3d;
        ** Set current posterior weights as priors for next management unit (dynamic case) **;
        pre_dynamic[1,i] = post_dynamic[1,i];
        pre_dynamic[2,i] = post_dynamic[2,i];
        pre_dynamic[3,i] = post_dynamic[3,i];
      end;
      run;

    ** Write out updated model weights based on MCMC likelihoods and successive updating **;
    data wt_out1;
      retain modelyear;
      set results (keep=postd3_1-postd3_6
                   rename=(postd3_1=&gtype.wt1
                           postd3_2=&gtype.wt2
                           postd3_3=&gtype.wt3
                           postd3_4=&gtype.wt4
                           postd3_5=tallwt5
                           postd3_6=tallwt6)) end=lastone;
      modelyear = &myear;
      keep modelyear &keepstring;
      if lastone then output;
      run;

    proc print data=results;
      title1 "Model weight updating for MCMC analysis (&mcmcreps samples)";
      title2 "Starting model priors = (&priors)";
      var name trtmnt def ys tran1 state1 tran2 state2 lik3_1-lik3_&models postd3_1-postd3_&models postf3_1-postf3_&models;
      format def defolfmt. state1 state2 statefmt.;
      run;

    ** These steps compute median likelihoods (from MCMC) over units,
    **   then updates model weights based on medians **;
    proc means data=results mean median;
      var lik3_1-lik3_6 postf3_1-postf3_6;
      output out=medliks mean(lik3_1-lik3_6)=mean1-mean6
                         median(lik3_1-lik3_6)=med1-med6;
      run;
    title1;

    ** Write out updated model weights based on median of MCMC likelihoods **;
    data wt_out2;
      set medliks;
      modelyear = &myear;
      array z[6] med1-med6;
      array pre[6] _temporary_ (&priors);
      array post[6] &gtype.wt1 &gtype.wt2 &gtype.wt3 &gtype.wt4 tallwt5 tallwt6;
      sum = 0;
      do i=1 to 6;
        if z[i]>0 then sum = sum + z[i]*pre[i];
      end;
      do i=1 to 6;
        post[i] = z[i]*pre[i]/sum;
      end;
      keep modelyear &keepstring;
      run;

    proc print data=wt_out1;
      title1 "Model weight updates for &gtype &year1-&myear based on MCMC successive-unit updating";
      run;

    proc print data=wt_out2;
      title1 "Model weight updates for &gtype &year1-&myear based on MCMC median-likelihood updating";
      run;

    title1;

    data modwtsout;
      merge modwts wt_out2;
      by modelyear;
      run;

    proc print data=modwtsout;  run;

    proc export data=modwtsout outtable="&mwtout" dbms=accesscs replace;
      database="&dbpath\&dbfile";  server="&server";  port=&port;
      run;

  %end;

%mend;

%read

quit;