/*
 Copyright (C) 2010 Marius Hofert and Martin Maechler

 This program is free software; you can redistribute it and/or modify it under
 the terms of the GNU General Public License as published by the Free Software
 Foundation; either version 3 of the License, or (at your option) any later 
 version.

 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 GNU General Public License for more
 details.

 You should have received a copy of the GNU General Public License along with
 this program; if not, see <http://www.gnu.org/licenses/>.
*/

#include <Rmath.h>

#include "nacopula.h"

/**< lang5 will probably be defined in future versions of R */
#ifndef lang5
static SEXP lang5(SEXP s, SEXP t, SEXP u, SEXP v, SEXP w)
{
    PROTECT(s);
    s = LCONS(s, list4(t, u, v, w));
    UNPROTECT(1);
    return s;
}
#endif

/**
 * Sample S0 ~ S(alpha, 1, (cos(alpha*pi/2))^{1/alpha}, I_{alpha == 1}; 1)
 * with Laplace-Stieltjes transform exp(-t^alpha) via the algorithm as described
 * in Chambers, Mallows, Stuck (1976) (S0 = S(alpha,1) in this reference), see
 * Nolan's book for the parametrization.
 * @param alpha parameter in (0,1]
 * @return S0
 * @author Marius Hofert, Martin Maechler
*/
double rstable0(double alpha){
	/**< S(1, 1, 0, 1; 1) with Laplace-Stieltjes transform exp(-t) */
	if(alpha == 1.) return 1.;

	/**< alpha in (0,1) */
	double U = unif_rand();
	double W;
	do { W = exp_rand(); } while (W == 0.);
	return pow(A_(M_PI*U,alpha)/pow(W,1.-alpha),1./alpha);
}

/**
 * Sample S ~ S(alpha, 1, 1, (1/cos(alpha*pi/2))*I_{alpha == 1}; 1)
 * with Laplace-Stieltjes transform exp(-(1/cos(alpha*pi/2))*t^alpha),
 * see Nolan's book for the parametrization.
 * @param alpha parameter in (0,1]
 * @return S
 * @author Marius Hofert, Martin Maechler
*/
double rstable(double alpha){
    if(alpha == 1.)
	return 1./0.; /**< Inf or NaN */
    else
	return pow(cos(M_PI_2*alpha), -1./alpha) * rstable0(alpha);
}

/**
 * Vectorize rstable. Generate a vector of random variates
 * S ~ S(alpha, 1, 1, (1/cos(alpha*pi/2))*I_{alpha == 1}; 1)
 * with Laplace-Stieltjes transform exp(-(1/cos(alpha*pi/2))*t^alpha),
 * see Nolan's book for the parametrization.
 * @param S vector of stable random variates (result)
 * @param n length of the vector S
 * @param alpha parameter in (0,1]
 * @return none
 * @author Marius Hofert, Martin Maechler
*/
void rstable_vec(double S[], const int n, const double alpha){
    if(n >= 1) {
	double cf = pow(cos(M_PI_2*alpha), -1./alpha);
	GetRNGstate();
	for(int i=0; i < n; i++)
	    S[i] = cf * rstable0(alpha);
	PutRNGstate();
    }
    return;
}

/**
 * Generate a vector of random variates
 * S ~ S(alpha, 1, 1, (1/cos(alpha*pi/2))*I_{alpha == 1}; 1)
 * with Laplace-Stieltjes transform exp(-(1/cos(alpha*pi/2))*t^alpha),
 * see Nolan's book for the parametrization.
 * @param n sample size
 * @param alpha parameter in (0,1]
 * @return vector of random variates S
 * @author Martin Maechler
*/
SEXP rstable_c(SEXP n, SEXP alpha)
{
    int nn = asInteger(n);
    SEXP res = PROTECT(allocVector(REALSXP, nn));

    rstable_vec(REAL(res), nn, asReal(alpha));
    UNPROTECT(1);
    return(res);
}

/**
 * Sample St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)*V_0)^{1/alpha},
 *			 V_0*I_{alpha = 1}, h*I_{alpha != 1}; 1)
 * with Laplace-Stieltjes transform exp(-V_0((h+t)^alpha-h^alpha)),
 * see Nolan's book for the parametrization, via the "fast rejection algorithm",
 * see Hofert (2010).
 * @param St vector of random variates (result)
 * @param V0 vector of random variates V0
 * @param h parameter in [0,infinity)
 * @param alpha parameter in (0,1]
 * @param n length of St
 * @return none
 * @author Marius Hofert, Martin Maechler
*/
void retstable_MH(double *St, const double V0[], double h, double alpha, int n){
    /**
     * alpha == 1 => St corresponds to a point mass at V0 with Laplace-Stieltjes
     * transform exp(-V0*t)
    */
    if(alpha == 1.){
	for(int i = 0; i < n; i++) {
		St[i] = V0[i];
	}
	return;
    }

    GetRNGstate();

    for(int i = 0; i < n; i++) { /**< for each of the n required variates */

	/**< find m := optimal number of summands, using asymptotic formula */
	int m;
	m = imax2(1, (int)round(V0[i] * pow(h,alpha)));

	/**< apply standard rejection m times and sum the variates St_k */
	double c = pow(V0[i]/m, 1./alpha);
	St[i] = 0.; /**< will be the result after the summation */
	for(int k = 0; k < m; k++){
	    double St_k, U;
	    /**
	     * apply standard rejection for sampling
	     * \tilde{S}(alpha, 1, (cos(alpha*pi/2)*V_0/m)^{1/alpha},
	     *           (V_0/m)*I_{alpha = 1}, h*I_{alpha != 1}; 1) with
	     * Laplace-Stieltjes transform exp(-(V_0/m)*((h+t)^alpha-h^alpha))
	    */
	    do {
		/**
		 * sample St_k ~ S(alpha, 1, (cos(alpha*pi/2)*V_0/m)^{1/alpha},
		 *		   (V_0/m)*I_{\alpha = 1}; 1) with
		 * Laplace-Stieltjes transform exp(-(V_0/m)*t^alpha)
		*/
		St_k = c * rstable0(alpha);
		U = unif_rand();
	    } while (U > exp(- h * St_k));
	    /**
	     * on acceptance, St_k ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)
	     *				       *V_0/m)^{1/alpha},
	     *				       (V_0/m)*I_{alpha = 1},
	     *				       h*I_{alpha != 1}; 1) with
	     * Laplace-Stieltjes transform exp(-(V_0/m)*((h+t)^alpha-h^alpha))
	    */
	    St[i] += St_k;
	}

    } /**< end for */

    PutRNGstate();
    return;
}

/**
 * Fast and accurate evaluation of sinc(x) := sin(x)/x, including the limit x=0.
 * @param x any (double precision) number
 * @return sinc(x)
 * @author Martin Maechler; 28 Apr 2010
*/
double sinc_MM(double x) {
    double ax = fabs(x);
    if(ax < 0.006) {
	if(x == 0.) return 1;
	double x2 = x*x;
	if(ax < 2e-4)
	     return 1. - x2/6.;
	else return 1. - x2/6.*(1 - x2/20.);
    }
    /**< else */
    return sin(x)/x;
}

/**< to be called from R */
SEXP sinc_c(SEXP x_) {
    int n = LENGTH(PROTECT(x_ = coerceVector(x_, REALSXP)));
    SEXP r_ = allocVector(REALSXP, n); /**< the result */
    double *x = REAL(x_), *r = REAL(r_);

    for(int i=0; i < n; i++)
	r[i] = sinc_MM(x[i]);

    UNPROTECT(1);
    return r_;
}

/**
 * Evaluation of Zolotarev's function, see Devroye (2009), to the power 1-alpha.
 * The 3-arg. version allows more precision for  alpha ~=~ 1
 * @param x argument
 * @param alpha parameter in (0,1]
 * @return sin(alpha*x)^alpha * sin((1-alpha)*x)^(1-alpha) / sin(x)
 * @author Martin Maechler; 28 Apr 2010
*/
#define _A_3(_x, _alpha_, _I_alpha)				\
    pow(_I_alpha* sinc_MM(_I_alpha*_x), _I_alpha) *		\
    pow(_alpha_ * sinc_MM(_alpha_ *_x), _alpha_ ) / sinc_MM(_x)

double A_(double x, double alpha) {
  double Ialpha = 1.-alpha;
  return _A_3(x, alpha, Ialpha);
}

/**< to be called from R---see experiments in ../tests/retstable-ex.R */
SEXP A__c(SEXP x_, SEXP alpha, SEXP I_alpha) {
    int n = LENGTH(PROTECT(x_ = coerceVector(x_, REALSXP)));
    double alp = asReal(alpha), I_alp = asReal(I_alpha);
    if(fabs(alp + I_alp - 1.) > 1e-12)
	error("'I_alpha' must be == 1 - alpha more accurately");
    SEXP r_ = allocVector(REALSXP, n); /**< the result */
    double *x = REAL(x_), *r = REAL(r_);

    for(int i=0; i < n; i++)
	r[i] = _A_3(x[i], alp, I_alp);

    UNPROTECT(1);
    return r_;
}

/**
 * Evaluation of B(x)/B(0), see Devroye (2009).
 * @param x argument
 * @param alpha parameter in (0,1]
 * @return sinc(x) / (sinc(alpha*x)^alpha * sinc((1-alpha)*x)^(1-alpha))
 * @author Martin Maechler; 28 Apr 2010
*/
double BdB0(double x,double alpha) {
    double Ialpha = 1.-alpha;
    double den = pow(sinc_MM(alpha*x),alpha) * pow(sinc_MM(Ialpha*x),Ialpha);
    return sinc_MM(x) / den;
}

/**
 * Sample St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)*V_0)^{1/alpha},
 *			 V_0*I_{alpha = 1}, h*I_{alpha != 1}; 1)
 * with Laplace-Stieltjes transform exp(-V_0((h+t)^alpha-h^alpha)),
 * see Nolan's book for the parametrization, via double rejection,
 * see Devroye (2009).
 * @param St vector of random variates (result)
 * @param V0 vector of random variates V0
 * @param h parameter in [0,infinity)
 * @param alpha parameter in (0,1]
 * @param n length of St
 * @return none
 * @author Marius Hofert, Martin Maechler
*/
void retstable_LD(double *St, const double V0[], double h, double alpha, int n)
{
    /**
     * alpha == 1 => St corresponds to a point mass at V0 with Laplace-Stieltjes
     * transform exp(-V0*t)
    */
    if(alpha == 1.){
	for(int i = 0; i < n; i++) {
		St[i] = V0[i];
	}
	return;
    }

    /**< compute variables not depending on V0 */
    const double c1 = sqrt(M_PI_2);
    const double c2 = 2.+c1;
    double b = (1.-alpha)/alpha;

    for(int i = 0; i < n; i++) { /**< for each of the n required variates */

	/**< set lambda for our parameterization */
	double V0alpha = pow(V0[i],1./alpha);
	double lambda = h*V0alpha;

	/**
	 * Apply the algorithm of Devroye (2009) to draw from
	 * \tilde{S}(alpha, 1, (cos(alpha*pi/2))^{1/alpha}, I_{alpha = 1},
	 * 	     lambda*I_{alpha != 1};1) with Laplace-Stieltjes transform
	 * exp(-((lambda+t)^alpha-lambda^alpha))
	*/
	double gamma = pow(lambda,alpha)*alpha*(1.-alpha);
	double sgamma = sqrt(gamma);
	double c3 = c2* sgamma;
	double xi = (1. + M_SQRT2 * c3)/M_PI; /**< according to John Lau */
	double psi = c3*exp(-gamma*M_PI*M_PI/8.)/M_SQRT_PI;
	double w1 = c1*xi/sgamma;
	double w2 = 2.*M_SQRT_PI * psi;
	double w3 = xi*M_PI;
	double X, c, E;
	do {
	    double U, z, Z;
	    do {
		double V = unif_rand();
		if(gamma >= 1) {
		    if(V < w1/(w1+w2)) U = fabs(norm_rand())/sgamma;
		    else{
			double W_ = unif_rand();
			U = M_PI*(1.-W_*W_);
		    }
		}
		else{
		    double W_ = unif_rand();
		    if(V < w3/(w2+w3)) U = M_PI*W_;
		    else U = M_PI*(1.-W_*W_);
		}
		double W = unif_rand();
		double zeta = sqrt(BdB0(U,alpha));
		double phi = pow(sgamma+alpha*zeta,1./alpha);
		z = phi/(phi-pow(sgamma,1./alpha));
		/**< compute rho */
		double rho = M_PI*exp(-pow(lambda,alpha)*(1.-1. \
							  /(zeta*zeta))) / \
							  ((1.+c1)*sgamma/zeta \
		   					  + z);
		double d = 0.;
		if(U >= 0 && gamma >= 1) d += xi*exp(-gamma*U*U/2.);
		if(U > 0 && U < M_PI) d += psi/sqrt(M_PI-U);
		if(U >= 0 && U <= M_PI && gamma < 1) d += xi;
		rho *= d;
		Z = W*rho;
	    } while( !(U < M_PI && Z <= 1.)); /**< check rejection condition */

	    double
		a = pow(A_(U,alpha), 1./(1.-alpha)),
		m = pow(b*lambda/a,alpha),
		delta = sqrt(m*alpha/a),
		a1 = delta*c1,
		a3 = z/a,
		s = a1+delta+a3;
	    double V_ = unif_rand(), N_ = 0., E_ = 0. /**< -Wall */;
	    if(V_ < a1/s) {
		N_ = norm_rand();
		X = m-delta*fabs(N_);
	    } else {
		if(V_ < (a1+delta)/s) /**< according to John Lau */
		    X = m+delta*unif_rand();
		else {
		    E_ = exp_rand();
		    X = m+delta+E_*a3;
		}
	    }
	    E = -log(Z);
	    /**< check rejection condition */
	    c = a*(X-m)+lambda*(pow(X,-b)-pow(m,-b));
	    if(X < m) c -= N_*N_/2.;
	    else if(X > m+delta) c -= E_;

	} while (!(X >= 0 && c <= E));
	/**
	 * Transform variates from the distribution corresponding to the
	 * Laplace-Stieltjes transform exp(-((lambda+t)^alpha-lambda^alpha))
	 * to those of the distribution corresponding to the Laplace-Stieltjes
	 * transform exp(-V_0((h+t)^alpha-h^alpha)).
	*/
	St[i] = V0alpha / pow(X,b);

    } /**< end for */
    return;
}

/**
 * Sample St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)*V_0)^{1/alpha},
 *			 V_0*I_{alpha = 1}, h*I_{alpha != 1}; 1)
 * with Laplace-Stieltjes transform exp(-V_0((h+t)^alpha-h^alpha)),
 * see Nolan's book for the parametrization, via the algorithms of
 * Devroye (2009) and Hofert (2010).
 * @param V0_ R vector of type numeric(n) containing the random variates V0
 * @param h R parameter of type numeric(1)
 * @param alpha R parameter in (0,1] of type numeric(1)
 * @param method R "string" of type character(1)
 * @return R vector of type numeric(n) containing the random variates
 * @author Martin Maechler
 */
SEXP retstable_c(SEXP V0_, SEXP h, SEXP alpha, SEXP method)
{
    int n = LENGTH(PROTECT(V0_ = coerceVector(V0_, REALSXP)));
    const char* meth_ch = CHAR(STRING_ELT(method, 0));
    enum method_kind { MH, LD
    } meth_typ = ((!strcmp(meth_ch, "MH")) ? MH :
		  ((!strcmp(meth_ch, "LD")) ? LD :
		   -1));
    SEXP St; /**< the result */
    PROTECT(St = allocVector(REALSXP, n));
    switch(meth_typ) {
    case MH:
	retstable_MH(REAL(St), REAL(V0_), asReal(h), asReal(alpha), n);
	break;
    case LD:
	retstable_LD(REAL(St), REAL(V0_), asReal(h), asReal(alpha), n);
	break;
    default:
	error(_("retstable_c(): invalid '%s'"), "method");
    }

    UNPROTECT(2);
    return(St);
}