Implement compression of ChemDB files and reactor's output files

docker/Dockerfile

@@ -1,5 +1,6 @@
 #FROM  ppernot1/reactorui:base
-FROM  ppernot1/test
+#FROM  ppernot1/test
+FROM  ppernot1/shiny_base
 
 WORKDIR /
 

renv.lock

@@ -82,6 +82,27 @@
       "Repository": "CRAN",
       "Hash": "08588806cba69f04797dab50627428ed"
     },
+    "R.methodsS3": {
+      "Package": "R.methodsS3",
+      "Version": "1.8.1",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "4bf6453323755202d5909697b6f7c109"
+    },
+    "R.oo": {
+      "Package": "R.oo",
+      "Version": "1.24.0",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "5709328352717e2f0a9c012be8a97554"
+    },
+    "R.utils": {
+      "Package": "R.utils",
+      "Version": "2.10.1",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "a9e316277ff12a43997266f2f6567780"
+    },
     "R6": {
       "Package": "R6",
       "Version": "2.4.1",

scripts/OneRun_Loc.sh

@@ -27,6 +27,7 @@ cd -
 ./reactor
 
 # Save results
-cp ./fracmol_out.dat ../MC_Output/fracmol_${NUM4}.dat
-cp ./rrates.out      ../MC_Output/reacs_rates_${NUM4}.dat
-cp ./phrates.out     ../MC_Output/photo_rates_${NUM4}.dat
+gzip -c9kf ./fracmol_out.dat  > ../MC_Output/fracmol_${NUM4}.dat.gz
+gzip -c9kf ./rrates.out       > ../MC_Output/reacs_rates_${NUM4}.dat.gz
+gzip -c9kf ./phrates.out      > ../MC_Output/photo_rates_${NUM4}.dat.gz
+gzip -c9kf ./integ_stats.out  > ../MC_Output/integ_stats_${NUM4}.dat.gz

server_files/analysis.R

@@ -154,6 +154,57 @@ getRates = function() {
     )
   )
 }
+getIntegStats = function() {
+  # Load integrator stats
+  files = list.files(
+    path    = paste0(ctrlPars$projectDir,'/MC_Output'),
+    pattern = 'integ_stats_',
+    full.names=TRUE)
+  nf = length(files)
+
+  ## Get MC values
+  x = as.data.frame(
+    data.table::fread(files[1], header=FALSE, skip=0)
+  )
+  time = x[,1]
+  nstat = ncol(x)-1
+  ntime = nrow(x)
+  stats = array(0,dim=c(nf,ntime,nstat))
+  for(i in seq_along(files) ) {
+    tab = as.data.frame(
+      data.table::fread(files[i], header=FALSE,
+                        skip=0, colClasses='numeric')
+    )
+    stats[i,1:ntime,1:nstat] = as.matrix(tab[,-1])
+  }
+
+  yMean = apply(stats,c(2,3),
+                function(x) mean(x,na.rm=TRUE))
+  ySd   = apply(stats,c(2,3),
+                function(x) sd(x,na.rm=TRUE))
+  yLow95= apply(stats,c(2,3),
+                function(x) quantile(x,probs = 0.025,na.rm=TRUE))
+  ySup95= apply(stats,c(2,3),
+                function(x) quantile(x,probs = 0.975,na.rm=TRUE))
+
+  statNames = c('NFE', 'NFI', 'NSTEPS', 'NACCPT',
+                'NREJCT', 'NFESIG', 'MAXM')
+
+  return(
+    list(
+      nfStat   = nf,
+      nt       = ntime,
+      nStat    = nstat,
+      time     = time,
+      stats    = stats,
+      statNames= statNames,
+      statMean = yMean,
+      statSd   = ySd,
+      statLow  = yLow95,
+      statSup  = ySup95
+    )
+  )
+}
 selectSpecies <- function(species, categs) {
   # Select species to plot from categsPlot
 
@@ -791,7 +842,106 @@ output$sensitivity <- renderPlot({
 
 })
 # Sanity ####
-output$sanity <- renderPrint({
+rangesSanityInteg <- reactiveValues(x = NULL, y = NULL)
+output$sanityInteg <- renderPlot({
+  if(is.null(concList()))
+    return(NULL)
+
+  if(is.null(statsList()))
+    statsList(getIntegStats())
+
+  # Extract stats
+  for (n in names(statsList()))
+    assign(n,rlist::list.extract(statsList(),n))
+
+  # Extract graphical params
+  for (n in names(gPars))
+    assign(n,rlist::list.extract(gPars,n))
+
+  sel = rep(TRUE, nStat)
+  sel[3] = FALSE # do not show NSTEPS (=NACCPT+NREJCT)
+
+  # Define zoom range
+  if (is.null(rangesSanityInteg$x)) {
+    xlim <- range(time)# * 100 # expand for labels
+  } else {
+    xlim <- rangesSanityInteg$x
+  }
+  if (is.null(rangesSanityInteg$y)) {
+    ylim = range(c(0, statSup[, sel]), na.rm = TRUE)
+  } else {
+    ylim <- rangesSanityInteg$y
+  }
+
+  # Generate colors per stat
+  colors = assignColorsToSpecies(
+    input$colSel, statNames, sel, nf=1,
+    cols, col_tr, col_tr2)
+
+  # Show uncertainty bands ?
+  showBands = (nfStat > 1)
+
+  # Plot
+  par(mfrow = c(1, 1),
+      cex = cex, cex.main = cex, mar = mar,
+      mgp = mgp, tcl = tcl, pty = pty, lwd = lwd)
+
+  plot(time, time,
+       type = 'n',
+       log  = 'x',
+       main = ifelse(!showBands,
+                     'Mean values',
+                     'Mean and 95 % Proba. Intervals'),
+       xlab = 'Time / s',
+       xlim = xlim,
+       xaxs = 'i',
+       ylab = 'Integrator statistics',
+       ylim = ylim,
+       yaxs = 'i')
+  grid()
+
+
+  # CI
+  if(showBands) {
+    for(isp in (1:nStat)[sel])
+      polygon(c(time     , rev(time      )),
+              c(statLow[,isp],rev(statSup[,isp])),
+              col = colors$col_trSp[isp],
+              border = NA)
+  }
+
+  # Mean concentrations
+  matplot(time, statMean[,sel], type='l', lty = 1,
+          col = colors$colsSp[sel],
+          lwd = 1.2*lwd, add = TRUE)
+
+  legend(
+    'topleft', bty = 'n',
+    legend = statNames[sel],
+    col    = if(showBands) colors$col_trSp[sel] else colors$colsSp[sel],
+    lty    = 1,
+    lwd    = ifelse(showBands,20,4)
+  )
+  box()
+
+})
+outputOptions(output, "sanityInteg",
+              suspendWhenHidden = FALSE)
+observeEvent(
+  input$sanity_dblclick,
+  {
+    brush <- input$sanity_brush
+    if (!is.null(brush)) {
+      rangesSanityInteg$x <- c(brush$xmin, brush$xmax)
+      rangesSanityInteg$y <- c(brush$ymin, brush$ymax)
+    } else {
+      rangesSanityInteg$x <- NULL
+      rangesSanityInteg$y <- NULL
+    }
+  }
+)
+
+output$sanityOutputs <- renderPrint({
   if(is.null(concList())   |
      is.null(ratesList())
   )
@@ -803,16 +953,25 @@ output$sanity <- renderPrint({
   for (n in names(ratesList()))
     assign(n,rlist::list.extract(ratesList(),n))
 
+  nMC = nf
+
   # Analyze final concentrations
   conc = conc[, nt, ]
-  lconc = ifelse(conc <= 0, NA, log10(conc))
-  sdc  = apply(lconc, 2, function(x) sd(x,na.rm=TRUE))
+  if(nMC == 1) {
+    conc = matrix(conc,nrow=1)
+    lconc = ifelse(conc <= 0, NA, log10(conc))
+    sdc  = rep(1,ncol(conc))
+  } else {
+    lconc = ifelse(conc <= 0, NA, log10(conc))
+    sdc  = apply(lconc, 2, function(x) sd(x,na.rm=TRUE))
+  }
   sel  = which(sdc == 0 | !is.finite(sdc))
 
+
   if(length(sel) != 0) {
     cat(' Species with suspicious final concentrations\n',
         '---------------------------------------------\n\n')
-    nMC = nrow(conc)
+
     sd0 = nzero = ninf = c()
     for (ii in 1:length(sel)) {
       i = sel[ii]

server_files/model.R

@@ -1104,26 +1104,15 @@ observeEvent(input$sampleChem, {
     neutralsSourceDir = file.path(chemDBDir(),input$neuVers)
     ionsSourceDir     = file.path(chemDBDir(),input$ionVers)
 
-    # # Generate data files for reactor code ###
-    # sp_aux =paste(
-    #   nbSpecies,'\n',
-    #   paste(species,collapse = ' '), '\n',
-    #   paste(mass,   collapse = ' '), '\n',
-    #   paste(rep(0,nbSpecies), collapse = ' ')
-    # )
-    # writeLines(
-    #   sp_aux,
-    #   file.path(outputDir,'Run','species_aux.dat')
-    # )
-
     # Split photodissociations and reactions
     photo = type == 'photo'
     Dphoto = matrix(D[ photo,],nrow=sum(photo), ncol=nbSpecies)
     D      = matrix(D[!photo,],nrow=sum(!photo),ncol=nbSpecies)
     Lphoto = matrix(L[ photo,],nrow=sum(photo), ncol=nbSpecies)
     L      = matrix(L[!photo,],nrow=sum(!photo),ncol=nbSpecies)
 
-    # Treat Reactions ###
+    # Treat Reactions ---
+    # Reduce full database and save to compressed file
 
     if(nMC > 1) {
       # Clean target dir
@@ -1139,20 +1128,20 @@ observeEvent(input$sampleChem, {
         sel = !photo
         origs = unique(unlist(orig[sel]))
         for (iMC in 0:nMC) {
+
           sampleFile = paste0('run_', sprintf('%04i', iMC), '.csv')
+
           data = ''
           for (io in origs) {
+
+            # Get full database
             sourceDir = io
             filename = file.path(sourceDir, sampleFile)
-            # scheme   = read.csv(file = filename,
-            #                    header = FALSE,
-            #                    sep = ';')
             scheme   = as.data.frame(
               data.table::fread(file = filename,
                                 header = FALSE,
                                 sep = ';')
             )
-
             scheme  = t(apply(scheme, 1, function(x)
               gsub(" ", "", x)))
             paramsLoc = list()
@@ -1162,6 +1151,8 @@ observeEvent(input$sampleChem, {
               terms = scheme[i, 8:ncol(scheme)]
               paramsLoc[[ii]] = terms[!is.na(terms) & terms != '']
             }
+
+            # Select and collate data for reduced scheme (locnum)
             selOrig = which (orig == io)
             paramsLoc = paramsLoc[unlist(locnum[selOrig])]
             for (j in 1:length(paramsLoc))
@@ -1170,17 +1161,25 @@ observeEvent(input$sampleChem, {
                 paste(paramsLoc[[j]], collapse = ' '),
                 '\n')
           }
+
+          # Save zipped data
           writeLines(
             data,
-            file.path(targetMCDir, 'Reactions', sampleFile)
+            gzfile(
+              file.path(targetMCDir, 'Reactions',
+                        paste0(sampleFile,'.gz'))
+            )
+            # file.path(targetMCDir, 'Reactions', sampleFile) # Unzipped
           )
 
           incProgress(1/nMC, detail = paste('Sample', iMC))
+
         }
       })
     }
 
-    # Treat Photo-processes ###
+    # Treat Photo-processes ---
+    # On a file copy basis (compression is managed beforehand)
 
     if (dim(Dphoto)[1]!=0) {
 
@@ -1196,40 +1195,53 @@ observeEvent(input$sampleChem, {
 
         withProgress(message = 'PhotoProcs', {
           for (iMC in 0:nMC) {
+
             prefix=paste0(sprintf('%04i',iMC),'_')
+
             for (i in (1:nbReac)[photo]) {
               sp=species[which(Lphoto[i,]!=0)]
-              file = paste0(prefix,'se',sp,'.dat')
+
+              # Cross-sections
+              file = paste0(prefix,'se',sp,'.dat.gz')
               fromFile = file.path(photoSourceDir, file)
               toFile   = file.path(targetMCDir, 'Photoprocs', file)
               file.copy(from = fromFile, to = toFile)
+
+              # Branching ratios
               if( params[[i]][1] != 0 ) {
-                file = paste0(prefix,'qy',sp,'_',params[[i]][1],'.dat')
+                file = paste0(prefix,'qy',sp,'_',params[[i]][1],'.dat.gz')
                 fromFile = file.path(photoSourceDir, file)
                 toFile   = file.path(targetMCDir, 'Photoprocs', file)
                 file.copy(from = fromFile, to = toFile)
               }
+
             }
+
             incProgress(1/nMC, detail = paste('Sample', iMC))
+
           }
         })
 
         # Copy of nominal values in ./Photo for standalone test run
         # + Remove run prefix
         iMC = 0
         prefix=paste0(sprintf('%04i',iMC),'_')
+
         for (i in (1:nbReac)[photo]) {
+
           sp=species[which(Lphoto[i,]!=0)]
-          file = paste0('se',sp,'.dat')
+          file = paste0('se',sp,'.dat.gz')
           fromFile = file.path(photoSourceDir, prefix, file)
           toFile   = file.path(outputDir, 'Run','Photo', file)
           file.copy(from = fromFile, to = toFile)
+
           if( params[[i]][1] != 0 ) {
-            file = paste0('qy',sp,'_',params[[i]][1],'.dat')
+            file = paste0('qy',sp,'_',params[[i]][1],'.dat.gz')
             fromFile = file.path(photoSourceDir, prefix, file)
             toFile   = file.path(outputDir, 'Run','Photo', file)
             file.copy(from = fromFile, to = toFile)
           }
+
         }
       }
     }

server_files/project.R

@@ -7,6 +7,7 @@ spectrumData = reactiveVal(NULL)
 reacScheme   = reactiveVal(NULL)
 concList     = reactiveVal(NULL)
 ratesList    = reactiveVal(NULL)
+statsList    = reactiveVal(NULL)
 fluxesList   = reactiveVal(NULL) # flMean, flSd...
 graphsList   = reactiveVal(NULL) # connectivity mats :linksR, linksR2
 stoechList   = reactiveVal(NULL) # stoechiometry mats: L, R