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

/*
 Copyright (C) 2025 William Day

 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/exercise.hpp>
#include <ql/pricingengines/barrier/analyticsoftbarrierengine.hpp>
#include <ql/instruments/barrieroption.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <utility>
#include <ql/termstructures/yield/flatforward.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
#include <ql/time/daycounters/actual360.hpp>
#include <ql/pricingengines/barrier/analyticbarrierengine.hpp>
#include <iostream>


namespace QuantLib {

    AnalyticSoftBarrierEngine::AnalyticSoftBarrierEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process)
    : process_(std::move(process)) {
        registerWith(process_);
    }


    void AnalyticSoftBarrierEngine::calculate() const {

        // Market data
        Real S = underlying();
        Real X = strike();
        Rate r = riskFreeRate();
        Rate q = dividendYield();
        Volatility sigma = volatility();

        // Barrier parameters
        Real U = barrierHi();
        Real L = barrierLo();
        Barrier::Type barrierType = arguments_.barrierType;
        
        // Stability tweak for r and q
        const Real epsilon = 1e-6;  
        if (std::abs(r - q) < 1e-10) {
            r = q + epsilon;  // Avoids mu = 0.5 singularity
        }

        // Option parameters
        Time T = residualTime();
        ext::shared_ptr<PlainVanillaPayoff> payoff = ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 
        Option::Type optionType = payoff->optionType();
        Integer eta = (optionType == Option::Call ? 1 : -1);
        Rate b = r - q; // cost of carry

        validateInputs(S, X, r, q, T, U, L, optionType, barrierType, sigma);

        bool isKnockedIn = (barrierType == Barrier::DownIn && S <= L) || 
                          (barrierType == Barrier::UpIn && S >= U);
        bool isKnockedOut = (barrierType == Barrier::DownOut && S <= L) || 
                          (barrierType == Barrier::UpOut && S >= U);

        bool isSingleBarrier = (std::fabs(U - L) < 1e-4);  
        

        // edge case 1: fully knocked in options should be priced as vanilla (there are no more barrier features to consider)
        if (isKnockedIn) {
              results_.value = vanillaEquivalent();
              return;
            }

        // edge case 2: knocked out options are worthless
        else if (isKnockedOut) {   
            results_.value = 0.0;  
            return;
            }
        
        // edge case 3: Haug formula breaks when U=L, use single barrier option formula instead
        if (isSingleBarrier) {
            results_.value = standardBarrierEquivalent();
            return;
        }

        // soft barrier pricing logic
        Real w = knockInValue(S, X, r, sigma, T, U, L, b, optionType,eta);
        results_.value = (barrierType == Barrier::DownIn || barrierType == Barrier::UpIn)
            ? w                     // knock in price
            : vanillaEquivalent() - w;  // knock out price
        }
    

    // Implements the formula to calculate 'w' from the Haug textbook, used in soft barrier pricing
    Real AnalyticSoftBarrierEngine::knockInValue(Real S, Real X, Rate r, Volatility sigma, Time T,
                                                Real U, Real L, Real b, Option::Type optionType,
                                                Integer eta) const {
        // constant terms                                              
        const Real mu = (b + 0.5 * sigma * sigma) / (sigma * sigma);
        const Real sqrtT = std::sqrt(T);
        const Real lambda1 = std::exp(-0.5 * sigma * sigma * T * (mu + 0.5) * (mu - 0.5));
        const Real lambda2 = std::exp(-0.5 * sigma * sigma * T * (mu - 0.5) * (mu - 1.5));
        const Real SX = S * X;
        const Real logU2_SX = std::log((U * U) / SX);
        const Real logL2_SX = std::log((L * L) / SX);

        // d and e terms
        const Real d1 = logU2_SX / (sigma * sqrtT) + mu * sigma * sqrtT;
        const Real d2 = d1 - (mu + 0.5) * sigma * sqrtT;
        const Real d3 = logU2_SX / (sigma * sqrtT) + (mu - 1) * sigma * sqrtT;
        const Real d4 = d3 - (mu - 0.5) * sigma * sqrtT;

        const Real e1 = logL2_SX / (sigma * sqrtT) + mu * sigma * sqrtT;
        const Real e2 = e1 - (mu + 0.5) * sigma * sqrtT;
        const Real e3 = logL2_SX / (sigma * sqrtT) + (mu - 1) * sigma * sqrtT;
        const Real e4 = e3 - (mu - 0.5) * sigma * sqrtT;

        const Real Nd1 = f_(eta * d1);
        const Real Nd2 = f_(eta * d2);
        const Real Nd3 = f_(eta * d3);
        const Real Nd4 = f_(eta * d4);
        const Real Ne1 = f_(eta * e1);
        const Real Ne2 = f_(eta * e2);
        const Real Ne3 = f_(eta * e3);
        const Real Ne4 = f_(eta * e4);


        // term 1
        Real term1 = eta * S * std::exp((b - r) * T) * std::pow(S, -2.0 * mu)
            * std::pow(SX, mu + 0.5) / (2.0 * (mu + 0.5));


        term1 *= std::pow(U * U / SX, mu + 0.5) * Nd1 - lambda1 * Nd2
            - std::pow(L * L / SX, mu + 0.5) * Ne1 + lambda1 * Ne2;


        // term 2
        Real term2 = eta * X * std::exp(-r * T) * std::pow(S, -2.0 * (mu - 1))
            * std::pow(SX, mu - 0.5) / (2.0 * (mu - 0.5));


        term2 *= std::pow(U * U / SX, mu - 0.5) * Nd3 - lambda2 * Nd4
            - std::pow(L * L / SX, mu - 0.5) * Ne3 + lambda2 * Ne4;


        // final result
        Real w = (1.0 / (U - L)) * (term1 - term2);
        return w;
    }


    // helper function to check inputs are reasonable
    void AnalyticSoftBarrierEngine::validateInputs(Real S, Real X, Rate r, Rate q, Time T, Real U, Real L,
                                                   Option::Type optionType, Barrier::Type barrierType,
                                                   Real sigma) const {
        // Core Parameter checks                                                
        QL_REQUIRE(S > 0.0, "Spot price must be > 0");
        QL_REQUIRE(X > 0.0, "Strike price must be > 0");
        QL_REQUIRE(T > 0.0, "Option must have time to maturity > 0");
        QL_REQUIRE(sigma > 0, "Volatility must be > 0");
        QL_REQUIRE(optionType == Option::Call || optionType == Option::Put, "Invalid option type");                                       
        QL_REQUIRE(r <= 1.0 && r >= -0.05, "Interest rate must be between -5% and 100%");
        QL_REQUIRE(q <= 1.0 && q >= -0.1, "Dividend yield must be between -10% and 100%");

        
        // Barrier type checks
        QL_REQUIRE(
          barrierType == Barrier::DownIn ||
          barrierType == Barrier::DownOut ||
          barrierType == Barrier::UpIn ||
          barrierType == Barrier::UpOut,
          "Invalid barrier type");
        QL_REQUIRE(L != Null<Real>(), "no low barrier given");
        QL_REQUIRE(U != Null<Real>(), "no high barrier given");
        QL_REQUIRE(U > 0.0 && L > 0.0, "Barrier levels must be positive");
        QL_REQUIRE(U >= L, "Upper barrier must be greater than or equal to lower barrier");
        }
    

    // helper functions 
    Real AnalyticSoftBarrierEngine::underlying() const {
        return process_->x0();
    }

    Real AnalyticSoftBarrierEngine::strike() const {
        ext::shared_ptr<PlainVanillaPayoff> payoff = ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);  
        QL_REQUIRE(payoff, "non-plain payoff given");
        return payoff->strike();
    }

    Time AnalyticSoftBarrierEngine::residualTime() const {
        return process_->time(arguments_.exercise->lastDate());
    }

    Volatility AnalyticSoftBarrierEngine::volatility() const {
        return process_->blackVolatility()->blackVol(residualTime(), strike());
    }

    Real AnalyticSoftBarrierEngine::stdDeviation() const {
        return volatility() * std::sqrt(residualTime());
    }

    Real AnalyticSoftBarrierEngine::barrierLo() const {
        return arguments_.barrier_lo;
    }

    Real AnalyticSoftBarrierEngine::barrierHi() const {
        return arguments_.barrier_hi;
    }

    Rate AnalyticSoftBarrierEngine::riskFreeRate() const {
        return process_->riskFreeRate()->zeroRate(residualTime(), Continuous, NoFrequency);
    }

    DiscountFactor AnalyticSoftBarrierEngine::riskFreeDiscount() const {
        return process_->riskFreeRate()->discount(residualTime());
    }

    Rate AnalyticSoftBarrierEngine::dividendYield() const {
        return process_->dividendYield()->zeroRate(residualTime(),Continuous, NoFrequency);
    }

    DiscountFactor AnalyticSoftBarrierEngine::dividendDiscount() const {
        return process_->dividendYield()->discount(residualTime());
    }
            

    Real AnalyticSoftBarrierEngine::vanillaEquivalent() const {
        ext::shared_ptr<StrikedTypePayoff> payoff =
            ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        Real forwardPrice = underlying() * dividendDiscount() / riskFreeDiscount();
        BlackCalculator black(payoff, forwardPrice, stdDeviation(), riskFreeDiscount());
        Real vanilla = black.value();
        return std::max(vanilla, 0.0);

    }       

    Real AnalyticSoftBarrierEngine::standardBarrierEquivalent() const {

    ext::shared_ptr<StrikedTypePayoff> payoff =
        ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
    QL_REQUIRE(payoff, "Payoff could not be cast to StrikedTypePayoff");

    BarrierOption tempOption(
        arguments_.barrierType,
        arguments_.barrier_hi,
        0.0,
        payoff,
        arguments_.exercise
    );

    Real spotVal = underlying();
    Real qVal = dividendYield();
    Real rVal = riskFreeRate();
    Volatility volVal = volatility();

    Handle<Quote> spot(ext::make_shared<SimpleQuote>(spotVal));
    Handle<YieldTermStructure> q(ext::make_shared<FlatForward>(0, TARGET(), qVal, Actual360()));
    Handle<YieldTermStructure> r(ext::make_shared<FlatForward>(0, TARGET(), rVal, Actual360()));
    Handle<BlackVolTermStructure> vol(ext::make_shared<BlackConstantVol>(0, TARGET(), volVal, Actual360()));


    ext::shared_ptr<GeneralizedBlackScholesProcess> process =
        ext::make_shared<GeneralizedBlackScholesProcess>(spot, q, r, vol);
    tempOption.setPricingEngine(ext::make_shared<AnalyticBarrierEngine>(process));

    Real npv = tempOption.NPV();
    Real result = std::max(npv, 0.0);
    return result;
}

}