#include <omp.h>
#include "lossdistrib.hpp"
#include <Rcpp.h>
#include <cstring>

#define USE_BLAS

using namespace Rcpp;

// [[Rcpp::plugins(cpp11, openmp)]]

// [[Rcpp::export]]
List GHquadCpp(int n) {
    // Setup for eigenvalue computations
    char JOBZ   = 'V'; // Compute eigenvalues & vectors
    int INFO;
    // Initialize array for workspace
    double * WORK   = new double[2*n-2];

    // Initialize array for eigenvectors
    double * V      = new double[n*n];

    std::vector<double> w(n);
    std::vector<double> Z(n);
    for(int i = 0; i<n-1; i++) {
      w[i] = sqrt((i+1.)/2);
    }

    // Run eigen decomposition
    dstev_(&JOBZ, &n, Z.data(), w.data(), V, &n, WORK, &INFO);

    for (int i=0; i<n; i++) {
      w[i] = V[i*n] * V[i*n];
      Z[i] *= sqrt(2);
    }

    // Deallocate temporary arrays
    delete[] WORK;
    delete[] V;
    return List::create(Named("Z") = Z,
                        Named("w") = w);
}

void lossdistrib(const std::vector<double>& p, const std::vector<double>& w,
                 const double S[], const int N, const int T, double q[],
                 bool defaultflag = false) {
    /* recursive algorithm with first order correction for computing
       the loss distribution.
       p: vector of default probabilities
       w: issuer weights
       S: vector of severities (should be same length as p)
       N: number of ticks in the grid
       defaultflag: if true compute the default distribution */

    int d1, d2;
    double d, p1, p2;
    double* qtemp = new double[T];
    q[0] = 1;
    int bound;
    double pbar;
    int one = 1;
    openblas_set_num_threads(1);
    double lu = 1./(N-1);

    int M = 1;
    for(size_t i = 0; i < p.size(); i++){
        d = defaultflag? w[i]/lu : S[i] * w[i]/ lu;
        d1 = floor(d);
        d2 = ceil(d);
        p1 = p[i] * (d2-d);
        p2 = p[i] - p1;
        std::memcpy(qtemp, q, std::min(M, T) * sizeof(double));
        pbar = 1-p[i];
        bound = std::min(M, T);
#ifdef USE_BLAS
        dscal_(&bound, &pbar, q, &one);
#else
        for(int j=0; j < bound; j++) {
            q[j] *= pbar * q[j];
        }
#endif
        bound = std::min(M, T-d2);
#ifdef USE_BLAS
        daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
        daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
#else
        for(int j=0; j < bound; j++) {
            q[d1+j] += p1 * qtemp[j];
            q[d2+j] += p2 * qtemp[j];
        }
#endif
        M += d2;
    }
    /* correction for weight loss */
    if(M > N && T==N){
        double sum = 0;
        for(int j=0; j < std::min(M, N); j++){
            sum += q[j];
        }
        q[N-1] += 1-sum;
    }
    delete[] qtemp;
}

// [[Rcpp::export]]
inline NumericVector  lossdistrib(const NumericVector& p, const NumericVector& w,
                                  const NumericVector& S, const int N,
                                  bool defaultflag = false) {
    double* q = new double[N];
    lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w),
                S.begin(), N, N, q, defaultflag);
    NumericVector r(N, q);
    delete[] q;
    return r;
}

inline std::vector<double> lossdistrib(const NumericVector& p, const NumericVector& w,
                                       const NumericVector& S, const int N, const int T,
                                       bool defaultflag = false) {
    std::vector<double> q(T);
    lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w),
                S.begin(), N, T, q.data(), defaultflag);
    return q;
}

// [[Rcpp::export]]
NumericMatrix lossdistrib_Z(const std::vector<double>& p, const std::vector<double>& w,
                            const NumericMatrix S, int N, const std::vector<double>& rho,
                            const std::vector<double>& Z, bool defaultflag=false){

    double* q = new double[N*Z.size()];
    int p_size = p.size();
    #pragma omp parallel for
    for(size_t i = 0; i < Z.size(); i++) {
        std::vector<double> pshocked(p.size());
        for(size_t j = 0; j < p_size; j++){
            pshocked[j] = shockprob(p[j], rho[j], Z[i], 0);
        }
        lossdistrib(pshocked, w, S(_,i).begin(), N, N, q+N*i, defaultflag);
    }
    NumericMatrix qmat(N, Z.size(), q);
    delete[] q;
    return qmat;
}

static inline double* posK(int T, double K, double lu){
    double* val = new double[T];
    for(int i = 0; i < T; i++){
        val[i] = K-lu*i;
    }
    return val;
}

// [[Rcpp::export]]
double exp_trunc(const NumericVector& p, const NumericVector& w, const NumericVector& S, int N, double K) {
    double lu = 1./(N+1);
    int T = (int) floor(K * N)+1;
    int one = 1;
    double* q  = new double[T];
    lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w),
                S.begin(), N, T, q);
    double* val = posK(T, K, lu);
    int r = ddot_(&T, val, &one, q, &one);
    delete[] q;
    delete[] val;
    return r;
}

void lossdistrib_joint(const std::vector<double>& p, const std::vector<double>& w,
                       const double S[], int N, matrix& q, bool defaultflag=false) {
    /* recursive algorithm with first order correction
       computes jointly the loss and recovery distribution
       p vector of default probabilities
       w vector of issuer weights (same length as p)
       S vector of severities (should be same length as p)
       N number of ticks in the grid
       defaultflag if true computes the default distribution
       returns the joint probability distribution */

    int i, j;
    double x, y;
    double alpha1, alpha2, beta1, beta2;
    double w1, w2, w3, w4;
    matrix qtemp(N, N);
    double lu = 1./(N-1);
    q(0,0) = 1;
    int Mx=1, My=1;
    int one = 1;
    for(size_t k = 0; k < p.size(); k++) {
        x = defaultflag? w[k] /lu : S[k] * w[k] / lu;
        y = (1-S[k]) * w[k] / lu;
        i = floor(x);
        j = floor(y);
        alpha1 = i + 1 - x;
        alpha2 = 1 - alpha1;
        beta1 = j + 1 - y;
        beta2 = 1 - beta1;
        w1 = alpha1 * beta1 * p[k];
        w2 = alpha1 * beta2 * p[k];
        w3 = alpha2 * beta2 * p[k];
        w4 = alpha2 * beta1 * p[k];

        for(int n=0; n<std::min(My, N); n++) {
            std::memcpy(qtemp(n), q(n), std::min(Mx, N) * sizeof(double));
        }
        int bound = std::min(Mx, N);
        double pbar = 1-p[k];
        for(int n=0; n < std::min(My, N); n++) {
#ifdef USE_BLAS
            dscal_(&bound, &pbar, q(n), &one);
#else
            for(int m=0; m < bound; m++){
                q(m, n) *= pbar;
            }
#endif
        }
        bound = std::min(Mx, N-i-1);
        for(int n=0; n < std::min(My, N-j-1); n++) {
#ifdef USE_BLAS
            daxpy_(&bound, &w1, qtemp(n), &one, &q(i,j+n), &one);
            daxpy_(&bound, &w2, qtemp(n), &one, &q(i,j+1+n), &one);
            daxpy_(&bound, &w3, qtemp(n), &one, &q(i+1,j+1+n), &one);
            daxpy_(&bound, &w4, qtemp(n), &one, &q(i+1,j+n), &one);
#else
            double temp;
            for(int m = 0; m < bound; m++) {
                temp = qtemp(m,n);
                q(i+m,j+n) += w1 * temp;
                q(i+m,j+1+n) += w2 * temp;
                q(i+1+m,j+1+n] += w3 * temp;
                q(i+1+m,j+n) += w4 * temp;
            }
#endif
        }
        Mx += i+1;
        My += j+1;
    }
    /* correction for weight loss */
    if(Mx>N || My>N){
        double sum = 0;
        for(int n=0; n < std::min(My, N); n++){
            for(int m=0; n < std::min(Mx, N); m++){
                sum += q(m,n);
            }
        }
        q[std::min(N, Mx) * std::min(N, My)-1] += 1 - sum;
    }
}

// [[Rcpp::export]]
NumericMatrix lossdistrib_joint(const std::vector<double>& p, const std::vector<double>& w,
                                const NumericVector S, int N, bool defaultflag=false) {
    matrix q(N,N);
    lossdistrib_joint(p, w, S.begin(), N, q, defaultflag);
    NumericMatrix M(N, N, q.data());
    return M;
}

// [[Rcpp::export]]
NumericVector recovdist(const NumericVector& dp, const NumericVector& pp,
                        const NumericVector& w, const NumericVector& S, int N) {
    /* recursive algorithm with first order correction for computing
       the recovery distribution in case of prepayment.
       dp vector of default probabilities
       pp vector of prepay probabilities
       S vector of severities (should be same length as p)
       w vector of weights
       N number of ticks in the grid
       returns the loss distribution */

    int d1l, d1u, d2l, d2u;
    double d1, d2, dp1, dp2, pp1, pp2;

    double lu = 1./(N - 1);
    NumericVector qtemp(N);
    NumericVector q(N);
    q[0] = 1;
    int M = 1;
    for(size_t i=0; i<dp.size(); i++) {
        d1 = w[i] * (1-S[i]) /lu;
        d2 = w[i]/lu;
        d1l = floor(d1);
        d1u = d1l + 1;
        d2l = floor(d2);
        d2u = d2l + 1;
        dp1 = dp[i] * (d1u - d1);
        dp2 = dp[i] - dp1;
        pp1 = pp[i] * (d2u - d2);
        pp2 = pp[i] - pp1;
        std::copy_n(q.begin(), std::min(M, N), qtemp.begin());
        for(int j = 0; j< std::min(M, N); j++) {
            q[j] = (1-dp[i]-pp[i]) * q[j];
        }
        for(int j=0; j < std::min(M, N-d2u); j++) {
            q[d1l+j] += dp1 * qtemp[j];
            q[d1u+j] += dp2 * qtemp[j];
            q[d2l+j] += pp1 * qtemp[j];
            q[d2u+j] += pp2 * qtemp[j];
        };
        M += d2u;
    }
    /* correction for weight loss */
    if(M > N) {
        double sum = 0;
        for(int j=0; j < std::min(M, N); j++) {
            sum += q[j];
        }
        q[N-1] += 1-sum;
    }
    return q;
}

NumericVector lossdistrib_prepay_joint(const NumericVector& dp, const NumericVector& pp,
                                       const NumericVector& w, const NumericVector& S,
                                       int N, bool defaultflag=false) {
    int i, j1, j2;
    double x, y1, y2;
    double alpha1, alpha2, beta1, beta2;
    double dpw1, dpw2, dpw3, dpw4;
    double ppw1, ppw2, ppw3;

    double lu = 1./(N-1);
    double* qtemp = new double[N*N];
    double* q = new double[N*N];
    q[0] = 1;

    int Mx = 1;
    int My = 1;
    for(size_t k = 0; k < dp.size(); k++) {
        y1 = (1-S[k]) * w[k]/lu;
        y2 = w[k]/lu;
        x = defaultflag? y2: y2-y1;
        i = floor(x);
        j1 = floor(y1);
        j2 = floor(y2);
        alpha1 = i + 1 - x;
        alpha2 = 1 - alpha1;
        beta1 = j1 + 1 - y1;
        beta2 = 1 - beta1;
        dpw1 = alpha1 * beta1 * dp[k];
        dpw2 = alpha1 * beta2 * dp[k];
        dpw3 = alpha2 * beta2 * dp[k];
        dpw4 = alpha2 * beta1 * dp[k];

        /* by default distribution, we mean names fractions of names that disappeared
           either because of default or prepayment */
        for(int n=0; n<std::min(My, N); n++){
            std::memcpy(qtemp+n*N, q+n*N, std::min(Mx, N) * sizeof(double));
        }
        int bound = std::min(Mx, N);
        double pbar = 1-dp[k]-pp[k];
        int one;
        for(int n = 0; n < std::min(My, N); n++) {
#ifdef USE_BLAS
            dscal_(&bound, &pbar, q+N*n, &one);
#else
            for(int m=0; m < bound; m++) {
                q[m+N*n] *= pbar;
            }
#endif
        }
        double temp;
        double* begin;
        if(defaultflag) {
            ppw1 = alpha1 * alpha1 * pp[k];
            ppw2 = alpha1 * alpha2 * pp[k];
            ppw3 = alpha2 * alpha2 * pp[k];
            bound = std::min(Mx, N-i-1);
            for(int n=0; n < std::min(My, N-j2-1); n++) {
                begin = qtemp+N*n;
#ifdef USE_BLAS

                daxpy_(&bound, &dpw1, begin, &one, q+i+N*(j1+n), &one);
                daxpy_(&bound, &dpw2, begin, &one, q+i+N*(j1+1+n), &one);
                daxpy_(&bound, &dpw3, begin, &one, q+i+1+N*(j1+1+n), &one);
                daxpy_(&bound, &dpw4, begin, &one, q+i+1+N*(j1+n), &one);

                daxpy_(&bound, &ppw1, begin, &one, q+i+N*(j2+n), &one);
                daxpy_(&bound, &ppw2, begin, &one, q+i+N*(j2+1+n), &one);
                daxpy_(&bound, &ppw3, begin, &one, q+i+1+N*(j2+1+n), &one);
                daxpy_(&bound, &ppw2, begin, &one, q+i+1+N*(j2+n), &one);
#else
                for(int m=0; m < bound; m++) {
                    temp = *(begin + m);
                    q[i+m+N*(j1+n)] += dpw1 * temp;
                    q[i+1+m+N*(j1+n)] += dpw4 * temp;
                    q[i+m+N*(j1+1+n)] += dpw2 * temp;
                    q[i+1+m+N*(j1+1+n)] += dpw3 * temp;

                    q[i+m+N*(j2+n)] += ppw1 * temp;
                    q[i+1+m+N*(j2+n)] += ppw2 * temp;
                    q[i+m+N*(j2+1+n)] += ppw2 * temp;
                    q[i+1+m+N*(j2+1+n)] += ppw3 * temp;

                }
#endif
            }
        } else {
            ppw1 = pp[k] * (j2+1-y2);
            ppw2 = pp[k] * (y2-j2);
            bound = std::min(Mx, N-i-1);
            for(int n = 0; n < std::min(My, N-j2-1); n++) {
                begin = qtemp + N*n;
#ifdef USE_BLAS

                daxpy_(&bound, &dpw1, begin, &one, q+i+N*(j1+n), &one);
                daxpy_(&bound, &dpw2, begin, &one, q+i+N*(j1+1+n), &one);
                daxpy_(&bound, &dpw3, begin, &one, q+i+1+N*(j1+1+n), &one);
                daxpy_(&bound, &dpw4, begin, &one, q+i+1+N*(j1+n), &one);
                daxpy_(&bound, &ppw1, begin, &one, q+N*(j2+n), &one);
                daxpy_(&bound, &ppw2, begin, &one, q+N*(j2+1+n), &one);
#else
                for(int m = 0; m < bound; m++) {
                    temp = *(begin+m);

                    q[i+m+N*(j1+n)] += dpw1 * temp;
                    q[i+1+m+N*(j1+n)] += dpw4 * temp;
                    q[i+m+N*(j1+1+n)] += dpw2 * temp;
                    q[i+1+m+N*(j1+1+n)] += dpw3 * temp;

                    q[m+N*(j2+n)] +=  ppw1 * temp;
                    q[m+N*(j2+1+n)] +=  ppw2 * temp;
                }
#endif
            }
        }
        Mx += i + 1;
        My += j2 + 1;
    }
    /* correction for weight loss */
    if(Mx > N || My > N){
        double sum = 0;
        for(int m = 0; m < std::min(Mx, N); m++) {
            for(int n=0; n < std::min(My, N); n++) {
                sum += q[m+n*N];
            }
        }
        q[std::min(N, Mx)*std::min(My,N)-1] += 1 - sum;
    }
    free(qtemp);
}

double shockprob(double p, double rho, double Z, bool give_log=false) {
    if(rho==1){
        return((double)(Z<=R::qnorm(p, 0, 1, 1, 0)));
    }else{
        return( R::pnorm( (R::qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
    }
}

static inline double dqnorm(double x){
    return 1/R::dnorm(R::qnorm(x, 0, 1, 1, 0), 0, 1, 0);
}

double dshockprob(double p, double rho, double Z){
    return( R::dnorm((R::qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) );
}

NumericMatrix shockprobvec2(double p, double rho, NumericVector Z){
    /* return a two column matrix with shockprob in the first column
       and dshockprob in the second column*/
    NumericMatrix q(Z.size(), 2);
    #pragma omp parallel for
    for(size_t i = 0; i < Z.size(); i++){
        q[i] = shockprob(p, rho, Z[i]);
        q[i + Z.size()] = dshockprob(p, rho, Z[i]);
    }
    return q;
}

double shockseverity(double S, double Z, double rho, double p){
    if(p==0){
        return 0;
    }else{
        return( exp(shockprob(S * p, rho, Z, true) - shockprob(p, rho, Z, true)) );
    }
}

double quantile(const NumericVector& Z, const NumericVector& w, double p0){
    double cumw = 0;
    size_t i;
    for(i=0; i < Z.size(); i++) {
        cumw += w[i];
        if(cumw > p0) {
            break;
        }
    }
    return( Z[i] );
}

// [[Rcpp::export]]
double fitprob(const NumericVector& Z, const NumericVector& w, double rho, double p0){
    if(p0 == 0){
        return 0.;
    }else if(rho == 1){
        return R::pnorm(quantile(Z, w, p0), 0, 1, 1, 0);
    }else{
        int one = 1;
        double eps = 1e-12;
        NumericMatrix q = shockprobvec2(p0, rho, Z);
        int nZ = Z.size();
        double dp = (ddot_(&nZ, q(_,0).begin(), &one, w.begin(), &one) - p0) /
            ddot_(&nZ, q(_,1).begin(), &one, w.begin(), &one);
        double p = p0;
        double phi = 1;
        while(fabs(dp) > eps){
            while( (p - phi * dp) < 0 || (p - phi * dp) > 1){
                phi *= 0.8;
            }
            p -= phi * dp;
            q = shockprobvec2(p, rho, Z);
            dp = (ddot_(&nZ, q(_,0).begin(), &one, w.begin(), &one) - p0) /
                ddot_(&nZ, q(_,1).begin(), &one, w.begin(), &one);
        }
        return p;
    }
}

// [[Rcpp::export]]
std::vector<double> stochasticrecov(double R, double Rtilde,
                                    const NumericVector& Z, const NumericVector& w,
                                    double rho, double porig, double pmod){
    if(porig==0){
        std::vector<double> q(Z.size(), R);
        return q;
    }else{
        double ptemp = (1 - R) / (1 - Rtilde) * porig;
        double ptilde = fitprob(Z, w, rho, ptemp);
        std::vector<double> q(Z.size());
        #pragma omp parallel for
        for(size_t i = 0; i < Z.size(); i++){
            q[i] = fabs(1 - (1 - Rtilde) * exp( shockprob(ptilde, rho, Z[i], true) -
                                                shockprob(pmod, rho, Z[i], true)));
        }
        return q;
    }
}

// void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
//                                 double *S, int *N, int *defaultflag, double *rho,
//                                 double *Z, double *wZ, int *nZ, double *q) {
//     int i, j;
//     double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
//     double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
//     int N2 = (*N) * (*N);
//     double* qmat = malloc(sizeof(double) * N2 * (*nZ));

//     double alpha = 1;
//     double beta = 0;
//     int one = 1;

//     #pragma omp parallel for private(j)
//     for(i = 0; i < *nZ; i++){
//         for(j = 0; j < *ndp; j++){
//             dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
//             ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
//         }
//         lossdistrib_prepay_joint_blas(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
//                                       w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
//     }

//     dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);

//     free(dpshocked);
//     free(ppshocked);
//     free(qmat);
// }

NumericMatrix lossdistrib_joint_Z(const std::vector<double>& dp, const std::vector<double>& w,
                                  const NumericMatrix S, int N, bool defaultflag, const std::vector<double>& rho,
                                  const std::vector<double>& Z, const std::vector<double>& wZ) {

    int m = Z.size(), ndp = dp.size();
    NumericMatrix q(N, N);
    int N2 = N * N;
    double* qmat = new double[N2 * m];
    double alpha = 1;
    double beta = 0;
    int one = 1;

    #pragma omp parallel for
    for(size_t i = 0; i < m; i++) {
        std::vector<double> dpshocked(ndp);
        matrix qZ(N, N, qmat+i);
        for(size_t j =0; j < ndp; j++) {
            dpshocked[j] = shockprob(dp[j], rho[j], Z[i], 0);
        }
        lossdistrib_joint(dpshocked, w, S(_,i).begin(), N, qZ, defaultflag);
    }
    dgemv_((char*)"n", &N2, &m, &alpha, qmat, &N2, wZ.data(), &one, &beta, q.begin(), &one);
    delete[] qmat;
    return q;
}

// [[Rcpp::export]]
List BCloss_recov_dist(NumericMatrix defaultprob,
                       const std::vector<double>& issuerweights, NumericVector recov,
                       NumericVector Z, NumericVector w,
                       NumericVector rho, int N, bool defaultflag = false) {
    /*
      computes the loss and recovery distribution over time with a flat gaussian
      correlation
      inputs:
      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
      and dim2 number of time steps
      issuerweights: vector of issuer weights (length dim1)
      recov: vector of recoveries (length dim1)
      Z: vector of factor values (length n)
      w: vector of factor weights (length n)
      rho: correlation beta vector (length dim1)
      N: number of ticks in the grid
      defaultflag: if true, computes the default distribution
      outputs:
      L: matrix of size (N, dim2)
      R: matrix of size (N, dim2)
    */
    int one = 1;
    double alpha = 1;
    double beta = 0;

    NumericMatrix L(N, defaultprob.ncol());
    NumericMatrix R(N, defaultprob.ncol());
    double* Lw = new double[N*Z.size()];
    double* Rw = new double[N*Z.size()];
    for(int t = 0; t < defaultprob.ncol(); t++) {
        #pragma omp parallel
        {
            std::vector<double> gshocked(defaultprob.nrow());
            double* Rshocked = new double[defaultprob.nrow()];
            double* Sshocked = new double[defaultprob.nrow()];
            double g;
            #pragma omp for
            for(int i=0; i < Z.size(); i++) {
                for(int j = 0; j < defaultprob.nrow(); j++) {
                    g = defaultprob(j, t);;
                    gshocked[j] = shockprob(g, rho[j], Z[i], 0);
                    Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
                    Rshocked[j] = 1 - Sshocked[j];
                }
                lossdistrib(gshocked, issuerweights, Sshocked, N, N,
                            Lw+i*N, defaultflag);
                lossdistrib(gshocked, issuerweights, Rshocked, N, N,
                            Rw+i*N, defaultflag);
            }
            delete[] Rshocked;
            delete[] Sshocked;
        }
        int n = Z.size();
        dgemv_((char*)"n", &N, &n, &alpha, Lw, &N, w.begin(), &one, &beta, L(_,t).begin(), &one);
        dgemv_((char*)"n", &N, &n, &alpha, Rw, &N, w.begin(), &one, &beta, R(_,t).begin(), &one);
    }
    delete[] Lw;
    delete[] Rw;
    return List::create(Named("L")=L, Named("R")=R);
}

// void BCloss_recov_trunc(double *defaultprob, int *dim1, int *dim2,
//                         double *issuerweights, double *recov, double *Z, double *w,
//                         int *n, double *rho, int *N, double * K, int *defaultflag,
//                         double *ELt, double *ERt) {
//     /*
//       computes EL_t = E[(K-L_t)^+] and ER_t = E[(K-(1-R_t))^+]
//        the the put options on loss and recovery over time
//        with a flat gaussian correlation.
//       inputs:
//       defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
//       and dim2 number of time steps
//       issuerweights: vector of issuer weights (length dim1)
//       recov: vector of recoveries (length dim1)
//       Z: vector of factor values (length n)
//       w: vector of factor weights (length n)
//       rho: correlation beta vector (length dim1)
//       N: number of ticks in the grid
//       K: the strike
//       defaultflag: if true, computes the default distribution
//       outputs:
//       ELt: vector of length dim2
//       ERt: vector of length dim2
//     */
//     int t, i, j;
//     double g, Ktilde;
//     int one = 1;
//     int T = (int) floor((*K) * (*N))+1;
//     int Ttilde;
//     double lu = 1./(*N+1);
//     double* valL = malloc( T * sizeof(double));
//     posK(T, *K, lu, valL);
//     double* EL = malloc( (*n) * sizeof(double));
//     double* ER = malloc( (*n) * sizeof(double));
//     double* Lw = malloc(T * (*n) * sizeof(double));
//     double alpha = 1;
//     double beta = 0;

//     for(t=0; t < (*dim2); t++) {
//         memset(Lw, 0, T * (*n) * sizeof(double));
//         #pragma omp parallel for private(j, g, Ktilde, Ttilde)
//         for(i=0; i < *n; i++){
//             double* Rw = NULL;
//             double* Rshocked = malloc((*dim1) * sizeof(double));
//             double* Sshocked = malloc((*dim1) * sizeof(double));
//             double* gshocked = malloc((*dim1) * sizeof(double));
//             double* gshockedbar = malloc((*dim1) * sizeof(double));
//             double* valR = NULL;

//             for(j=0; j < (*dim1); j++){
//                 g = defaultprob[j + (*dim1) * t];
//                 gshocked[j] = shockprob(g, rho[j], Z[i], 0);
//                 Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
//                 gshockedbar[j] = 1 - gshocked[j];
//                 Rshocked[j] = 1 - Sshocked[j];
//             }

//             lossdistrib_truncated(gshocked, dim1, issuerweights, Sshocked,
//                                   N, &T, defaultflag, Lw + i * T);
//             ER[i] = 0;
//             Ktilde = *K - ddot_(dim1, issuerweights, &one, Sshocked, &one);
//             if(Ktilde > 0){
//                 Ttilde = (int) floor(Ktilde * (*N))+1;
//                 Rw = calloc(Ttilde, sizeof(double));
//                 lossdistrib_truncated(gshockedbar, dim1, issuerweights, Rshocked,
//                                       N, &Ttilde, defaultflag, Rw);
//                 valR = malloc(Ttilde * sizeof(double));
//                 posK(Ttilde, Ktilde, lu, valR);
//                 ER[i] = ddot_(&Ttilde, Rw, &one, valR, &one);
//             }
//             if(Rw != NULL){
//                 free(Rw);
//             }
//             if(valR != NULL){
//                 free(valR);
//             }
//             free(Rshocked);
//             free(Sshocked);
//             free(gshocked);
//             free(gshockedbar);
//         }
//         dgemv_("t", &T, n, &alpha, Lw, &T, valL, &one, &beta, EL, &one);
//         ELt[t] = ddot_(n, EL, &one, w, &one);
//         ERt[t] = ddot_(n, ER, &one, w, &one);
//     }
//     free(Lw);
//     free(EL);
//     free(ER);
//     free(valL);
// }

// void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
//                  double *issuerweights, double *recov, double *Z, double *w,
//                  int *n, double *rho, int *N, int *defaultflag,
//                  double *L) {
//     /*
//       computes the loss distribution over time with a flat gaussian
//       correlation
//       inputs:
//       defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
//       and dim2 number of time steps
//       issuerweights: vector of issuer weights (length dim1)
//       recov: vector of recoveries (length dim1)
//       Z: vector of factor values (length n)
//       w: vector of factor weights (length n)
//       rho: correlation beta vector (length dim1)
//       N: number of ticks in the grid
//       defaultflag: if true, computes the default distribution
//       outputs:
//       L: matrix of size (N, dim2)
//     */
//     int t, i, j;
//     double g;
//     double *gshocked, *Sshocked, *Lw;
//     int one = 1;
//     double alpha = 1;
//     double beta = 0;

//     gshocked = malloc((*dim1) * (*n) * sizeof(double));
//     Sshocked = malloc((*dim1) * (*n) * sizeof(double));
//     Lw = malloc((*N) * (*n) * sizeof(double));

//     for(t=0; t < (*dim2); t++) {
//         memset(Lw, 0, (*N) * (*n) * sizeof(double));
//         #pragma omp parallel for private(j, g)
//         for(i=0; i < *n; i++){
//             for(j=0; j < (*dim1); j++){
//                 g = defaultprob[j + (*dim1) * t];
//                 gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
//                 Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
//             }
//             lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
//                              Lw + i * (*N));
//         }
//         dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
//     }
//     free(gshocked);
//     free(Sshocked);
//     free(Lw);
// }
//