R/toLatex_marssMODEL.R

toLatex.marssMODEL = function(object, ..., file=NULL, digits=2, greek=TRUE, orientation="landscape", math.sty="amsmath", output=c("pdf","tex"),replace=TRUE,simplify=TRUE){
  #by default uses geometry package and amsmath
  require(Hmisc)
  require(stringr)
  x=object
  #set the environment of the subfunction to this function since I want them to have access to the passed in variables
  #environment(build.eqn.tex)=environment()
  #environment(get.mat.tex)=environment()
  #environment(parameters.time.varying)=environment()
  
  #set up the headers and footers
  eqn.head="\\begin{gather*}"
  eqn.foot="\\end{gather*}"
  mat.head="\\begin{bmatrix}"
  mat.foot="\\end{bmatrix}"
  
  marssMODEL=x #used below
  #get the model from x (x is a marssMODEL object)
  fixed=x$fixed
  free=x$free
  #btw, this is not robust to name changes
  m=dim(fixed$x0)[1]; n=dim(x$data)[1]; TT=dim(x$data)[2]
  model.dims=attr(x,"model.dims")
  model.el=attr(x,"par.names")
  
  #read in the equation from the attributes.  This will be used to construct the tex
  #needs to look like x=A*x+B*C*D+E+..+w, w~MVN(E+F+..,G+H+...)\n next eqn
  eqns=attr(x,"equation")
  eqns=strsplit(eqns,"\n")[[1]]
  
  #Build a vector of all the special names in all the equations; lhs of = and ~
  #Do some testing too and determine if time will be expanded
  eqn.special=c()
  # the default behavior is no time expansion
  eqn.ts = 1
  time.expand = FALSE     
  for(i in 1:length(eqns)){
    eqn=eqns[i]
    if(!str_detect(eqn,"=")) stop("toLatex.marssMODEL: the atttibute, equation, must have the form x=... or y=....")
    eqn.name=str_trim(strsplit(eqn,"=")[[1]][1])
    #the part to the left of , is the deterministic part + random bit in last term
    #the eqn for the errors is after the ,
    eqn.stoc=str_trim(strsplit(eqn,";")[[1]][2])
    #the name is before the ~
    err.name=str_trim(strsplit(eqn.stoc,"~")[[1]][1])
    #the eqn.name (like x or y) and the err.name (like w or v) are treated differently
    eqn.special=c(eqn.special,str_trim(strsplit(eqn.name,"_")[[1]][1]),str_trim(strsplit(err.name,"_")[[1]][1]))

    if(!str_detect(eqn,"=")) stop("toLatex.marssMODEL: the atttibute, equation, must have the form x=... or y=....")
    eqn.name=str_trim(strsplit(eqn,"=")[[1]][1])
    #the part to the left of , is the deterministic part + random bit in last term
    eqn.det=str_trim(strsplit(eqn,";")[[1]][1])
    #the eqn for the errors is after the ,
    eqn.stoc=str_trim(strsplit(eqn,";")[[1]][2])
    #the name is before the ~
    err.name=str_trim(strsplit(eqn.stoc,"~")[[1]][1])
    #check that the name of the errors matches before and after the ,
    tmp=strsplit(eqn.det,"[+]")[[1]] #tmp[length(tmp)] should be error part
    tmp=strsplit(tmp[length(tmp)],"[*]")[[1]] #tmp[length(tmp)] should be err.name
    if(err.name != str_trim(tmp[length(tmp)])) stop("toLatex.marssMODEL: in the attribute, equation, the name of the errors after the , does not match the name before the ,")
    #these 2 lines get the bit between ( and ,
    eqn.stoc.mean=strsplit(eqn.stoc,",")[[1]][1]
    eqn.stoc.mean=str_trim(strsplit(eqn.stoc.mean,"[(]")[[1]][2])
    #these 2 lines get the bit between , and )
    eqn.stoc.var=strsplit(eqn.stoc,"[)]")[[1]][1]
    eqn.stoc.var=str_trim(strsplit(eqn.stoc.var,",")[[1]][2])
    #If any of the equation is time-varying then will have 1 equation for each time
    #Determine if that is necessary
    if(!simplify & (parameters.time.varying(eqn.det, eqn.special, marssMODEL) |
      parameters.time.varying(eqn.stoc.mean, eqn.special, marssMODEL) |
      parameters.time.varying(eqn.stoc.var, eqn.special, marssMODEL)) ){
      eqn.ts=TT
      time.expand=TRUE
    }
  }
  
  all.t.eqn.tex=c() #this will be all the equations
  for(t in 1:eqn.ts){ #one eqn for each time step; default is to simplify and not time expand
  eqn.tex=c() #this is just for time=t    
  #eqns is from the equation attribute of marssMODEL
  for(i in 1:length(eqns)){
    eqn=eqns[i]
    ##############################################################
    # First read in and parse the equation attribute of marssMODEL
    ##############################################################
    eqn.name=str_trim(strsplit(eqn,"=")[[1]][1])
    #the part to the left of , is the deterministic part + random bit in last term
    eqn.det=str_trim(strsplit(eqn,";")[[1]][1])
    #the eqn for the errors is after the ,
    eqn.stoc=str_trim(strsplit(eqn,";")[[1]][2])
    #the name is before the ~
    err.name=str_trim(strsplit(eqn.stoc,"~")[[1]][1])
    #check that the name of the errors matches before and after the ,
    tmp=strsplit(eqn.det,"[+]")[[1]] #tmp[length(tmp)] should be error part
    tmp=strsplit(tmp[length(tmp)],"[*]")[[1]] #tmp[length(tmp)] should be err.name
    #these 2 lines get the bit between ( and ,
    eqn.stoc.mean=strsplit(eqn.stoc,",")[[1]][1]
    eqn.stoc.mean=str_trim(strsplit(eqn.stoc.mean,"[(]")[[1]][2])
    #these 2 lines get the bit between , and )
    eqn.stoc.var=strsplit(eqn.stoc,"[)]")[[1]][1]
    eqn.stoc.var=str_trim(strsplit(eqn.stoc.var,",")[[1]][2])
    
    #Set up and test the time designation for the equation term 
    if(!str_detect(eqn.name,"_")) stop(paste("toLatex.marssMODEL: in the attribute, equation, ", eqn.name, " is missing the time designation.  should look like _{t}.",sep=""))
    eqn.name.time=str_trim(strsplit(eqn.name,"_")[[1]][2])
    eqn.name=str_trim(strsplit(eqn.name,"_")[[1]][1]) 
    eqn.labs=attr(x,paste(toupper(eqn.name),".names",sep=""))
    #get rid of the bit before the =
    eqn.det=str_trim(strsplit(eqn.det,"=")[[1]][2])
    
    #Set up and test the time designation for error term
    if(!str_detect(err.name,"_")) stop(paste("toLatex.marssMODEL: in the attribute, equation, ", err.name, " is missing the time designation.  should look like _{t}.",sep=""))
    err.name.time=str_trim(strsplit(err.name,"_")[[1]][2])
    err.name=str_trim(strsplit(err.name,"_")[[1]][1])    
    
    
    #############################################################
    # Build the latex for equation
    #############################################################
      #########################################
      # start with the eqn name (left of =) 
      #########################################
      if(!time.expand) time.sub=eqn.name.time
      else time.sub=paste("{",as.character(t),"}",sep="")
      mat.body=paste(paste(eqn.name,"_{",eqn.labs,"}",sep=""),collapse="\\\\\n")
      mat.tex=paste(mat.head,"\n",mat.body,"\n",mat.foot,"_",time.sub,"=",sep="")
      eqn.tex=c(eqn.tex,mat.tex)
      headfoot = list(eqn.name=eqn.name, mat.head=mat.head, mat.foot=mat.foot)
    
      #############################################################
      # send the deterministic part to right of = as char string like Bx_t+U 
      # and build the latex for the equation
      ##############################################################
      if(time.expand){
        eqn.tex=c(eqn.tex, build.eqn.tex(eqn.det, eqn.special, marssMODEL, greek, digits, simplify, headfoot, t=t)) 
      }else{
        eqn.tex=c(eqn.tex, build.eqn.tex(eqn.det, eqn.special, marssMODEL, greek, digits, simplify, headfoot)) 
      }
      #put the stochastic part on it's own line
      eqn.tex=c(eqn.tex,"\\\\")
      
      ##############################################################
      # build the stochastic part of the equation
      ##############################################################
      #start with the eqn name (left of ~) 
      if(!time.expand) time.sub=err.name.time
      else time.sub=paste("{",as.character(t),"}",sep="")
      mat.body=paste(paste(err.name,"_{",eqn.labs,"}",sep=""),collapse="\\\\\n")
      mat.tex=paste(mat.head,"\n",mat.body,"\n",mat.foot,"_",time.sub," \\sim",sep="")
      eqn.tex=c(eqn.tex,mat.tex)
      #add the distribution designation. will be in front of (
      tmp=strsplit(eqn.stoc,"[(]")[[1]][1]; tmp=str_trim(strsplit(tmp,"~")[[1]][2])
      #wrap in pmatrix
      eqn.tex=c(eqn.tex,tmp,"\\begin{pmatrix}")
      
      #add the mean
      if(!time.expand)
        eqn.tex=c(eqn.tex,build.eqn.tex(eqn.stoc.mean, eqn.special, marssMODEL, greek, digits,simplify, headfoot))
      else
        eqn.tex=c(eqn.tex,build.eqn.tex(eqn.stoc.mean, eqn.special, marssMODEL, greek, digits,simplify, headfoot, t=t))
      
      #separate from variance with a ,
      eqn.tex=c(eqn.tex,", ")
      #add the variance
      if(!time.expand)
        eqn.tex=c(eqn.tex,build.eqn.tex(eqn.stoc.var, eqn.special, marssMODEL, greek, digits,simplify, headfoot))
      else
        eqn.tex=c(eqn.tex,build.eqn.tex(eqn.stoc.var, eqn.special, marssMODEL, greek, digits,simplify, headfoot, t=t))
      
      #close
      eqn.tex=c(eqn.tex,"\\end{pmatrix}")
      
      #end of equation i
      eqn.tex=c(eqn.tex,"\\\\\n \\\\")
  } #end of equation (like x=..., w ~ MVN(...)
  #get rid of trailing two \\\\ after the last eqn
  eqn.tex=eqn.tex[-1*length(eqn.tex)]
  #add the header and footer for the equation
  all.t.eqn.tex=c(all.t.eqn.tex,eqn.head,eqn.tex,eqn.foot)
  } #over all t for which we should print
  
  
  ##############################################################
  # Create the latex file
  ##############################################################
  if (length(math.sty)) 
    sty = paste("\\usepackage{", math.sty, "}", sep = "")
  sty = c(sty, paste("\\usepackage[",orientation,"]{geometry}", sep = ""))
  if("Rnw" %in% output)
    sty = c(sty, "\\usepackage{Sweave}")
  
  #Set up the file name
  if(missing(file) & replace) tmp="marssMODEL"
  if(missing(file) & !replace) tmp=tempfile("marssMODEL",tmpdir=".")
  if(!missing(file)) tmp=file
  #If the user put on a extension of .tex, .dvi or .pdf remove that
  if(tolower(str_sub(tmp,-4))==".tex") tmp=str_sub(tmp,-4)
  if(tolower(str_sub(tmp,-4))==".pdf") tmp=str_sub(tmp,-4)
  if(tolower(str_sub(tmp,-4))==".dvi") tmp=str_sub(tmp,-4)
  
  tmptex.filename <- paste(tmp, "tex", sep = ".")
  cat("\\documentclass{report}", sty, "\\begin{document}\\pagestyle{empty}",
      all.t.eqn.tex, "\\end{document}\n", file = tmptex.filename, sep = "\n")
  
  if("pdf" %in% output) latexcmd="pdflatex"
  if("dvi" %in% output) latexcmd="latex"
  if("pdf" %in% output | "dvi" %in% output){
    #sys(paste(latexcmd,  "-interaction=scrollmode", shQuote(tmp)), output = FALSE)
    sys(paste(latexcmd, shQuote(tmp)), output = FALSE)
    sys(paste("rm", paste(str_replace(tmp,"[\\]","/"),"aux",sep=".")))
    sys(paste("rm", paste(str_replace(tmp,"[\\]","/"),"log",sep=".")))
  }
  if("Rnw" %in% output)
    sys(paste("cp", paste(str_replace(tmp,"[\\]","/"),"tex",sep="."), paste(str_replace(tmp,"[\\]","/"),"Rnw",sep=".")))
  if(!("tex" %in% output))
    sys(paste("rm", paste(str_replace(tmp,"[\\]","/"),"tex",sep=".")))
} #end of toLatex.marssMODEL

get.mat.tex=function(x,greek,digits){
  #called by build.eqn.tex
  #x is a matrix for the parameter
  #use Hmisc util to convert to tabular form
  mat=array(list(),dim=dim(x))
  tmp.fun=function(x,digits=2,greek=greek){
    ifelse(is.numeric(x[[1]]),round(x[[1]],digits=digits),latexTranslate(x[[1]],greek=greek))
  }
  #this will be a list matrix whether or not x is a list matrix
  mat[,]=lapply(x,tmp.fun,digits=digits,greek=greek)
  #now deal with any +,*,-1* in the matrices
  tmp.fun.2=function(x){
    if(is.character(x[[1]])){
      x=str_replace_all(x,"[+][-]","-")
      x=str_replace_all(x,"[-]1[*]","-")
      x=str_replace_all(x,"[+]1[*]","+")
      x=str_replace_all(x,"[*]","$\\\times$")
      x=str_replace_all(x,"[+]","$+$")
      x=str_replace_all(x,"[-]","$-$")
      #get rid of all the $ because the matrix will be wrapped in equation
      x = str_replace_all(x,"[$]"," ") 
    }else{ x[[1]] }
  }
  #mat is still a list matrix
  mat[,]=lapply(mat,tmp.fun.2)
  
  #now use the Hmisc util to translate to tabular form
  tmp.tex=latexTabular(mat)[[1]]
  #latexTabular is replacing _ with \\\\_
  tmp.tex=str_replace_all(tmp.tex,"[\\][\\]_","_")
  tmp.tex=strsplit(tmp.tex,"\n")[[1]]
  tmp.tex=tmp.tex[-1]
  tmp.tex=tmp.tex[-1*length(tmp.tex)]
  tmp.tex=tmp.tex[-1*length(tmp.tex)]
  tmp.tex[length(tmp.tex)]=str_replace_all(tmp.tex[length(tmp.tex)],"[\\]","")
  
  tmp.body=paste(tmp.tex,collapse="\n")
  return(tmp.body)
}

build.eqn.tex = function(eqn, eqn.special, x, greek, digits, simplify, headfoot, t=NULL){
  #if simplify is FALSE, then will need t if parameters are time-varying
  #eqn.special are the terms like x, y, v, w
  #greek means whether to turn things like alpha into \alpha
  tmp.eqn.tex=c()
  #get the model from x (x is a marssMODEL object)
  fixed=x$fixed
  free=x$free
  model.dims=attr(x,"model.dims")
  m=model.dims[["x"]][1]; n=model.dims[["y"]][1]  
  model.el=attr(x,"par.names")
  eqn.name = headfoot$eqn.name
  mat.head = headfoot$mat.head
  mat.foot = headfoot$mat.foot
  #this takes a text form of an equation and builds the tex
  for(el in strsplit(eqn,"[+]")[[1]]){
    show.el=TRUE
    for(el2 in strsplit(el,"[*]")[[1]]){
      show.el2=TRUE
      #if el2 has a time subscript, get the time (usually t or t-1 but might be anything)
      #el2.time is the time subscript and el2 is the parameter name
      if(str_detect(el2,"_")){
        if(missing(t)) #then not time.expand
          el2.time=str_trim(strsplit(el2,"_")[[1]][2])
        else el2.time=paste("{",as.character(t),"}",sep="")
        el2=strsplit(el2,"_")[[1]][1]
      }else{
        el2.time=NULL
      }
      el2=str_trim(el2) #strip any extra while space
      #Detect if el2 is a number, like 0, and print as a number
      options(warn=-1) #don't print warnings if el2 cannot be converted to numeric
      if(!is.na(as.numeric(el2))) el2=as.numeric(el2)
      options(warn=0)
      
      if(!is.numeric(el2) & !(el2 %in% c(model.el,eqn.special))) stop("toLatex.marssMODEL: in the attribute, equation, one of the names does not match a name in fixed or model.dims.")
      #if els is a number, just print as a number
      if(is.numeric(el2)){ mat.tex=as.character(el2)
      }else{
        if(el2 %in% eqn.special){ 
          ######################################
          #if el2 is something like x, y, v or w
          ######################################
          #if el2 is a eqn.name then use its names, otherwise it is a error name
          if(el2 %in% eqn.special[seq(1,length(eqn.special),2)]){ el2.names = el2
          }else{ el2.names = eqn.special[which(eqn.special==el2)-1] }
          labs=attr(x,paste(toupper(el2.names),".names",sep=""))
          mat.body=paste(paste(el2,"_{",labs,"}",sep=""),collapse="\\\\\n")
          if(is.null(el2.time)) stop(paste("toLatex.marssMODEL: in the attribute, equation, ", el2, " is missing the time designation.  should look like _{t}.",sep=""))
          mat.tex=paste(mat.head,"\n",mat.body,"\n",mat.foot,"_",el2.time,sep="")
        }else{ 
          ######################################
          # el2 is a parameter name
          ######################################
          if(!is.null(colnames(free[[el2]]))) #this adds the parameter name to the element name, like B.1
            colnames(free[[el2]])=paste(el2,colnames(free[[el2]]),sep=".")
          if(model.dims[[el2]][3]==1){
            mat= fixed.free.to.formula(fixed[[el2]],free[[el2]],model.dims[[el2]][1:2])
          }else{
            mat = array(NA, dim=model.dims[[el2]])
            for(i in 1:model.dims[[el2]][3])
              mat[,,i]=fixed.free.to.formula(fixed[[el2]],free[[el2]],model.dims[[el2]][1:2])
          }
          #if simplify=TRUE only show mats that are not all 0 and that are not identity
          if(simplify) show.el=show.el&!(all(sapply(mat,identical,0)))
          if(simplify){
            if(length(dim(mat))==3) show.el2=!all(apply(mat,3,is.identity))
            if(length(dim(mat))==2) show.el2=!is.identity(mat)
          }
          if(show.el2 & show.el){
            #if mat is 2D then it is not time-varying so no time subscript
            if(!is.matrix(mat)){ #then it is 3D
               #if mat is 3D but all the mats are equal, then it is not time-varying
              #The following tests for that
              mat2D=array(mat,dim=c(dim(mat)[1]*dim(mat)[2],dim(mat)[3])) #reform to be 2D
              #test equality within row; will return TRUE if dim 3 = 1
              if(all(apply(mat,1,all.equal.vector))){ #not time-varying 
                mat=sub3D(mat,t=1) #not time-varying so just use the mat at t=1
              }else{ #is time-varying
                if(!simplify){ #then there will be an equation for each parameter at time t
                  if(missing(t)) stop("build.eqn.tex: if simplify is FALSE and parameters are time-varying, need to specify t")
                  mat=sub3D(mat,t=t)
                }else{ #do simplify
                  #step 1: figure out which parameter elements are time-varying
                  mat2D=array(mat,dim=c(dim(mat)[1]*dim(mat)[2],dim(mat)[3])) #reform to be 2D, matrix(vec(mat),t)
                  tv.elem=!apply(mat2D,1,all.equal.vector) #which elements are time varying
                  #step 2: replace that element with the name el2(i)_t
                  tv.name = paste(el2,"(",1:dim(mat2D)[1],")_",el2.time,sep="")
                  #step 3: set up a mat to be t=1; will replace time=varying ones
                  mat=sub3D(mat,t=1)
                  #replace any time-varying elements with tv.name
                  mat[tv.elem]=tv.name[tv.elem]
                }
              }
            }
            mat.body=get.mat.tex(mat,greek,digits)
            mat.tex=paste(mat.head,"\n",mat.body,"\n",mat.foot,sep="")
          }
        }
      }
      if(show.el & show.el2) tmp.eqn.tex=c(tmp.eqn.tex,mat.tex)
    }
    if(show.el) tmp.eqn.tex=c(tmp.eqn.tex,"+")
  } 
  tmp.eqn.tex=tmp.eqn.tex[-1*length(tmp.eqn.tex)]
  return(tmp.eqn.tex)
}

parameters.time.varying = function(eqn, eqn.special, x){
  #this will search through the equation parameters and determine if any parameters are time-varying
  #get the model from x (x is a marssMODEL object)
  fixed=x$fixed
  free=x$free
  if(!is.null(attr(x,"model.dims"))){ model.dims=attr(x,"model.dims")
  }else{ model.dims=x$model.dims }
  model.el=names(fixed)
  
  #this finds the bits between + in the equation
  for(el in strsplit(eqn,"[+]")[[1]]){
    for(el2 in strsplit(el,"[*]")[[1]]){ #this finds the bits between * in el
      if(str_detect(el2,"_")){
        el2.time=str_trim(strsplit(el2,"_")[[1]][2])
        el2=strsplit(el2,"_")[[1]][1]
      }else{
        el2.time=NULL  #no time subscript
        #the parameter should not be time-varying then....
      }
      el2=str_trim(el2) #strip any extra while space
      #Detect if el2 is a number, like 0, and print as a number
      options(warn=-1) #don't print warnings if el2 cannot be converted to numeric
      if(!is.na(as.numeric(el2))){
        el2=as.numeric(el2)
        if(!is.null(el2.time)) stop("toLatex.marssMODEL: equation has a number with a time-subscript.  That doesn't make sense.\nCheck the equation attribute for the marssMODEL.\n")
      }
      options(warn=0)
      
      if(!is.numeric(el2) & !(el2 %in% c(model.el,eqn.special))) stop("toLatex.marssMODEL: in the attribute, equation, one of the names does not match a name in fixed or model.dims.")
      #if els is a number, just print as a number
      if(is.numeric(el2)){ 
        next #go to next parameter in model
      }
      if(el2 %in% eqn.special){ 
        next #ignore the specials.  They always are time-varying
      } 
      if(model.dims[[el2]][3]==1) next # it is not time-varying
      #Even if 3D, it might not be time-varying if all elements are the same
      testtv=FALSE
      if(!all(apply(fixed[[el2]],1,all.equal.vector))) testtv=TRUE #time-varying
      mat2D=array(free[[el2]],dim=c(dim(free[[el2]])[1]*dim(free[[el2]])[2],dim(free[[el2]])[3]))
      if(!all(apply(mat2D,1,all.equal.vector))) testtv=TRUE #time-varying
      if(!testtv){ #not time-varying 
        next
      }else{
        if(is.null(el2.time)) stop(paste("toLatex.marssMODEL: something is wrong. ",el2," in the equation attributes has no time subscript, yet it is time-varying in the marssMODEL object.\nCannot make latex because the correct time-subscript label is unknown.\n",sep=""))
        return(TRUE)
      }
      
    } #el2 the bits between * for el's
  }# el the bits between +
  return(FALSE)
} #end of function
gragusa/MARSS documentation built on May 17, 2019, 8:18 a.m.