/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.analytics.financial.model.option.pricing.fourier;

import static com.opengamma.analytics.math.ComplexMathUtils.add;
import static com.opengamma.analytics.math.ComplexMathUtils.divide;
import static com.opengamma.analytics.math.ComplexMathUtils.log;
import static com.opengamma.analytics.math.ComplexMathUtils.mod;
import static com.opengamma.analytics.math.ComplexMathUtils.multiply;
import static com.opengamma.analytics.math.ComplexMathUtils.sqrt;
import static com.opengamma.analytics.math.ComplexMathUtils.subtract;
import static com.opengamma.analytics.math.number.ComplexNumber.I;

import org.apache.commons.lang.NotImplementedException;

import com.opengamma.analytics.math.TrigonometricFunctionUtils;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.number.ComplexNumber;

/**
 * The Cox-Ingersoll-Ross process is a mean-reverting positive process, with SDE:
 * $$
 * \begin{align*}
 * dy_t = \kappa(\theta - y_t)dt + \lambda\sqrt{y_t}dW_t
 * \end{align*}
 * $$
 * and characteristic exponent
 * $$
 * \begin{align*}
 * \psi(u, t; \kappa, \theta, \lambda) &= \frac{2\kappa\theta}{\lambda^2}\left[ 
 * \frac{\kappa t}{2} - \ln\left(\cosh\left(\frac{\gamma t}{2}\right) + \frac{\kappa}{\gamma}\sinh\left(\frac{\gamma t}{2}\right)\right)
 * + \frac{2iu}{\kappa + \gamma \coth\left(\frac{\gamma t}{2}\right)}\right]\\
 * \text{where}\\
 * \gamma &= \sqrt{\kappa^2 - 2 \lambda^2 iu}
 * \end{align*}
 * $$
 */
public class IntegratedCIRTimeChangeCharacteristicExponent implements StocasticClockCharcteristicExponent {
  private final double _kappa;
  private final double _theta;
  private final double _lambda;
  private final double _alphaMax;

  /**
   * 
   * @param kappa The mean-reverting speed
   * @param theta The mean 
   * @param lambda The volatility
   */
  public IntegratedCIRTimeChangeCharacteristicExponent(final double kappa, final double theta, final double lambda) {
    _kappa = kappa;
    _theta = theta;
    _lambda = lambda;
    _alphaMax = _kappa * _kappa / 2 / _lambda / _lambda;
  }

  @Override
  public Function1D<ComplexNumber, ComplexNumber> getFunction(final double t) {
    return new Function1D<ComplexNumber, ComplexNumber>() {
      @Override
      public ComplexNumber evaluate(final ComplexNumber u) {
        return getValue(u, t);
      }
    };
  }
  
  @Override
  public ComplexNumber getValue(ComplexNumber u, double t) {
    if (u.getReal() == 0.0 && u.getImaginary() == 0.0) {
      return new ComplexNumber(0.0);
    }

    final ComplexNumber ui = multiply(I, u);

    //handle small lambda properly 
    if (2 * mod(u) * _lambda * _lambda / _kappa / _kappa < 1e-6) {
      final double d = _theta * t + (1 - _theta) * (1 - Math.exp(-_kappa * t)) / _kappa;
      return multiply(d, ui);
    }

    ComplexNumber temp = subtract(_kappa * _kappa, multiply(2 * _lambda * _lambda, ui));
    final ComplexNumber gamma = sqrt(temp);
    final ComplexNumber gammaHalfT = multiply(gamma, t / 2.0);
    temp = divide(multiply(2, ui), add(_kappa, divide(gamma, TrigonometricFunctionUtils.tanh(gammaHalfT))));
    final ComplexNumber kappaOverGamma = divide(_kappa, gamma);
    final double power = 2 * _kappa * _theta / _lambda / _lambda;
    final ComplexNumber res = add(multiply(power, subtract(_kappa * t / 2, getLogCoshSinh(gammaHalfT, kappaOverGamma))), temp);
    return res;
  }

  // ln(cosh(a) + bsinh(a)
  private ComplexNumber getLogCoshSinh(final ComplexNumber a, final ComplexNumber b) {
    final ComplexNumber temp = add(TrigonometricFunctionUtils.cosh(a), multiply(b, TrigonometricFunctionUtils.sinh(a)));
    return log(temp);
  }

  /**
   * 
   * @return $\frac{\kappa^2}{2\lambda^2}$
   */
  @Override
  public double getLargestAlpha() {
    return _alphaMax;
  }

  /**
   * 
   * @return $-\infty$
   */
  @Override
  public double getSmallestAlpha() {
    return Double.NEGATIVE_INFINITY;
  }

  /**
   * Gets the mean-reverting speed.
   * @return kappa
   */
  public double getKappa() {
    return _kappa;
  }

  /**
   * Gets the mean.
   * @return theta
   */
  public double getTheta() {
    return _theta;
  }

  /**
   * Gets the volatility.
   * @return lambda
   */
  public double getLambda() {
    return _lambda;
  }

  @Override
  public int hashCode() {
    final int prime = 31;
    int result = 1;
    long temp;
    temp = Double.doubleToLongBits(_kappa);
    result = prime * result + (int) (temp ^ (temp >>> 32));
    temp = Double.doubleToLongBits(_lambda);
    result = prime * result + (int) (temp ^ (temp >>> 32));
    temp = Double.doubleToLongBits(_theta);
    result = prime * result + (int) (temp ^ (temp >>> 32));
    return result;
  }

  @Override
  public boolean equals(final Object obj) {
    if (this == obj) {
      return true;
    }
    if (obj == null) {
      return false;
    }
    if (getClass() != obj.getClass()) {
      return false;
    }
    final IntegratedCIRTimeChangeCharacteristicExponent other = (IntegratedCIRTimeChangeCharacteristicExponent) obj;
    if (Double.doubleToLongBits(_kappa) != Double.doubleToLongBits(other._kappa)) {
      return false;
    }
    if (Double.doubleToLongBits(_lambda) != Double.doubleToLongBits(other._lambda)) {
      return false;
    }
    return Double.doubleToLongBits(_theta) == Double.doubleToLongBits(other._theta);
  }

  @Override
  public ComplexNumber[] getCharacteristicExponentAdjoint(ComplexNumber u, double t) {
    throw new NotImplementedException();
  }

  @Override
  public Function1D<ComplexNumber, ComplexNumber[]> getAdjointFunction(double t) {
    throw new NotImplementedException();
  }

 

}