Nothing
R/fits.R
In BBEST: Bayesian Estimation of Incoherent Neutron Scattering Backgrounds
Defines functions
do.iter
do.fit.banks
do.fit
set.DEoptim.control
set.control
Documented in
do.fit
do.fit.banks
do.iter
set.control
set.DEoptim.control
#####################################################################################
#
#
# FUNCTIONS TO PERFORM BACKGROUND FIT VIA DIFFERENTIAL EVOLUTION ALGORITHM
#
#
#####################################################################################
# set.control(CR, F, NP, itermax, parallelType)
#
# returns: wrapper to set.DEoptim.control
# arguments
# CR: crossover probability from interval [0,1]
# F: differential weighting factor from interval [0,2]
# NP: number of population members
# itermax: the maximum iteration (population generation) allowed
# parallelType: defines the type of parallelization to employ.
set.control <- function(CR=.85, F=.7, NP=300, itermax=2000, parallelType=1){
control <- list()
control$CR <- CR
control$F <- F
control$NP <- NP
control$itermax <- itermax
control$parallelType <- parallelType
return(control)
}
#####################################################################################
# set.DEoptim.control(control, knots.n)
#
# returns: DEoptim control parameters
# arguments
# control: return value of set.control
# knots.n: knots number
set.DEoptim.control <- function(control=list(), knots.n){
if(is.null(control$CR))
CR <- 0.85
else
CR <- control$CR
if(is.null(control$F))
F <- 0.7
else
F <- control$F
if(is.null(control$NP))
NP <- round(2.5*knots.n*10)
else
NP <- control$NP
if(is.null(control$itermax))
itermax <- 2000
else
itermax <- control$itermax
if(is.null(control$parallelType))
parallelType <- 1
else
parallelType <- control$parallelType
DE.control <- DEoptim.control(CR=CR, F=F, NP=NP, itermax=itermax, parallelType=parallelType,
# packages = list("PerformanceAnalytics"),
parVar=list("basisMatrix", "basisSpline",
"logPosterior", "logPriorF", "DMatrix", "logLikelihoodBkg", "logLikelihoodSignal",
"get.bkg", "logProbabilityBkgR", "logLikelihoodGrGauss", "logLikelihoodGrCorr",
"logPriorBkgRSmooth", "logPriorBkgRGP", "Dx", "invert.order", "set.SB",
"logPosteriorAnalyt", "bkg.analyt", "logLikelihoodBkgAnalyt", "logLikelihoodSignalAnalyt"))
DE.control
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, scale, knots.x, knots.n, stdev, control, p.bkg, save.to)
#
# returns: Performs evolutionary global background optimization via DEoptim
# arguments
# data: object of type data. Contatins experimental data and fit parameters
# bounds.lower, lower and upeer boundaries for background
# bounds.upper:
# scale: scale factor
# knots.x: the spline knot x-values.
# knots.n: number of knots.
# stdev: Whether to calculate uncertainty in background estimation
# control: he return value of set.control
# p.bkg: the probability that a single pixel contains "only" background
# save.to: the name of the file where the results will be saved
do.fit <- function(data, bounds.lower, bounds.upper, scale=c(1,1), knots.x=NA,
knots.n=NA, analytical=FALSE, stdev=TRUE, control=list(), p.bkg=.5, save.to=""){
# 1. Prepare data to fit...
cat("Preparing data... \n")
if( (scale[1]==1) && (scale[2]==1) ) scale <- NA
cov.r <- ADP <- knots.y <- pars <- NA
alpha <- 1
if(is.null(data$Gr))
Gr <- NA
else{
Gr <- data$Gr
data$Gr <- NA
}
if(is.null(data$SB))
data$SB <- rep(0, length(data$x))
else if (is.na(data$SB[1]))
data$SB <- rep(0, length(data$x))
# 2. prepare fit params...
if(is.na(knots.x[1]) && !is.na(knots.n))
knots.x <- seq(min(data$x), max(data$x), length=knots.n)
if(is.na(knots.n))
knots.n <- length(knots.x)
if(analytical==TRUE){
iter.0 <- as.relistable(list(pars=rep(0, 6)))
bounds.lower <- c(pars=c(-abs(bounds.lower[1]*10), -5, -5, -20, 0, 0))
bounds.upper <- c(pars=c(abs(bounds.upper[1]*10), 10, 5, 2000,
max(data$x), max(data$x)*2))
knots.n <- 6
}
else{
iter.0 <- as.relistable(list(knots.y=rep(0, knots.n)))
bounds.lower <- c(knots.y=rep(bounds.lower[1], knots.n))
bounds.upper <- c(knots.y=rep(bounds.upper[1], knots.n))
}
if(!is.na(scale[1])){
iter.0$alpha <- 0
bounds.lower <- c(bounds.lower, alpha=scale[1])
bounds.upper <- c(bounds.upper, alpha=scale[2])
}
if(!is.null(data$fitADP)){
if(data$fitADP$oneADP==TRUE)
nn <- 1
else
nn <- length(data$fitADP$n.atoms)
iter.0$ADP <- rep(0, nn)
bounds.lower <- c(bounds.lower, ADP=rep(data$fitADP$ADP.lim[1], nn))
bounds.upper <- c(bounds.upper, ADP=rep(data$fitADP$ADP.lim[2], nn))
}
DE.control <- set.DEoptim.control(control, knots.n)
# 3. Starting fit!
#
# initial guess
# cc <- matrix(nrow=DE.control$NP, ncol=length(knots.x))
# Phi <- basisMatrix(data$x, knots.x)
# for(ii in 1:DE.control$NP){
# spar=(0.9-0.5)/DE.control$NP*ii+0.5
# lowpass.spline <- smooth.spline(data$x,data$y-data$SB, spar = spar) ## Control spar for amount of smoothing
# cc[ii,] <- c(solve(t(Phi)%*%Phi) %*% t(Phi) %*% (predict(lowpass.spline, data$x)$y) )
# }
# DE.control$initialpop=cc
cat("Starting DifEv algorithm... \n")
doPbkgIter <- FALSE
if(p.bkg==-1){
doPbkgIter <- TRUE
p.bkg <- 0.02
}
DEoptim.fit <- DEoptim(get.posterior, lower = bounds.lower, upper = bounds.upper,
control= DE.control, skel=iter.0, data=data,
knots.x=knots.x, Gr=Gr, p.bkg=p.bkg)
if(doPbkgIter){
cat("\n\n P.bkg iteration... \n\n")
if(analytical==TRUE){
pars <- DEoptim.fit$optim$bestmem[1:6]
bkg <- bkg.analyt(pars=pars, x=data$x)
}
else{
knots.y <- DEoptim.fit$optim$bestmem[1:knots.n]
bkg <- get.bkg(x=data$x, knots.x=knots.x, knots.y=knots.y)
}
dev.norm <- (data$y-bkg-data$SB)/data$sigma
p.bkg <- 0.5*50^((1-dev.norm^2)/8)
DEoptim.fit <- DEoptim(get.posterior, lower = bounds.lower, upper = bounds.upper,
control= DE.control, skel=iter.0, data=data,
knots.x=knots.x, Gr=Gr, p.bkg=p.bkg)
}
if(analytical==TRUE){
pars <- DEoptim.fit$optim$bestmem[1:6]
bkg <- bkg.analyt(pars=pars, x=data$x)
}
else{
knots.y <- DEoptim.fit$optim$bestmem[1:knots.n]
bkg <- get.bkg(x=data$x, knots.x=knots.x, knots.y=knots.y)
}
if(!is.null(data$fitADP)){
ADP <- tail(DEoptim.fit$optim$bestmem, nn)
data <- set.SB(data, n.atoms=data$fitADP$n.atoms,
scatter.length=data$fitADP$scatter.length, ADP=ADP)
}
data$Gr <- Gr
if(!is.na(scale[1])){ # recalculating signal and params...
alpha <- DEoptim.fit$optim$bestmem[knots.n+1]
data$y <- (data$y-bkg-data$SB)*alpha + data$SB + bkg
data$lambda <- data$lambda*alpha
data$sigma <- data$sigma*alpha
if(!is.na(Gr[1]))
data$Gr$sigma.r <- Gr$sigma.r * alpha
}
cat("Background estimation complete! \n")
if(stdev==TRUE && analytical==FALSE){
cat("Calculating uncertainty in background... \n")
stdev <- get.hess(data, knots.x, knots.y, Gr=data$Gr, r=seq(0, 2, 0.01), p.bkg=p.bkg)
}
curves <- list(y=data$y, bkg=bkg, SB=data$SB)
knots <- list(x=knots.x, y=knots.y)
fit.details <- list(lambda=data$lambda, sigma=data$sigma, knots.n=knots.n,
control=control, Gr=data$Gr, n.atoms=data$fitADP$n.atoms,
scatter.length=data$fitADP$scatter.length, id=data$id,
bounds.lower=bounds.lower['knots.y1'],
bounds.upper=bounds.upper['knots.y1'])
fit.results <- list(x=data$x, curves=curves, uncrt=stdev, knots=knots,
pars=pars, scale=alpha, ADP=ADP, fit.details=fit.details)
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.results, file=save.to)
}
cat("...done! \n")
return(fit.results)
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, knots.n.left, knots.n.right, x.boundary, control, save.to)
#
# returns: Performs evolutionary global background optimization via DEoptim
# (wrapper to do.fit) for several banks.
# arguments
# data: object of type data. Contatins experimental data and fit parameters
# bounds.lower, lower and upeer boundaries for background
# bounds.upper:
# knots.n.left, specify knot positions. knots.n.left and knots.n.right knots are created on the
# knots.n.right, left and on the right of x.boundary point, respectively
# x.boundary:
# control: he return value of set.control
# save.to: the name of the file where the results will be saved
do.fit.banks <- function(data, bounds.lower, bounds.upper, knots.n.left=NA,
knots.n.right=NA, x.boundary=NA, analytical=FALSE, control, save.to=""){
N <- length(data)
fit.res <- list()
knots.x <- NA
for(i in 1:N){
cat("\n\n ===================================\n\n")
cat("Fitting bank # ",i,"\n\n")
x <- data[[i]]$x
if(analytical==FALSE){
dx <- ((max(x)-x.boundary)/(knots.n.right-1) + x.boundary/(knots.n.left-1) )/2
knots.x <- seq(0, x.boundary, length=knots.n.left)
knots.x <- c(knots.x, seq(x.boundary+dx, max(x), length=knots.n.right))
}
fit.res[[i]] <- do.fit(data[[i]], bounds.lower, bounds.upper, knots.x=knots.x,
analytical=analytical, stdev=FALSE, control=control, save.to="")
fit.res[[i]]$fit.details$id <- data[[i]]$id
}
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.res, file=save.to)
}
fit.res
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, knots.n.left, knots.n.right, x.boundary, control, save.to)
#
do.iter <- function(fit.results, local=TRUE, eps=1e-4, n.iter=10000, save.to=""){
dat <- list(x=fit.results$x, y=fit.results$curves$y, SB=fit.results$curves$SB,
lambda=fit.results$fit.details$lambda, sigma=fit.results$fit.details$sigma)
knots.x <- fit.results$knots$x
knots.y <- fit.results$knots$y
cat("Adjusting baseline...", save.to, "\n")
if(local)
cc <- grad.descent(data=dat, knots.x, knots.y, Gr=NA, p.bkg=0.5, eps=eps, N=n.iter)
else{
control <- fit.results$fit.details$control
bounds.lower <- fit.results$fit.details$bounds.lower
bounds.upper <- fit.results$fit.details$bounds.upper
cat("\n\n Starting DifEv to find bkg with no low-r constraints... \n\n")
ff <- do.fit(dat, bounds.lower, bounds.upper, knots.x=knots.x,
stdev=FALSE, control=control)
cc <- ff$knots$y
}
if(any(is.na(cc))) stop("perhaps adjust parameters")
l <- fit.results$curves$bkg
l.no.r <- get.bkg(dat$x, knots.x, cc)
r <- seq(0, 1.0, 0.005)
gr <- sineFT(f.Q=dat$y-l-1, Q=dat$x, r=r)
gr.no.r <- sineFT(f.Q=dat$y-l.no.r-1, Q=dat$x, r=r)
d <- gr.no.r- gr # < 0!!! Difference between fits with and without Gr info
l.corr <- 0
dr <- r[2]-r[1]
for(i in 1:length(dat$x)){
l.corr[i] <- sum(d*sin(dat$x[i]*r)*dr/dat$x[i])
}
dat$SB <- dat$SB - l.corr
cat("\n\n Estimating background for the corrected data... \n\n")
if(local)
cc2 <- grad.descent(data=dat, knots.x, cc, Gr=fit.results$fit.details$Gr, p.bkg=0.5, eps=eps, N=n.iter)
else{
Gr <- fit.results$fit.details$Gr
dat <- set.Gr(dat, r1=Gr$r1, r2=Gr$r2, rho.0=Gr$rho.0,
type1=Gr$type1, type2=Gr$type2)
ff2 <- do.fit(dat, bounds.lower, bounds.upper, knots.x=knots.x,
stdev=FALSE, control=control)
cc2 <- ff2$knots$y
}
if(any(is.na(cc2))) stop("perhaps adjust parameters")
l2 <- get.bkg(dat$x, knots.x, cc2)
fit.results$curves$bkg <- l2
fit.results$knots$y <- cc2
fit.results$curves$corr <- l.corr
if(length(fit.results$uncrt)>1){
cat("Calculating uncertainty in background... \n")
stdev <- get.hess(dat, knots.x, knots.y, Gr=fit.results$fit.details$Gr, r=seq(0, 2, 0.01), p.bkg=.5)
}
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.results, file=save.to)
}
return(fit.results)
}
# MATMUL
#matmul <- function(A, B){
# if(!is.loaded("matmul"))
# dyn.load("./matmul.dll")
# matrix(.C("matmul_matmul_R",
# heightA=as.integer(nrow(A)),
# widthA=as.integer(ncol(A)),
# widthB=as.integer(ncol(B)),
# A=as.double(c(A)),
# rstrideA=as.integer(ncol(A)),
# B=as.double(c(B)),
# rstrideB=as.integer(ncol(B)),
# C=as.double(rep(0, nrow(A)*ncol(B))),
# rstrideC=as.integer(ncol(B)))$C,
# nrow=nrow(A), ncol=ncol(B))
#}
Try the BBEST package in your browser
Any scripts or data that you put into this service are public.
BBEST documentation built on Jan. 8, 2021, 2:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions do.iter do.fit.banks do.fit set.DEoptim.control set.control
Documented in do.fit do.fit.banks do.iter set.control set.DEoptim.control
#####################################################################################
#
#
# FUNCTIONS TO PERFORM BACKGROUND FIT VIA DIFFERENTIAL EVOLUTION ALGORITHM
#
#
#####################################################################################
# set.control(CR, F, NP, itermax, parallelType)
#
# returns: wrapper to set.DEoptim.control
# arguments
# CR: crossover probability from interval [0,1]
# F: differential weighting factor from interval [0,2]
# NP: number of population members
# itermax: the maximum iteration (population generation) allowed
# parallelType: defines the type of parallelization to employ.
set.control <- function(CR=.85, F=.7, NP=300, itermax=2000, parallelType=1){
control <- list()
control$CR <- CR
control$F <- F
control$NP <- NP
control$itermax <- itermax
control$parallelType <- parallelType
return(control)
}
#####################################################################################
# set.DEoptim.control(control, knots.n)
#
# returns: DEoptim control parameters
# arguments
# control: return value of set.control
# knots.n: knots number
set.DEoptim.control <- function(control=list(), knots.n){
if(is.null(control$CR))
CR <- 0.85
else
CR <- control$CR
if(is.null(control$F))
F <- 0.7
else
F <- control$F
if(is.null(control$NP))
NP <- round(2.5*knots.n*10)
else
NP <- control$NP
if(is.null(control$itermax))
itermax <- 2000
else
itermax <- control$itermax
if(is.null(control$parallelType))
parallelType <- 1
else
parallelType <- control$parallelType
DE.control <- DEoptim.control(CR=CR, F=F, NP=NP, itermax=itermax, parallelType=parallelType,
# packages = list("PerformanceAnalytics"),
parVar=list("basisMatrix", "basisSpline",
"logPosterior", "logPriorF", "DMatrix", "logLikelihoodBkg", "logLikelihoodSignal",
"get.bkg", "logProbabilityBkgR", "logLikelihoodGrGauss", "logLikelihoodGrCorr",
"logPriorBkgRSmooth", "logPriorBkgRGP", "Dx", "invert.order", "set.SB",
"logPosteriorAnalyt", "bkg.analyt", "logLikelihoodBkgAnalyt", "logLikelihoodSignalAnalyt"))
DE.control
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, scale, knots.x, knots.n, stdev, control, p.bkg, save.to)
#
# returns: Performs evolutionary global background optimization via DEoptim
# arguments
# data: object of type data. Contatins experimental data and fit parameters
# bounds.lower, lower and upeer boundaries for background
# bounds.upper:
# scale: scale factor
# knots.x: the spline knot x-values.
# knots.n: number of knots.
# stdev: Whether to calculate uncertainty in background estimation
# control: he return value of set.control
# p.bkg: the probability that a single pixel contains "only" background
# save.to: the name of the file where the results will be saved
do.fit <- function(data, bounds.lower, bounds.upper, scale=c(1,1), knots.x=NA,
knots.n=NA, analytical=FALSE, stdev=TRUE, control=list(), p.bkg=.5, save.to=""){
# 1. Prepare data to fit...
cat("Preparing data... \n")
if( (scale[1]==1) && (scale[2]==1) ) scale <- NA
cov.r <- ADP <- knots.y <- pars <- NA
alpha <- 1
if(is.null(data$Gr))
Gr <- NA
else{
Gr <- data$Gr
data$Gr <- NA
}
if(is.null(data$SB))
data$SB <- rep(0, length(data$x))
else if (is.na(data$SB[1]))
data$SB <- rep(0, length(data$x))
# 2. prepare fit params...
if(is.na(knots.x[1]) && !is.na(knots.n))
knots.x <- seq(min(data$x), max(data$x), length=knots.n)
if(is.na(knots.n))
knots.n <- length(knots.x)
if(analytical==TRUE){
iter.0 <- as.relistable(list(pars=rep(0, 6)))
bounds.lower <- c(pars=c(-abs(bounds.lower[1]*10), -5, -5, -20, 0, 0))
bounds.upper <- c(pars=c(abs(bounds.upper[1]*10), 10, 5, 2000,
max(data$x), max(data$x)*2))
knots.n <- 6
}
else{
iter.0 <- as.relistable(list(knots.y=rep(0, knots.n)))
bounds.lower <- c(knots.y=rep(bounds.lower[1], knots.n))
bounds.upper <- c(knots.y=rep(bounds.upper[1], knots.n))
}
if(!is.na(scale[1])){
iter.0$alpha <- 0
bounds.lower <- c(bounds.lower, alpha=scale[1])
bounds.upper <- c(bounds.upper, alpha=scale[2])
}
if(!is.null(data$fitADP)){
if(data$fitADP$oneADP==TRUE)
nn <- 1
else
nn <- length(data$fitADP$n.atoms)
iter.0$ADP <- rep(0, nn)
bounds.lower <- c(bounds.lower, ADP=rep(data$fitADP$ADP.lim[1], nn))
bounds.upper <- c(bounds.upper, ADP=rep(data$fitADP$ADP.lim[2], nn))
}
DE.control <- set.DEoptim.control(control, knots.n)
# 3. Starting fit!
#
# initial guess
# cc <- matrix(nrow=DE.control$NP, ncol=length(knots.x))
# Phi <- basisMatrix(data$x, knots.x)
# for(ii in 1:DE.control$NP){
# spar=(0.9-0.5)/DE.control$NP*ii+0.5
# lowpass.spline <- smooth.spline(data$x,data$y-data$SB, spar = spar) ## Control spar for amount of smoothing
# cc[ii,] <- c(solve(t(Phi)%*%Phi) %*% t(Phi) %*% (predict(lowpass.spline, data$x)$y) )
# }
# DE.control$initialpop=cc
cat("Starting DifEv algorithm... \n")
doPbkgIter <- FALSE
if(p.bkg==-1){
doPbkgIter <- TRUE
p.bkg <- 0.02
}
DEoptim.fit <- DEoptim(get.posterior, lower = bounds.lower, upper = bounds.upper,
control= DE.control, skel=iter.0, data=data,
knots.x=knots.x, Gr=Gr, p.bkg=p.bkg)
if(doPbkgIter){
cat("\n\n P.bkg iteration... \n\n")
if(analytical==TRUE){
pars <- DEoptim.fit$optim$bestmem[1:6]
bkg <- bkg.analyt(pars=pars, x=data$x)
}
else{
knots.y <- DEoptim.fit$optim$bestmem[1:knots.n]
bkg <- get.bkg(x=data$x, knots.x=knots.x, knots.y=knots.y)
}
dev.norm <- (data$y-bkg-data$SB)/data$sigma
p.bkg <- 0.5*50^((1-dev.norm^2)/8)
DEoptim.fit <- DEoptim(get.posterior, lower = bounds.lower, upper = bounds.upper,
control= DE.control, skel=iter.0, data=data,
knots.x=knots.x, Gr=Gr, p.bkg=p.bkg)
}
if(analytical==TRUE){
pars <- DEoptim.fit$optim$bestmem[1:6]
bkg <- bkg.analyt(pars=pars, x=data$x)
}
else{
knots.y <- DEoptim.fit$optim$bestmem[1:knots.n]
bkg <- get.bkg(x=data$x, knots.x=knots.x, knots.y=knots.y)
}
if(!is.null(data$fitADP)){
ADP <- tail(DEoptim.fit$optim$bestmem, nn)
data <- set.SB(data, n.atoms=data$fitADP$n.atoms,
scatter.length=data$fitADP$scatter.length, ADP=ADP)
}
data$Gr <- Gr
if(!is.na(scale[1])){ # recalculating signal and params...
alpha <- DEoptim.fit$optim$bestmem[knots.n+1]
data$y <- (data$y-bkg-data$SB)*alpha + data$SB + bkg
data$lambda <- data$lambda*alpha
data$sigma <- data$sigma*alpha
if(!is.na(Gr[1]))
data$Gr$sigma.r <- Gr$sigma.r * alpha
}
cat("Background estimation complete! \n")
if(stdev==TRUE && analytical==FALSE){
cat("Calculating uncertainty in background... \n")
stdev <- get.hess(data, knots.x, knots.y, Gr=data$Gr, r=seq(0, 2, 0.01), p.bkg=p.bkg)
}
curves <- list(y=data$y, bkg=bkg, SB=data$SB)
knots <- list(x=knots.x, y=knots.y)
fit.details <- list(lambda=data$lambda, sigma=data$sigma, knots.n=knots.n,
control=control, Gr=data$Gr, n.atoms=data$fitADP$n.atoms,
scatter.length=data$fitADP$scatter.length, id=data$id,
bounds.lower=bounds.lower['knots.y1'],
bounds.upper=bounds.upper['knots.y1'])
fit.results <- list(x=data$x, curves=curves, uncrt=stdev, knots=knots,
pars=pars, scale=alpha, ADP=ADP, fit.details=fit.details)
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.results, file=save.to)
}
cat("...done! \n")
return(fit.results)
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, knots.n.left, knots.n.right, x.boundary, control, save.to)
#
# returns: Performs evolutionary global background optimization via DEoptim
# (wrapper to do.fit) for several banks.
# arguments
# data: object of type data. Contatins experimental data and fit parameters
# bounds.lower, lower and upeer boundaries for background
# bounds.upper:
# knots.n.left, specify knot positions. knots.n.left and knots.n.right knots are created on the
# knots.n.right, left and on the right of x.boundary point, respectively
# x.boundary:
# control: he return value of set.control
# save.to: the name of the file where the results will be saved
do.fit.banks <- function(data, bounds.lower, bounds.upper, knots.n.left=NA,
knots.n.right=NA, x.boundary=NA, analytical=FALSE, control, save.to=""){
N <- length(data)
fit.res <- list()
knots.x <- NA
for(i in 1:N){
cat("\n\n ===================================\n\n")
cat("Fitting bank # ",i,"\n\n")
x <- data[[i]]$x
if(analytical==FALSE){
dx <- ((max(x)-x.boundary)/(knots.n.right-1) + x.boundary/(knots.n.left-1) )/2
knots.x <- seq(0, x.boundary, length=knots.n.left)
knots.x <- c(knots.x, seq(x.boundary+dx, max(x), length=knots.n.right))
}
fit.res[[i]] <- do.fit(data[[i]], bounds.lower, bounds.upper, knots.x=knots.x,
analytical=analytical, stdev=FALSE, control=control, save.to="")
fit.res[[i]]$fit.details$id <- data[[i]]$id
}
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.res, file=save.to)
}
fit.res
}
#####################################################################################
# do.fit(data, bounds.lower, bounds.upper, knots.n.left, knots.n.right, x.boundary, control, save.to)
#
do.iter <- function(fit.results, local=TRUE, eps=1e-4, n.iter=10000, save.to=""){
dat <- list(x=fit.results$x, y=fit.results$curves$y, SB=fit.results$curves$SB,
lambda=fit.results$fit.details$lambda, sigma=fit.results$fit.details$sigma)
knots.x <- fit.results$knots$x
knots.y <- fit.results$knots$y
cat("Adjusting baseline...", save.to, "\n")
if(local)
cc <- grad.descent(data=dat, knots.x, knots.y, Gr=NA, p.bkg=0.5, eps=eps, N=n.iter)
else{
control <- fit.results$fit.details$control
bounds.lower <- fit.results$fit.details$bounds.lower
bounds.upper <- fit.results$fit.details$bounds.upper
cat("\n\n Starting DifEv to find bkg with no low-r constraints... \n\n")
ff <- do.fit(dat, bounds.lower, bounds.upper, knots.x=knots.x,
stdev=FALSE, control=control)
cc <- ff$knots$y
}
if(any(is.na(cc))) stop("perhaps adjust parameters")
l <- fit.results$curves$bkg
l.no.r <- get.bkg(dat$x, knots.x, cc)
r <- seq(0, 1.0, 0.005)
gr <- sineFT(f.Q=dat$y-l-1, Q=dat$x, r=r)
gr.no.r <- sineFT(f.Q=dat$y-l.no.r-1, Q=dat$x, r=r)
d <- gr.no.r- gr # < 0!!! Difference between fits with and without Gr info
l.corr <- 0
dr <- r[2]-r[1]
for(i in 1:length(dat$x)){
l.corr[i] <- sum(d*sin(dat$x[i]*r)*dr/dat$x[i])
}
dat$SB <- dat$SB - l.corr
cat("\n\n Estimating background for the corrected data... \n\n")
if(local)
cc2 <- grad.descent(data=dat, knots.x, cc, Gr=fit.results$fit.details$Gr, p.bkg=0.5, eps=eps, N=n.iter)
else{
Gr <- fit.results$fit.details$Gr
dat <- set.Gr(dat, r1=Gr$r1, r2=Gr$r2, rho.0=Gr$rho.0,
type1=Gr$type1, type2=Gr$type2)
ff2 <- do.fit(dat, bounds.lower, bounds.upper, knots.x=knots.x,
stdev=FALSE, control=control)
cc2 <- ff2$knots$y
}
if(any(is.na(cc2))) stop("perhaps adjust parameters")
l2 <- get.bkg(dat$x, knots.x, cc2)
fit.results$curves$bkg <- l2
fit.results$knots$y <- cc2
fit.results$curves$corr <- l.corr
if(length(fit.results$uncrt)>1){
cat("Calculating uncertainty in background... \n")
stdev <- get.hess(dat, knots.x, knots.y, Gr=fit.results$fit.details$Gr, r=seq(0, 2, 0.01), p.bkg=.5)
}
if(save.to!=""){
cat("Saving results to file ", save.to, "\n")
save(fit.results, file=save.to)
}
return(fit.results)
}
# MATMUL
#matmul <- function(A, B){
# if(!is.loaded("matmul"))
# dyn.load("./matmul.dll")
# matrix(.C("matmul_matmul_R",
# heightA=as.integer(nrow(A)),
# widthA=as.integer(ncol(A)),
# widthB=as.integer(ncol(B)),
# A=as.double(c(A)),
# rstrideA=as.integer(ncol(A)),
# B=as.double(c(B)),
# rstrideB=as.integer(ncol(B)),
# C=as.double(rep(0, nrow(A)*ncol(B))),
# rstrideC=as.integer(ncol(B)))$C,
# nrow=nrow(A), ncol=ncol(B))
#}
Try the BBEST package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.