Integrate/Differentiation

#include<stdio.h>

#include<stdlib.h>

#include<time.h>

#include<math.h>

#define RND_MAX 2147480000.0

int main()

{

time_t seconds;

time(&seconds);

srand((unsigned int) seconds);

int i;

float x,y;

for(i=-10;i<=10;i=i+1)

{

/*

if(i%2)

printf("%d 1\n%d 2\n%d 3\n",i,i+1,i+2);

//printf("1\n1\n%d\n",rand()%4);

else

printf("%d 4\n%d 3\n%d 2\n",i,i+1,i+2);

//printf("4\n3\n%d\n",rand()%15);

*/

//printf("%d\n",i*i*i+3*i*i+2);

//printf("%d\n",3*i*i+6*i);

printf("%d\n",6*i+6);

}

return 0;

}

i<- -5:5

x<-i*i*i+3*i*i+5*i-50

fx_spline <- splinefun(x,i)

plot(i,x, type='l' )

#lines(i,fx_spline(x, deriv=1), type='l', col='green')

#lines(i, fx_spline(x, deriv=2), type='l', col='red')

lines(i,diff(x), type='l', col='green')

#lines(i, fx_spline(x, deriv=2), type='l', col='red')

lines(i,3*i*i+6*i+5, type='l', col='blue')

lines(i,6*i+6, type='l', col='orange')

# numerical integration in R

# example based on: http://tolstoy.newcastle.edu.au/R/help/04/10/6138.html

# sample data: FTIR spectra

x <- read.csv(url('http://casoilresource.lawr.ucdavis.edu/drupal/files/fresh_li_material.CSV'), header=FALSE)[100:400,]

names(x) <- c('wavenumber','intensity')

# fit a piece-wise linear function

fx.linear <- approxfun(x$wavenumber, x$intensity)

# integrate this function, over the original limits of x

Fx.linear <- integrate(fx.linear, min(x$wavenumber), max(x$wavenumber))

# fit a smooth spline, and return a function describing it

fx.spline <- splinefun(x$wavenumber, x$intensity)

# integrate this function, over the original limits of x

Fx.spline <- integrate(fx.spline, min(x$wavenumber), max(x$wavenumber))

# visual check, linear and spline fits shifted up for clarity

par(mar=c(0,0,0,0))

plot(x, type = "p", las=1, axes=FALSE, cex=0.5, ylim=c(0,0.12))

lines(x$wavenumber, fx.linear(x$wavenumber) + 0.01, col=2)

lines(x$wavenumber, fx.spline(x$wavenumber) + 0.02, col=3)

grid(nx=10, col=grey(0.5))

legend(x=615, y=0.11, legend=c('original','linear','spline'), col=1:3, pch=c(1,NA,NA), lty=c(NA, 1, 1), bg='white')

# results are pretty close

data.frame(method=c('linear', 'spline'), area=c(Fx.linear$value, Fx.spline$value),error=c(Fx.linear$abs.error,Fx.spline$abs.error))

method area error

1 linear 27.71536 0.0005727738

2 spline 27.71527 0.0030796529

par(mar=c(0,0,0,0), mfcol=c(2,1))

plot(x, type = "l", lwd=2, axes=FALSE)

grid(nx=10, col=grey(0.5))

plot(x$wavenumber, fx.spline(x$wavenumber, deriv=1), type='l', axes=FALSE)

lines(x$wavenumber, fx.spline(x$wavenumber, deriv=2), col='red')

grid(nx=10, col=grey(0.5))

abline(h=0, lty=2)

legend('topright', legend=c('1st derivative','2nd derivative'), lty=1, col=1:2, bg='white')

source : http://casoilresource.lawr.ucdavis.edu/drupal/node/896

features <- function(x, ymat, smoother=c("smooth.spline", "glkerns"), fits.return=FALSE, control=list( ), ...)

# Extracting features of a noisy and discretely sampled function.

# Author: Ravi Varadhan, October, 2010

#

# Function arguments:

# x = vector of independent variable, e.g. time

# ymat = matrix of time-series to be smoothed

# smoother = technique to use for automatic smoothing; only two options are allowed

# fits.return = whether to return the smoothed function and its first and second derivatives

#

# control = a list of control variables describe below

# npts = number of points to use in computing features

# plot.it = a logical variable indicating whether the smoothed results should be plotted

# c.outlier = a constant denoting number of standard deviations away from smooth fit for determining whether a point is an outlier

# decim.out = number of decimals to display in the features output

#

# Extracted features:

# frms = Root-mean-squared function value over the range of data

# fmean = Mean function value over the range of data

# fsd = square-root of the average variance of the function around its mean

# noise = standard deviation of residuals after smoothing

# snr = signal-to-noise ratio = fsd/noise

# fwiggle = root-mean-squared second-derivative of function; a measure of "wiggliness" in the function

# deriv.range = minimum and maximum of first derivative of smoothed function scaled by fsd

# critical.pts = X-locations where the first derivative is zero

# curvature = second derivative of smoothed function at critical points ( > 0 for minimum, and < 0 for maximum)

# scaled by fsd

# outliers = points that are potential outliers

{

ctrl <- list(npts=max(100, 2*length(x)), plot.it=FALSE, c.outlier=3, decim.out=2)

namc <- names(control)

if (!all(namc %in% names(ctrl)))

stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])

ctrl[namc] <- control

npts <- ctrl$npts

plot.it <- ctrl$plot.it

c.outlier <- ctrl$c.outlier

decim.out <- ctrl$decim.out

trapezoid <- function(x,y) sum(diff(x)*(y[-1] + y[-length(y)]))/2

smoother <- match.arg(smoother, c("smooth.spline", "glkerns"))

ipkgs <- installed.packages()

if (smoother == "glkerns") {

if ("lokern" %in% ipkgs[,1]) require(lokern, quietly=TRUE)

else stop("Install package `lokern'", call.=FALSE)

deriv1 <- function(z) glkerns(x, y, deriv=1, x.out=z, hetero=TRUE)$est

deriv2 <- function(z) glkerns(x, y, deriv=2, x.out=z, hetero=TRUE)$est

}

if (smoother == "smooth.spline") {

deriv1 <- function(z) predict(fit, deriv=1, x=z)$y

deriv2 <- function(z) predict(fit, deriv=2, x=z)$y

}

x.out <- seq(min(x), max(x), length=npts)

if (is.null(dim(ymat))) ymat <- matrix(ymat, nrow=1, ncol=length(ymat))

nr <- nrow(ymat)

nc <- ncol(ymat)

fmat <- matrix(NA, nr, 10)

cpmat <- matrix(NA, nr, max(ceiling(nc/10), 100))

cvmat <- matrix(NA, nr, max(ceiling(nc/10), 100))

olmat <- matrix(NA, nr, max(ceiling(nc/10), 100))

if (plot.it & nr > 1) par(ask=TRUE)

if (fits.return) fits <- fits1 <- fits2 <- matrix(NA, nr, nc)

for (k in 1:nr) {

y <- as.numeric(ymat[k, ])

if (smoother == "glkerns") {

fit <- glkerns(x, y, deriv=0, x.out=x, hetero=FALSE)

fit0 <- glkerns(x, y, deriv=0, x.out=x.out, hetero=TRUE)$est

fit1 <- glkerns(x, y, deriv=1, x.out=x.out, hetero=TRUE)$est

fit2 <- glkerns(x, y, deriv=2, x.out=x.out, hetero=TRUE)$est

resid <- y - fit$est

resid.scaled <- abs(scale(resid))

if (fits.return) {

fits[k, ] <- fit$est

fits1[k, ] <- glkerns(x, y, deriv=1, x.out=x, hetero=TRUE)$est

fits2[k, ] <- glkerns(x, y, deriv=2, x.out=x, hetero=TRUE)$est

}

}

if (smoother == "smooth.spline") {

fit <- smooth.spline(x, y, cv=FALSE)

fit0 <- predict(fit, deriv=0, x=x.out)$y

fit1 <- predict(fit, deriv=1, x=x.out)$y

fit2 <- predict(fit, deriv=2, x=x.out)$y

resid <- y - predict(fit)$y

resid.scaled <- abs(scale(resid))

if (fits.return) {

fits[k, ] <- predict(fit)$y

fits1[k, ] <- predict(fit, deriv=1)$y

fits2[k, ] <- predict(fit, deriv=2)$y

}

}

fmean <- trapezoid(x.out, fit0) / diff(range(x.out))

fstar <- sqrt(trapezoid(x.out, fit0^2) / diff(range(x.out)))

fsd <- sqrt(fstar^2 - fmean^2)

noise <- attr(resid.scaled, "scaled:scale")

snr <- fsd/noise

fwiggle <- sqrt(trapezoid(x.out, fit2^2) / diff(range(x.out)))

d0 <- (fit1[-1] * fit1[-npts]) < 0

ncrt <- sum(d0)

if (ncrt == 0) crtpts <- curv <- NULL else {

crtpts <- curv <- rep(NA, ncrt)

ind0 <- (1:npts)[d0]

for (i in 1:ncrt) {

temp <- try(uniroot(interval=c(x.out[ind0[i]], x.out[1 + ind0[i]]), f=deriv1), silent=TRUE)

if (class(temp) != "try-error") {

crtpts[i] <- temp$root

curv[i] <- deriv2(temp$root)

}

}

}

rm(fit)

outl <- resid.scaled > c.outlier

if (sum(outl) == 0) outliers <- NULL else outliers <- x[outl]

#plot

if (plot.it) {

plot(x, y, type="p", ...)

lines(x.out, fit0, col=2)

}

if (!is.null(crtpts) ) {

cpmat[k, 1:length(crtpts)] <- crtpts

cvmat[k, 1:length(crtpts)] <- curv

}

if (!is.null(outliers) ) olmat[k, 1:length(outliers)] <- outliers

fmat[k, ] <- as.numeric(c(fmean, as.numeric(range(fit0)), fsd, noise, snr, as.numeric(range(fit1)), fwiggle, length(crtpts)))

}

fmat <- as.data.frame(fmat)

names(fmat) <- c("fmean", "fmin", "fmax", "fsd", "noise", "snr", "d1min", "d1max", "fwiggle", "ncpts")

ncpmat <- max(apply(cpmat, 1, function(x) sum(!is.na(x))))

ncvmat <- max(apply(cvmat, 1, function(x) sum(!is.na(x))))

nolmat <- max(apply(olmat, 1, function(x) sum(!is.na(x))))

if (ncpmat==0) cpmat <- NULL else cpmat <- round(cpmat[, 1:ncpmat, drop=FALSE], decim.out)

if (ncvmat==0) cvmat <- NULL else cvmat <- round(cvmat[, 1:ncvmat, drop=FALSE], decim.out)

if (nolmat==0) olmat <- NULL else olmat <- olmat[, 1:nolmat, drop=FALSE]

ret.obj <- list(fmat=round(fmat, decim.out), cptmat=cpmat, curvmat=cvmat, olmat=olmat)

if (fits.return) attr(ret.obj, "fits") <- list(x=x, fits0=fits, fits1=fits1, fits2=fits2)

return(ret.obj)

}

---------------------------------------------------------------------------------------------

op <- par(mfrow = c(3, 1))

source("features.R")

#x <- sort(runif(n))

#y <- exp(-0.2 * sin(2*pi*x)) + rnorm(n, sd=0.05)

a<-read.table("WEBM1")

x<- seq(1, 2000, by = 1)

y<-a$V1

#y<-a$V2

#x<- seq(-5, 5, by = 0.01)

#y<- x*x*x

ans1 <- features(x, y, fits.return=TRUE, control=list(plot.it=FALSE))

#ans <- features(x, y, fits.return=TRUE, control=list(plot.it=TRUE))

fits1 <- attr(ans1, "fits")

#yexact <- -0.2 * 2*pi * cos(2*pi*x) * exp(-0.2 * sin(2*pi*x))

#lines(x, yexact, col=3)

plot(x,y,type="l",col=2,ylab="Original data",main="Video Trace")

grid(col=grey(0.3))

plot(x, fits1$fits1, type="l", xlab="x",col=3, ylab="First derivative")

#plot(fits1$x, fits1$fits1, type="l", xlab="x",col=3, ylab="First derivative")

grid(col=grey(0.3))

ans2 <- features(x, fits1$fits1, fits.return=TRUE, control=list(plot.it=FALSE))

fits2 <- attr(ans2, "fits")

plot(x, fits2$fits1, type="l", xlab="x",col=3, ylab="2nd derivative")

grid(col=grey(0.3))

ans3 <- features(x, fits2$fits1, fits.return=TRUE, control=list(plot.it=FALSE))

fits3 <- attr(ans3, "fits")

#plot(x, fits3$fits1, type="l", xlab="x",col=3, ylab="3rd derivative")

#grid(col=grey(0.3))

par(op)

-------------------------------------------------------------------------------------

op <- par(mfrow = c(2, 1))

x<-scan("a15_m4_mpeg")

y<- seq(1, 299, by = 1)

source("features.R")

ans1 <- features(x, y, fits.return=TRUE, control=list(plot.it=FALSE))

plot(x,y,type="l",col=2,ylab="Original data")

grid(col=grey(0.3))

par(op)

-------------------------------------------------------------------------------------

op <- par(mfrow = c(3, 1))

x <- seq(0, 9, by = 0.1)

# u is the true function, u1-u3 are its symbolic derivatives wrt x

# u1 = f', u2 = f'', u3 = f'''

#u <- expression(x*x*x+3*x*x+5*x-50)

u <- expression(sin(pi * (x - 0.5)))

u1 <- as.expression(D(u, 'x'))

u2 <- as.expression(D(u1, 'x'))

u3 <- as.expression(D(u2, 'x'))

#plot( eval(u), type = 'l',main = expression(f(x) == sin(pi * (x - 0.5))))

plot(eval(u), type = 'l', ylim = c(-30, 30),ylab = "", xlab = "",main = expression(f(x) == sin(pi * (x - 0.5))))

lines( eval(u1), type = 'l', col = 'red')

lines(eval(u2), type = 'l', col = 'blue')

lines( eval(u3), type = 'l', col = 'green')

grid(col=grey(0.3))

#legend('topright', legend = c("f", "f'", "f''", "f'''"),col = c('black', 'red', 'blue', 'green'), lwd = 2)

# y2 is an evaluation of u at each x point

y2 <- sin(pi * (x - 0.5))

#y2 <- x*x*x+3*x*x+5*x-50

# Construct the spline function approximation to u

f <- splinefun(x, y2)

# Plot the interpolation function and its first three derivatives

curve(f(x), 1, 10, col = "black", lwd = 1.5, ylim = c(-30, 30), xlim = c(2,11),ylab = "", main = 'Spline interpolation of f')

curve(f(x, deriv=1), 1, 10, col = 'red', lwd = 1.5, add = TRUE)

curve(f(x, deriv=2), 1, 10, col = 'blue', lwd = 1.5, add = TRUE)

curve(f(x, deriv=3), 1, 10, col = 'green', lwd = 1.5, add = TRUE)

grid(col=grey(0.3))

x<- seq(0, 9, by = 0.01)

plot(sin(pi * (x - 0.5)), type = 'l', ylim = c(-30, 30),ylab = "", xlab = "",main ='Direct Derivatives')

lines((cos(pi * (x - 0.5))*pi), type='l', col='red')

lines((-1)*(sin(pi * (x - 0.5))*(pi*pi)), type='l', col='blue')

lines((-1)*(cos(pi * (x - 0.5))*(pi*pi*pi)), type='l', col='green')

grid(col=grey(0.3))

par(op)