# loading source code
source('vrstFunctions.R')

# loading simple cell data (responses with a mean of 5 spikes per image)
show('loading data ...')
data <- loadData(xFilename='xSC.dat', yFilename='ySC_5spi.dat')

# first portion of the data is for estimating the models
xTrain <- data$x[1:5000,]
yTrain <- data$y[1:5000]

# try Volterra model of order 1 to 5
orders <- 1:5

# use 90% of the train data to estimate the models
quiet <- FALSE
show('Be patient please. For this data set, PPR takes around 9 minutes on a Pentium III, 750MHz, 750MB')
estimationRes <- estimateVolterraModels(x=xTrain[1:4500,], y=yTrain[1:4500],
                                                           nroRDs=2, maxTerms=6,
                                                           orders=orders,
                                                           quiet=quiet)

# use remaining 10% of the train data to select the order of the Volterra model
validationRes <- predictVolterraModels(models=estimationRes$volterraModels,
                                        rds=estimationRes$rds,
                                        x=xTrain[4501:5000,], 
                                        y=yTrain[4501:5000])$correlations

# second portion of the data is for testing the predictive power of the model
# that best achieved best predictive power in the previous step
order <- 3
volterraModel <- estimationRes$volterraModels[[order]]
xTest <- data$x[5001:7000,]
yTest <- data$y[5001:7000]
prediction_3rdOrder <- predictVolterraModel(coefs=volterraModel$coefs, 
                                             order=order, 
                                             rds=estimationRes$rds, 
                                             x=xTest, 
                                             y=yTest)$correlation

# show results
layout(matrix(c(1,1,2,3), 2, 2, byrow=TRUE))

# plot validation results
plot(orders, validationRes, xlab='order', ylab='correlation coefficients')
grid()
title('Simple Cell')

# plot 1st PPR relevant dimension 
image(matrix(estimationRes$rds[,1], 10, 10), axes=FALSE, asp=1)
contour(matrix(estimationRes$rds[,1], 10, 10), add = TRUE, drawlabels = FALSE)
title('1st rel. dim.')

# plot 2nd PPR relevant dimension 
image(matrix(estimationRes$rds[,2], 10, 10), axes=FALSE, asp=1)
contour(matrix(estimationRes$rds[,2], 10, 10), add = TRUE, drawlabels = FALSE)
title('2nd rel. dim.')

# display the prediction with the test data
show(sprintf('correlation coef., order=%d, test data: %f', as.integer(order),
             prediction_3rdOrder))