/*  This program computes the covariance matrix of the Conditionalized      */
/*  Random Midpoint Displacement (RMDmn) method for simulating              */
/*  fractional Gaussian noise.                                              */
/*  The C-packages Ranlib and Meschach are used, both available             */
/*  via Netlib (http://www.netlib.org).                                     */
/*  Part of the code is taken from a file containing the RMDmn algorithm,   */
/*  written by Petteri Mannersalo and Ilkka Norros of the Technical         */
/*  Research Centre of Finland (VTT).                                       */

/*  References:                                                             */
/*  1) I. Norros, P. Mannersalo, and J.L. Wang (1999),                      */
/*     Simulation of fractional Brownian motion with Conditionalized Random */
/*     Midpoint Displacement,                                               */
/*     Advances in Performance Analysis, Vol. 2, pp. 77--101.               */
/*  2) A.B. Dieker (2002),                                                  */
/*     Simulation of fractional Brownian motion,                            */
/*     Master's thesis, Vrije Universiteit Amsterdam.                       */

/*  Copyright Ton Dieker                                                    */
/*  Centre of Mathematics and Computer Science (CWI) Amsterdam              */
/*  April 2002                                                              */

/*  ton@cwi.nl                                                              */


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "matrix.h"
#include "matrix2.h"

#define Power(x, y)	(pow((double)(x), (double)(y)))

MAT *A,*D;
VEC *c, *d;
VEC *cD;
double *cDc;


double covar(long i1, long j1, long i2, long j2, double h) {
  /* covariance of the increments of the fractional brownian motion */
  /* multiplied by factor 2^(2*i1*h). Non-intersecting intervals    */
  /* [jx 2^(-ix), (jx+1) 2^(-ix)].                                  */ 
  double cov;

  cov =      Power(( j2   *Power(2,i1-i2)-(j1-1)), 2*h);
  cov =  cov+Power(((j2-1)*Power(2,i1-i2)- j1   ), 2*h);
  cov =  cov-Power(((j2-1)*Power(2,i1-i2)-(j1-1)), 2*h);
  cov = (cov-Power(( j2   *Power(2,i1-i2)- j1   ), 2*h))/2.0;

  return(cov);
}


MAT *covmat(long left, long right, double h) {
  /* scaled covariance matrix for fractional Brownian motion. */
  /* number of intervals counted is left+right.               */  
  MAT *A;
  long dim,i,j;
  double foo;

  dim = left+right;
  A = m_get(dim,dim);

  for(i=0; i<left; i++) {
    A->me[i][i] = 1;
    for(j=i+1; j<left; j++) {
	foo = covar(0,-j,0,-i,h); 
	A->me[left-1-i][left-j-1] = foo;
    }
    for(j=left;j<dim;j++) {
      A->me[j][left-1-i] = covar(0,-i,-1,j-left+1,h); 
    }
  }

  for(i=left; i<dim; i++) {
    A->me[i][i] = Power(2,2*h);
    for(j=i+1; j<dim; j++) {
      A->me[j][i] = covar(-1,i-left+1,-1,j-left+1,h)*Power(2,2*h); 
    }
  }

  for(i=1; i<dim; i++) {
    for(j=0; j<i; j++) {
      A->me[j][i] = A->me[i][j];
    }
  }
  return(A);
}


VEC *covvec(long left, long right, double h) {
  /* scaled covariance vector cov(x0,xi), where i varies from */
  /* x0-left to x0+right, (xi not x0).                        */
  VEC *x;
  long dim,i,j;

  dim = left+right;
  x = v_get(dim);
  
  for(i=0;i<left;i++) {
    x->ve[left-i-1] = covar(0,-i,0,1,h);
  }
  if(right > 0) {
    x->ve[left] = Power(2,2*h-1);
    for(i=left+1; i<dim; i++) {
      x->ve[i] = covar(0,1,-1 ,i-left+1,h);
    }
  }
  return(x);
}


void init(long left, long right, double h) {
  long i,j;
  /* compute the covariance vectors cov(x0,xi) and  */
  /* inverses of the covariance matrices cov(xi,xj) */
  for(i=0; i<=left; i++) {
    for(j=0; j<=right; j++) {
      if(i>0 || j>0) {
	A=covmat(i,j,h);
	D=m_get(A->m,A->n);
	c=covvec(i,j,h);
	d=v_get(c->dim);
	D=CHfactor(covmat(i,j,h));
	CHsolve(D,covvec(i,j,h),d);
	cD[j+i*(right+1)]=*d;
	cDc[j+i*(right+1)]=in_prod(d,c);
      }     
    }
  }
}


void end(long left, long right) {
  /* time to clean up a bit... */
  long i,j;
  m_free(A);
  m_free(D);
  v_free(c);

  /* the following code caused memory problems on my computer. */
  /* you may try if it works for you...                        */
  /* for(i=0; i<(left+1)*(right+1); i++) {                     */
  /*   v_free(&cD[i]);                                         */
  /* }                                                         */
  /* v_free(cD);                                               */
  /* free(cDc);                                                */

  /* this is how I partially solved these problems. Please     */
  /* notify me if you know why the above code didn't work.     */
  for(i=0; i<=left; i++) {
    for(j=0; j<=right; j++) {
	d = &cD[j+i*(right+1)];
        if(i>0 || j>0) free((char *) d->ve);
    }
  }
  v_free(cD);
  free(cDc);
}


MAT *fbmcov(long left, long right, long N, double h) {
  long rind,lind;
  long generator,indj;
  long i,j,m,k,l;
  double coe,foo,aux;
  MAT *y1,*y2;
  
  y1 = m_get(1,1);
  y1->me[0][0] = 1;

  for(i=1; i<=N; i++) {
    y2 = m_get(pow(2,i),pow(2,i));
    for (j=0; j<pow(2,i-1); j++) {
      for (k=0; k<pow(2,i-1); k++) {
	y2->me[2*j][2*k] = y2->me[2*j+1][2*k+1] = y1->me[j][k];
	y2->me[2*j+1][2*k] = y2->me[2*j][2*k+1] = y1->me[j][k];
      }
    }
    coe=Power(2.0,-2*i*h);
    indj=-1; 
    
    /* left hand side */ 
    for(j=1; j<=min(left,pow(2,i)); j=j+2) {
      lind = j-1;
      rind = min(right,pow(2,i-1)-(j-1)/2);    
      d = &cD[rind+lind*(right+1)];
      foo = coe*(1-cDc[rind+lind*(right+1)]);
      y2->me[j-1][j-1] += 1./foo;
      for(k=0; k<lind; k++) {
	for (l=0; l<lind; l++) {
	  y2->me[k][l] += d->ve[k]*d->ve[l]/foo;
	}
	for (l=0; l<rind; l++) {
	  y2->me[k][lind+2*l] += d->ve[k]*d->ve[lind+l]/foo;
	  y2->me[k][lind+2*l+1] += d->ve[k]*d->ve[lind+l]/foo;
	}
      }
      for(k=0; k<rind; k++) {
	for (l=0; l<lind; l++){
	  y2->me[lind+2*k][l] += d->ve[lind+k]*d->ve[l]/foo;
	  y2->me[lind+2*k+1][l] += d->ve[lind+k]*d->ve[l]/foo;
	}
	for (l=0; l<rind; l++) {
	  aux = d->ve[lind+k]*d->ve[lind+l]/foo;
	  y2->me[lind+2*k][lind+2*l] += aux;
	  y2->me[lind+2*k+1][lind+2*l+1] += aux;
	  y2->me[lind+2*k][lind+2*l+1] += aux;
	  y2->me[lind+2*k+1][lind+2*l] += aux;
	}
      }

      for (k=0;k<lind;k++) {
	y2->me[lind][k] -= d->ve[k]/foo;
	y2->me[k][lind] -= d->ve[k]/foo;
      }
      for (k=0;k<rind;k++) {
	y2->me[lind][lind+2*k] -= d->ve[lind+k]/foo;
	y2->me[lind][lind+2*k+1] -= d->ve[lind+k]/foo;
	y2->me[lind+2*k][lind] -= d->ve[lind+k]/foo;
	y2->me[lind+2*k+1][lind] -= d->ve[lind+k]/foo;
      }
      
      indj=j;
    }
    
    /* middle */ 
    if(indj < pow(2,i)-2*right) {
      d = &cD[(left+1)*(right+1)-1];
      foo = coe*(1-cDc[(left+1)*(right+1)-1]);   
      for(j=indj+2; j<pow(2,i)+2-2*right; j=j+2) {
	y2->me[j-1][j-1] += 1./foo;    
	for(k=0; k<left; k++) {
	  for(l=0; l<left; l++) {
	    y2->me[j-1-left+k][j-1-left+l] += d->ve[l]*d->ve[k]/foo;
	  }
	  for(l=0; l<right; l++) {
	    y2->me[j-1-left+k][j-1+2*l] += d->ve[left+l]*d->ve[k]/foo;
	    y2->me[j-1-left+k][j+2*l] += d->ve[left+l]*d->ve[k]/foo;
	  }
	}
	for (k=0; k<right; k++) {
	  for(l=0; l<left; l++) {
	    y2->me[j-1+2*k][j-1-left+l] += d->ve[l]*d->ve[left+k]/foo;
	    y2->me[j+2*k][j-1-left+l] += d->ve[l]*d->ve[left+k]/foo;
	  }
	  for(l=0; l<right; l++) {
	    aux = d->ve[left+l]*d->ve[left+k]/foo;
	    y2->me[j-1+2*k][j-1+2*l] += aux;
	    y2->me[j-1+2*k][j+2*l] += aux;
	    y2->me[j+2*k][j-1+2*l] += aux;
	    y2->me[j+2*k][j+2*l] += aux;
	  }
	}
	
	for (k=0; k<left; k++) {
	  y2->me[j-1][j-1-left+k] -= d->ve[k]/foo;
	  y2->me[j-1-left+k][j-1] -= d->ve[k]/foo;
	}
	for (k=0; k<right; k++) {
	  y2->me[j-1][j-1+2*k] -= d->ve[left+k]/foo;
	  y2->me[j-1][j+2*k] -= d->ve[left+k]/foo;
	  y2->me[j-1+2*k][j-1] -= d->ve[left+k]/foo;
	  y2->me[j+2*k][j-1] -= d->ve[left+k]/foo;
	}
	
	indj=j;
      }
    }
    
    /* right hand side */ 
    for(j=indj+2; j<pow(2,i); j=j+2) {
      lind = min(left,j-1);
      rind = min(rind,(pow(2,i-1)-(j-1)/2));
      
      d = &cD[rind+lind*(right+1)];
      foo = coe*(1-cDc[rind+lind*(right+1)]);
      y2->me[j-1][j-1] += 1./foo;    
      for(k=0; k<lind; k++) {
	for(l=0; l<lind; l++) {
	  y2->me[j-1-lind+k][j-1-lind+l] += d->ve[l]*d->ve[k]/foo;
	}
	for(l=0; l<rind; l++) {
	  y2->me[j-1-lind+k][j-1+2*l] += d->ve[lind+l]*d->ve[k]/foo;
	  y2->me[j-1-lind+k][j+2*l] += d->ve[lind+l]*d->ve[k]/foo;
	}
      }
      for (k=0; k<rind; k++) {
	for(l=0; l<lind; l++) {
	  y2->me[j-1+2*k][j-1-lind+l] += d->ve[l]*d->ve[lind+k]/foo;
	  y2->me[j+2*k][j-1-lind+l] += d->ve[l]*d->ve[lind+k]/foo;
	}
	for(l=0; l<rind; l++) {
	  aux = d->ve[lind+l]*d->ve[lind+k]/foo;
	  y2->me[j-1+2*k][j-1+2*l] += aux;
	  y2->me[j-1+2*k][j+2*l] += aux;
	  y2->me[j+2*k][j-1+2*l] += aux;
	  y2->me[j+2*k][j+2*l] += aux;
	}
      }
      
      for (k=0; k<lind; k++) {
	y2->me[j-1][j-1-lind+k] -= d->ve[k]/foo;
	y2->me[j-1-lind+k][j-1] -= d->ve[k]/foo;
      }
      for (k=0; k<rind; k++) {
	y2->me[j-1][j-1+2*k] -= d->ve[lind+k]/foo;
	y2->me[j-1][j+2*k] -= d->ve[lind+k]/foo;
	y2->me[j-1+2*k][j-1] -= d->ve[lind+k]/foo;
	y2->me[j+2*k][j-1] -= d->ve[lind+k]/foo;
      }
    }
    m_free(y1);
    y1 = y2;
  }
  m_inverse(y1,y1);
  return(y1);
}  


void rmdcov(long *N, long *m, long *n, double *h, double *out) {
  /* function that computes the covariance matrix of the RMD method      */
  /* by Norros et al.                                                    */
  /* Input:  *N      the resulting matrix is 2^(*N)x2^(*N)               */
  /*         *m      the truncation parameter on the left                */
  /*         *n      the truncation parameter on the right               */
  /*         *h      the Hurst parameter                                 */
  /* Output: *out    the resulting sample is stored in this array        */
  long i,j,aux;
  long left = *m;
  long right = *n;
  long S = *N;
  double H= *h;
  MAT *cov;

  aux = Power(2,S);
  cD = (VEC *) calloc((left+1)*(right+1),sizeof(VEC));
  cDc = (double *) malloc((left+1)*(right+1)*sizeof(double));
  init(left,right,H);

  cov = fbmcov(left, right, S, H);
  sm_mlt(Power(2,2*S*H),cov,cov);

  /* put the result in array `out' for exporting to R */
  for (i=0; i<aux; i++)
    for (j=0; j<aux; j++)
      out[i+j*aux] = cov->me[i][j];
  end(left,right);

  m_free(cov);
}