Function to calculate the RSS between the simulated PDF or the total scattering structure function and target data

Share:

Description

Given a nanoparticle model and PDF or total scattering structure function data, this function calculates the PDF or total scattering structure function associated with the model. It then determines the residual sum of squares (RSS) between the simulated signals and the data.

Usage

1
2
3
4
5
rss(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)					 

Arguments

par

numeric vector of parameter values. These should be given along with the appropriate argument skel. See details and examples.

dataG

numeric vector, which, if not NA, determines the PDF data.

dataS

numeric vector, which, if not NA, determines the total scattering structure function data.

dataSAS

numeric vector, which, if not NA, determines the total scattering structure function data in SAS region.

part

either NA or the atomic positions in a nanoparticle given as a numeric matrix in which each row gives the (3D) position coordinates. A matrix of this form is the return value of the function simPart. Default value of NA means that the particle is to be re-simulated in each call to rss (necessary for stochastic model functions).

type

character indicating type of scattering. Either "X-ray" or "neutron".

simPar

arguments are passed to the function simPart to simulate the nanoparticles if these values are not to be optimized; any of the following arguments may be specified: atoms, atomsShell, sym,symShell, latticep, latticepShell, r, rcore, box, ellipse, shell, rcenter, move, rotShell, rcenterShell and distr. The latter should be set up to "normal" or "lognormal" to simulate polydisperse particles with normal or log-normal distribution, respectively. See simPart for the other parameters notation details.

PDF.fixed

arguments are passed to the functions calcPDF and broadPDF. Any of the following arguments may be specified: dr, minR, maxR, scatterLength, scatterFactor, delta, n, Qmax, termRip, maxRTermRip, N1, N2, and Qdamp; see calcPDF and GrSAS. If termRip is TRUE termination ripples are simulated with parameters dr, Qmax and maxRTermRip (see termRip). If Qdamp is not NA the PDF is multiplied by the function exp(-Q_{damp}^2 r^2/2).

TotalScatt.fixed

arguments are passed to the functions calcTotalScatt; any of the following arguments may be specified: dQ, minQ, maxQ, minQ_SAS, maxQ_SAS, dQ_SAS, scatterLength, scatterFactor, delta, n, dr, del, eps, kind, N1, N2, f2, f2, convolution, Qdamp, SASscale, and paramSASQ; see calcTotalScatt. If convolution is TRUE function gaussConvol is called with Qdamp parameter. If paramSASQ is TRUE the total scattering function in SAS region is calculated using functions IqSAS or IqSASP (with f1 and f2 parameters indicating average scattering lengths in the particle core and shell). Parameter SASscale sets weighting scheme for total scattering function in SAS region residual. Should be "normal", "log" (for logarithmic scale) or "weighted" (to normalize the residual at each grid point with the function value).

parscale

either NA or a numeric vector of the same length as par indicating values by which the values in par should be divided for the purpose of parameter scaling.

skel

an object of class relistable. First the parameters to be optimized should be written as a named list of form parameters<-list(a=1,b=2,c=3). Then the argument par can be given as par=unlist(parameters) and skel can be given as skel=as.relistable(parameters).

con

logical indicating whether to reset r to be equal to rcore if rcore > r.

oneDW

logical indicating whether to use similar sigma parameter for all atom types within the particle.

punish

if the inequalities described above for con are violated, return a large RSS value (10e15).

gammaR

either "param" or a numeric vector of gamma baseline term data. If "param" the baseline is simulated with functions GrSAS and GrSASCS.

rvector

numeric vector of grid points at which the PDF was evaluated.

fixed

list; if not NA indicates parameter names.

wG, wSAS

weights for the PDF and total scattering function in SAS region residuals.

avRes

numeric value; if the model for the data is stochastic, he number of particles to simulate; the model PDF or total scattering structure function is calculated for each particle and then averaged.

pareto

logical; if TRUE the result value is a vector of residuals for the given data sets that can be used for the subsequent Pareto analysis.

Details

The following parameters may be fitted in rss (i.e., specified in par or skel): r, rsigma (sets standard deviation in the normal or log-normal distribution with mean value r), box, ellipse, bkg (background line that should be extracted from the SAS data), shell, rcore, scaleSq (scale factor for the total scattering function or the PDF), scaleSASq (scale factor for the total scattering function in SAS region or the gamma baseline term), pStack, pDimer, delta, sigma, latticep and latticepShell. For wurtzite structure to fit u parameter (second atom z-coordinate) third element in the vectors latticep and latticepShell may be specified.

Value

If pareto=FALSE the R-factor: the square root of the sum of squared distances between the target data and simulated signal divided by a the sum of squared target data values. If pareto=TRUE a vector with R-factors as components.

Note

Functions

simPartPar(sym = "fcc", symShell=NA, latticep = 4.08, r=10,

atoms=list(), atomsShell=list(), distr = "lognormal",

rcenter=FALSE, latticepShell=NA, rcore=NA, shell=NA, move=TRUE,

box=NA, ellipse=NA, rotShell=FALSE, rcenterShell=FALSE),

PDFPar(dr=.01, minR=1, maxR=20, termRip=FALSE, Qmax=30, maxRTermRip=20,

scatterLength=NA, n=2, delta=NA, Qdamp=NA, scatterFactor=NA, N1=4, N2=4),

and

TotalScattPar(dQ=.01, minQ=0.771, maxQ=35, minQ_SAS=0.001, maxQ_SAS=0.771,

dQ_SAS=0.005, scatterLength=NA, scatterFactor=NA, dr = 0.001,

del = 0.01, eps=1e-3, kind="fastHist", N1=4, N2=4, f1=5, f2=5,

SASscale="normal", convolution = FALSE, Qdamp=0.0457, delta=NA, n=2, paramSASQ=FALSE)

return values suitable for simPar, PDF.fixed and TotalScatt.fixed parameters, respectively.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
part <- res_gQ <- res_gQ_WAS <- res_gQ_SAS <- list() 

## prepare polydisperse sample consisting of 30 particles:
size <- sort(rlnorm(20, meanlog = log(10), sdlog = log(1.1)))
base_Cu <- getBase("Cu")
Au <- createAtom("Au", base=base_Cu)

for(i in 1:20) {
  cat("r = ", size[i], "\n")
  part[[i]] <- simPart(atoms=list(Au), r=size[i], latticep=4.08, 
      rcenter=TRUE) 
  res_gQ_WAS[[i]] <- calcTotalScatt(part[[i]],minQ=.771, maxQ=20, 
      dr=0.02, dQ=0.02, sigma=c(0.02))							
  res_gQ_SAS[[i]] <- calcTotalScatt(part[[i]],minQ=.001, 
      dr=0.02, maxQ=.771, dQ=0.01, sigma=c(0.02))
  cat(i,"\n")
}

## calculate average values:
gQ_av_WAS <- 0
for(i in 1:length(res_gQ_WAS)) {
  gQ_av_WAS <- res_gQ_WAS[[i]]$gQ + gQ_av_WAS
}
gQ_av_WAS <- gQ_av_WAS/length(res_gQ_WAS)

gQ_av_SAS <- 0
for(i in 1:length(res_gQ_SAS)) {
  gQ_av_SAS <- res_gQ_SAS[[i]]$gQ + gQ_av_SAS
}
gQ_av_SAS <- gQ_av_SAS/length(res_gQ_SAS)

## calculate PDF and gamma baseline term
resSAS  <- calcQDepPDF(minR=0, maxR=30, dr=0.02, minQ=.001,
    maxQ=.771, verbose=100, 
    preTotalScat=list(Q=res_gQ_SAS[[1]]$Q,gQ=gQ_av_SAS))
  
mr <- which(res_gQ_WAS[[1]]$Q > 17)[1]
mxr <- which(res_gQ_WAS[[1]]$Q > 19)[1]
cuto <- res_gQ_WAS[[1]]$Q[mr:mxr][which(abs(gQ_av_WAS[mr:mxr])
    ==min(abs(gQ_av_WAS[mr:mxr]) ))[1] ]
resWAS <- calcQDepPDF(minR=0, maxR=30, dr=0.02, minQ=.771, 
    maxQ=cuto, verbose=100, 
    preTotalScat=list(Q=res_gQ_WAS[[1]]$Q,gQ=gQ_av_WAS))

## set boundaries for fitting procedure
iter_0 <- as.relistable(list(latticep=0, r=0, sigma=c(0), rsigma=0))
boundsL <- c(latticep=3.5,  r=10.0, sigma=c(.01), rsigma=1.1)
boundsU <- c(latticep=4.5,  r=14.0,  sigma=c(.03), rsigma=1.3)

## in order to estimate the parameters that were used to
## simulate the particles, the DEoptim package may be
## used. Install it, remove the comment symbols '#' below,
## and use a call like:
#library(DEoptim)

#resDE <- DEoptim(rss, lower = boundsL, upper = boundsU, 
#    oneDW=FALSE, type="neutron",                   
#    control=DEoptim.control(CR=0.85, F=0.7, NP=40, itermax=30, 
#        parallelType=1, packages = list("PerformanceAnalytics"),
#        parVar=list("rss", "simPart", "calcPDF", 
#        "IqSASP", "GrSAS", "GrSASCS", "broadPDF", "termRip",
#        "getSymEl", "IqSAS", "displacePart", "calcTotalScatt")),       
#    dataS=NA, dataSAS=gQ_av_SAS, dataG=resWAS$gr, 
#    simPar = simPartPar(atoms=list(Au), rcenter=TRUE, 
#        move=TRUE, rot=FALSE),
#    gammaR = resSAS$gr, rvector = resSAS$r, skel=iter_0, 
#    PDF.fixed=PDFPar(minR=0, maxR=30, dr=.01, 
#        scatterLength=Au$scatterLength),
#    TotalScatt.fixed=TotalScattPar(minQ=0.771, maxQ=20, 
#        dQ=0.02, dQ_SAS=0.01, minQ_SAS=.001, maxQ_SAS=.771, 
#        scatterLength=Au$scatterLength, kind="fastHist", 
#        SASscale="normal", convolution = FALSE),
#    wG=1.0, wSAS=0.05, avRes=10, pareto=FALSE) 
#
## now resDE$optim contains estimates for the lattice parameter,
## particle radius, displacement variance sigma, and standard 
## deviation rsigma.
## show results:
#resDE$optim$bestmem

## package mco may be used to construct a Pareto front of 
## solutions. Note that calculations may take considerable time!

#library(mco)
#resMCO <- nsga2(rss, 4, 2, 
#    dataS=NA, dataSAS=gQ_av_SAS, dataG=resWAS$gr,  oneDW=FALSE,
#    simPar = simPartPar(rcenter=TRUE, move=TRUE, rot=FALSE),
#    gammaR = resSAS$gr, rvector = resSAS$r, skel=iter_0,
#    PDF.fixed=PDFPar(minR=0, maxR=30, dr=.01, 
#        scatterLength=scatterLength),
#    TotalScatt.fixed=TotalScattPar(minQ=0.771, maxQ=20, dQ_SAS=0.01, 
#        dQ=.02, minQ_SAS=.001, maxQ_SAS=.771, scatterLength=scatterLength,
#        kind="fastHist", SASscale="normal", convolution = FALSE),
#    wG=1.0, wSAS=0.2, avRes=10, pareto=TRUE,							 
#    constraints=NULL, cdim=0, popsize=40, generations=c(40), 
#    cprob=0.85, lower.bounds=boundsL, upper.bounds=boundsU)

## show results
#plot(resMCO, xlab="SAS(Q) residual", ylab="G(r) residual")			
#paretoSet(resMCO)