/*
 Copyright (C) 2003 Niels Elken Sønderby

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it under the
 terms of the QuantLib license.  You should have received a copy of the
 license along with this program; if not, please email ferdinando@ametrano.net
 The license is also available online at http://quantlib.org/html/license.html

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/*
 NesQuant is an extension to QuantLib
 http://www.nielses.dk/quantlib/nesquant
*/

#include <nq/Math/polynomialfit.hpp>

namespace NesQuant {

    const Array PolynomialFit::value() {
        perform();
        return value_;
    }

    const Array PolynomialFit::parameters() {
        perform();
        return parameters_;
    }

    void PolynomialFit::perform() {
        if (!isPerformed_) {
            try {
                // SVD can only handle mxn matrices where m>=n
                bool useTranspose = false;
                if (data_.rows() < fitFunctions_)
                    useTranspose = true;

                // perform the fit using singular value decomposition

                // make design matrix M
                Matrix designMat(data_.rows(), fitFunctions_);

                for (Size i = 0; i < data_.rows(); ++i) {
                    designMat[i][0] = 1;
                    Size j = 1;
                    for (Size k = 0; k < data_.columns() - 1; ++k) {
                        double powx = data_[i][k];
                        designMat[i][j] = powx;
                        j++;
                        for (Size l = 2; l <= degree_; ++l) {
                            powx *= data_[i][k];
                            designMat[i][j] = powx;
                            j++;
                        }
                    }
                }

                if (useTranspose)
                    designMat = transpose(designMat);

                // do singular value decomposition
                SVD svd(designMat);
                Matrix U, UT, V;
                Array s;
                svd.getU(U);
                UT = transpose(U);
                svd.getV(V);
                svd.getSingularValues(s);

                // create diagonalmatrix S = diag(1/s)
                Matrix S(s.size(), s.size(), 0.0);
                for (i = 0; i < s.size(); ++i) {
                    if (QL_FABS(s[i]) > QL_EPSILON)
                        S[i][i] = 1 / s[i];
                    else
                        S[i][i] = 0.0;
                }

                Array f_(data_.rows());
                for (i = 0; i < f_.size(); ++i) {
                    f_[i] = data_[i][data_.columns() - 1];
                }

                if (!useTranspose) {
                    // calculate parameters a = V diag(1/s) UT y
                    parameters_ = V * S * UT * f_;
                    // calculate value v = M a
                    value_ = designMat * parameters_;
                } else {
                    parameters_ = U * S * transpose(V) * f_;
                    value_ = transpose(designMat) * parameters_;
                }

                isPerformed_ = true;
            } catch (...) {
                // if an error occurs unset flag and rethrow error
                isPerformed_ = false;
                throw;
            }
        }
    }

}