Nothing
R/rss.R
In nanop: Tools for Nanoparticle Simulation and Calculation of PDF and Total Scattering Structure Function
###################################################################################
## rss
## Objective function
## this function can be minimized to find the parameter values that
## minimize the RSS for a given model and given data
## In our experience this is best done with the package DEoptim.
## for pareto task use package mco
# dataSAS: SAS part of S(Q)
# paramSASQ: if SAS(Q) should be calculated parametrically
# gammaR: SAS part of G(r)
# = "param" => parametric envelope is calculated
# avRes: averages results over avRes simulations
# wG: weight for G(r)
# wSAS: weight for SAS(Q)
# pareto: if TRRUE then fuction value is a vector of residuals
# kind: "exact", "fast", "fastHist", "fast_av" (for polydisperse particles)
# SASscale: "normal", "log", weighted"
# scaleSq: calculated PDF or S(Q) function should be multiplied by scaleSq
# scaleSASq: calculated SAS(Q) should be multiplied by scaleSASq
# if gammaR is given it is divided by scaleSASq
rss <- function(par, dataG=NA, dataS=NA, dataSAS=NA,
part=NA, type="neutron", simPar=NA,
PDF.fixed=list(), TotalScatt.fixed=list(),
parscale=NA, skel=NA, con=TRUE, oneDW=FALSE,
punish=FALSE,
gammaR=NA, rvector=NA, fixed=NA,
wG=0.03, wSAS=0.03, avRes=1, pareto=FALSE){
if(!is.na(parscale[1]))
par2 <- relist( par / parscale, skel)
else
par2 <- relist(par, skel)
if(is.list(fixed))
par2[names(fixed)] <- fixed
if(con && !(is.null(par2$r))) {
if(length(c(par2$rcore)) == 1) ## test for non-param
if(par2$rcore > par2$r)
par2$rcore <- par2$r
}
else if(punish) { ## only makes sense now for core shell
if(length(c(par2$rcore,par2$r)) == 2)
if(par2$rcore>par$r)
return(10e15)
}
r <- if(is.null(par2$r)) NA else par2$r
rsigma <- if(is.null(par2$rsigma)) NA else par2$rsigma
bkg <- if(is.null(par2$bkg)) 0 else par2$bkg
shell <- if(is.null(par2$shell)) NA else par2$shell
rcore <- if(is.null(par2$rcore)) NA else par2$rcore
scaleSq <- if(is.null(par2$scaleSq)) 1 else par2$scaleSq
scaleSASq <- if(is.null(par2$scaleSASq)) 1 else par2$scaleSASq
pStack <- if(is.null(par2$pStack)) 0 else par2$pStack
pDimer <- if(is.null(par2$pDimer)) 0 else par2$pDimer
delta <- if(is.null(par2$delta)) 0 else par2$delta
box <- if(is.null(par2$box)) NA else par2$box
ellipse <- if(is.null(par2$ellipse)) NA else par2$ellipse
if(delta==0 && !is.na(PDF.fixed$delta))
delta <- PDF.fixed$delta
if(delta==0 && !is.na(TotalScatt.fixed$delta))
delta <- TotalScatt.fixed$delta
if(!is.null(par2$sigma))
sig <- par2$sigma
else{
sig <- 0
for(i in 1:length(simPar$atoms))
sig <- append(sig, simPar$atoms[[i]]$sigma)
if(length(simPar$atomsShell)>0){
for(i in 1:length(simPar$atomsShell))
sig <- append(sig, simPar$atomsShell[[i]]$sigma)
}
sig <- sig[-1]
}
if(is.na(r))
r<- simPar$r
if(is.na(rsigma))
rsigma<- simPar$rsigma
rIn <- r
if(length(par2$latticep)>0)
latticep <- par2$latticep
else
latticep <- simPar$latticep
if(length(latticep)==3){ #works for wurtzite structure
simPar$atoms[[2]]$base[1,3] <- latticep[3]
simPar$atoms[[2]]$base[2,3] <- 0.5-latticep[3]
}
if(length(par2$latticepShell)>0)
latticepShell <- par2$latticepShell
else
latticepShell <- simPar$latticepShell
if(length(latticepShell)==3){ #works for wurtzite structure
simPar$atomsShell[[2]]$base[1,3] <- latticepShell[3]
simPar$atomsShell[[2]]$base[2,3] <- 0.5-latticepShell[3]
}
if(!is.na(shell) && !is.na(rsigma)){
if(!is.na(rcore)){
if(simPar$distr=="lognormal")
core <- rlnorm(avRes, meanlog = log(rcore), sdlog = log(rsigma))
else
core <- exp(rnorm(avRes, log(rcore), log(rsigma)))
core <- sort(core)
r <- core + shell
rcore <- core
}
if(!is.na(r)){
if(simPar$distr=="lognormal")
r <- rlnorm(avRes, meanlog = log(r), sdlog = log(rsigma))
else
r <- exp(rnorm(avRes, log(r), log(rsigma)))
r <- sort(r)
rcore <- r-shell
}
}
if(is.na(shell) && !is.na(rsigma)){
if(simPar$distr=="lognormal")
r <- rlnorm(avRes, meanlog = log(r), sdlog = log(rsigma))
else
r <- exp(rnorm(avRes, log(r), log(rsigma)))
r <- sort(r)
}
kind <- TotalScatt.fixed$kind
dr <- TotalScatt.fixed$dr
del <- TotalScatt.fixed$del
eps <- TotalScatt.fixed$eps
SASscale <- TotalScatt.fixed$SASscale
paramSASQ <- TotalScatt.fixed$paramSASQ
gQ_SAS_av <- gQ_av <- mod_av <- gQ_SAS <- gQ <- mod <- 0
simulate_particle <- TRUE
####################
#
# !MAIN CYCLE OVER AVRES!
#
####################
for(j in 1:avRes){
p <- 0
if(kind !="fast_av" && !is.na(rsigma)){
r_tmp <- r[j]
rcore_tmp <- rcore[j]
}
else{
r_tmp <- r
rcore_tmp <- rcore
}
if(!is.na(part[1]) && j==1) # if on the first step we already have particle
simulate_particle <- FALSE # we don't have to simulate new ones...
if(simulate_particle)
part <- simPart(atoms=simPar$atoms, sym=simPar$sym, latticep=latticep, r=r_tmp,
atomsShell=simPar$atomsShell, symShell=simPar$symShell,
latticepShell=latticepShell, rcore=rcore_tmp,
rcenter=simPar$rcenter, box=box, ellipse=ellipse,
move=simPar$move, rotShell=simPar$rotShell, rcenterShell=simPar$rcenterShell,
pDimer=pDimer, pStack=pStack)
if(oneDW && j==1)
sig <- rep(sig, attributes(part)$nAtomTypes)
#========================================================
if(!is.na(dataS[1])){
cat("Calculating Scattering function... \n")
gQ <- calcTotalScatt(part, type=type, kind=kind,
dQ=TotalScatt.fixed$dQ,
minQ=TotalScatt.fixed$minQ,
maxQ=TotalScatt.fixed$maxQ,
scatterLength=TotalScatt.fixed$scatterLength,
scatterFactor=TotalScatt.fixed$scatterFactor,
sigma=sig,
delta=delta, n=TotalScatt.fixed$n,
dr=dr, del=del, eps=eps)$gQ
}
#========================================================
if(!is.na(dataG[1])){
cat("Calculating PDF... \n")
mod <- calcPDF(part,type=type,
dr=PDF.fixed$dr,
minR=PDF.fixed$minR,
maxR=PDF.fixed$maxR,
scatterLength=PDF.fixed$scatterLength,
scatterFactor=PDF.fixed$scatterFactor,
)
mod <- broadPDF(mod, sigma=sig, delta=delta, n=PDF.fixed$n, nAtomTypes = attributes(part)$nAtomTypes)$gr
mod[which(is.na(mod))] <- 0
if(any(is.na(mod)))
mod <- rep(10e15, length(rvector))
mod[which(is.infinite(mod))] <- 0
}
TotalScattSAS=list(dQ=TotalScatt.fixed$dQ_SAS, minQ=TotalScatt.fixed$minQ_SAS, maxQ=TotalScatt.fixed$maxQ_SAS, scatterLength=TotalScatt.fixed$scatterLength)
#========================================================
if(!is.na(dataSAS[1]) && (!paramSASQ)){
cat("Calculating SAS ... \n")
gQ_SAS <- calcTotalScatt(part,type=type, kind=kind,
dQ=TotalScattSAS$dQ,
minQ=TotalScattSAS$minQ,
maxQ=TotalScattSAS$maxQ,
scatterLength=TotalScatt.fixed$scatterLength,
scatterFactor=TotalScatt.fixed$scatterFactor,
sigma=sig,
delta=delta, n=TotalScatt.fixed$n,
dr=dr, del=del, eps=eps)$gQ
}
#=============================================================
gQ_SAS_av <- gQ_SAS_av + gQ_SAS
gQ_av <- gQ_av + gQ
mod_av <- mod_av + mod
}
gQ_SAS <- gQ_SAS_av/avRes
gQ <- gQ_av/avRes
mod <- mod_av/avRes
if(!is.na(dataG[1])){
if(!any(is.na(rvector)))
mod <- 4*pi*rvector*mod
if(!any(is.na(gammaR))){
if(gammaR[1]=="param"){
if(!is.na(shell))
sasvector <- GrSASCS(rvector, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
N1=PDF.fixed$sym, N2=PDF.fixed$sym, sym=simPar$sym, symShell=simPar$symShell)
else
sasvector <- GrSAS(rvector, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
N1=PDF.fixed$sym, N2=PDF.fixed$sym, sym=simPar$sym, symShell=simPar$symShell)
mod <- mod - sasvector
}
else
mod <- mod - gammaR/scaleSASq
}
if(PDF.fixed$termRip && !any(is.na(mod)))
mod <- termRip(mod, Qmax=PDF.fixed$Qmax, dr=PDF.fixed$dr, maxR=PDF.fixed$maxR, maxRTermRip=PDF.fixed$maxRTermRip)
if(!is.na(PDF.fixed$Qdamp))
mod <- mod*exp(-(PDF.fixed$Qdamp*rvector)^2/2)
}
if(!is.na(dataS[1]) && TotalScatt.fixed$convolution && !any(is.na(gQ)))
gQ <- gaussConvol(SQ=gQ, Q=seq(TotalScatt.fixed$minQ, TotalScatt.fixed$maxQ, TotalScatt.fixed$dQ), Qdamp=TotalScatt.fixed$Qdamp)
if(!is.na(dataSAS[1]) && paramSASQ){
cat("Calculating parametric SAS(Q)... \n")
Q=seq(TotalScatt.fixed$minQ_SAS, TotalScatt.fixed$maxQ_SAS, TotalScatt.fixed$dQ_SAS)
if(!is.na(rsigma))
gQ_SAS <- IqSASP(Q=Q, shell=shell, Rpart=rIn, latticep=latticep, latticepShell=latticepShell,
scatterLength=c(TotalScatt.fixed$f1, TotalScatt.fixed$f2), pDimer=pDimer,
N1=TotalScatt.fixed$N1, N2=TotalScatt.fixed$N2,
sym=simPar$sym, symShell=simPar$symShell, rsigma=rsigma)
else
gQ_SAS <- IqSAS(Q=Q, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
scatterLength=c(TotalScatt.fixed$f1, TotalScatt.fixed$f2),
N1=TotalScatt.fixed$N1, N2=TotalScatt.fixed$N2,
pDimer=pDimer, sym=simPar$sym, symShell=simPar$symShell)
}
##
mod <- scaleSq*mod
if(!is.na(dataG[1]))
rssG <- sqrt(sum((mod-dataG)^2)/( sum(dataG^2))) * wG
else
rssG <- 0
##
gQ <- scaleSq*gQ
if(!is.na(dataS[1]))
rssS <- sqrt(sum((gQ-dataS)^2)/(sum(dataS^2)))
else
rssS <- 0
##
if(!is.na(dataSAS[1])){
dataSAS <- dataSAS - bkg
gQ_SAS <- scaleSASq*gQ_SAS
if(SASscale=="normal")
rssSAS <- sqrt(sum((gQ_SAS-dataSAS)^2)/(sum(dataSAS^2))) * wSAS
if(SASscale=="log")
rssSAS <- sqrt(sum((log(abs(gQ_SAS))-log(abs(dataSAS)))^2)/sum( (log(abs(dataSAS)))^2 )) * wSAS
if(SASscale=="weighted")
rssSAS <- sqrt(sum((gQ_SAS/dataSAS-1)^2)/(length(dataSAS))) * wSAS
}
else
rssSAS <- 0
if(is.na(rssS)) rssS <- 1e15
if(is.na(rssSAS)) rssSAS <- 1e15
if(is.na(rssG)) rssG <- 1e15
include <- logical(3)
include[1] <- !is.na(dataSAS[1])
include[2] <- !is.na(dataG[1])
include[3] <- !is.na(dataS[1])
res <- numeric(3)
if(include[1])
res[1] <- rssSAS
else
res[1] <- 0
if(include[2])
res[2] <- rssG
else
res[2] <- 0
if(include[3])
res[3] <- rssS
else
res[3] <- 0
if(pareto)
res[which(include==TRUE)]
else
sqrt(res[1]^2 + res[2]^2 + res[3]^2)
}
Try the nanop package in your browser
Any scripts or data that you put into this service are public.
nanop documentation built on May 2, 2019, 3:39 p.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.
###################################################################################
## rss
## Objective function
## this function can be minimized to find the parameter values that
## minimize the RSS for a given model and given data
## In our experience this is best done with the package DEoptim.
## for pareto task use package mco
# dataSAS: SAS part of S(Q)
# paramSASQ: if SAS(Q) should be calculated parametrically
# gammaR: SAS part of G(r)
# = "param" => parametric envelope is calculated
# avRes: averages results over avRes simulations
# wG: weight for G(r)
# wSAS: weight for SAS(Q)
# pareto: if TRRUE then fuction value is a vector of residuals
# kind: "exact", "fast", "fastHist", "fast_av" (for polydisperse particles)
# SASscale: "normal", "log", weighted"
# scaleSq: calculated PDF or S(Q) function should be multiplied by scaleSq
# scaleSASq: calculated SAS(Q) should be multiplied by scaleSASq
# if gammaR is given it is divided by scaleSASq
rss <- function(par, dataG=NA, dataS=NA, dataSAS=NA,
part=NA, type="neutron", simPar=NA,
PDF.fixed=list(), TotalScatt.fixed=list(),
parscale=NA, skel=NA, con=TRUE, oneDW=FALSE,
punish=FALSE,
gammaR=NA, rvector=NA, fixed=NA,
wG=0.03, wSAS=0.03, avRes=1, pareto=FALSE){
if(!is.na(parscale[1]))
par2 <- relist( par / parscale, skel)
else
par2 <- relist(par, skel)
if(is.list(fixed))
par2[names(fixed)] <- fixed
if(con && !(is.null(par2$r))) {
if(length(c(par2$rcore)) == 1) ## test for non-param
if(par2$rcore > par2$r)
par2$rcore <- par2$r
}
else if(punish) { ## only makes sense now for core shell
if(length(c(par2$rcore,par2$r)) == 2)
if(par2$rcore>par$r)
return(10e15)
}
r <- if(is.null(par2$r)) NA else par2$r
rsigma <- if(is.null(par2$rsigma)) NA else par2$rsigma
bkg <- if(is.null(par2$bkg)) 0 else par2$bkg
shell <- if(is.null(par2$shell)) NA else par2$shell
rcore <- if(is.null(par2$rcore)) NA else par2$rcore
scaleSq <- if(is.null(par2$scaleSq)) 1 else par2$scaleSq
scaleSASq <- if(is.null(par2$scaleSASq)) 1 else par2$scaleSASq
pStack <- if(is.null(par2$pStack)) 0 else par2$pStack
pDimer <- if(is.null(par2$pDimer)) 0 else par2$pDimer
delta <- if(is.null(par2$delta)) 0 else par2$delta
box <- if(is.null(par2$box)) NA else par2$box
ellipse <- if(is.null(par2$ellipse)) NA else par2$ellipse
if(delta==0 && !is.na(PDF.fixed$delta))
delta <- PDF.fixed$delta
if(delta==0 && !is.na(TotalScatt.fixed$delta))
delta <- TotalScatt.fixed$delta
if(!is.null(par2$sigma))
sig <- par2$sigma
else{
sig <- 0
for(i in 1:length(simPar$atoms))
sig <- append(sig, simPar$atoms[[i]]$sigma)
if(length(simPar$atomsShell)>0){
for(i in 1:length(simPar$atomsShell))
sig <- append(sig, simPar$atomsShell[[i]]$sigma)
}
sig <- sig[-1]
}
if(is.na(r))
r<- simPar$r
if(is.na(rsigma))
rsigma<- simPar$rsigma
rIn <- r
if(length(par2$latticep)>0)
latticep <- par2$latticep
else
latticep <- simPar$latticep
if(length(latticep)==3){ #works for wurtzite structure
simPar$atoms[[2]]$base[1,3] <- latticep[3]
simPar$atoms[[2]]$base[2,3] <- 0.5-latticep[3]
}
if(length(par2$latticepShell)>0)
latticepShell <- par2$latticepShell
else
latticepShell <- simPar$latticepShell
if(length(latticepShell)==3){ #works for wurtzite structure
simPar$atomsShell[[2]]$base[1,3] <- latticepShell[3]
simPar$atomsShell[[2]]$base[2,3] <- 0.5-latticepShell[3]
}
if(!is.na(shell) && !is.na(rsigma)){
if(!is.na(rcore)){
if(simPar$distr=="lognormal")
core <- rlnorm(avRes, meanlog = log(rcore), sdlog = log(rsigma))
else
core <- exp(rnorm(avRes, log(rcore), log(rsigma)))
core <- sort(core)
r <- core + shell
rcore <- core
}
if(!is.na(r)){
if(simPar$distr=="lognormal")
r <- rlnorm(avRes, meanlog = log(r), sdlog = log(rsigma))
else
r <- exp(rnorm(avRes, log(r), log(rsigma)))
r <- sort(r)
rcore <- r-shell
}
}
if(is.na(shell) && !is.na(rsigma)){
if(simPar$distr=="lognormal")
r <- rlnorm(avRes, meanlog = log(r), sdlog = log(rsigma))
else
r <- exp(rnorm(avRes, log(r), log(rsigma)))
r <- sort(r)
}
kind <- TotalScatt.fixed$kind
dr <- TotalScatt.fixed$dr
del <- TotalScatt.fixed$del
eps <- TotalScatt.fixed$eps
SASscale <- TotalScatt.fixed$SASscale
paramSASQ <- TotalScatt.fixed$paramSASQ
gQ_SAS_av <- gQ_av <- mod_av <- gQ_SAS <- gQ <- mod <- 0
simulate_particle <- TRUE
####################
#
# !MAIN CYCLE OVER AVRES!
#
####################
for(j in 1:avRes){
p <- 0
if(kind !="fast_av" && !is.na(rsigma)){
r_tmp <- r[j]
rcore_tmp <- rcore[j]
}
else{
r_tmp <- r
rcore_tmp <- rcore
}
if(!is.na(part[1]) && j==1) # if on the first step we already have particle
simulate_particle <- FALSE # we don't have to simulate new ones...
if(simulate_particle)
part <- simPart(atoms=simPar$atoms, sym=simPar$sym, latticep=latticep, r=r_tmp,
atomsShell=simPar$atomsShell, symShell=simPar$symShell,
latticepShell=latticepShell, rcore=rcore_tmp,
rcenter=simPar$rcenter, box=box, ellipse=ellipse,
move=simPar$move, rotShell=simPar$rotShell, rcenterShell=simPar$rcenterShell,
pDimer=pDimer, pStack=pStack)
if(oneDW && j==1)
sig <- rep(sig, attributes(part)$nAtomTypes)
#========================================================
if(!is.na(dataS[1])){
cat("Calculating Scattering function... \n")
gQ <- calcTotalScatt(part, type=type, kind=kind,
dQ=TotalScatt.fixed$dQ,
minQ=TotalScatt.fixed$minQ,
maxQ=TotalScatt.fixed$maxQ,
scatterLength=TotalScatt.fixed$scatterLength,
scatterFactor=TotalScatt.fixed$scatterFactor,
sigma=sig,
delta=delta, n=TotalScatt.fixed$n,
dr=dr, del=del, eps=eps)$gQ
}
#========================================================
if(!is.na(dataG[1])){
cat("Calculating PDF... \n")
mod <- calcPDF(part,type=type,
dr=PDF.fixed$dr,
minR=PDF.fixed$minR,
maxR=PDF.fixed$maxR,
scatterLength=PDF.fixed$scatterLength,
scatterFactor=PDF.fixed$scatterFactor,
)
mod <- broadPDF(mod, sigma=sig, delta=delta, n=PDF.fixed$n, nAtomTypes = attributes(part)$nAtomTypes)$gr
mod[which(is.na(mod))] <- 0
if(any(is.na(mod)))
mod <- rep(10e15, length(rvector))
mod[which(is.infinite(mod))] <- 0
}
TotalScattSAS=list(dQ=TotalScatt.fixed$dQ_SAS, minQ=TotalScatt.fixed$minQ_SAS, maxQ=TotalScatt.fixed$maxQ_SAS, scatterLength=TotalScatt.fixed$scatterLength)
#========================================================
if(!is.na(dataSAS[1]) && (!paramSASQ)){
cat("Calculating SAS ... \n")
gQ_SAS <- calcTotalScatt(part,type=type, kind=kind,
dQ=TotalScattSAS$dQ,
minQ=TotalScattSAS$minQ,
maxQ=TotalScattSAS$maxQ,
scatterLength=TotalScatt.fixed$scatterLength,
scatterFactor=TotalScatt.fixed$scatterFactor,
sigma=sig,
delta=delta, n=TotalScatt.fixed$n,
dr=dr, del=del, eps=eps)$gQ
}
#=============================================================
gQ_SAS_av <- gQ_SAS_av + gQ_SAS
gQ_av <- gQ_av + gQ
mod_av <- mod_av + mod
}
gQ_SAS <- gQ_SAS_av/avRes
gQ <- gQ_av/avRes
mod <- mod_av/avRes
if(!is.na(dataG[1])){
if(!any(is.na(rvector)))
mod <- 4*pi*rvector*mod
if(!any(is.na(gammaR))){
if(gammaR[1]=="param"){
if(!is.na(shell))
sasvector <- GrSASCS(rvector, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
N1=PDF.fixed$sym, N2=PDF.fixed$sym, sym=simPar$sym, symShell=simPar$symShell)
else
sasvector <- GrSAS(rvector, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
N1=PDF.fixed$sym, N2=PDF.fixed$sym, sym=simPar$sym, symShell=simPar$symShell)
mod <- mod - sasvector
}
else
mod <- mod - gammaR/scaleSASq
}
if(PDF.fixed$termRip && !any(is.na(mod)))
mod <- termRip(mod, Qmax=PDF.fixed$Qmax, dr=PDF.fixed$dr, maxR=PDF.fixed$maxR, maxRTermRip=PDF.fixed$maxRTermRip)
if(!is.na(PDF.fixed$Qdamp))
mod <- mod*exp(-(PDF.fixed$Qdamp*rvector)^2/2)
}
if(!is.na(dataS[1]) && TotalScatt.fixed$convolution && !any(is.na(gQ)))
gQ <- gaussConvol(SQ=gQ, Q=seq(TotalScatt.fixed$minQ, TotalScatt.fixed$maxQ, TotalScatt.fixed$dQ), Qdamp=TotalScatt.fixed$Qdamp)
if(!is.na(dataSAS[1]) && paramSASQ){
cat("Calculating parametric SAS(Q)... \n")
Q=seq(TotalScatt.fixed$minQ_SAS, TotalScatt.fixed$maxQ_SAS, TotalScatt.fixed$dQ_SAS)
if(!is.na(rsigma))
gQ_SAS <- IqSASP(Q=Q, shell=shell, Rpart=rIn, latticep=latticep, latticepShell=latticepShell,
scatterLength=c(TotalScatt.fixed$f1, TotalScatt.fixed$f2), pDimer=pDimer,
N1=TotalScatt.fixed$N1, N2=TotalScatt.fixed$N2,
sym=simPar$sym, symShell=simPar$symShell, rsigma=rsigma)
else
gQ_SAS <- IqSAS(Q=Q, Rcore=rcore, Rpart=r, latticep=latticep, latticepShell=latticepShell,
scatterLength=c(TotalScatt.fixed$f1, TotalScatt.fixed$f2),
N1=TotalScatt.fixed$N1, N2=TotalScatt.fixed$N2,
pDimer=pDimer, sym=simPar$sym, symShell=simPar$symShell)
}
##
mod <- scaleSq*mod
if(!is.na(dataG[1]))
rssG <- sqrt(sum((mod-dataG)^2)/( sum(dataG^2))) * wG
else
rssG <- 0
##
gQ <- scaleSq*gQ
if(!is.na(dataS[1]))
rssS <- sqrt(sum((gQ-dataS)^2)/(sum(dataS^2)))
else
rssS <- 0
##
if(!is.na(dataSAS[1])){
dataSAS <- dataSAS - bkg
gQ_SAS <- scaleSASq*gQ_SAS
if(SASscale=="normal")
rssSAS <- sqrt(sum((gQ_SAS-dataSAS)^2)/(sum(dataSAS^2))) * wSAS
if(SASscale=="log")
rssSAS <- sqrt(sum((log(abs(gQ_SAS))-log(abs(dataSAS)))^2)/sum( (log(abs(dataSAS)))^2 )) * wSAS
if(SASscale=="weighted")
rssSAS <- sqrt(sum((gQ_SAS/dataSAS-1)^2)/(length(dataSAS))) * wSAS
}
else
rssSAS <- 0
if(is.na(rssS)) rssS <- 1e15
if(is.na(rssSAS)) rssSAS <- 1e15
if(is.na(rssG)) rssG <- 1e15
include <- logical(3)
include[1] <- !is.na(dataSAS[1])
include[2] <- !is.na(dataG[1])
include[3] <- !is.na(dataS[1])
res <- numeric(3)
if(include[1])
res[1] <- rssSAS
else
res[1] <- 0
if(include[2])
res[2] <- rssG
else
res[2] <- 0
if(include[3])
res[3] <- rssS
else
res[3] <- 0
if(pareto)
res[which(include==TRUE)]
else
sqrt(res[1]^2 + res[2]^2 + res[3]^2)
}
Try the nanop 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.