New function for erosion calculation

R/solvers.R

@@ -110,6 +110,118 @@ solv_conc_lag <- function(t,z,C0,Psp0,Pmu0,lambda,S,final=FALSE){
 
 
 
+#######################################################################
+#' Compute steady state denudation
+#' Eulerian formulation
+#'
+#' @param Cobs Measured concentration (at/g)
+#' @param p production and decay parameters (4 elements vector)
+#' p[1] -> unscaled spallation production rate (at/g/a)
+#' p[2] -> unscaled stopped muons production rate (at/g/a)
+#' p[3] -> unscaled fast muons production rate (at/g/a)
+#' p[4] -> decay constant (1/a)
+#' @param S scaling factors (2 elements vector)
+#' S[1] -> scaling factor for spallation
+#' S[2] -> scaling factor for muons
+#' @param L attenuation length (3 elements vector)
+#' L[1] -> neutrons
+#' L[2] -> stopped muons
+#' L[3] -> fast muons
+#' @param Cobs Uncertainty on measured concentration (at/g), optional but required ffor method MC and MCMC
+#' @param method One of single (default), MC, MCMC
+#' @n number of draws for MC and MCMC methods (default 10000)
+#'
+#' @return
+#' @export
+#'
+#' @examples
+solv_ss_eros_eul <- function(Cobs,p,S,L,Cobs_e=0,method="single",n=10000){
+  Cmax = solv_conc_eul(0,0,Inf,0,p,S,L)
+  if (Cobs>Cmax){stop("Observed concentration higher than theoretical maximum")}
+  if (!(method %in% c("single","MC","MCMC"))){stop("Parameter method should be one of single, MC or MCMC")}
+  if (method!="single" & Cobs_e==0){stop("Need to define the error on the observed concentration")}
+  if (method == "single") {
+    res = uniroot(fun_solv_ss_eros_eul,c(0,1),Cobs,p,S,L,tol=1e-5)
+  } else if (method == "MC") {
+    C0 = rtruncnorm(n,a=0,mean=Cobs,sd=Cobs_e)
+    ero = rep(NA,n)
+    val = rep(NA,n)
+    iter = rep(NA,n)
+    prec = rep(NA,n)
+    for (i in 1:n){
+      a = uniroot(fun_solv_ss_eros_eul,c(0,1),C0[i],p,S,L,tol=1e-6) # attention à la borne sup de l'interval que le taux d'érosion max soit assez rapide pour les cocentrations faibles
+      ero[i] = a$root
+      val[i] = a$f.root
+      iter[i] = a$iter
+      prec[i] = a$estim.prec
+    }
+    res = data.frame(C,ero,val,iter,prec)
+  } else if (method == "MCMC") {
+    ero0 = uniroot(fun_solv_ss_eros_eul,c(0,0.1),Cobs,p,S,L)$root
+    binf = max(0,ero0*(1-10*Cobs_e/Cobs))
+    bsup = ero0*(1+10*Cobs_e/Cobs)
+    res = run_mcmc1(p,S,L,Cobs,Cobs_e,binf,bsup,n)
+
+  }
+  return(res)
+}
+
+#' Function to be solved by solv_ss_eros_eul when method="single"
+#' @param ero
+#' @param Cobs
+#' @param p
+#' @param S
+#' @param L
+#'
+#'
+fun_solv_ss_eros_eul <- function(ero,Cobs,p,S,L){
+  C = solv_conc_eul(0,ero,4.55e9,0,p,S,L)
+  return((C-Cobs)/Cobs)
+}
+
+#' Likelihood function : P(obs|params) for solv_ss_eros_eul when method="MCMC"
+likelihood1 <- function(ero,p,S,L,Cobs,Cobs_e){
+  C = solv_conc_euler(0,ero,4.55e9,0,p,S,L)
+  chi2 = (C-Cobs)^2/Cobs_e^2
+  res = (-1*chi2/2) + log(1/(sqrt(2*pi)*Cobs_e))
+  return(res)
+}
+#' Prior distribution  P(params) for solv_ss_eros_eul when method="MCMC"
+prior1 <- function(ero,binf,bsup){
+  return(dunif(ero, min=binf, max=bsup, log = T)) # log prior
+}
+#' Un-normalized posterior : numerator of Bayes formula : P(obs|params)*P(params) for solv_ss_eros_eul when method="MCMC"
+##->  likelihood*prior (but we work with log)
+posterior1 <- function(ero,p,S,L,Cobs,Cobs_e,binf,bsup){
+  return (likelihood1(ero,p,S,L,Cobs,Cobs_e) + prior1(ero,binf,bsup)) # log posterior
+}
+#' Proposal function for solv_ss_eros_eul when method="MCMC"
+proposal1 <- function(ero,binf,bsup){
+  return(rtruncnorm(1,a=binf,b=bsup,mean=ero,sd=(bsup-binf)/20))
+}
+#' RUN MCMC function for solv_ss_eros_eul when method="MCMC"
+run_mcmc1 <- function(p,S,L,Cobs,Cobs_e,binf,bsup,n){
+  # start chain
+  chain = as.data.frame(matrix(NA,nrow=n,ncol=3))
+  colnames(chain) <- c("ero","posterior","probab")
+  chain$ero[1] = proposal1(runif(1,binf,bsup),binf,bsup)
+  chain$posterior[1] = posterior1(chain$ero[1],p,S,L,Cobs,Cobs_e,binf,bsup)
+  for (i in 1:(n-1)){
+    #    if (i%%(round(iter*1/100)+1)==0) {cat(round(i/iter*100),"%  -")}
+    proposal = proposal1(chain$ero[i],binf,bsup)
+    pst = posterior1(proposal,p,S,L,Cobs,Cobs_e,binf,bsup)
+    probab = exp(pst - chain$posterior[i]) # ratio of posterior probability proposal/current
+    chain$probab[i]=probab
+    if (runif(1) < probab){
+      chain$ero[i+1] = proposal
+      chain$posterior[i+1] = pst
+    }else{
+      chain$ero[i+1] = chain$ero[i]
+      chain$posterior[i+1] = chain$posterior[i]
+    }
+  } # end loop
+  return(chain)
+} # end function run_mcmc