# I use the function knn.reg to evaluate the numerical solution at a fixed grid of points
#
if( !require('FNN') ){
    install.packages('FNN', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(FNN)

source('../Chapter2/euler.R')

t0 = 0
y0 = 1
t1 = 1.0

dy.dt = function(t, y){ cos(5*pi*t) }

h=0.2
res = euler(dy.dt, h=h, start=t0, y0=y0, end=t1)
t_grid_1 = seq( t0, 0.6, by=h )
e_yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( t_grid_1 ), k=1, algorithm='brute' )
print('Part (b):')
print(e_yhat_h1$pred[-1])


h=0.1
res = euler(dy.dt, h=h, start=t0, y0=y0, end=t1)
t_grid_2 = seq( t0, 0.4, by=h )
e_yhat_h2 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( t_grid_2 ), k=1, algorithm='brute' )
print('Part (c):')
print(e_yhat_h2$pred[-1])


# Truth:
#
t_truth = seq( 0, 1, length.out=100 )
y_truth = sin(5*pi*t_truth)/(5*pi) + 1
y_truth_ds = knn.reg( as.matrix( t_truth ), y=y_truth, test=as.matrix(t_grid_2), k=1, algorithm='brute' )


#postscript("../../WriteUp/Graphics/Chapter8/chap_8_sect_1_prob_22_plot.eps", onefile=FALSE, horizontal=FALSE)
ymin = min( c(y_truth, e_yhat_h1$pred, e_yhat_h2$pred) )
ymax = max( c(y_truth, e_yhat_h1$pred, e_yhat_h2$pred) )

plot(t_truth, y_truth, type='l', xlab='t', ylab='y(t)', ylim=c(ymin, ymax), col='green')
lines(t_grid_1, e_yhat_h1$pred, type='b', col='blue', pch=19)
lines(t_grid_2, e_yhat_h2$pred, type='b', col='red', pch=19)
legend( 'topright', c('Truth', 'Euler (h=0.2)', 'Euler (h=0.1)'), lty=c(1, 1, 1), col=c('green', 'blue', 'red') )
grid()
#dev.off()