/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*

Copyright (C) 2007, 2008 Mark Joshi

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

<quantlib-dev@lists.sf.net>. The license is also available online at

<https://www.quantlib.org/license.shtml>.

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.

*/

#include <ql/math/matrixutilities/basisincompleteordered.hpp>

#include <algorithm>

namespace QuantLib {

    BasisIncompleteOrdered::BasisIncompleteOrdered(Size euclideanDimension)

        : euclideanDimension_(euclideanDimension) {}

    bool BasisIncompleteOrdered::addVector(const Array& newVector1) {

        QL_REQUIRE(newVector1.size() == euclideanDimension_,

            "missized vector passed to "

            "BasisIncompleteOrdered::addVector");

        newVector_ = newVector1;

        if (currentBasis_.size()==euclideanDimension_)

            return false;

        for (auto& currentBasi : currentBasis_) {

            Real innerProd =

                std::inner_product(newVector_.begin(), newVector_.end(), currentBasi.begin(), Real(0.0));

            for (Size k=0; k<euclideanDimension_; ++k)

                newVector_[k] -= innerProd * currentBasi[k];

        }

        Real norm = std::sqrt(std::inner_product(newVector_.begin(),

            newVector_.end(),

            newVector_.begin(), Real(0.0)));

        if (norm<1e-12) // maybe this should be a tolerance

            return false;

        for (Size l=0; l<euclideanDimension_; ++l)

            newVector_[l]/=norm;

        currentBasis_.push_back(newVector_);

        return true;

    }

    Size BasisIncompleteOrdered::basisSize() const {

        return currentBasis_.size();

    }

    Size BasisIncompleteOrdered::euclideanDimension() const {

        return euclideanDimension_;

    }

    Matrix BasisIncompleteOrdered::getBasisAsRowsInMatrix() const {

        Matrix basis(currentBasis_.size(), euclideanDimension_);

        for (Size i=0; i<basis.rows(); ++i)

            for (Size j=0; j<basis.columns(); ++j)

                basis[i][j] = currentBasis_[i][j];

        return basis;

    }

    namespace

    {

        Real normSquared(const Matrix& v, Size row)

        {

            Real x=0.0;

            for (Size i=0; i < v.columns(); ++i)

                x += v[row][i]*v[row][i];

            return x;

        }

        Real norm(const Matrix& v, Size row)

        {

            return std::sqrt(normSquared( v,  row));

        }

        Real innerProduct(const Matrix& v, Size row1, const Matrix& w, Size row2)

        {

            Real x=0.0;

            for (Size i=0; i < v.columns(); ++i)

                x += v[row1][i]*w[row2][i];

            return x;

        }

    }

    OrthogonalProjections::OrthogonalProjections(const Matrix& originalVectors,

                                                 Real multiplierCutoff,

                                                 Real tolerance)

    : originalVectors_(originalVectors),

      multiplierCutoff_(multiplierCutoff),

      numberVectors_(originalVectors.rows()),

      dimension_(originalVectors.columns()),

      validVectors_(true,originalVectors.rows()), // opposite way round from vector constructor

      orthoNormalizedVectors_(originalVectors.rows(),

                              originalVectors.columns())

    {

        std::vector<Real> currentVector(dimension_);

        for (Size j=0; j < numberVectors_; ++j)

        {

            if (validVectors_[j])

            {

                for (Size k=0; k< numberVectors_; ++k) // create an orthormal basis not containing j

                {

                    for (Size m=0; m < dimension_; ++m)

                        orthoNormalizedVectors_[k][m] = originalVectors_[k][m];

                    if ( k !=j && validVectors_[k])

                    {

                        for (Size l=0; l < k; ++l)

                        {

                            if (validVectors_[l] && l !=j)

                            {

                                Real dotProduct = innerProduct(orthoNormalizedVectors_, k, orthoNormalizedVectors_,l);

                                for (Size n=0; n < dimension_; ++n)

                                    orthoNormalizedVectors_[k][n] -=  dotProduct*orthoNormalizedVectors_[l][n];

                            }

                        }

                        Real normBeforeScaling= norm(orthoNormalizedVectors_,k);

                        if (normBeforeScaling < tolerance)

                        {

                            validVectors_[k] = false;

                        }

                        else

                        {

                            Real normBeforeScalingRecip = 1.0/normBeforeScaling;

                            for (Size m=0; m < dimension_; ++m)

                                orthoNormalizedVectors_[k][m] *= normBeforeScalingRecip;

                        } // end of else (norm < tolerance)

                    } // end of if k !=j && validVectors_[k])

                }// end of  for (Size k=0; k< numberVectors_; ++k)

                // we now have an o.n. basis for everything except  j

                Real prevNormSquared = normSquared(originalVectors_, j);

                for (Size r=0; r < numberVectors_; ++r)

                    if (validVectors_[r] && r != j)

                    {

                        Real dotProduct = innerProduct(orthoNormalizedVectors_, j, orthoNormalizedVectors_,r);

                        for (Size s=0; s < dimension_; ++s)

                           orthoNormalizedVectors_[j][s] -= dotProduct*orthoNormalizedVectors_[r][s];

                    }

               Real projectionOnOriginalDirection = innerProduct(originalVectors_,j,orthoNormalizedVectors_,j);

               Real sizeMultiplier = prevNormSquared/projectionOnOriginalDirection;

               if (std::fabs(sizeMultiplier) < multiplierCutoff_)

               {

                    for (Size t=0; t < dimension_; ++t)

                        currentVector[t] = orthoNormalizedVectors_[j][t]*sizeMultiplier;

               }

               else

                   validVectors_[j] = false;

            } // end of  if (validVectors_[j])

            projectedVectors_.push_back(currentVector);

        } //end of j loop

        numberValidVectors_ =0;

        for (Size i=0; i < numberVectors_; ++i)

            numberValidVectors_ +=  validVectors_[i] ? 1 : 0;

    } // end of constructor

    const std::valarray<bool>& OrthogonalProjections::validVectors() const

    {

        return validVectors_;

    }

    const std::vector<Real>& OrthogonalProjections::GetVector(Size index) const

    {

        return projectedVectors_[index];

    }

  Size OrthogonalProjections::numberValidVectors() const

  {

        return numberValidVectors_;

  }

}