R/fraping.R

Defines functions compareFit compareMean compareParam plotFit plotRecover tableFit simFit newFit newDataFrap is.dataFrap.class

is.dataFrap.class <- function(x) list.compareItem(x,"class","frap.000000uFpZ4K5kgF0kXZ8RaD5PknOYGvqg99t")
#is.Fit.class <- function(x) list.compareItem(x,"class","fit.000000Wy43veDItMRutEj1wG3qbLRSxFhPF6TtunN0b")

newDataFrap <- function(name, bleach){
  fmin=1-bleach/100
  Color<-newVec(color())
  Seconds<-newVec()
  Photo  <-newVec()
  Time0  <-newVec()
  Time  <-newList()
  IntDen<-newList()
  Area  <-newList()
  Rlist<-newList()
  Tablist<-newList()
  newAdded<-newVec(FALSE)
  
  addFile <- function(seconds, photo, fileR, fileC, fileB=""){
    newAdded$set(TRUE)
    tabR <- read.csv(fileR)
    tabC <- read.csv(fileC)
    tabB <- if(fileB!="") {
      read.csv(fileB)
    }else{
      datf<-data.frame(rep(1, nrow(tabR)), rep(0, nrow(tabR)))
      colnames(datf)<-c("Area", "IntDen")
      datf
    }
    time<-seq(0, seconds, seconds/(nrow(tabR)-1))
    Time$add(time)
    intden<-newList()
    intden$set(tabR$IntDen, "R")
    intden$set(tabC$IntDen, "C")
    intden$set(tabB$IntDen, "B")
    IntDen$add(intden)
    area<-newList()
    area$set(tabR$Area[1], "R")
    area$set(tabC$Area[1], "C")
    area$set(tabB$Area[1], "B")
    Area$add(area)
    
    Seconds$add(seconds)
    Photo  $add(photo)
    Time0  $add(time[photo])
  }
  addDir <- function(seconds, photo, dirR, dirC, dirB=""){
    DIRR <- paste(dirR, list.files(dirR), sep="/")
    DIRC <- paste(dirC, list.files(dirC), sep="/")
    DIRB <- paste(dirB, list.files(dirB), sep="/")
    for (i in 1:length(DIRR))
      addFile(seconds, photo,
              DIRR[i], DIRC[i],if(dirB=="") "" else DIRB[i])
  }
  area <- function(index=NA, type="R") {
    if(Time$length()==0) return(NULL)
    if(is.na(index[1])) index<-1:Time$length()
    if(length(index)>1) return(funVec(area, index, type=type))
    Area$get(index)$get(type)
  }
  recover <- function(index=NA, type="RCB", area=T, stand=T, AB=F){
    if(Time$length()==0) return(NULL)
    if(is.na(index[1])) index<-1:Time$length()
    if(length(index)>1) 
      return(funList(recover, index, type=type, area=area, stand=stand, AB=AB))
    code<-join(type, area, stand, AB, index, sep=".")
    recover<-Rlist$get(code)
    if(!is.null(recover)) return(recover)
    getIntDen <- function(index, type="R", area=T) {
      IntDen$get(index)$get(type)/(if(area) area(index, type) else 1)
    }
    recover<-if(type=="R" | type=="C" | type=="B"){
      getIntDen(index, type, area)
    }else if(type=="RC"){
      getIntDen(index, "R",area)/getIntDen(index, "C", area)
    }else if(type=="RB"){
      getIntDen(index, "R",area)-getIntDen(index, "B", area)
    }else if(type=="CB"){
      getIntDen(index, "C",area)-getIntDen(index, "B", area)
    }else if(type=="RCB"){
      (getIntDen(index, "R",area)-getIntDen(index, "B", area))/
        (getIntDen(index, "C",area)-getIntDen(index, "B", area))
    }else {
      NULL  
    }
    if(stand)
      recover <- reparam(recover, c(recover[Photo$get(index)], recover[1]), c(fmin, 1))
    if(AB)
      recover <- recover[Photo$get(index):length(Time$get(index))]
    Rlist$set(recover, code)
    return(recover)
  }
  time <- function(index=NA, AB=F){
    if(Time$length()==0) return(NULL)
    if(is.na(index[1])) index<-1:Time$length()
    if(length(index)>1) return(funList(time, index, AB=AB))
    if(!AB) return(Time$get(index))
    Time$get(index)[1:(length(Time$get(index))-Photo$get(index)+1)]
  }
  getColor<-function() Color$get()
  
  setColor<-function(col) Color$set(col)
  
  table<- function(type="RCB", area=T, stand=T, AB=F){
    code<-join(type, area, stand, AB, sep=".")
    tab<-Tablist$get(code)
    if(newAdded$get() | is.null(tab)){
      tab<-tableLinesInter(
        time(AB=AB), recover(type=type,area=area,stand=stand,AB=AB))
      Tablist$set(tab, code)
      newAdded$set(FALSE)
    }
    return(tab)
  }
  
  list(addFile=addFile, addDir=addDir,
       recover=recover, time=time, area=area, seconds=Seconds$get, 
       photo=Photo$get, time0=Time0$get, getColor=getColor, 
       setColor=setColor, table=table, length=Time$length, 
       name=name, bleach=bleach, fmin=fmin,
       class="frap.000000uFpZ4K5kgF0kXZ8RaD5PknOYGvqg99t")
}

newFit <- function(name, fun, param, interval){
  fn<-newFunction("fun")
  fn$add("fit<-function(",join(param, sep=", "),", alpha, fmin, x)"," fmin+alpha*(1-fmin)*fun(x, ",join(pst(param,"=",param), sep=", "),")")
  fn$add("function(",join(param, sep=", "),", alpha, fmin, x, y) sum(abs(y-fit(",join(param, sep=", "),", alpha, fmin, x)))")
  error<-fn$run(fun)
  param<-c(param, "alpha")
  interval<-c(interval, list(c(0,1)))
  n<-length(param)
  function(data, type="RCB", area=T, stand=T){#########No perder de vista
    code<- join("\\fit",data$name, name, param, round(listToVec(interval,1:length(interval)),10), type, area, stand, sep="_")
    op<-globalVar(code)
    if(!base::is.null(op)) return(op)
    table<-data.frame()
    for (i in 1:data$length()) {
      x<-data$time(i, AB=T)
      y<-data$recover(i, type=type, area=area, stand=stand, AB=T)
      op<-TRY(optima(error, param, interval, fmin=data$fmin, x=x, y=y),
              {print("Optimization failed");rep(NA, n+1)})
      table<-rbind(table, op)
    }
    names(table) <-c(param,"error")
    error<-table$error
    table$error<-NULL
    table$fmax<-data$fmin+table$alpha*(1-data$fmin)
    table$UF<-(1-table$fmax)/(1-data$fmin)
    table$MF<-1-table$UF
    table$error<-error
    evaluate<-function(table) {
      fn<-newFunction("fun")
      fn$add("fit<-list()")
      fn$add("prob<-list()")
      for (i in 1:nrow(table)){
        fn$add("fit[[",i,"]]<-function(x) ",data$fmin,"+",table[i, n]*(1-data$fmin),"*fun(x, ",
               join(pst(param[1:(n-1)],"=",table[i, 1:(n-1)]), sep=", "),")")
        fn$add("prob[[",i,"]]<-function(x) fun(x, ", join(pst(param[1:(n-1)],"=",table[i, 1:(n-1)]), sep=", "),")")  
      }    
      fn$add("return(list(fit=fit, prob=prob))")
      write(fn$getCode(),"out.R")
      fn$run(fun)
    }
    eva<-evaluate(table)
    op<-list(name=name, data=data, evaluate=evaluate, table=table,
             fit=eva$fit, prob=eva$prob, saveable=TRUE, simulated=FALSE, code=code)
    globalVar(code, op)
    return(op)
  }
}

simFit <- function(fit, Nsim=50, seed=NA, saveable=F){
  code=NULL
  if(!is.na(seed) & saveable) {
    saveable<-TRUE
    set.seed(seed)
    code <- join("\\simfit",fit$data$name, fit$name, Nsim, seed, fit$code, sep="_")
    sft <-globalVar(code)
    if(!base::is.null(sft)) return(sft)
  }
  n<-ncol(fit$table)-5
  tab<-fit$table[,1:(n+1)]
  fn<-newFunction("tab")
  fn$add("lm(tab$alpha~",join(pst("tab$",colnames(tab)[1:n]),sep="+"),")$coefficients")
  coef<-fn$run(tab)
  tab$alpha<-NULL
  sim<-simulator(tab)
  table<-NULL
  while (IF(is.null(table), TRUE, nrow(table)<Nsim)) {
    sm<-sim(1)
    alpha<-sum(c(1, sm)*coef)
    table<-rbind(table, c(sm,alpha))
  }
  table<-as.data.frame(table)
  colnames(table)<-c(colnames(tab),"alpha")
  table$fmax <- fit$data$fmin+table$alpha*(1-fit$data$fmin)
  table$UF<-(1-table$fmax)/(1-fit$data$fmin)
  table$MF<-1-table$UF
  eva<-fit$evaluate(table)
  sft<-list(name=pst(fit$name,"Sim"), data=fit$data, table=table, 
            fit=eva$fit, prob=eva$prob, saveable=saveable, simulated=TRUE, code=code)
  if(saveable) globalVar(code, sft)
  return(sft)
}

tableFit <- function(fit, timelim=NULL, npoints=100, displacement=F) {
  code<-NULL
  if(fit$saveable){
    code <- join("\\fitable",fit$data$name, fit$name, timelim, npoints, displacement, fit$code, sep="_")
    ft <-globalVar(code)
    if(!base::is.null(ft)) return(ft)
  }
  n<-fit$data$length()
  if(!is.null(timelim)){
    minTime<-rep(timelim[1],n)
    maxTime<-rep(timelim[2],n)
  }else{
    minTime<-rep(0,n)
    maxTime<-IF(displacement,
                fit$data$seconds()-fit$data$time0(),
                fit$data$seconds())
  }
  DISP <- IF(displacement, fit$data$time0(), rep(0,n))
  TIM<-funList(function(i) seq(minTime[i], maxTime[i],(maxTime[i]-minTime[i])/(npoints-1))+DISP[i], 1:n)
  FIT<-funList(function(i) (fit$fit[[i]])(TIM[[i]]-DISP[i]), 1:n)
  ft <- list(TIM=TIM, FIT=FIT, table=tableLinesInter(TIM, FIT), code=code)
  if(fit$saveable) globalVar(code, ft)
  return(ft)
}

plotRecover <- function(..., type="RCB", area=T, stand=T, AB=F, #########No perder de vista
                        #Plot
                        new.plot = T, plot.lines=F, plot.points=F,
                        plot.shadow=F, plot.mean=F, col=NULL, getGroup=F) {
  type <- iterVector(type, cyclic=T)
  area <- iterVector(area, cyclic=T)
  stand <- iterVector(stand, cyclic=T)
  AB <- iterVector(AB, cyclic=T)
  groupRecover <- function(data){
    list(X =data$time(AB=AB$This()),
         Y =data$recover(type=type$This(), area=area$This(), stand=stand$This(), AB=AB$This()),
         table =data$table(type=type$Next(), area=area$Next(), stand=stand$Next(), AB=AB$Next()),
         col=fixCol(data$getColor()), class=group.class())
  }
  if(getGroup) return(groupRecover)
  ret<-substrae.LIST.DAT(...,
                         CONDITION=is.dataFrap.class,
                         TRANSFORM=function(x) groupRecover(x))
  plotGroup(
    lista.123=ret$lista.123, data.123=ret$data.123,
    new.plot=new.plot, plot.lines=plot.lines, plot.points=plot.points, 
    plot.shadow=plot.shadow, plot.mean=plot.mean, col=col)
}

plotFit <- function(..., fit, type="RCB", area=T, stand=T, #########No perder de vista
                    simulated=F, Nsim=50, seed=NA,
                    npoints=100, displacement=F,
                    new.plot = T, plot.lines=F, plot.points=F,
                    plot.shadow=F, plot.mean=F, col=NULL) {
  type <- iterVector(type, cyclic=T)
  area <- iterVector(area, cyclic=T)
  stand <- iterVector(stand, cyclic=T)
  simulated <- iterVector(simulated, cyclic=T)
  Nsim <- iterVector(Nsim, cyclic=T)
  seed <- iterVector(seed, cyclic=T)
  fit <- iterList(IF(is.list(fit), fit, list(fit)), cyclic=T)
  groupFit <- function(data) {
    type<-type$Next()
    area<-area$Next()
    stand<-stand$Next()
    simulated<-simulated$Next()
    Nsim<-Nsim$Next()
    seed<-seed$Next()
    col<-data$getColor()
    fit<-fit$Next()
    if(is.null(fit)){
      return(plotRecover(type=type, area=area, stand=stand, AB=T, col=col, getGroup=T)(data))
    }
    ft<-fit(data, type=type, area=area, stand=stand)
    if(simulated)ft<-simFit(ft, Nsim=Nsim, seed=seed)
    tab<-tableFit(ft, timelim=list(...)$xlim, npoints=npoints, displacement=displacement)
    list(X=tab$TIM, Y=tab$FIT, table=tab$table, col=fixCol(col), class=group.class())
  }
  ret<-substrae.LIST.DAT(...,
                         CONDITION=is.dataFrap.class,
                         TRANSFORM=function(x) groupFit(x))
  plotGroup(
    lista.123=ret$lista.123, data.123=ret$data.123,
    new.plot=new.plot, plot.lines=plot.lines, plot.points=plot.points, 
    plot.shadow=plot.shadow, plot.mean=plot.mean, col=col)
}

compareParam <- function(data1, data2, fit, param,
                         type="RCB", area=T, stand=T,
                         simulated=F, Nsim=50, seed=NA,
                         conf.level=0.95, alternative="two.sided", view=F,
                         new.plot=T, plot.lines=T, plot.points=F, col=NULL,
                         ...){
  lst<-list(...)
  fit<-iterList(IF(is.list(fit), fit, list(fit)), cyclic=T)
  type<-iterVector(type, cyclic=T)
  area<-iterVector(area, cyclic=T)
  stand<-iterVector(stand, cyclic=T)
  fit1<-fit$Next()(data1, type=type$Next(), area=area$Next(), stand=stand$Next())
  fit2<-fit$Next()(data2, type=type$Next(), area=area$Next(), stand=stand$Next())
  simulated<-iterVector(simulated, cyclic=T)
  Nsim<-iterVector(Nsim, cyclic=T)
  seed<-iterVector(seed, cyclic=T)
  if(simulated$Next()) fit1<-simFit(fit1, Nsim=Nsim$Next(), seed=seed$Next())
  if(simulated$Next()) fit2<-simFit(fit2, Nsim=Nsim$Next(), seed=seed$Next())
  p<-(0:100)/100
  n1<-sum((names(fit1$table)==param)*(1:ncol(fit1$table)))
  n2<-sum((names(fit2$table)==param)*(1:ncol(fit2$table)))
  param1<-fit1$table[,n1]
  param2<-fit2$table[,n2]
  name1<-join(names(fit1$table)[n1],fit1$name,fit1$data$name,sep="_")
  name2<-join(names(fit2$table)[n2],fit2$name,fit2$data$name,sep="_")
  if(name1==name2){
    name1<-pst(name1, "1",sep="_")
    name2<-pst(name2, "2",sep="_")
  }
  fn<-newFunction(name1,name2)
  fn$add("t.test(",name1,", ",name2,", conf.level=",conf.level,", alternative=\"",alternative,"\")")
  fn$getCode()
  #####Plot
  name1<-join(names(fit1$table)[n1],fit1$name,fit1$data$name,sep=" ")
  name2<-join(names(fit2$table)[n2],fit2$name,fit2$data$name,sep=" ")
  if(name1==name2){
    name1<-pst(name1, "1",sep=" ")
    name2<-pst(name2, "2",sep=" ")
  }
  q1<-quantile(param1, p)
  q2<-quantile(param2, p)
  if(new.plot){
    lim<-c(min(q1, q2), max(q1, q2))
    if(is.null(lst$xlim)) lst$xlim<-lim
    if(is.null(lst$ylim)) lst$ylim<-lim
    if(is.null(lst$xlab)) lst$xlab<-name1
    if(is.null(lst$ylab)) lst$ylab<-name2
    ref<-c(max(lst$xlim), min(lst$ylim))
    
    G1<-list(X=ref,Y=ref, col=color("#333333"),class=group.class())
    G2<-list(X=q1,Y=q2, col=color("blue"),class=group.class())
    plotGroup(data.123=newList(G1, G2), lista.123=lst,
              plot.lines=plot.lines, plot.points=plot.points,
              plot.shadow=F, plot.mean=F, col=col)
  }else{
    G2<-list(X=q1,Y=q2, col=color("blue"), class=group.class())
    linesGroup(data.123=newList(G2), lista.123=lst,
              plot.lines=plot.lines, plot.points=plot.points,
              plot.shadow=F, plot.mean=F, col=col)
  }
  if(view) fn$run(param1, param2)
}

compareMean <- function(data1, data2, fit,
                        type="RCB", area=T, stand=T,
                        simulated=F, Nsim=50, seed=NA,
                        conf.level=0.95, alternative="two.sided",
                        view=F, p.value=T, npoints=100,
                        new.plot=T, plot.lines=T, plot.points=F, col=NULL,
                        ...){
  lst<-list(...)
  fit<-iterList(IF(is.list(fit), fit, list(fit)), cyclic=T)
  type<-iterVector(type, cyclic=T)
  area<-iterVector(area, cyclic=T)
  stand<-iterVector(stand, cyclic=T)
  fit1<-fit$Next()(data1, type=type$Next(), area=area$Next(), stand=stand$Next())
  fit2<-fit$Next()(data2, type=type$Next(), area=area$Next(), stand=stand$Next())
  simulated<-iterVector(simulated, cyclic=T)
  Nsim<-iterVector(Nsim, cyclic=T)
  seed<-iterVector(seed, cyclic=T)
  simulated$Restart()
  if(simulated$Next()) fit1<-simFit(fit1, Nsim=Nsim$Next(), seed=seed$Next())
  if(simulated$Next()) fit2<-simFit(fit2, Nsim=Nsim$Next(), seed=seed$Next())
  if(is.null(lst$xlim)) lst$xlim<-c(0, min(fit1$data$seconds(), fit2$data$seconds()))
  table1<-tableFit(fit1, timelim=lst$xlim, npoints=npoints, displacement=F)
  table2<-tableFit(fit2, timelim=lst$xlim, npoints=npoints, displacement=F)
  if(p.value){
    pval<-NULL
    ttest<-NULL
    for (i in 1:npoints) {
      group1<-table1$table$tab[i,]
      group2<-table2$table$tab[i,]
      ttest[[i]]<-TRY(t.test(group1, group2, conf.level=conf.level, alternative=alternative), list(p.value=NA))
      pval[i] <- ttest[[i]]$p.value
    }
    if(new.plot){
      refy<-rep(1-conf.level,2)
      line<-pval
      if(is.null(lst$ylim)) lst$ylim<-c(0, max(refy, pval, na.rm=T))
      if(is.null(lst$xlab)) lst$main<-"P-Value"
      if(is.null(lst$xlab)) lst$xlab<-"Time"
      if(is.null(lst$ylab)) lst$ylab<-"Probability"
    }
  }else{
    if(new.plot){
      line<-table1$table$mean-table2$table$mean
      if(is.null(lst$ylim)) lst$ylim<-c(min(line, na.rm=T), max(line, na.rm=T))
      if(is.null(lst$ylab)) lst$main<-"Difference of means"
      if(is.null(lst$xlab)) lst$xlab<-"Time"
      if(is.null(lst$ylab)) lst$ylab<-"Difference"
      refy<-rep(0, 2)
    }
  }
  if(new.plot){
    G1<-list(X=lst$xlim,Y=refy, col=color("#333333"),class=group.class())
    G2<-list(X=table1$table$z, Y=line, col=color("blue"), class=group.class())
    plotGroup(data.123=newList(G1, G2), lista.123=lst,
              plot.lines=plot.lines, plot.points=plot.points,
              plot.shadow=F, plot.mean=F, col=col)
  }else{
    G2<-list(X=table1$table$z, Y=line, col=color("blue"),class=group.class())
    linesGroup(data.123=newList(G2), lista.123=lst,
               plot.lines=plot.lines, plot.points=plot.points,
               plot.shadow=F, plot.mean=F, col=col)
  }
  if(view) if(p.value) return(ttest) else return(line)
}

compareFit <- function(data, fit){
  myhist <- function(x, digits=5, ...){
    lst<-list(...)
    mlst<-Lines("\\get.list")
    mlst$alp.lines<-isnt.null(lst$alp.lines,mlst$alp.lines)
    mlst$col.lines<-isnt.null(lst$col.lines,mlst$col.lines)
    mlst$lty.lines<-isnt.null(lst$lty.lines,mlst$lty.lines)
    mlst$lwd.lines<-isnt.null(lst$lwd.lines,mlst$lwd.lines)
    
    lst$alp.abline<-mlst$alp.lines
    lst$col.abline<-mlst$col.lines
    lst$lty.abline<-mlst$lty.lines
    lst$lwd.abline<-mlst$lwd.lines
    lst$col.abline.auxiliar<-mlst$col.lines.auxiliar
    lst$alp.lin<-mlst$alp.lines
    lst$col.lin<-mlst$col.lines
    lst$lty.lin<-mlst$lty.lines
    lst$lwd.lin<-mlst$lwd.lines
    lst$col.lin.auxiliar<-mlst$col.lines.auxiliar
    
    mhist(x, param=lst)
    mn<-round(mean(x), digits=digits)
    mabline(v=mn, sld.abline=F, param=lst)
    mlegend(join("Mean = ", mn), param=lst)
  }
  n<- length(fit)
  par(mfrow=c(n,n))
  for (i in 1:n) {
    for (j in 1:n) {
      if(i==j){
        myhist(fit[[i]](data)$table$error,cex.legend=.6)
      }else{
        compareParam(data, data, fit=c(fit[[i]],fit[[j]]), param="error")
      }    
    } 
  }
  par(mfrow=c(1,1))
  
}
artitzco/prueba documentation built on June 2, 2022, 10:28 p.m.