myprojects-hg-reborn

Fork

rm(list=ls()) #first I want to clear the workspace by getting rid of 
              #the previous variables which may interfere with the computations
              #I am going to carry out
library(minpack.lm)

f <- function(x, A1,A2,mu1,mu2,myvar1,myvar2) {
    expr <- expression(
log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
         +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2))
         )
    eval(expr) ### expression I want to use for my NLS
}



j <- function(x, A1,A2,mu1,mu2,myvar1,myvar2) {
    expr <- expression(
log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
         +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2))
         )

    c(eval(D(expr, "A1")),   
      eval(D(expr, "A2" )),
eval(D(expr, "mu1" )),
eval(D(expr, "mu2" )),
eval(D(expr, "myvar1" )),
eval(D(expr, "myvar2" )))
}

#the function above contains the Jacobian of the function I want to use to carry out my fittings

fcn  <- function(p, x, N, N.Err, fcall, jcall)
    (N - do.call("fcall", c(list(x = x), as.list(p))))/N.Err



## Now I read the experimental data
luisa<-read.csv("volume-distribution.csv",header=FALSE)
luisadata<-as.matrix(luisa)

vecdim<-dim(luisadata)
Niter<-vecdim[2]-1

results<-matrix(nrow=Niter,ncol=9)

for (p in 1:Niter)

{

x<-luisadata[ ,1]
N<-luisadata[ ,(p+1)] ###vector with the experimental data

fcn     <- function(p, x, N, N.Err, fcall, jcall)
    (N - do.call("fcall", c(list(x = x), as.list(p))))/N.Err
  # Now I define the function I want
  # to minimize as the difference between 
  # the simulated data and the function
  #I plane to use to model them (eventually). 

fcn.jac <- function(p, x, N, N.Err, fcall, jcall) {
    N.Err <- rep(N.Err, length(p))                   
    -do.call("jcall", c(list(x = x), as.list(p)))/N.Err   
}

#Above I also do the same for the Jacobian

guess <- c("A1" = 85, "A2"=23,"mu1"=430,"mu2"=1670,
         "myvar1"=1.59,"myvar2"=1.5)

#the guess vector contains the initial ansatz for the fitting parameters

#I consider the fact that N.Err has to be related to the sqrt(N) since we
# are dealing with least-squares optimization.

myN.Err<-seq(1:length(x))
for (i in 1:length(x))
{
if (N[i] !=0)

{
myN.Err[i]<-sqrt(N[i]) 
}
else
{
myN.Err[i]<-1
}
}

 myN.Err<-1



# just a test to find out if I can modify the file or I have troubles
# kate is slow to save the changes but it does not crash ---at least so far

#initial ftol=0.5

out <- nls.lm(par = guess, fn = fcn,   jac = fcn.jac,
              fcall = f,   jcall = j,
              x = x, N = N, N.Err = myN.Err,
              control = list(nprint = 3,maxfew=200,factor=100,ftol=0.5, diag = numeric()))
            #   control = list(maxfew=200,factor=100,ftol=0.00001, diag = numeric()))

print("p is ")
print(p)
myvec<-out$par
results[p,1:6 ]<-myvec
results[p,7]<-myvec[1]/(myvec[1]+myvec[2])
results[p,8]<-myvec[2]/(myvec[1]+myvec[2])
results[p,9]<-myvec[1]+myvec[2]
}


# now I plot some results

T2<-seq(1:84)
for (i in 1:84)
{T2[i]<-i/2}
pdf("Mean-first-mode.pdf")
plot(T2,results[ ,3],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="Mu_1 [nm]",main="Fitted Mean for the First Mode")
dev.off()
pdf("Mean-second-mode.pdf")
plot(T2,results[ ,4],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="Mu_2 [nm]",main="Fitted Mean for the Second Mode")
dev.off()
pdf("Sigma-first-mode.pdf")
plot(T2,results[ ,5],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="Sigma_1 [nm]",main="Fitted Sigma for the First Mode")
dev.off()
pdf("Sigma-second-mode.pdf")
plot(T2,results[ ,6],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="Sigma_2 [nm]",main="Fitted Sigma for the Second Mode")
dev.off()
pdf("Weight-first-mode.pdf")
plot(T2,results[ ,7],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="w1",main="Fitted Weight for the First Mode")
dev.off()
pdf("Weight-second-mode.pdf")
plot(T2,results[ ,8],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="w2",main="Fitted Weight for the Second Mode")
dev.off()
pdf("Total-particle-number.pdf")
plot(T2,log(10)*results[ ,9],"l",col=300,lwd=2.5,xlab="Time [hours]",ylab="Particle Number",main="Total Particle Number")
dev.off()

#now I choose to print a selected fitting

mycol<-47
x<-luisadata[ ,1]
newx<-seq(min(x),max(x),len=1000)
g<-function(x,A1,A2,mu1,mu2,myvar1,myvar2)
{(log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
 +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2)))
        }
#mypar<-c(85,40.9,409,933,1.81,2.39)
mypar<-results[mycol,1:6]
N<-luisadata[ ,(mycol+1)]
pdf("Two-modes-nls-fitting-anomalous.pdf")
plot(x, N, bg = "black", pch = 21, cex = 1.5,xlab="Dp",ylab="dV/dlog(Dp)")
M <- do.call("g", c(list(x = newx), mypar))
lines(newx, M, col="black", lwd=1.)

h <- function(x, A1,mu1,myvar1) {
    expr <- expression(
log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1)))
    eval(expr) ### expression I want to use for my NLS
}

mypar1<-c(mypar[1],mypar[3],mypar[5])

M1 <- do.call("h", c(list(x = newx), mypar1))
lines(newx, M1, col="orange", lwd=1.)

mypar2<-c(mypar[2],mypar[4],mypar[6])
M2 <- do.call("h", c(list(x = newx), mypar2))
lines(newx, M2, col="green", lwd=1.)
dev.off()


print("So far so good")