RxpsG_2.3-1: Processing Tool for X-ray Photoelectron Spectroscopy Data

Documented in XPSmodFit

## ===========================================================
## XPSfit function of Object = XPS coreline
## call the XPSdofit AND plot the result
## ===========================================================

#' @title XPSmodFit performs the fit of a Core Line
#' @description Provides a userfriendly interface to select a fitting algorithms to
#'   fit an XPSCoreline. Fitting algorithms are:
#'   (I) Classic algorithms: Levenberg-Marquardt, Newton and Port based on the minimization
#'     of the Squares of the Differences and handled in a more robust way (although slower)
#'     compared to (\code{nlsLM}) function.
#'   (II) Conjugate-Gradient Algorithms: General-purpose optimization based on Nelder-Mead,
#'     quasi-Newton and conjugate-gradient algorithms.
#'   (III) Pseudo Algorithms: random based algorithm (rather slow) but ensure the convergence.
#' @param Object XPSCoreLine object
#' @param plt logical, to enable plot XPSCoreLine and fit. By default is set to TRUE
#' @param \dots  further parameters to the fitting function
#' @return The Object slot \code{XPSmodFit} will be filled with the result of the
#'  calculation. All the values of the components will be update and the result
#'  will be displyed.
#' @seealso \link{nlsLM}
#' @examples
#' \dontrun{
#' 	SampData[["C1s"]] <- XPSmodFit(SampData[["C1s"]])
#' }
#' @export
#'

XPSmodFit <- function(Object, plt=TRUE, ...) {
       import::from(FME, modFit)
       import::from(rootSolve,gradient)    #cannot use import::here because rootSolve not listed in DESCRIPTION
       assign("gradient", gradient, envir=.GlobalEnv)  #renders the function 'gradient' globally available

# ==============================================================
# showListParam: formatted print of a list of parameters
# ==============================================================

   showListParam <- function(lista, decimals, ParamNames){
       maxLen <- 15  #max length of the list
       ParamNames <- unlist(ParamNames)
       ParamNames <- sapply(ParamNames, function(x) encodeString(x, width=15, justify="right")) #max element width=15 and align text of list elements on the rigth
       cat("\n", ParamNames)
#--- 17 = max length of list element
       lista <- lapply(lista, function(x) round(x, digits=decimals)) #rounds the numbers of decimals
       lista <- lapply(lista, function(x) as.character(x)) # transform each list element in string
       lista <- lapply(lista, function(x) encodeString(x, width=17, justify="right")) #max element width=17 and align text of list elements on the rigth
       NN <- length(lista)
       for (ii in 1:NN){
           lista[[ii]] <- substr(lista[[ii]], 4, nchar(lista[[ii]]))  #drop the first element of the list (component identifier)
           xx <- paste(lista[[ii]], sep="", collapse="")
           xx <- paste("C", ii, ": ", xx, sep="")    #reformat component name
           cat("\n", xx)
       }
   }

# ===========================================================
# XPSdofit: performs the FIT on the original spectrum - background
#
#--- FitFunct() estimates the fitting function using the fit parameters. It is needed for the calculation of the residuals
#      FitFunct(XX, Param, FitExpr)
#       XX = vector of abscissas
#       Param = string containing fit parameters such as: "h1 <- 12569.234; mu1 <- 284.27;, sigma1 <- 1.56....
#       FitExpr = string describing the fitting function such as: Gauss(x, h1, mu1, sigma1)+Lorentz(x, h2, mu2, sigma2)+...
#--- FitResiduals() usinf the actual fitting parameters computes the difference between fitting function and original data
#       FitResiduals(Parms)             FitResiduals transform numerical parameters in string needed to evaluate the fitting function with FitFunct()
#       Parms = numeirc array numeric of fitting parametrs
#--- modFit() minimized the residuals between fitting function and original data modifying the fitting parameters
#       modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Pseudo")
#       f= vectors of residuals generated by FitResiduals()
#       p= Parms vector containing the fitting parameters
#       lower= lower limit for the fitting parameters
#       upper= upper limit for the fitting parameters
#       method= string containinf the selected fitting method
#
# ===========================================================

   XPSdofit <- function(Object, method, ...) {
#--- data to fit: curve - baseline
 	     datafit <- data.frame(x = Object@RegionToFit$x,
 				  y = Object@RegionToFit$y - Object@Baseline$y)
#--- FitExpr generates the fit expression for each of the fitting components
#--- generates the fitting parameters
	      FitExpr <- sapply(names(Object@Components), function(x) {
		           number <- unlist(strsplit(x, "C"))[2]
		           FuncName <- slot(Object@Components[[x]],"funcName")
		           Fargs <- formalArgs(FuncName)
		           lenvar <- length(Fargs)-1
#--- add index to each of the fitting parameters
		           VarNames <- c("", rep.int(number, lenvar))
		           Fparm <- paste(Fargs,VarNames, sep="")
#--- fitting formula is the ensamble of the fitting component functions
		           funct <- paste(paste(FuncName, "(", sep=""), paste(Fparm, collapse=","),")", sep="")

#--- modify fitting functions if links on fitting parameters are present
		           if (length(Object@Components[[x]]@link) ) {
                    for (idx in seq_along(slot(Object@Components[[x]],"link"))) {
   			                  funct <- sub(Object@Components[[x]]@link[[idx]]$variable, Object@Components[[x]]@link[[idx]]$expr, funct)
                    }
   		        }
		           return(funct)
	      })

       FitExpr <- paste(FitExpr, collapse= "+")

#--- parameters values
	      Parms <- lapply(Object@Components, function(x) getParam(x, parameter="start"))
	      Ubounds <- lapply(Object@Components, function(x) getParam(x, parameter="max"))
	      Lbounds <- lapply(Object@Components, function(x) getParam(x, parameter="min"))
#--- extract names of the fitting parameters
	      ParamNames <- lapply(seq_along(Object@Components), function(x,y) paste(rownames(y[[x]]@param),x,sep=""), y = Object@Components )

       cat("\n-------------------------------------------------")
       cat("\n\n START:")
       showListParam(Parms, 3, ParamNames[[1]])
       cat("\n\n UpLim:")
       showListParam(Ubounds, 3, ParamNames[[1]])
       cat("\n\n LwLim:")
       showListParam(Lbounds, 3, ParamNames[[1]])
       cat("\n-------------------------------------------------")
       cat("\n\n")

#--- #controls if a FIX constraint is set on fitting parameters
       Ncomp <- length(Object@Components)
       for(ii in 1:Ncomp){
   	      Npar <- length(Parms[[ii]])  #Componants may have different fitting functions with different number of parameters
          for(jj in 1:Npar){
              #if FIX present => Ubound==Lbound but in modFit this does not work
              if ((Ubounds[[ii]][jj] == Lbounds[[ii]][jj]) && (jj != 3)) { #here we run on all the parameters except sigma (jj=3)
                  Ubounds[[ii]][jj] <- Parms[[ii]][jj]+Parms[[ii]][jj]/10000 #Ubound~=LBound
                  Lbounds[[ii]][jj] <- Parms[[ii]][jj]-Parms[[ii]][jj]/10000
              }
          }
       }

#--- create a sequence containing all the fitting parameters
       Parms <- unlist(Parms)
       Ubounds <- unlist(Ubounds)
       Lbounds <- unlist(Lbounds)
#--- associate the correspondent parameter names
       ParamNames <- unlist(ParamNames)
	      names(Parms) <- ParamNames
	      names(Ubounds) <- ParamNames
	      names(Lbounds) <- ParamNames
#--- create a backup of the initial parameter values valori di start: questo vettore viene poi riempito
#--- start values of parameters are updated during fitting procedure
	      tmpfit <- Parms

#--- index = vector identifying the components where a link is present on a fit paramenetr
	      index <- NULL
	      for(number in seq_along(Object@Components)) {
	          if (length(Object@Components[[number]]@link)) {
	        	     for(idx in seq_along(slot(Object@Components[[number]],"link"))) {
		                index <- c(index,Object@Components[[number]]@link[[idx]]$position)
	        	     }
	          }
	      }
#--- eliminates the parameters which are linked
	      if (! is.null(index) ) {
		        Parms <- Parms[-c(index)]
		        Ubounds <- Ubounds[-c(index)]
		        Lbounds <- Lbounds[-c(index)]
	      }

#--- ParamName: mu1 = mu of FitComp 1,   sigma3 = sigma of FitComp 3,  h5 = h of FitComp 5.
#--- generates a string of type: mu1 <- 285.34;   sigma3 <- 1.4    ecc..
	     Param <- paste(ParamNames, " <- ", Parms, " ;", sep="")

      LL <- length(datafit$x)
#spectral data decimation: extract indexes of original data to use for fitting. If mesh=2 one value is taken one is drop
      indexes <- seq(from=1, to=LL, by=mesh)
      XX <- datafit$x[indexes]
      YY <- datafit$y[indexes]


#--- evaluate the fitting function using the fitting Parameters
      FitFunct <- function(XX, Param, FitExpr) {
                          eval(parse(text=Param))  #eval creates in the .GlobalEnv the fitting parameters having the correspondent values
                          eval(parse(text="x <- XX")) # FitExpr = fit function(x) needs the whole vector XX be available
                          estimate <- eval(parse(text=FitExpr))
                          return(as.vector(estimate))
                       }

#--- estimate the residuals as difference of the FitFunction and the original spectral data
      FitResiduals <- function(Parms) { #function to be minimized by modFit
                              Param <- paste(names(Parms), " <- ", Parms, " ;", sep="") #String containing FitParam evaluated by FitFunction
                              Fitting <- FitFunct(XX, Param, FitExpr) #fitting evaluates the fitting functon using the Param
                              residuals <- (YY - Fitting)
                              return(residuals)
                       }

#--- FITTING ALGORITHMS
     if (method=="Marq") {
        if (is.na(MaxIter)) MaxIter <<- 10000
        ctrl <- list(ptol= Tolerance, maxiter=MaxIter, nprint=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Marq", control=ctrl)
     }
     if (method=="Newton") {
        if (is.na(MaxIter)) MaxIter <<- 10000
        ctrl <- list(ptol= Tolerance, maxiter=MaxIter, nprint=traceFit)
        cat("\n Iterations are starting please wait")
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Newton", control=ctrl)
     }
     if (method=="Port") {
        if (is.na(MaxIter)) MaxIter <<- 200
        ctrl <- list(rel.tol=Tolerance, eval.max=200, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Port", hessian=TRUE, control=ctrl)
     }
     if (method=="Nelder-Mead") {
        if (is.na(MaxIter)) MaxIter <<- 200
        ctrl <- list(reltol=Tolerance, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Nelder-Mead", control=ctrl)
     }
     if (method=="CG") {
        if (is.na(MaxIter)) MaxIter <<- 200
        ctrl <- list(reltol=Tolerance, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="CG", control=ctrl)
     }
     if (method=="BFGS") {
        if (is.na(MaxIter)) MaxIter <<- 200
        ctrl <- list(reltol=Tolerance, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="BFGS", control=ctrl)
     }
     if (method=="L-BFGS-B") {
        if (is.na(MaxIter)) MaxIter <<- 200
        ctrl <- list(reltol=Tolerance, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="L-BFGS-B", control=ctrl)
     }
     if (method=="SANN") {
        if (is.na(MaxIter)) MaxIter <<- 200
        if (Verbose) traceFit=1
        ctrl <- list(reltol=Tolerance, iter.max=MaxIter, trace=traceFit)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="SANN", control=ctrl)
     }
     if (method=="Pseudo") {
        if (is.na(MaxIter)) MaxIter <<- 5000
        ctrl <- list(varleft=Tolerance, numiter=MaxIter, verbose=Verbose)
        FitEstimation <- modFit(f = FitResiduals, p = Parms, lower=Lbounds, upper=Ubounds, method="Pseudo", control=ctrl)
     }

     #--- SUMMARY of fit model and PLOT FIT RESULT
     cat("\n ----Best Fit Param.----\n")
     print(FitEstimation$par)
     cat("\n -----------------------\n")

     Param <- paste(names(FitEstimation$par), " <- ", FitEstimation$par, " ;", sep="") # qui metto i parametri in una stringa che valutero' in fitFunct
     Fit <- FitFunct(datafit$x, Param, FitExpr)

	    FitParam <- FitEstimation$par
#--- Update the FitParam values in the tmpfit vector
	    for (nomi in names(FitParam) ){
		       num <- grep(nomi, ParamNames)
		       tmpfit[num] <- FitParam[nomi]
	    }
#--- if there are links modify correspondent values
	    for (idx in seq_along(Object@Components) ) {
	       for ( jj in seq_along(Object@Components[[idx]]@link) ) {
		         FUN <- Object@Components[[idx]]@link[[jj]]$FUN
		         if ( ! is.na(FUN) ) {
		            origname <- Object@Components[[idx]]@link[[jj]]$newvar # mu1
		            num <- grep(origname, names(FitParam))
		            origvalue <- FitParam[num]
		            newvalue <- Object@Components[[idx]]@link[[jj]]$value
		            value <- sapply(origvalue,FUN,newvalue)
		            names(value) <- Object@Components[[idx]]@link[[jj]]$variable
		         } else {
		            origname <- Object@Components[[idx]]@link[[jj]]$expr # sigma1
  		          num <- grep(origname, names(FitParam))
              value <- FitParam[num]
		            names(value) <- Object@Components[[idx]]@link[[jj]]$variable
		         }
#--- link[[jj]]$position indicates the position of the actual parameter in FitParam
           index <- Object@Components[[idx]]@link[[jj]]$position
           tmpfit[index] <- value
	       }
	    }
     for ( idx in seq_along(Object@Components) ) {
        num1 <- grep(idx, names(tmpfit)) # Attention: this greep for idx=1 gets h1, mu1, sigma1... but also: h10, mu10, sigma10..., h11, mu11, sigma11... ecc.
        LL=length(num1)
        num2 <- gsub("[^0-9]", "", names(tmpfit))  #separates the characters form numbers in names(tmpfit)

        for (jj in LL:1) { #run on all the elements of num1 extracted by greep and selects only those equal to idx
            if (idx!= num2[num1[jj]]) { num1 <- num1[-jj] }  #drop all the elements of num1 not equal to idx
        }

        Object@Components[[idx]] <- setParam(Object@Components[[idx]],
                                             parameter="start",
                                             value=tmpfit[num1]) # start
#--- Fit components computation
	       Object@Components[[idx]] <- Ycomponent(Object@Components[[idx]], x=Object@RegionToFit$x, y=Object@Baseline$y)
	  }
#--- update slot Fit
        Object@Fit$y <- Fit
        Object@Fit$fit <- FitEstimation
        if (plt == TRUE){
            plot(Object)
            lines(datafit$x, Fit+Object@Baseline$y, lwd = 1, col = "red")
        }
        return(Object)
   }


#===== variables =====
   FittedObject <- Object  #creates a copy of original data
   mesh <- 1
   MaxIter <- NA
   Tolerance <- 1e-9
   Verbose <- TRUE
   traceFit <- 1


#=====  MAIN =========

   Fwin <- gwindow("XPS MODEL FITTING", parent=c(50, 10), visible=FALSE)
   FgroupMain <- ggroup(horizontal=FALSE, container = Fwin)

   Fgroup1 <- ggroup(horizontal=TRUE, container = FgroupMain)
   Fframe1 <- gframe("Data Decimation", spacing=5, container=Fgroup1)
   Fobj1 <- gcombobox(c("default","1","2","3","4","5","6","7","8","9","10"), selected=1, editable=FALSE, handler=function(h, ...){
                              mesh <<- (svalue(Fobj1))
                              if (mesh == "default") mesh <<- 1
                              if (mesh != "default") mesh <<- as.numeric(mesh)
                              LL <- length(Object@RegionToFit$x)
                              indices <- seq(from=1, to=LL, by=mesh)
                              XX <- Object@RegionToFit$x[indices] #if mesh > 1 data decimation
                              YY <- Object@RegionToFit$y[indices]
#                              BL <- Object@Baseline$y[indices]
                              if (plt==TRUE){
                                  plot(Object)
                                  points(XX, YY, type="p", pch=16, cex=0.6, col="limegreen")
                              }
                        }, container=Fframe1)

   Fframe2 <- gframe("Max Iterations", spacing=5, container=Fgroup1)
   Fobj2 <- gcombobox(c("default","100","200","500","1000","3000","5000","10000"), selected=1, editable=FALSE, handler=function(h, ...){
                              MaxIter <<- svalue(Fobj2)
                              if (MaxIter == "default") MaxIter <<- NA
                              if (MaxIter != "default") MaxIter <<- as.numeric(MaxIter)
                        }, container=Fframe2)

   Fframe3 <- gframe("Tolerance", spacing=5, container=Fgroup1)
   Fobj3 <- gcombobox(c("default", "1e-3","1e-4","1e-5","1e-6","1e-7","1e-8", "1e-9"), selected=1, editable=FALSE, handler=function(h, ...){
                              Tolerance <<- svalue(Fobj3)
                              if (Tolerance == "default") Tolerance <<- 1e-9
                              if (Tolerance != "default") Tolerance <<- as.numeric(Tolerance)
                        }, container=Fframe3)

   Fframe4 <- gframe("Verbose", spacing=5, container=Fgroup1)
   Fobj4 <- gcombobox(c("Yes", "No"), selected=1, editable=FALSE, handler=function(h, ...){
                              Verbose <- svalue(Fobj4)
                              if (Verbose == "Yes"){
                                 Verbose <<- TRUE
                                 traceFit <<- 1
                              } else {
                                 Verbose <<- FALSE
                                 traceFit <<- -1
                              }
                        }, container=Fframe4)

   mesh <<- MaxIter <<- Tolerance <<- NA

   glabel("Decimation accelerates fits but decreases precision;", container=FgroupMain)
   glabel("MaxIter: 10000 for classic Algorithms;", container=FgroupMain)
   glabel("         150 for Nelder-Mead and CG, 500 for the other Conjugate-gradient Algorithms;", container=FgroupMain)
   glabel("         500 for SANN, 10000 for Pseudo Algorithm;", container=FgroupMain)
   glabel("Tolerance: 1e-8 is the default;", container=FgroupMain)


   Fgroup5 <- ggroup(container = FgroupMain)
   Fframe5 <- gframe("Classic Algorithms", horizontal=TRUE, spacing=5, container=Fgroup5)
   gbutton("Lev.Marq.", handler=function(...){
                              FitModel <- "Marq"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe5)
   gbutton("Newton", handler=function(...){
                              FitModel <- "Newton"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe5)
   gbutton("Port", handler=function(...){
                              FitModel <- "Port"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe5)

   Fgroup6 <- ggroup(container = FgroupMain)
   Fframe6 <- gframe("Conjugate-gradient Algorithms", horizontal=TRUE, spacing=5, container=Fgroup6)
   gbutton("Nelder-Mead", handler=function(...){
                              FitModel <- "Nelder-Mead"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe6)
   gbutton("ConjGrad", handler=function(...){
                              FitModel <- "CG"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe6)
   gbutton("BFGS", handler=function(...){
                              FitModel <- "BFGS"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe6)
   gbutton("L-BFGS-B", handler=function(...){
                              FitModel <- "L-BFGS-B"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe6)

   Fgroup7 <- ggroup(container = FgroupMain)
   Fframe7 <- gframe("Random Algorithms", horizontal=TRUE, spacing=5, container=Fgroup7)
   gbutton("SANN", handler=function(...){
                              FitModel <- "SANN"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe7)
   gbutton("Pseudo", handler=function(...){
                              FitModel <- "Pseudo"
                              FittedObject <<- XPSdofit(Object, FitModel,  ...)
                        },container=Fframe7)

   Fgroup8 <- ggroup(horizontal=TRUE, container = FgroupMain)
   gbutton("SAVE & EXIT", handler=function(...){
                              SpectIndx <- get("activeSpectIndx", envir=.GlobalEnv)
                              activeFName <- get("activeFName", envir=.GlobalEnv)
                              FName <- get(activeFName, envir=.GlobalEnv)
                              FName[[SpectIndx]] <- FittedObject
                              assign(activeFName, FName, envir=.GlobalEnv)
                              if (plt == TRUE){
                                  plot(FittedObject)
                              }
                              Object <<- FittedObject
                              dispose(Fwin)
                              XPSSaveRetrieveBkp("save")
                        },container=Fgroup8)

   gbutton("RESET", handler=function(...){
                              mesh <<- 1   #the default values of variables is set in the ctrl list for all the algorithms
                              MaxIter <<- NA
                              Tolerance <<- 1e-9
                              traceFit <<- 1
                              svalue(Fobj1) <- "default"
                              svalue(Fobj2) <- "default"
                              svalue(Fobj3) <- "default"
                              svalue(Fobj4) <- "Yes"
                              SpectIndx <- get("activeSpectIndx", envir=.GlobalEnv)
                              activeFName <- get("activeFName", envir=.GlobalEnv)
                              FName <- get(activeFName, envir=.GlobalEnv)
                              Object <<- FName[[SpectIndx]]
                              FittedObject <<- NULL
                              if (plt == TRUE){
                                  plot(Object)
                              }
                        },container=Fgroup8)

   gbutton("EXIT", handler=function(...){
                              dispose(Fwin)
                              XPSSaveRetrieveBkp("save")
                        },container=Fgroup8)

   visible(Fwin) <- TRUE
   Fwin$set_modal(TRUE)  #set_modal == TRUE: no other operation allowed until modFit active

   return(Object)
}