estimatePPRRDs <- function(x, y, nTerms, maxTerms=6, quiet=FALSE) {
    timeRes <- system.time(rds <- ppr(x=x, y=y, nterms=nTerms,
                                           max.terms=maxTerms)$alpha)
    if(!quiet) show(sprintf('ppr used %.4f seconds of user CPU time',
                    timeRes[1]))
    rds <- orthogonalizeMatrix(rds)
    return(rds)
}
estimateVolterraModel <- function(x, y, rds, order) {
    # estimate coefficients of the low-dimensional Volterra problem
    coefs <- estimateCoefsTruncatedPInv(x=x, y=y, rds=rds, order=order)

    # reconstruct the first and second order kernels
    k1 <- getK1FromCoefsAndRDs(coefs=coefs, rds=rds)
    k2 <- getK2FromCoefsAndRDs(coefs=coefs, rds=rds)

    return(list(coefs=coefs, k1=k1, k2=k2))
}
estimateCoefsTruncatedPInv <- function(x, y, rds, order) {
    # project stimuli in relevant dimensions
    px <- x%*%rds

    # build data matrix of Volterra linear equation
    dataMatrix <- buildDataMatrix(px=px, order=order, nRDs=ncol(rds))

    # solve the Volterra linear equation 
    coefs <- solveTruncatedPInv(x=dataMatrix, y=y)

    return(coefs)
}
predictVolterraModel <- function(coefs, order, rds, x, y) {
    px <- x%*%rds
    dataMatrix <- buildDataMatrix(px=px, rds=rds, order=order, nRDs=ncol(rds))
    predictions <- dataMatrix%*%coefs
    correlation <- cor(y, predictions)
    mse <- mean((y-predictions)^2)
    return(list(predictions=predictions, correlation=correlation, mse=mse))
}
buildDataMatrix <- function(px, rds, order, nRDs) {
    totalNroCoefs <- getTotalNroCoefs(order=order, nRDs=nRDs)
    dataMatrix <- matrix(nrow=nrow(px), ncol=totalNroCoefs)
    index <- 1
    # zeroth order term
    dataMatrix[, index] <- rep(1, nrow(dataMatrix))
    index <- index+1

    # first order terms
    for(j in 1:nRDs) {
        dataMatrix[, index] <- px[, j]
        index <- index+1
    }
    if(order<2) {
        return(dataMatrix)
    }

    # second order terms
    items <- array(dim=2)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
            index <- index+1
        } 
    }
    if(order<3) {
        return(dataMatrix)
    }
        
    # third order terms
    items <- array(dim=3)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                index <- index+1
            }
        }
    } 
    if(order<4) {
        return(dataMatrix)
    }
        
    # fourth order terms
    items <- array(dim=4)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                    index <- index+1
                }
            }
        }
    } 
    if(order<5) {
        return(dataMatrix)
    }
        
    # fifth order terms
    items <- array(dim=5)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    for(j5 in j4:nRDs) {
                        items[5] <- j5
                        dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                        index <- index+1
                    }
                }
            }
        }
    } 
    if(order<6) {
        return(dataMatrix)
    }
        
    # sixth order terms
    items <- array(dim=6)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    for(j5 in j4:nRDs) {
                        items[5] <- j5
                        for(j6 in j5:nRDs) {
                            items[6] <- j6
                            dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                            index <- index+1
                        }
                    }
                }
            }
        }
    } 
    if(order<7) {
        return(dataMatrix)
    }
    
    # seventh order terms
    items <- array(dim=7)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
        for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    for(j5 in j4:nRDs) {
                        items[5] <- j5
                        for(j6 in j5:nRDs) {
                            items[6] <- j6
                            for(j7 in j6:nRDs) {
                                items[7] <- j7
                                dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                                   index <- index+1
                            }
                        }
                    }
                }
            }
        }
    } 
    if(order<8) {
        return(dataMatrix)
    }
    
    # eight order terms
    items <- array(dim=8)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    for(j5 in j4:nRDs) {
                        items[5] <- j5
                        for(j6 in j5:nRDs) {
                            items[6] <- j6
                            for(j7 in j6:nRDs) {
                                items[7] <- j7
                                for(j8 in j7:nRDs) {
                                    items[8] <- j8
                                    dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                                    index <- index+1
                                }
                            }
                        }
                    }
                }
            }
        }
    } 
    if(order<9) {
        return(dataMatrix)
    }
    
    # ninth order terms
    items <- array(dim=9)
    for(j1 in 1:nRDs) {
        items[1] <- j1
        for(j2 in j1:nRDs) {
            items[2] <- j2
            for(j3 in j2:nRDs) {
                items[3] <- j3
                for(j4 in j3:nRDs) {
                    items[4] <- j4
                    for(j5 in j4:nRDs) {
                        items[5] <- j5
                        for(j6 in j5:nRDs) {
                            items[6] <- j6
                            for(j7 in j6:nRDs) {
                                items[7] <- j7
                                for(j8 in j7:nRDs) {
                                    items[8] <- j8
                                    for(j9 in j8:nRDs) {
                                        items[9] <- j9
                                        dataMatrix[, index] <-
                    computeSymmetriesCoef(items)*
                    computeProductCols(px, items)
                                        index <- index+1
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    } 
    if(order>9) {
        stop('order <= 9')
    }
    return(dataMatrix)
}
getK1FromCoefsAndRDs <- function(coefs, rds) {
    q1 <- getQ1(nrds=ncol(rds), coefs=coefs)
    k1 <- getK1(rds=rds, q1=q1)
}
getK2FromCoefsAndRDs <- function(coefs, rds) {
    q2 <- getQ2(nrds=ncol(rds), coefs=coefs)
    kernel <- rds%*%q2%*%t(rds)
    return(kernel)
}
solveTruncatedPInv <- function(x=x, y=y, percentageLargerSV=1e-10) {
    pInv <- truncatedPInv(x, percentageLargerSV)
    coefs <- pInv%*%y
    return(coefs)
}
estimateVolterraModels <- function(x, y, nroRDs, maxTerms=6, orders=1:5,
                                      quiet=FALSE) {

    if(!quiet) {
        show('estimating relevant dimensions ...')
    }
    rds <- estimatePPRRDs(x=x, y=y, nTerms=nroRDs, maxTerms=maxTerms,
                               quiet=quiet)
    if(!quiet) {
        show('estimating Volterra models ...')
    }
    volterraModels <- estimateVolterraModelsWithRDs(x=x, y=y, rds=rds, 
                                                         orders=orders,
                                                         quiet=quiet)
    return(list(rds=rds, volterraModels=volterraModels))
}
estimateVolterraModelsWithRDs <- function(x, y, rds, orders, quiet=FALSE) {
    result <- list()
    i <- 1
    for(order in orders) {
        if(!quiet) {
            show(sprintf('estimating volterra model %d ...', as.integer(i)))
            i <- i+1
        }
        result <- c(result, list(c(order=order, 
                                    estimateVolterraModel(x=x, y=y, rds=rds, 
                                                               order=order))))
    }
    return(result)
}
getK1 <- function(rds, q1) {
    kernel <- rds%*%q1
    return(kernel)
}
getNroCoefsForTerm <- function(order, nRDs) {
    return(choose(nRDs+order-1, order))
}
getQ1 <- function(nrds, coefs) {
    nCoefs0thOrder <- 1
    nCoefs1rstOrder <- nrds
    q1 <- coefs[(nCoefs0thOrder+1):(nCoefs0thOrder+1+nCoefs1rstOrder-1)]
    return(q1)
}
getQ2 <- function(nrds, coefs) {
    matrix <- matrix(nrow=nrds, ncol=nrds)
    index <- 1
    nCoefs0thOrder <- 1
    nCoefs1rstOrder <- nrds
    nCoefs2ndOrder <- (nrds+1)*nrds/2
    coefs2ndOrder <- coefs[(nCoefs0thOrder+nCoefs1rstOrder+1):
                           (nCoefs0thOrder+nCoefs1rstOrder+1+nCoefs2ndOrder-1)]

    if(nrds>1) {
        for(i1 in 1:(nrds-1)) {
            matrix[i1, i1] <- coefs2ndOrder[index]
            index <- index+1
            for(i2 in (i1+1):nrds) {
                matrix[i1, i2] <- coefs2ndOrder[index]
                matrix[i2, i1] <- coefs2ndOrder[index]
                index <- index+1
            }
        }
    }
    matrix[nrds, nrds] <- coefs2ndOrder[index]
    return(matrix)
}
getTotalNroCoefs <- function(order, nRDs) {
    nroUnknowns <- 1
    for(termOrder in 1:order) {
        nroUnknowns <- nroUnknowns + getNroCoefsForTerm(termOrder, nRDs)
    }
    return(nroUnknowns)
}
computeSymmetriesCoef <- function(indices) {
    multiplicities <- as.vector(1)
    lastDifferentIndex <- indices[1]

    for(i in 2:length(indices)) {
        if(indices[i]==lastDifferentIndex) {
            multiplicities[length(multiplicities)] =
                multiplicities[length(multiplicities)] + 1
        } else {
            lastDifferentIndex <- indices[i]
            multiplicities = cbind(multiplicities, as.vector(1))
        }
    }
    coef <- factorial(length(indices))
    for(i in 1:length(multiplicities)) {
        coef <- coef/factorial(multiplicities[i])
    }
    return(coef)
}
predictVolterraModels <- function(models, rds, x, y) {
    predictions <- vector()
    correlations <- vector()
    mses <- vector()
    for(model in models) {
        res <- predictVolterraModel(model$coefs, model$order, rds, x, y)
        predictions <- cbind(predictions, res$predictions)
        correlations <- c(correlations, res$correlation)
        mses <- c(mses, res$mse)
    }
    return(list(predictions=predictions, correlations=correlations, mses=mses))
}
truncatedPInv <- function(x, percentageLargerSV) {
    svdRes <- svd(x)
    dInv <- rep(0, length(svdRes$d))
    nroSingularValues <-
        length(svdRes$d[svdRes$d>(svdRes$d[1]*percentageLargerSV)])
    for(i in 1:nroSingularValues) {
        dInv[i] <- 1/svdRes$d[i]
    }
    truncatedPInv <- svdRes$v%*%diag(dInv, nrow=length(dInv))%*%t(svdRes$u)
    return(truncatedPInv)
}
computeProductCols <- function(matrix, indices) {
    i <- 1
    product <- matrix[,indices[i]]
    while(i < length(indices)) {
        i <- i+1
        product <- product*matrix[,indices[i]]
    }
    return(product)
}
orthogonalizeMatrix <- function(m) {
    qrRes <- qr(m)
    return(qr.Q(qrRes))
}
loadData <- function(xFilename='xSC.dat', yFilename='ySC_5spi.dat') {
    x <- as.matrix(read.table(xFilename))
    y <- as.vector(as.matrix(read.table(yFilename)))
    return(list(x=x, y=y))
}