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

/*
 Copyright (C) 2006 StatPro Italia srl

 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
 <http://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.
*/

// ===========================================================================
// NOTE: The following copyright notice applies to the original code,
//
// Copyright (C) 2002 Peter Jäckel "Monte Carlo Methods in Finance".
// All rights reserved.
//
// Permission to use, copy, modify, and distribute this software is freely
// granted, provided that this notice is preserved.
// ===========================================================================

#include <ql/methods/montecarlo/brownianbridge.hpp>

namespace QuantLib {

    BrownianBridge::BrownianBridge(Size steps)
    : size_(steps), t_(size_), sqrtdt_(size_),
      bridgeIndex_(size_), leftIndex_(size_), rightIndex_(size_),
      leftWeight_(size_), rightWeight_(size_), stdDev_(size_) {
        for (Size i=0; i<size_; ++i)
            t_[i] = static_cast<Time>(i+1);
        initialize();
    }

    BrownianBridge::BrownianBridge(const std::vector<Time>& times)
    : size_(times.size()), t_(times), sqrtdt_(size_),
      bridgeIndex_(size_), leftIndex_(size_), rightIndex_(size_),
      leftWeight_(size_), rightWeight_(size_), stdDev_(size_) {
        initialize();
    }

    BrownianBridge::BrownianBridge(const TimeGrid& timeGrid)
    : size_(timeGrid.size()-1), t_(size_), sqrtdt_(size_),
      bridgeIndex_(size_), leftIndex_(size_), rightIndex_(size_),
      leftWeight_(size_), rightWeight_(size_), stdDev_(size_) {
        for (Size i=0; i<size_; ++i)
            t_[i] = timeGrid[i+1];
        initialize();
    }


    void BrownianBridge::initialize() {

        sqrtdt_[0] = std::sqrt(t_[0]);
        for (Size i=1; i<size_; ++i)
            sqrtdt_[i] = std::sqrt(t_[i]-t_[i-1]);

        // map is used to indicate which points are already constructed.
        // If map[i] is zero, path point i is yet unconstructed.
        // map[i]-1 is the index of the variate that constructs
        // the path point # i.
        std::vector<Size> map(size_, 0);

        //  The first point in the construction is the global step.
        map[size_-1] = 1;
        //  The global step is constructed from the first variate.
        bridgeIndex_[0] = size_-1;
        //  The variance of the global step
        stdDev_[0] = std::sqrt(t_[size_-1]);
        //  The global step to the last point in time is special.
        leftWeight_[0] = rightWeight_[0] = 0.0;
        for (Size j=0, i=1; i<size_; ++i) {
            // Find the next unpopulated entry in the map.
            while (map[j])
                ++j;
            Size k = j;
            // Find the next populated entry in the map from there.
            while (!map[k])
                ++k;
            // l-1 is now the index of the point to be constructed next.
            Size l = j + ((k-1-j)>>1);
            map[l] = i;
            // The i-th Gaussian variate will be used to set point l-1.
            bridgeIndex_[i] = l;
            leftIndex_[i]   = j;
            rightIndex_[i]  = k;
            if (j != 0) {
                leftWeight_[i]= (t_[k]-t_[l])/(t_[k]-t_[j-1]);
                rightWeight_[i] = (t_[l]-t_[j-1])/(t_[k]-t_[j-1]);
                stdDev_[i] =
                    std::sqrt(((t_[l]-t_[j-1])*(t_[k]-t_[l]))
                              /(t_[k]-t_[j-1]));
            } else {
                leftWeight_[i]  = (t_[k]-t_[l])/t_[k];
                rightWeight_[i] =  t_[l]/t_[k];
                stdDev_[i] = std::sqrt(t_[l]*(t_[k]-t_[l])/t_[k]);
            }
            j=k+1;
            if (j>=size_)
                j=0;    //  wrap around
        }
    }

}