### R code from vignette source 'vignettes/IPPD/inst/doc/IPPD.Rnw'

###################################################
### code chunk number 1: setup
###################################################

options(prompt = "R> ", continue = " ", warn = -1, width = 90)                                       



###################################################
### code chunk number 2: loadmyo500
###################################################

library(IPPD) 
data(myo500)
x <- myo500[,"mz"]
y <- myo500[,"intensities"]
y <- y[x <= 2500]
x <- x[x <= 2500]


###################################################
### code chunk number 3: plotmyo500first1000
###################################################

layout(matrix(c(1,2), 1, 2))

plot(x[1:1000], y[1:1000], xlab = expression(x[1]~~ldots~~x[1000]), 
     cex.lab = 1.5, cex.axis = 1.25, ylab = expression(y))

plot(x[x >= 804 & x <= 807], y[x >= 804 & x <= 807], 
     xlab = "x: 804 <= x <= 807", 
     cex.lab = 1.5, cex.axis = 1.25, ylab = expression(y), type = "b")

layout(matrix(1))



###################################################
### code chunk number 4: fitGauss
###################################################

fitGauss <- fitModelParameters(mz = x, intensities = y,
                   model = "Gaussian", fitting = "model", formula.sigma = formula(~mz),
                   control = list(window = 6, threshold = 200))


###################################################
### code chunk number 5: fitEMG
###################################################

fitEMG <- fitModelParameters(mz = x, intensities = y,
                   model = "EMG", fitting = "model",
                           formula.alpha = formula(~mz),
                           formula.sigma = formula(~mz),
                           formula.mu = formula(~1),
                           control = list(window = 6, threshold = 200))


###################################################
### code chunk number 6: assessfit
###################################################

show(fitEMG)
mse.EMG <- data.frame(mse = slot(fitEMG,"peakfitresults")[,"rss"] 
                      / slot(fitEMG,"peakfitresults")[,"datapoints"], 
                      peakshape = rep("EMG", nrow( slot(fitEMG,"peakfitresults"))))
mse.Gauss <- data.frame(mse = slot(fitGauss,"peakfitresults")[,"rss"] 
                        / slot(fitGauss,"peakfitresults")[,"datapoints"], 
                        peakshape = rep("Gaussian", nrow( slot(fitGauss,"peakfitresults"))))
mses <- rbind(mse.EMG, mse.Gauss)
with(mses, boxplot(mse ~ peakshape, cex.axis = 1.5, cex.lab = 1.5, ylab = "MSE"))



###################################################
### code chunk number 7: assessfitEMGpeak
###################################################
visualize(fitEMG, type = "peak", cex.lab = 1.5, cex.axis = 1.25)


###################################################
### code chunk number 8: assessfit2
###################################################

visualize(fitEMG, type = "model", modelfit = TRUE, 
          parameters = c("sigma", "alpha"), 
          cex.lab = 1.5, cex.axis = 1.25)



###################################################
### code chunk number 9: getlist
###################################################

EMGlist <- getPeaklist(mz = x, intensities = y, model = "EMG",
model.parameters = fitEMG, 
loss = "L2", trace = FALSE,
control.localnoise = list(factor.place = 2),
control.basis = list(charges = c(1, 2)),
control.postprocessing = list(ppm = 200))

show(EMGlist)



###################################################
### code chunk number 10: showlist
###################################################

threshold(EMGlist, threshold = 3, refit = TRUE, trace = FALSE)



###################################################
### code chunk number 11: plot1
###################################################

visualize(EMGlist, x, y, lower= 963, upper = 973,
           fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
           localnoise = TRUE, quantile = 0.5,
           cutoff.functions = 3)



###################################################
### code chunk number 12: plot2
###################################################


visualize(EMGlist, x, y, lower= 1502, upper = 1510,
           fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
           localnoise = TRUE, quantile = 0.5,
           cutoff.functions = 2)




###################################################
### code chunk number 13: plot3
###################################################

visualize(EMGlist, x, y, lower= 1360, upper = 1364,
          fit = FALSE, fittedfunction = TRUE, fittedfunction.cut = TRUE,
          localnoise = TRUE, quantile = 0.5, 
          cutoff.functions = 2)