GLUE first working example

RBBGCMuso/DESCRIPTION

@@ -24,7 +24,8 @@ Imports:
 	tibble,
 	tidyr,
 	tcltk,
-	digest
+	digest,
+	data.table
 LinkingTo: Rcpp
 SystemRequirements: C++11
 Maintainer: Roland Hollo's <hollorol@gmail.com>

RBBGCMuso/NAMESPACE

@@ -18,10 +18,12 @@ export(musoRand)
 export(musoRandomizer)
 export(musoSensi)
 export(normalMuso)
+export(optiMuso)
 export(paramSweep)
 export(plotMuso)
 export(plotMusoWithData)
 export(randEpc)
+export(readMeasuredMuso)
 export(runMuso)
 export(rungetMuso)
 export(saveAllMusoPlots)
@@ -32,6 +34,7 @@ export(updateMusoMapping)
 import(ggplot2)
 import(utils)
 importFrom(Rcpp,evalCpp)
+importFrom(data.table,fread)
 importFrom(digest,digest)
 importFrom(dplyr,'%>%')
 importFrom(dplyr,filter)
@@ -40,6 +43,10 @@ importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(dplyr,summarize)
 importFrom(dplyr,tbl_df)
+importFrom(ggplot,aes_string)
+importFrom(ggplot,geom_point)
+importFrom(ggplot,ggplot)
+importFrom(ggplot,ggsave)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,aes_string)
 importFrom(ggplot2,element_blank)
@@ -57,7 +64,9 @@ importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_classic)
 importFrom(ggplot2,xlab)
 importFrom(ggplot2,ylab)
+importFrom(gridExtra,grid.arrange)
 importFrom(limSolve,xsample)
+importFrom(lubridate,leap_year)
 importFrom(magrittr,'%<>%')
 importFrom(magrittr,'%>%')
 importFrom(rmarkdown,pandoc_version)

RBBGCMuso/R/assistantFunctions.R

@@ -124,3 +124,51 @@ dynRound <- function(x,y,seqLen){
     return(round(a,digitNum))
 }
 
+
+readValuesFromFile  <- function(epc, linums){
+	epcFile <- readLines(epc)
+	rows <- numeric(2)
+	values <- sapply(linums, function(x){
+	    rows[1] <- as.integer(x)
+	    rows[2] <- as.integer(round(100*x)) %% 10 + 1 
+	    epcFile <- readLines(epc)
+	    selRow <- unlist(strsplit(epcFile[rows[1]], split= "[\t ]"))
+	    selRow <- selRow[selRow!=""]
+	    return(as.numeric(selRow[rows[2]]))
+
+	})
+
+	return(values)
+}
+#' readMeasuredMuso
+#'
+#' MuSo data reader
+#' @importFrom data.table fread
+#' @export
+
+readMeasuredMuso <- function(inFile,
+                             naString = getOption("datatable.na.strings","NA"), sep = ",",
+                             leapYearHandling = TRUE,
+                             convert.var = NULL,
+                             convert.scalar = 1,
+                             convert.fun = (function (x) {x * convert.scalar}),
+                             convert.file = NULL,
+                             filterCol = NULL,
+                             filterVal = 1,
+                             selVar = NULL
+                             ){
+    
+    baseData <- fread(file = inFile, na.strings = as.character(naString), sep=sep)
+    baseData <- as.data.frame(baseData)
+    if(!is.null(filterCol)){
+        filterVar<- colnames(baseData)[filterCol]
+        baseData[(baseData[,filterVar] == filterVal),selVar] <- NA
+    }
+    head(baseData)
+    if(!is.null(selVar)){
+        baseData <- cbind.data.frame(baseData,convert.fun(baseData[,selVar]))
+        colnames(baseData)[ncol(baseData)]<- paste0("M",selVar) 
+    }
+    
+    return(data.table(baseData))
+}

RBBGCMuso/R/calibration.R

@@ -1,2 +1,147 @@
-## GPP_mes <- read.csv("hhs_GPP_measured.csv", stringAsFactors)
-## head(GPP_mes$year)
+#' optiMuso
+#'
+#' This function calculates the -users specified- likelihood for random model input.
+#'
+#' @author Roland HOLLOS
+#' @param measuredDataFile a
+#' @param parameters b
+#' @param sep c
+#' @param startDate d
+#' @param endDate e
+#' @param formatString a
+#' @param filterCol a
+#' @param filterVal b
+#' @param selVar c
+#' @param outLoc c
+#' @param pretag a
+#' @param calPar a
+#' @param skipSpinup a
+#' @param iterations c
+#' @param constrains d
+#' @param likelihood d
+#' @param settings e
+#' @param leapYear b
+#' @param plotName u
+#' @importFrom ggplot ggplot aes_string geom_point ggsave
+#' @importFrom magrittr '%>%'
+#' @importFrom gridExtra grid.arrange
+#' @export 
+optiMuso <- function(measuredDataFile, parameters = NULL,
+                     sep = ",", startDate,
+                     endDate, formatString = "%Y-%m-%d",
+                     naString = NULL, leapYear = TRUE,
+                     filterCol = NULL, filterVal = 1,
+                     selVar, outLoc = "./calib",
+                     preTag = "cal-",
+                     settings =  NULL,
+                     outVars = NULL,
+                     iterations = 30,
+                     skipSpinup = TRUE,
+                     constrains = NULL,
+                     plotName = "calib.jpg",
+		     likelihood = function(x, y){
+			 exp(-sqrt(mean((x-y)^2)))
+		     },
+		     calPar = 3009)
+{
+ measuredData <- readMeasuredMuso(inFile = measuredDataFile, sep = sep, selVar = selVar,filterCol = filterCol, filterVal = filterVal)
+
+ if(is.null(parameters)){
+     parameters <- tryCatch(read.csv("parameters.csv", stringsAsFactor=FALSE), error = function (e) {
+         stop("You need to specify a path for the parameters.csv, or a matrix.")
+     })
+ } else {
+     if((!is.list(parameters)) & (!is.matrix(parameters))){
+         parameters <- tryCatch(read.csv(parameters, stringsAsFactor=FALSE), error = function (e){
+             stop("Cannot find neither parameters file neither the parameters matrix")
+         })
+     }}
+    
+ outLoc <- normalizePath(outLoc)
+ outLocPlain <- basename(outLoc)
+ currDir <- getwd()
+ 
+ if(!dir.exists(outLoc)){
+     dir.create(outLoc)
+     warning(paste(outLoc," is not exists, so it was created"))
+ }
+ 
+ outLoc <- normalizePath(outLoc)
+ 
+ if(is.null(settings)){
+     settings <- setupMuso()
+ }
+ 
+ parameterNames <- parameters[,1]
+ pretag <- file.path(outLoc,preTag)
+ npar <- length(settings$calibrationPar)
+
+    ##reading the original epc file at the specified
+ ## row numbers
+ if(iterations < 3000){
+     randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = 3000)
+     randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),]
+ } else {
+     randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = iterations)
+ }
+ 
+ origEpc <- readValuesFromFile(settings$epc[2],parameters[,2])
+ 
+ ## Prepare the preservedCalib matrix for the faster
+ ##  run.
+ 
+ pretag <- file.path(outLoc,preTag)
+ 
+ ## Creating function for generating separate
+ ## csv files for each run
+ 
+ progBar <- txtProgressBar(1,iterations,style=3)
+ colNumb <- which(settings$dailyVarCodes == calPar)
+ settings$iniInput[2] %>%
+     (function(x) paste0(dirname(x),"/",tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))) %>%
+     unlink
+ randValues <- randVals[[2]]
+ settings$calibrationPar <- randVals[[1]]
+ list2env(alignData(measuredData,dataCol = 8,modellSettings = settings,startDate = startDate,endDate = endDate,leapYear = FALSE),envir=environment())
+ 
+ modellOut <- numeric(iterations + 1) # single variable solution
+ origModellOut <- calibMuso(settings=settings,silent=TRUE)
+ write.csv(x=origModellOut, file=paste0(pretag,1,".csv"))
+ modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb])
+ for(i in 2:(iterations+1)){
+     tmp <- tryCatch(calibMuso(settings = settings,
+                               parameters = randValues[(i-1),],
+                               silent= TRUE,
+                               skipSpinup = skipSpinup)[modIndex,colNumb], error = function (e) NA)
+     
+     modellOut[i]<- likelihood(measuredData,tmp)
+     write.csv(x=tmp, file=paste0(pretag,(i+1),".csv"))
+     setTxtProgressBar(progBar,i)
+ }
+ modellOut
+ paramLines <- parameters[,2]
+ paramLines <- order(paramLines)
+ randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])]
+ randInd <- order(randInd)
+ 
+        
+
+ epcStrip <- rbind(origEpc[order(parameters[,2])],
+                   randValues[,randVals[[1]] %in% parameters[,2]][,randInd])
+ 
+ 
+ preservedCalib <- cbind(epcStrip,
+                         modellOut)
+ colnames(preservedCalib) <- c(parameterNames[paramLines], "likelihood")
+ p<-list()
+ 
+ for(i in seq_along(colnames(preservedCalib)[-ncol(preservedCalib)])){
+     p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood"))+geom_point(size=0.9)
+ }
+
+ ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))),dpi = 600)
+ 
+ return(preservedCalib[preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"]),])
+}
+
+

RBBGCMuso/R/musoTime.R

@@ -25,7 +25,10 @@ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leap
     dates <- seq(from = as.Date(paste0(startYear,"01","01"),format = "%Y%m%d"), to =  as.Date(paste0(endYear,"12","31"),format = "%Y%m%d"), by = "day")
     if(leapYearHandling){
         if(prettyOut){
-            return(cbind(format(dates,"%d.%m.%Y"),as.numeric(format(dates,"%d")),as.numeric(format(dates,"%m")),as.numeric(format(dates,"%Y")))   )
+            return(cbind(format(dates,"%d.%m.%Y"),
+                         as.numeric(format(dates,"%d")),
+                         as.numeric(format(dates,"%m")),
+                         as.numeric(format(dates,"%Y")))   )
         }
         
         if(combined == FALSE){
@@ -35,10 +38,14 @@ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leap
         }
 
     } else {
+         dates <- dates[format(dates,"%m%d")!="0229"]
         if(prettyOut){
-            return(cbind(format(dates,"%d.%m.%Y"),as.numeric(format(dates,"%d")),as.numeric(format(dates,"%m")),as.numeric(format(dates,"%Y")))   )
+            return(data.frame(date = format(dates,"%d.%m.%Y"),
+                              day = as.numeric(format(dates,"%d")),
+                              month = as.numeric(format(dates,"%m")),
+                              year = as.numeric(format(dates,"%Y"))))
         }
-        dates <- dates[format(dates,"%m%d")!="0229"]
+       
 
         if(combined == FALSE){
             return(cbind(format(dates,"%d"),format(dates,"%m"),format(dates,"%Y")))
@@ -48,29 +55,48 @@ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leap
     }
     
 }
+#' alignData
+#'
+#' This function align the data to the model and the model to the data
+#' @importFrom lubridate leap_year
+#' @keywords internal
+alignData  <- function(mdata, dataCol, modellSettings = NULL, startDate, endDate, formatString = "%Y-%m-%d", leapYear = TRUE){
 
-corrigLeapYear  <- function(data, dataCol, modellSettings = NULL, startYear, fromDate = NULL,toDate = NULL,formatString = "%Y-%m-%d"){
-    data <- as.data.frame(data) 
-    numDays <- nrow(data)
-    dates <- seq(as.Date(paste0(startYear,"01","01"),format = "%Y%m%d"), by= "day", length = numDays)
-    goodInd <- which(!(leap_year(dates)&
-                       (format(date,"%m") == "12")&
-                       (format(date,"%d") == "31")))
-    realDate <- musoDate(startYear = startYear, numYears = numYears)
+    startDate <- as.Date(startDate, format = formatString)
+    endDate <- as.Date(endDate, format = formatString)
+    mdata <- as.data.frame(mdata)
     
-    data <- cbind.data.frame(real,data[goodInd])
-    
-    modellDates <- musoDate(startYear = settings$startYear,numYears = settings$numYears)
-    
-    
-
     if(is.null(modellSettings)){
         modellSettings <- setupMuso()
     }
-
-   
-}
-
-alignDataWithModelIndex  <- function(){
     
-} 
+    dates <- seq(startDate, to = endDate, by= "day")
+    if(!leapYear){
+      dates <- dates[which(format(dates,"%m%d") != "0229")]
+    }
+    mdata <- mdata[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d"),]
+    dates <- dates[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d")]
+    goodInd <- which(!(leap_year(dates)&
+                       (format(dates,"%m") == "12")&
+                       (format(dates,"%d") == "31")))
+    if(leapYear){
+        goodInd <- which(!(leap_year(dates)&
+                       (format(dates,"%m") == "12")&
+                       (format(dates,"%d") == "31")))
+    } else {
+        goodInd <-seq_along(dates)
+    }
+    realDate <- dates[which(format(dates,"%m%d") != "0229")]
+    if(leapYear){
+        mdata <- cbind.data.frame(realDate,mdata)        
+    } else {
+        mdata <- cbind.data.frame(dates,mdata)
+    }
+    modellDates <- as.Date(musoDate(startYear = modellSettings$startYear,numYears = modellSettings$numYears), format = "%d.%m.%Y")
+    mdata <- mdata[mdata[,1] %in% modellDates,]
+    nonEmpty <- which(!is.na(mdata[,dataCol+1]))
+    mdata <- mdata[nonEmpty,]
+    modIndex <- which(modellDates %in% mdata[,1])
+
+    list(measuredData = mdata[,dataCol +1], modIndex = modIndex)
+}

RBBGCMuso/man/alignData.Rd

@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/musoTime.R
+\name{alignData}
+\alias{alignData}
+\title{alignData}
+\usage{
+alignData(mdata, dataCol, modellSettings = NULL, startDate, endDate,
+  formatString = "\%Y-\%m-\%d", leapYear = TRUE)
+}
+\description{
+This function align the data to the model and the model to the data
+}
+\keyword{internal}

RBBGCMuso/man/calibMuso.Rd

@@ -26,8 +26,6 @@ keepEpc=FALSE, export=FALSE, silent=FALSE, aggressive=FALSE, leapYear=FALSE)
 
 \item{aggressive}{It deletes every possible modell-outputs from the previous modell runs.}
 
-\item{leapYear}{Should the function do a leapyear correction on the outputdata? If TRUE, then the 31.12 day will be doubled.}
-
 \item{keepBinary}{In default RBBGCMuso to keep  working area as clean as possible, deletes all the regular output files. The results are directly printed to the standard output, but you can redirect it, and save it to a variable, or you can export your results to the desired destination in a desired format. Whith this variable you can enable to keep the binary output files. If you want to set the location of the binary output, please take a look at the binaryPlace argument.}
 
 \item{binaryPlace}{The place of the binary output files.}
@@ -37,6 +35,8 @@ keepEpc=FALSE, export=FALSE, silent=FALSE, aggressive=FALSE, leapYear=FALSE)
 \item{skipSpinup}{If TRUE, calibMuso wont do spinup simulation}
 
 \item{prettyOut}{date ad Date type, separate year, month, day vectors}
+
+\item{leapYear}{Should the function do a leapyear correction on the outputdata? If TRUE, then the 31.12 day will be doubled.}
 }
 \value{
 No return, outputs are written to file

RBBGCMuso/man/dayOfMonths.Rd

@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/musoTime.R
-\name{dayOfMonths}
-\alias{dayOfMonths}
-\title{dayOfMonths}
-\usage{
-dayOfMonths(year, corrigated=TRUE)
-}
-\arguments{
-\item{year}{a year}
-
-\item{corrigated}{Do you want to handle the leapyears, if yes choose TRUE}
-}
-\value{
-vector with 12 element. First is January, the last is December. All of the vector element represents the number of the days in that specific month
-}
-\description{
-This function gives as a vector which contains the number of the days per each month
-}
-\keyword{internal}

RBBGCMuso/man/dayOfYears.Rd

@@ -1,15 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/musoTime.R
-\name{dayOfYears}
-\alias{dayOfYears}
-\title{This function tells us how many days are in the given year.}
-\usage{
-dayOfYears(year, corrigated = TRUE)
-}
-\description{
-This function tells us how many days are in the given year.
-}
-\author{
-Roland Hollos
-}
-\keyword{internal}

RBBGCMuso/man/isLeapyear.Rd

@@ -1,18 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/musoTime.R
-\name{isLeapyear}
-\alias{isLeapyear}
-\title{isLeapyear}
-\usage{
-isLeapyear(year)
-}
-\arguments{
-\item{year}{a year}
-}
-\value{
-TRUE, if leapyear, FALSE if dont.
-}
-\description{
-This function tells us if its argument a leapyear or not.
-}
-\keyword{internal}

RBBGCMuso/man/musoDate.Rd

@@ -2,29 +2,14 @@
 % Please edit documentation in R/musoTime.R
 \name{musoDate}
 \alias{musoDate}
-\title{It generates BiomeBGC-MuSo dates}
+\title{musoDate}
 \usage{
-musoDate(startYear, numYears, timestep = "d", combined = TRUE,
-  corrigated = TRUE, format = "en")
-}
-\arguments{
-\item{timestep}{timestep, which can be daily ("d"), monthly ("m"), yearly("y")}
-
-\item{combined}{If FALSE the output is a vector of 3 string: day, month year, if true, these strings will be concatenated.}
-
-\item{corrigated}{If True it counts with leapyears, else dont.}
-
-\item{format}{This is the format of the date. It can be "en" (dd.mm.yyyy), or "hu" (yyyy.mm.dd)}
-
-\item{settings}{You have to run the setupMuso function before musoDate. It is its output which contains all of the necessary system variables. It sets the whole environment.}
-}
-\value{
-The exact date-vectors for the BioBGC-MuSo output. You can use this for labelling purpose for example.
+musoDate(startYear, endYears = NULL, numYears, combined = TRUE,
+  leapYearHandling = FALSE, prettyOut = FALSE)
 }
 \description{
-It generates all of the day-dates which are between the start and endyear of BiomeBGC-MuSo run.
-How many days are from the given date and given period length(periodlen)?
+This function generates MuSo compatibla dates for the data
 }
 \author{
-Roland Hollos
+Roland HOLLOS
 }

RBBGCMuso/man/musoLeapYears.Rd

@@ -1,15 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/musoTime.R
-\name{musoLeapYears}
-\alias{musoLeapYears}
-\title{Musoleapyear}
-\usage{
-musoLeapYears(settings)
-}
-\description{
-How many days are from the given date and given period length(periodlen)?
-}
-\author{
-Roland Hollos
-}
-\keyword{internal}

RBBGCMuso/man/optiMuso.Rd

@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calibration.R
+\name{optiMuso}
+\alias{optiMuso}
+\title{optiMuso}
+\usage{
+optiMuso(measuredDataFile, parameters = NULL, sep = ",", startDate,
+  endDate, formatString, naString = NULL, leapYear = TRUE,
+  filterCol = NULL, filterVal = 1, selVar, outLoc = "./calib",
+  preTag = "cal-", settings = NULL, outVars = NULL,
+  iterations = 30, skipSpinup = TRUE, constrains = NULL,
+  plotName = "calib.jpg", likelihood = function(x, y) {    
+  exp(-sqrt(mean((x - y)^2))) }, calPar = 3009)
+}
+\arguments{
+\item{measuredDataFile}{a}
+
+\item{parameters}{b}
+
+\item{sep}{c}
+
+\item{startDate}{d}
+
+\item{endDate}{e}
+
+\item{formatString}{a}
+
+\item{leapYear}{b}
+
+\item{filterCol}{a}
+
+\item{filterVal}{b}
+
+\item{selVar}{c}
+
+\item{outLoc}{c}
+
+\item{settings}{e}
+
+\item{iterations}{c}
+
+\item{skipSpinup}{a}
+
+\item{constrains}{d}
+
+\item{plotName}{u}
+
+\item{likelihood}{d}
+
+\item{calPar}{a}
+
+\item{pretag}{a}
+}
+\description{
+This function calculates the -users specified- likelihood for random model input.
+}
+\author{
+Roland HOLLOS
+}

RBBGCMuso/man/readMeasuredMuso.Rd

@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/assistantFunctions.R
+\name{readMeasuredMuso}
+\alias{readMeasuredMuso}
+\title{readMeasuredMuso}
+\usage{
+readMeasuredMuso(inFile, naString = getOption("datatable.na.strings",
+  "NA"), sep = ",", leapYearHandling = TRUE, convert.var = NULL,
+  convert.scalar = 1, convert.fun = (function(x) {     x *
+  convert.scalar }), convert.file = NULL, filterCol = NULL,
+  filterVal = 1, selVar = NULL)
+}
+\description{
+MuSo data reader
+}

RBBGCMuso/man/sumDaysOfPeriod.Rd

@@ -1,15 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/musoTime.R
-\name{sumDaysOfPeriod}
-\alias{sumDaysOfPeriod}
-\title{How many days are from the given date and given period length(periodlen)?}
-\usage{
-sumDaysOfPeriod(year, periodlen, corrigated = TRUE)
-}
-\description{
-How many days are from the given date and given period length(periodlen)?
-}
-\author{
-Roland Hollos
-}
-\keyword{internal}