fixed adding of data

demo/demo-kin-surr.R

+## Licence: LGPL version 2.1
+## Time-stamp: "Last modified 2018-05-06 21:29:53 delucia"
+
+
+library(RedModRphree)
+db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat", package="RedModRphree"), is.db=TRUE)
+
+base <- c("SOLUTION 1", "units mol/kgw", "temp 25.0", "water 1",
+          "pH 9.91 charge", "pe 4.0", "C   1.2279E-04",
+          "Ca  1.2279E-04", "Mg 0.001", "Cl 0.002", "PURE 1",
+          "O2g -0.1675 10", "KINETICS 1", "-steps 100",
+          "-bad_step_max 1000", "Calcite", "-m      0.000207",
+          "-parms  3.20", "Dolomite", "-m      0.0", "-parms  0.32",
+          "END")
+
+
+initsim <- c("SELECTED_OUTPUT", "-high_precision true",
+             "-reset false", "USER_PUNCH",
+             "-head C Ca Cl Mg H e O2g Calcite Dolomite",
+             "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), ACT(\"H+\"), ACT(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
+             "SOLUTION 1", "units mol/kgw", "temp 25.0", "water 1",
+             "pH 9.91 charge", "pe 4.0", "C   1.2279E-04",
+             "Ca  1.2279E-04", "Cl 1E-12", "Mg 1E-12", "PURE 1",
+             "O2g     -0.6788 10.0", "Calcite  0.0 2.07E-3",
+             "Dolomite 0.0 0.0", "END")
+ 
+
+selout <- c("SELECTED_OUTPUT", "-high_precision true", "-reset false",
+            "-time", "-soln", "-temperature true", "-water true",
+            "-pH", "-pe", "-totals C Ca Cl Mg",
+            "-kinetic_reactants Calcite Dolomite", "-equilibrium O2g")
+
+prop <- c("C","Ca","Cl","Mg","pH","pe","O2g", "Calcite","Dolomite")
+maxsim <- 20         ## minimum number of single chemical simulations in each time step to spawn parallelization
+grid   <- 50         ## grid discretization (number of elements)
+L      <- 0.5        ## total system length [m]
+U      <- 9.375e-6   ## imposed (constant) Darcy velocity [m3/s]
+Cf     <- 0.9        ## Courant number to add some numerical dispersion
+
+## boundary conditions (concentration of the injected fluid entering into the 1D column)
+bound <- c(H=1E-7, e=1E-4, Mg=0.001, Cl=0.002)
+
+exper <- c(50,100)
+expsim <- vector(mode="list",length=length(exper))
+timexpsim <- vector(mode="list",length=length(exper))
+
+for (i in seq_along(exper)) {
+    setup <- list(n=exper[i], len=L, U=U, base=base, first=selout, Cf=Cf, bound=bound, dt=210,
+                  prop=prop, immobile=c(7,8,9), kin=c(8,9), ann=list(O2g=-0.1675), initsim=initsim)
+    
+    timexpsim[[i]] <- system.time({
+        expsim[[i]] <- ReactTranspBalanceKin(setup=setup, db=db, step="fix_dt", maxtime=21000,
+                                             procs=4, maxsim=maxsim, writeout=FALSE, 
+                                             reduce=FALSE, surrogate=FALSE)
+    })
+}
+
+    
+par(mfrow=c(2,1))
+MatplotSingle(sim=expsim[[1]], xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]", main="Kinetic simulation, grid 50")
+MatplotSingle(sim=expsim[[2]], xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]", main="Kinetic simulation, grid 100")
+
+sam <- ExtractSamples(expsim,
+                      var_design=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"),
+                      var_result=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"))
+
+
+
+
+## library(data.table)
+
+## data.table::fwrite(sam$design, file = "RedModRphree/inst/extdata/demo_kinetic_design.csv", append = FALSE, quote = "auto",
+##   sep = ",",  eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE)
+
+## data.table::fwrite(sam$result, file = "RedModRphree/inst/extdata/demo_kinetic_result.csv", append = FALSE, quote = "auto",
+##   sep = ",",  eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE)
+
+set.seed(3560)
+
+library(caret)
+fitControl <- trainControl(method = "none",
+                           allowParallel = TRUE,
+                           verboseIter=TRUE,
+                           returnData = FALSE,
+                           returnResamp="none",
+                           savePredictions="none",
+                           predictionBounds=rep(TRUE,2))
+
+
+modk <- parallel::mclapply(colnames(sam$result), function (x) 
+    return(train(x=sam$design, y=sam$result[,x],trControl = fitControl,  ##tuneGrid = tg, 
+                 method="parRF")), mc.cores=8)
+names(modk) <- colnames(sam$result)
+
+saveRDS(modk, "models_kinetics_RF.rds")
+
+PlotModsInSample(modk, sam$result, sam$design)
+
+
+baleq <- list(Calcite=c(C=1, Ca=1), Dolomite=c(C=2, Ca=1, Mg=1))
+
+balance <- CheckBalance(baleq, sam$design, sam$result)
+
+
+MySurrKin <- function(state, model){
+    excl <- which(colnames(state)=="O2g")
+    remember_state_names <- colnames(state)
+    des <- state[,-excl]
+    if (NA %in% des) {
+        print(des)
+        stop("Mysurr: NA in data")
+    }
+    if (NaN %in% des) {
+        print(des)
+        stop("Mysurr: NaN in data")
+    }    
+    order_for_surrogate <- c("C","Ca","Cl","Mg","pH","pe","Calcite","Dolomite")
+    rem_attr <- attr(state,"immobile")
+    ss <- des[, order_for_surrogate]
+    pred <- lapply(names(model), function(x) as.numeric(predict.train(model[[x]], ss)))
+    out <- cbind(as.data.frame(pred), state[,excl])
+    colnames(out) <- c(names(model),"O2g")
+
+    ## order back to how it is expected in "state"
+    out <- out[,remember_state_names]
+    attr(out,"immobile") <- rem_attr 
+    return(out)
+}
+
+exper <- expand.grid(grid=c(50, 100), tol=c(1E-6, 1E-5))
+simsur    <- vector(mode="list",length=nrow(exper))
+timesur   <- vector(mode="list",length=nrow(exper))
+outbalsur <- vector(mode="list",length=nrow(exper))
+
+for (i in seq(nrow(exper))) {
+
+    setup <- list(n=exper$grid[i], len=L, U=U, base=base, first=selout, Cf=Cf, bound=bound, dt=210,
+                  prop=prop, immobile=c(7,8,9), kin=c(8,9), ann=list(O2g=-0.1675), initsim=initsim)
+
+    timesur[[i]] <- system.time( simsur[[i]] <- ReactTranspBalanceKin(setup=setup, db=db, step="fix_dt", maxtime=21000,
+                                                                      procs=4, maxsim=maxsim, writeout=FALSE, 
+                                                                      reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurrKin, 
+                                                                      model=modk, baleq=baleq, tol=exper$tol[i], call_pqc=TRUE) )
+    outbalsur[[i]] <- out_bal
+}
+
+
+par(mfrow=c(2, 2))
+inds <- rep(c(1,2),2)
+for (i in seq(nrow(exper))) {
+    Matplot(sur=simsur[[i]], sim=expsim[[ inds[i]  ]],
+            xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]",
+            main=paste("RandomForest surrogate, tol=", exper$tol[i]))
+}