#include <R.h>
#include <Rmath.h>
#include <string.h>

#define MIN(x, y) (((x) < (y)) ? (x) : (y))

void lossdistrib(double *p, int *np, double *w, double *S, int *N, int* defaultflag, double *q) {
  /* recursive algorithm with first order correction for computing
     the loss distribution.
     p vector of default probabilities
     np length of p
     S vector of severities (should be same length as p)
     N number of ticks in the grid
     defaultflat if true compute the default distribution
     q the loss distribution */

  int i, j, d1, d2, M;
  double lu, d, p1, p2, sum;
  double *qtemp;

  lu = 1./(*N-1);
  qtemp = Calloc(*N, double);
  q[0] = 1;
  M = 1;
  for(i=0; i<(*np); 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;
    memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
    for(j=0; j < MIN(M, *N); j++){
      q[j] = (1-p[i]) * q[j];
    }
    for(j=0; j < MIN(M, *N-d2); j++){
      q[d1+j] += p1 * qtemp[j];
      q[d2+j] += p2 * qtemp[j];
    };
    M+=d2;
  }

  /* correction for weight loss */
  if(M>*N){
    sum = 0;
    for(j=0; j<MIN(M, *N); j++){
      sum += q[j];
    }
    q[*N-1] += 1-sum;
  }
  Free(qtemp);
}

void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
			   int *T, int *defaultflag, double *q) {
  /* recursive algorithm with first order correction for computing
     the loss distribution.
     p vector of default probabilities
     np length of p
     S vector of severities (should be same length as p)
     N number of ticks in the grid
     T where to truncate the distribution
     defaultflat if true computes the default distribution
     q the loss distribution */

  int i, j, d1, d2, M;
  double lu, d, p1, p2, sum;
  double *q1, *q2;

  lu = 1./(*N-1);
  q1 = Calloc( *T, double);
  q2 = Calloc( *T, double);
  q[0] = 1;
  M = 1;
  for(i=0; i<(*np); 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;
    for(j=0; j < MIN(M, *T); j++){
      q1[j] = p1 * q[j];
      q2[j] = p2 * q[j];
      q[j] = (1-p[i]) * q[j];
    }
    for(j=0; j < MIN(M, *T-d1); j++){
      q[d1+j] += q1[j];
    };
    for(j=0; j < MIN(M, *T-d2); j++){
      q[d2+j] += q2[j];
    };
    M += d2;
  }
  Free(q1);
  Free(q2);
}

void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
  /* recursive algorithm with first order correction
     computes jointly the loss and recovery distribution
     p vector of default probabilities
     np length of p
     w vector of issuer weights (length np)
     S vector of severities (should be same length as p)
     N number of ticks in the grid
     defaultflag if true computes the default distribution
     q the joint probability distribution */

  int i, j, k, m, n;
  int Mx, My;
  double lu, x, y, sum;
  double alpha1, alpha2, beta1, beta2;
  double w1, w2, w3, w4;
  double *qtemp;

  lu = 1./(*N-1);
  qtemp = Calloc( (*N) * (*N), double);
  q[0] = 1;
  Mx=1;
  My=1;
  for(k=0; k<(*np); 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;
    w2 = alpha1 * beta2;
    w3 = alpha2 * beta2;
    w4 = alpha2 * beta1;

    for(n=0; n<MIN(My, *N); n++){
      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
    }
    for(n=0; n<MIN(My, *N); n++){
      for(m=0; m<MIN(Mx, *N); m++){
	q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
      }
    }
    for(n=0; n < MIN(My, *N-j-1); n++){
      for(m=0; m < MIN(Mx, *N-i-1); m++){
	q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
	q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
	q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
	q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
      }
    }
    Mx += i+1;
    My += j+1;
  }
  /* correction for weight loss */
  if(Mx>*N || My>*N){
    sum = 0;
    for(m=0; m < MIN(Mx, *N); m++){
      for(n=0; n < MIN(My, *N); n++){
  	sum += q[m+n*(*N)];
      }
    }
    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
  }
  Free(qtemp);
}

void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
  /* 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
     n length of p
     S vector of severities (should be same length as p)
     w vector of weights
     N number of ticks in the grid
     q the loss distribution */

  int i, j, d1l, d1u, d2l, d2u;
  int M;
  double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
  double *qtemp;

  lu = 1./(*N - 1);
  qtemp = Calloc( (*N), double);
  q[0] = 1;
  M=1;
  for(i=0; i<(*n); 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;
    memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
    for(j = 0; j< MIN(M, *N); j++){
      q[j] = (1-dp[i]-pp[i]) * q[j];
    }
    for(j=0; j < 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){
    sum = 0;
    for(j=0; j<MIN(M, *N); j++){
      sum += q[j];
    }
    q[*N-1] += 1-sum;
  }
  Free(qtemp);
}

void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
			      double *S, int *N, int *defaultflag, double *q) {
  int i, j1, j2, k, m, n;
  double lu, x, y1, y2, sum;
  double alpha1, alpha2, beta1, beta2;
  double dpw1, dpw2, dpw3, dpw4;
  double ppw1, ppw2, ppw3;
  double temp;
  double *qtemp;
  int Mx, My;

  lu = 1./(*N-1);
  qtemp = Calloc((*N) * (*N), double);
  q[0] = 1;
  Mx=1;
  My=1;
  for(k=0; k<(*ndp); 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(n=0; n<MIN(My, *N); n++){
      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
    }
    for(n=0; n<MIN(My, *N); n++){
      for(m=0; m<MIN(Mx, *N); m++){
	q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
      }
    }
    if(*defaultflag){
      ppw1 = alpha1 * alpha1 * pp[k];
      ppw2 = alpha1 * alpha2 * pp[k];
      ppw3 = alpha2 * alpha2 * pp[k];
      for(n=0; n < MIN(My, *N-j2-1); n++){
	for(m=0; m < MIN(Mx, *N-i-1); m++){
	  q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
	  q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
	  q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
	  q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];

	  q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
	  temp = ppw2 * qtemp[m+(*N)*n];
	  q[i+m+(*N)*(j2+1+n)] += temp;
	  q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
	  q[i+1+m+(*N)*(j2+n)] += temp;
	}
      }
    }else{
      for(n=0; n < MIN(My, *N-j2-1); n++){
	for(m=0; m < MIN(Mx, *N-i-1); m++){
	  q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
	  q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
	  q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
	  q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
	  q[m+(*N)*(j2+n)] +=  pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
	  q[m+(*N)*(j2+1+n)] +=  pp[k]*(y2-j2) * qtemp[m+(*N)*n];
	}
      }
    }
    Mx += i + 1;
    My += j2 + 1;
  }
  /* correction for weight loss */
  if(Mx>*N || My>*N){
    sum = 0;
    for(m=0; m < MIN(Mx, *N); m++){
      for(n=0; n < MIN(My, *N); n++){
	sum += q[m+n*(*N)];
      }
    }
    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
  }
  Free(qtemp);
}

double shockprob(double p, double rho, double Z, int give_log){
  return( pnorm( (qnorm(p, 0, 1, 1, 0) -sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
}

double shockseverity(double S, double Z, double rho, double p){
  return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );

}

void addandmultiply(double *X, double alpha, double *Y, int n) {
  int i;
  for(i = 0; i<n; i++){
    Y[i] += alpha*X[i];
  }
}

void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
		double *recov, double *Z, double *w, int *n, double *rho, int *N,
		int *defaultflag, double *L, double *R) {
  /*
    computes the loss and recovery distribution over time with a flat gaussiancorrelation
    inputs:
      Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers and dim2 number of time steps
      recov: vector of recoveries (length dim1)
      issuerweights: vector of issuer weights (length dim2)
      Z: vector of factor values (length n)
      w: vector of factor weights (length n)
      rho: correlation beta
      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 t, i, j;
  double g;
  double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;

  gshocked = Calloc((*dim1), double);
  Rshocked = Calloc((*dim1), double);
  Sshocked = Calloc((*dim1), double);
  Lw = malloc((*N) * sizeof(double));
  Rw = malloc((*N) * sizeof(double));

  for(t=0; t < (*dim2); t++) {
    for(i=0; i < *n; i++){
      for(j=0; j < (*dim1); j++){
	g = 1 - SurvProb[j + (*dim1) * t];
	gshocked[j] = shockprob(g, *rho, Z[i], 0);
	Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
	Rshocked[j] = 1 - Sshocked[j];
      }
      /* reset Lw and Rw to 0 */
      memset(Lw, 0, *N * sizeof(double));
      memset(Rw, 0, *N * sizeof(double));
      lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag, Lw);
      lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, defaultflag, Rw);
      addandmultiply(Lw, w[i], L + t * (*N), *N);
      addandmultiply(Rw, w[i], R + t * (*N), *N);
    }
  }
  Free(gshocked);
  Free(Rshocked);
  Free(Sshocked);
  free(Lw);
  free(Rw);
}