/*
 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"

/**
 * Sample V ~ F with F(n) = 1-1/(n*B(n,1-alpha)), n in IN, with 
 * Laplace-Stieltjes transform 1-(1-exp(-t))^alpha via the algorithm of 
 * Hofert (2010).
 * Note: The caller of this function must use GetRNGstate() and PutRNGstate().
 * @param alpha parameter theta0/theta1 in (0,1]
 * @param iAlpha 1-alpha 
 * @param gamma_1_a Gamma(1-alpha)
 * @return a random variate from F
 * @author Marius Hofert, Martin Maechler
*/
double rFJoe(double alpha, double iAlpha /**< := 1 - alpha */, 
	     double gamma_1_a /**< == Gamma(1 - alpha) == Gamma(iALpha) */){ 
    double U, I_al = 1./alpha;

    /**< FIXME(MM): (for alpha not too close to 1): re-express using 1-U */
    U = unif_rand();
    if(U <= alpha)
	return 1.;
    else { /**< alpha < U < 1 */
	double Ginv = pow((1-U)*gamma_1_a, -I_al);
	double fGinv = floor(Ginv);
	if(1-U < 1./(fGinv*beta(fGinv, iAlpha)))
	    return ceil(Ginv);
	else return fGinv;
    }
}

/** 
 * Vectorize rFJoe. Generate a vector of variates 
 * V ~ F with F(n) = 1-1/(n*B(n,1-alpha)), n in IN, with Laplace-Stieltjes 
 * transform 1-(1-exp(-t))^alpha.
 * @param V vector of random variates from F (result)
 * @param n length of the vector V
 * @param alpha parameter theta0/theta1 in (0,1]
 * @param iAlpha 1-alpha 
 * @return none
 * @author Marius Hofert, Martin Maechler
*/
void rFJoe_vec(double V[], const int n,
	       const double alpha, const double iAlpha /**< := 1 - alpha */){
    if(n >= 1) {
	double G1_a = gammafn(iAlpha);
	GetRNGstate();

	for(int i=0; i < n; i++)
	    V[i] = rFJoe(alpha, iAlpha, G1_a);

	PutRNGstate();
    }
    return;
}

/**
 * Generate a vector of variates 
 * V ~ F with F(n) = 1-1/(n*B(n,1-alpha)), n in IN, with Laplace-Stieltjes 
 * transform 1-(1-exp(-t))^alpha.
 * Note: Should be fast as it is used as a building block in different places.
 * @param n sample size
 * @param alpha parameter theta0/theta1 in (0,1]
 * @return vector of random variates V
 * @author Martin Maechler
*/
SEXP rFJoe_c(SEXP n, SEXP alpha)
{
    int nn = asInteger(n);
    double alp = asReal(alpha);
    SEXP res = PROTECT(allocVector(REALSXP, nn));

    rFJoe_vec(REAL(res), nn, alp, 1. - alp);

    UNPROTECT(1);
    return(res);
}