libs/dtgf/fahlman_ulrych_code_C99.h

// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/dtgf/fahlman_ulrych_code_C99.h,v 1.1 2013/04/28 07:46:58 tplarson Exp $

#include "ctypes_def.h"

#if TYPE == FLOAT 
   #define GAPSTRUCTURE sgapstructure
   #define BURG smulti_burg
   #define FAHLMAN_ULRYCH sfahlman_ulrych
   #define FILL_GAPS sfill_gaps
   #define TOEPLITZ_MULT stoeplitz_mult
   #define PRECONDITION sprecondition
   #define PCG spcg
   #define LEVINSON slevinson
   #define EPSILON FLT_EPSILON
#elif TYPE == DOUBLE 
   #define GAPSTRUCTURE dgapstructure
   #define BURG dmulti_burg
   #define FAHLMAN_ULRYCH dfahlman_ulrych
   #define FILL_GAPS dfill_gaps
   #define TOEPLITZ_MULT dtoeplitz_mult
   #define PRECONDITION dprecondition
   #define PCG dpcg
   #define LEVINSON dlevinson
   #define EPSILON DBL_EPSILON
#elif TYPE == COMPLEXFLOAT
   #define GAPSTRUCTURE cgapstructure
   #define BURG cmulti_burg
   #define FAHLMAN_ULRYCH cfahlman_ulrych
   #define FILL_GAPS cfill_gaps
   #define TOEPLITZ_MULT ctoeplitz_mult
   #define PRECONDITION cprecondition
   #define PCG cpcg
   #define LEVINSON clevinson
   #define EPSILON FLT_EPSILON
#elif TYPE == COMPLEXDOUBLE
   #define GAPSTRUCTURE zgapstructure
   #define BURG zmulti_burg
   #define FAHLMAN_ULRYCH zfahlman_ulrych
   #define FILL_GAPS zfill_gaps
   #define TOEPLITZ_MULT ztoeplitz_mult
   #define PRECONDITION zprecondition
   #define PCG zpcg
   #define LEVINSON zlevinson
   #define EPSILON DBL_EPSILON
#endif

#define PCG_MAXIT n_gap
#define PCG_TOL 1e-6

static int GAPSTRUCTURE(int n, CTYPE *data, int *isgood, 
                        gapped_timeseries *ts, int **data_length, 
                        CTYPE ***good_data, int minpercentage);
static void FILL_GAPS(int n, CTYPE *data, int *isgood, int order, CTYPE *ar_coeff);
static void TOEPLITZ_MULT(int n, CTYPE *x, CTYPE *y, void **data);
static void PRECONDITION(int n, CTYPE *x, CTYPE *y, void **data);

//static void identity(int n, CTYPE *x, CTYPE *y, void **data);


int FAHLMAN_ULRYCH(int n, CTYPE *data_in, int *isgood_in, int minpercentage, 
                   int maxorder, int iterations, int padends, 
                   int *order, CTYPE *ar_coeff_in)
{
  int i,j, first, last, iter, effective_order;
  RTYPE *err    ;
  gapped_timeseries ts,tsm;
  CTYPE **good_data;
  int *isgood;
  CTYPE *data, nold;
  int *data_length;
  int *model_gaps, oldorder;
  int model_order;
  CTYPE *ar_coeff;

  
  if (ar_coeff_in)
    ar_coeff = ar_coeff_in;
  else
    ar_coeff = malloc(maxorder*sizeof(CTYPE));

// the following check will be skipped in production executables since they will be built with NDEBUG defined.
// it doesn't matter because the check is already done in jtsfiddle.c which calls this function.
  assert(minpercentage>0 && minpercentage<=100);
  // Execute the Fahlman-Ulrych algorithm: 
  // FOR i=1 TO iterations DO
  //   Step 1. Compute AR model using known sample values.
  //   Step 2. For the unknown samples, compute predicted values that
  //           minimize the forward and backward prediction errors in 
  //           the least squares sense.
  // END FOR

  if (padends)
  {
    nold = n;
    n = nold + 2*maxorder;
    data = calloc(n,sizeof(CTYPE));
    isgood = calloc(n,sizeof(int));
    memcpy(&data[maxorder],data_in,nold*sizeof(CTYPE));
    memcpy(&isgood[maxorder],isgood_in,nold*sizeof(int));
  }
  else
  {
    isgood = isgood_in;
    data = data_in;
  }

  oldorder = -1;
  model_order = GAPSTRUCTURE(n, data, isgood, &ts, &data_length, 
                             &good_data, minpercentage);
  tsm = ts;
  model_gaps = malloc(n*sizeof(int));
  memcpy(model_gaps, isgood, n*sizeof(int));
  for (iter=0; iter<iterations; iter++)
  {
    // Analyze gapstructure.
    if (iter>0)
      model_order = GAPSTRUCTURE(n, data, model_gaps, &tsm, &data_length, 
                            &good_data, minpercentage);


    model_order = min( model_order,maxorder);

    // Sort gaps by length and determine the largest AR model order
    // such that "minpercentage" percent of the data will be included
    // in the estimation of the AR coefficients.
    //    model_order =  maxorder(&ts,data_length,minpercentage);    
    err  = (RTYPE *) malloc((model_order+1)*sizeof(RTYPE));

#ifdef DEBUGOUT
    {
      printf("Determining AR(%d) model...",model_order);  
      fflush(stdout);
    }
#endif

    effective_order = BURG(tsm.m_data, data_length, good_data, model_order, 
                           ar_coeff, err);

#ifdef DEBUGOUT
    {
      printf("done\n");
      printf("Specified order = %d, Effective order = %d.\n", model_order, effective_order);
      printf("Filling gaps...\n");  
      fflush(stdout);
    }
#endif

    model_order = effective_order;
    if (ts.first_is_good)
    {
      i = 1; 
      first = max(0,ts.data_int[0].last-(model_order)+1); 
      last=-1;
    }
    else    
    {
      i = 0;
      first = 0; 
      last=-1;
    }
    memcpy(model_gaps, isgood, n*sizeof(int));
    while ( i<ts.m_data && last<n-1)
    {
      // move the right-hand end of the interval until we find a
      // data section longer than the model order, i.e. which decouples
      // the unknowns in the gaps on either side.
      while (i<ts.m_data && 
             (ts.data_int[i].last-ts.data_int[i].first+1)<(model_order) &&
             (i==0 || (ts.data_int[i].first-ts.data_int[i-1].last-1)<model_order)
             )
      { i++; }


      if (i!=0 && (ts.data_int[i].first-ts.data_int[i-1].last)>=model_order)
        last = ts.data_int[i-1].last+(model_order)-1;
      else if ( i>=ts.m_data-1 )
        last = n-1;
      else
        last = ts.data_int[i].first+(model_order)-1;
      
#ifdef DEBUGOUT
      printf("Filling missing data in [ %d : %d ].\n",first,last);
#endif

      FILL_GAPS(last-first+1, &data[first], &isgood[first], 
                (model_order), ar_coeff);
    
      first = ts.data_int[i].last-(model_order)+1;
      i++;
    }
    // Mark filled gaps shorter than effective order as good.
    for(i=0; i<ts.m_gap; i++)
    {
      if ( ts.gap_int[i].last-ts.gap_int[i].first+1 <= model_order)
        for (j = ts.gap_int[i].first; j<=ts.gap_int[i].last; j++)
          model_gaps[j] = 1;
    }
#ifdef DEBUGOUT
    printf("done\n");
#endif
    if (iter>0)
    {
      free(tsm.gap_int);  free(tsm.data_int);
    }
    free(data_length);  free(good_data);
    free(err); 

#ifdef DEBUGFILE
char name[20];
FILE *tsfile, *gapfile;
sprintf(name,"tsiter%d.out",iter);
tsfile=fopen(name,"w");
sprintf(name,"gapsiter%d.out",iter);
gapfile=fopen(name,"w");
int ii;
for (ii=0;ii<n;ii++)
{
fprintf(tsfile,"%e %e\n",REAL(data[ii]),IMAG(data[ii]));
fprintf(gapfile,"%d\n",model_gaps[ii]);
}
fclose(tsfile);
fclose(gapfile);
#endif

    if (oldorder==model_order || tsm.m_gap==0)
      break;

    oldorder = model_order;
  }
  free(ts.gap_int);  free(ts.data_int);
  memcpy(isgood,model_gaps,n*sizeof(int));
  for (i=0;i<n;i++)
    if (model_gaps[i]==0)
      data[i] = 0;

  free(model_gaps);  


  if (padends)
  {
    memcpy(data_in, &data[maxorder],nold*sizeof(CTYPE));
    memcpy(isgood_in,&isgood[maxorder],nold*sizeof(int));
    free(data);
    free(isgood);
  }

  if (ar_coeff_in==NULL)
    free(ar_coeff);

  if (order)
    *order = model_order;

  return 0;
}


static void FILL_GAPS(int n, CTYPE *data, int *isgood, int order, 
                      CTYPE *ar_coeff)
{
  int iter, idx, i, shift, *gap_idx, *gap_idx_inv, *block;
  int n_gap, n_data, n_blocks, maxblock;
  CTYPE *rhs, *rhs2, *data_gap, *gamma2, *scratch;
  RTYPE rnorm, tol;
  void *aarg[3], *marg[3];
  CTYPE *gamma;


  //  for (i=0; i<10; i++)
  //  printf("data[%d] = %f%+fi.\n",i,creal(data[i]),cimag(data[i]));

  // Compute gamma containing the first row of T^H*T.
  gamma = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
  rhs = (CTYPE *)malloc(2*(n-order)*sizeof(CTYPE));
  for (shift = 0; shift<=order; shift++)
  {
    gamma[shift] = 0;
    for (i=0; (i+shift)<=order; i++) 
      gamma[shift] += 2*CONJ(ar_coeff[i+shift])*ar_coeff[i];
    //     printf("gamma[%d] = %f%+fi\n",shift,creal(gamma[shift]),cimag(gamma[shift]));
  }
  // for (i=0; i<=order; i++)
  //  printf("ar_coeff[%d] = %f%+fi.\n",i,creal(ar_coeff[i]),cimag(ar_coeff[i]));
  

  // Find gaps.
  gap_idx = (int*)calloc(n,sizeof(int));
  gap_idx_inv = (int*)calloc(n,sizeof(int));
  block = (int*)calloc(n,sizeof(int));
  n_data = 0;
  n_gap = 0;
  n_blocks = 0;
  maxblock = 0;
  for (i=0; i<n; i++)
  {
    if (isgood[i]) 
    {
      n_data++;
      if (i>0 && !isgood[i-1])
        n_blocks++;
    }
    else
    {
      // printf("gap at %4d\n",i);
      gap_idx[i] = n_gap;
      gap_idx_inv[n_gap] = i;
      n_gap++;
      block[n_blocks]++;
      if (block[n_blocks] > maxblock)
        maxblock = block[n_blocks];
    }
  }

  /* Set up right-hand side for block Toeplitz system. */
  for (shift=0; (shift+order)<n; shift++)
  {
    rhs[shift] = 0;
    rhs[shift+(n-order)] = 0;
    for (i=0; i<=order; i++)
    {
      if (isgood[i+shift])
      {
        rhs[shift] -= ar_coeff[order-i]*data[i+shift];
        rhs[shift+(n-order)] -= CONJ(ar_coeff[i])*data[i+shift];        
      }
    }
  }
  
  rhs2 = (CTYPE *)calloc(n_gap,sizeof(CTYPE));   
  for (shift=0; (shift+order)<n; shift++)
  {
    for (i=0; i<=order; i++)
    {
      if ( !isgood[i+shift] )
      {
        idx = gap_idx[i+shift];
        //      assert(idx>=0 && idx < n_gap);
        rhs2[idx] += CONJ(ar_coeff[order-i])*rhs[shift];
        rhs2[idx] += ar_coeff[i]*rhs[shift+(n-order)];
        
      }
    }
  }
  
  //  printf("n_gap = %d\n",n_gap);
  //for (i=0; i<n_gap; i++)
  //  printf("rhs2[%d] = %f%+fi.\n",i,creal(rhs2[i]),cimag(rhs2[i]));
  

  data_gap = (CTYPE*)calloc(n_gap,sizeof(CTYPE));
  if (n_gap == 1)
    data_gap[0] = rhs2[0] / gamma[0];
  else
  {
    // Prepare PCG inputs defining the matrix.
    aarg[0] = (void *) &order;
    aarg[1] = (void *) gap_idx_inv;
    aarg[2] = (void *) gamma;

    // Prepare PCG inputs defining the preconditioner.
    gamma2 = (CTYPE *) calloc(maxblock,sizeof(CTYPE));
    scratch = (CTYPE *) calloc(maxblock,sizeof(CTYPE));
    for (i=0; i<=min(order,maxblock-1); i++)
      gamma2[i] = gamma[i];
    marg[0] = (void *) block;
    marg[1] = (void *) gamma2;
    marg[2] = (void *) scratch;
  
    // Solve normal equations using PCG with a block Toeplitz preconditioner.
    memset(data_gap, 0, n_gap*sizeof(CTYPE));
    tol = sqrt(order)*PCG_TOL;
    iter = PCG(n_gap, PCG_MAXIT+2, tol, TOEPLITZ_MULT, 
               PRECONDITION, rhs2, data_gap, &rnorm, aarg, marg);
    if (rnorm>tol)
      fprintf(stderr, "Warning: PCG did not converge. "
              "After %d iterations (maxit=%d) the relative"
              " residual was %e (tol=%e)\n",
              iter, PCG_MAXIT+2, rnorm,tol);
 
    /*    printf("n_gap=%3d, maxblock=%3d, iter=%2d, rnorm=%e\n", 
          n_gap, maxblock, iter,rnorm); */
    //    for (i=0; i<n_gap; i++)
    //printf("x[%d] = %f%+fi\n",i,creal(data_gap[i]),cimag(data_gap[i]));
    free(scratch);  free(gamma2); 
  }
  // Fill gap the original time series with extrapolated data.
  for (i=0; i<n_gap; i++)
    data[gap_idx_inv[i]] = data_gap[i];

  free(block); free(gap_idx); free(gap_idx_inv);  
  free(data_gap); free(gamma); free(rhs); free(rhs2);
}


static int GAPSTRUCTURE(int n, CTYPE *data, int *isgood, gapped_timeseries *ts,
                  int **data_length, CTYPE ***good_data, int minpercentage)
{
  int i,j,first;

 // Count number of good and bad points and
  // number of good and bad intervals.
  ts->n_data = 0;  ts->m_data = 0;
  ts->n_gap = 0;   ts->m_gap = 0;
  for (i=0; i<n-1; i++)
  {
//    assert(isgood[i]==1 || isgood[i]==0);
    if (isgood[i]==1)
    {
      ts->n_data++;
      if (i==0)
      {
        ts->first_is_good = 1;
        ts->m_data++;
      }
      if (isgood[i+1]==0)
        ts->m_gap++;
    }
    else if (isgood[i]==0)
    {
      ts->n_gap++;
      if (i==0)
      {
        ts->first_is_good = 0;
        ts->m_gap++;
      }
      if (isgood[i+1]==1)
        ts->m_data++;
    }
    else
    {
      fprintf(stderr, "ERROR: problem in routine gapstructure: value of window function must be either 0 or 1\n");
      exit(1);
    }
  }
  if (isgood[n-1]==1)
      ts->n_data++;
  else    
      ts->n_gap++;


#ifdef DEBUGOUT
    printf("n_data = %d, m_data = %d\nn_gap = %d, m_gap = %d\n",
           ts->n_data,ts->m_data,ts->n_gap,ts->m_gap);
#endif
/*
  assert((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) ||
         (isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) ||
         (isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data));
*/

  if  (!((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) ||
         (isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) ||
         (isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data)))
  {
    fprintf(stderr, "ERROR: problem in routine gapstructure\n");
    exit(1);
  }


  // Allocate space for data pointers and AR coefficients. 
  *good_data  = (CTYPE **) malloc(ts->m_data*sizeof(CTYPE *));
  *data_length = (int *) malloc(ts->m_data*sizeof(int));
  ts->gap_int = (interval *) malloc(ts->m_gap*sizeof(interval));
  ts->data_int = (interval *) malloc(ts->m_data*sizeof(interval));
  

  // Scan again and set up data structures describing gap structure. 
  j = 0;
  first=0;
  if (isgood[0]==1)
  {
    first=0;
    ts->data_int[0].first = 0;
    (*good_data)[0] = &data[0];
    j++;
  }
  for (i=0; i<n-1; i++)
  {
    // first of an interval
    if (isgood[i]==0 && isgood[i+1]==1)
    {
      first = i+1;
      ts->data_int[j].first = first;
      (*good_data)[j] = &(data[i+1]);
      //      printf("first[%d] = %d\n",j,i+1);
      j++;
    }
    // end of an interval
    if (isgood[i]==1 && isgood[i+1]==0)
    {
      (*data_length)[j-1] = i-first+1;
      ts->data_int[j-1].last = i;
    }
  }
  if (isgood[n-1]==1)
  {
    (*data_length)[j-1] = n-1-first+1;
    ts->data_int[j-1].last = n-1;
  }
  if (ts->first_is_good)
  {
    for(j=0; j<ts->m_data;j++)
    {
      if (j<ts->m_gap)
          ts->gap_int[j].first = ts->data_int[j].last+1;
      if (j>0)
          ts->gap_int[j-1].last = ts->data_int[j].first-1;
    }
    if (ts->m_data==ts->m_gap)
      ts->gap_int[ts->m_gap-1].last = n-1;
  }
  else
  {
    ts->gap_int[0].first = 0;
    for(j=0; j<ts->m_data;j++)
    {
      if (j<ts->m_gap)
          ts->gap_int[j].last = ts->data_int[j].first-1;
      if (j<ts->m_gap-1)
          ts->gap_int[j+1].first = ts->data_int[j].last+1;
    }
    if (ts->m_gap==ts->m_data+1)
      ts->gap_int[ts->m_gap-1].last = n-1;
  }
  return maxorder(ts, *data_length, minpercentage);
}


// Matrix-vector multiply with symmetric Toeplitz matrix.
// y = T(r) x

static void TOEPLITZ_MULT(int n, CTYPE x[n], CTYPE y[n], void **data)
{
  int i,j,order,idx,ii;
  int *gap_idx_inv;
  CTYPE *r;

  order = *((int *) data[0]);
  gap_idx_inv = (int *) data[1];
  r = (CTYPE *) data[2];

  for (i=0; i<n; i++)
    y[i] = r[0]*x[i];

  for (i=0; i<n; i++)
  {
    ii = gap_idx_inv[i];
    for (j=i+1; j<n; j++)
    {
      idx = ii-gap_idx_inv[j];
      if (idx <= order && idx >= -order)
      {
        if (idx>=0)
        {
          y[i] += CONJ(r[idx])*x[j];
          y[j] += r[idx]*x[i];
        }
        else
        {
          y[i] += r[-idx]*x[j];
          y[j] += CONJ(r[-idx])*x[i];
        }
      }
    }
  }
}

// Solve M y = x;
static void PRECONDITION(int n, CTYPE *b, CTYPE *x, void **data)
{
  int i,first,*block;
  CTYPE *gamma, *scratch;
  block = (int *) data[0];
  gamma = (CTYPE *) data[1];
  scratch = (CTYPE *) data[2];
  

  for(i=0, first = 0; first<n; first+=block[i], i++)
  {
    memcpy(scratch,&b[first],block[i]*sizeof(CTYPE));
    LEVINSON(block[i],gamma,scratch,&x[first]);
  }
}


#undef GAPSTRUCTURE
#undef FAHLMAN_ULRYCH 
#undef BURG
#undef FILL_GAPS
#undef TOEPLITZ_MULT
#undef PRECONDITION
#undef PCG
#undef LEVINSON
#undef EPSILON
#include "ctypes_undef.h"
Revision:	1.2
Committed:	Fri Feb 13 19:40:03 2015 UTC (8 years, 7 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_8-7, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, globalhs_version_8, globalhs_version_9, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-10, Ver_9-1, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.1:	+18 -1 lines
Log Message:	add printing of debugging files when DEBUGFILE is defined
#	Content
1	// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/dtgf/fahlman_ulrych_code_C99.h,v 1.1 2013/04/28 07:46:58 tplarson Exp $
2
3	#include "ctypes_def.h"
4
5	#if TYPE == FLOAT
6	#define GAPSTRUCTURE sgapstructure
7	#define BURG smulti_burg
8	#define FAHLMAN_ULRYCH sfahlman_ulrych
9	#define FILL_GAPS sfill_gaps
10	#define TOEPLITZ_MULT stoeplitz_mult
11	#define PRECONDITION sprecondition
12	#define PCG spcg
13	#define LEVINSON slevinson
14	#define EPSILON FLT_EPSILON
15	#elif TYPE == DOUBLE
16	#define GAPSTRUCTURE dgapstructure
17	#define BURG dmulti_burg
18	#define FAHLMAN_ULRYCH dfahlman_ulrych
19	#define FILL_GAPS dfill_gaps
20	#define TOEPLITZ_MULT dtoeplitz_mult
21	#define PRECONDITION dprecondition
22	#define PCG dpcg
23	#define LEVINSON dlevinson
24	#define EPSILON DBL_EPSILON
25	#elif TYPE == COMPLEXFLOAT
26	#define GAPSTRUCTURE cgapstructure
27	#define BURG cmulti_burg
28	#define FAHLMAN_ULRYCH cfahlman_ulrych
29	#define FILL_GAPS cfill_gaps
30	#define TOEPLITZ_MULT ctoeplitz_mult
31	#define PRECONDITION cprecondition
32	#define PCG cpcg
33	#define LEVINSON clevinson
34	#define EPSILON FLT_EPSILON
35	#elif TYPE == COMPLEXDOUBLE
36	#define GAPSTRUCTURE zgapstructure
37	#define BURG zmulti_burg
38	#define FAHLMAN_ULRYCH zfahlman_ulrych
39	#define FILL_GAPS zfill_gaps
40	#define TOEPLITZ_MULT ztoeplitz_mult
41	#define PRECONDITION zprecondition
42	#define PCG zpcg
43	#define LEVINSON zlevinson
44	#define EPSILON DBL_EPSILON
45	#endif
46
47	#define PCG_MAXIT n_gap
48	#define PCG_TOL 1e-6
49
50	static int GAPSTRUCTURE(int n, CTYPE data, int isgood,
51	gapped_timeseries ts, int *data_length,
52	CTYPE ***good_data, int minpercentage);
53	static void FILL_GAPS(int n, CTYPE data, int isgood, int order, CTYPE *ar_coeff);
54	static void TOEPLITZ_MULT(int n, CTYPE x, CTYPE y, void **data);
55	static void PRECONDITION(int n, CTYPE x, CTYPE y, void **data);
56
57	//static void identity(int n, CTYPE x, CTYPE y, void **data);
58
59
60
61	int FAHLMAN_ULRYCH(int n, CTYPE data_in, int isgood_in, int minpercentage,
62	int maxorder, int iterations, int padends,
63	int order, CTYPE ar_coeff_in)
64	{
65	int i,j, first, last, iter, effective_order;
66	RTYPE *err ;
67	gapped_timeseries ts,tsm;
68	CTYPE **good_data;
69	int *isgood;
70	CTYPE *data, nold;
71	int *data_length;
72	int *model_gaps, oldorder;
73	int model_order;
74	CTYPE *ar_coeff;
75
76
77	if (ar_coeff_in)
78	ar_coeff = ar_coeff_in;
79	else
80	ar_coeff = malloc(maxorder*sizeof(CTYPE));
81
82	// the following check will be skipped in production executables since they will be built with NDEBUG defined.
83	// it doesn't matter because the check is already done in jtsfiddle.c which calls this function.
84	assert(minpercentage>0 && minpercentage<=100);
85	// Execute the Fahlman-Ulrych algorithm:
86	// FOR i=1 TO iterations DO
87	// Step 1. Compute AR model using known sample values.
88	// Step 2. For the unknown samples, compute predicted values that
89	// minimize the forward and backward prediction errors in
90	// the least squares sense.
91	// END FOR
92
93	if (padends)
94	{
95	nold = n;
96	n = nold + 2*maxorder;
97	data = calloc(n,sizeof(CTYPE));
98	isgood = calloc(n,sizeof(int));
99	memcpy(&data[maxorder],data_in,nold*sizeof(CTYPE));
100	memcpy(&isgood[maxorder],isgood_in,nold*sizeof(int));
101	}
102	else
103	{
104	isgood = isgood_in;
105	data = data_in;
106	}
107
108	oldorder = -1;
109	model_order = GAPSTRUCTURE(n, data, isgood, &ts, &data_length,
110	&good_data, minpercentage);
111	tsm = ts;
112	model_gaps = malloc(n*sizeof(int));
113	memcpy(model_gaps, isgood, n*sizeof(int));
114	for (iter=0; iter<iterations; iter++)
115	{
116	// Analyze gapstructure.
117	if (iter>0)
118	model_order = GAPSTRUCTURE(n, data, model_gaps, &tsm, &data_length,
119	&good_data, minpercentage);
120
121
122	model_order = min( model_order,maxorder);
123
124	// Sort gaps by length and determine the largest AR model order
125	// such that "minpercentage" percent of the data will be included
126	// in the estimation of the AR coefficients.
127	// model_order = maxorder(&ts,data_length,minpercentage);
128	err = (RTYPE ) malloc((model_order+1)sizeof(RTYPE));
129
130	#ifdef DEBUGOUT
131	{
132	printf("Determining AR(%d) model...",model_order);
133	fflush(stdout);
134	}
135	#endif
136
137	effective_order = BURG(tsm.m_data, data_length, good_data, model_order,
138	ar_coeff, err);
139
140	#ifdef DEBUGOUT
141	{
142	printf("done\n");
143	printf("Specified order = %d, Effective order = %d.\n", model_order, effective_order);
144	printf("Filling gaps...\n");
145	fflush(stdout);
146	}
147	#endif
148
149	model_order = effective_order;
150	if (ts.first_is_good)
151	{
152	i = 1;
153	first = max(0,ts.data_int[0].last-(model_order)+1);
154	last=-1;
155	}
156	else
157	{
158	i = 0;
159	first = 0;
160	last=-1;
161	}
162	memcpy(model_gaps, isgood, n*sizeof(int));
163	while ( i<ts.m_data && last<n-1)
164	{
165	// move the right-hand end of the interval until we find a
166	// data section longer than the model order, i.e. which decouples
167	// the unknowns in the gaps on either side.
168	while (i<ts.m_data &&
169	(ts.data_int[i].last-ts.data_int[i].first+1)<(model_order) &&
170	(i==0 \|\| (ts.data_int[i].first-ts.data_int[i-1].last-1)<model_order)
171	)
172	{ i++; }
173
174
175	if (i!=0 && (ts.data_int[i].first-ts.data_int[i-1].last)>=model_order)
176	last = ts.data_int[i-1].last+(model_order)-1;
177	else if ( i>=ts.m_data-1 )
178	last = n-1;
179	else
180	last = ts.data_int[i].first+(model_order)-1;
181
182	#ifdef DEBUGOUT
183	printf("Filling missing data in [ %d : %d ].\n",first,last);
184	#endif
185
186	FILL_GAPS(last-first+1, &data[first], &isgood[first],
187	(model_order), ar_coeff);
188
189	first = ts.data_int[i].last-(model_order)+1;
190	i++;
191	}
192	// Mark filled gaps shorter than effective order as good.
193	for(i=0; i<ts.m_gap; i++)
194	{
195	if ( ts.gap_int[i].last-ts.gap_int[i].first+1 <= model_order)
196	for (j = ts.gap_int[i].first; j<=ts.gap_int[i].last; j++)
197	model_gaps[j] = 1;
198	}
199	#ifdef DEBUGOUT
200	printf("done\n");
201	#endif
202	if (iter>0)
203	{
204	free(tsm.gap_int); free(tsm.data_int);
205	}
206	free(data_length); free(good_data);
207	free(err);
208
209	#ifdef DEBUGFILE
210	char name[20];
211	FILE tsfile, gapfile;
212	sprintf(name,"tsiter%d.out",iter);
213	tsfile=fopen(name,"w");
214	sprintf(name,"gapsiter%d.out",iter);
215	gapfile=fopen(name,"w");
216	int ii;
217	for (ii=0;ii<n;ii++)
218	{
219	fprintf(tsfile,"%e %e\n",REAL(data[ii]),IMAG(data[ii]));
220	fprintf(gapfile,"%d\n",model_gaps[ii]);
221	}
222	fclose(tsfile);
223	fclose(gapfile);
224	#endif
225
226	if (oldorder==model_order \|\| tsm.m_gap==0)
227	break;
228
229	oldorder = model_order;
230	}
231	free(ts.gap_int); free(ts.data_int);
232	memcpy(isgood,model_gaps,n*sizeof(int));
233	for (i=0;i<n;i++)
234	if (model_gaps[i]==0)
235	data[i] = 0;
236
237	free(model_gaps);
238
239
240	if (padends)
241	{
242	memcpy(data_in, &data[maxorder],nold*sizeof(CTYPE));
243	memcpy(isgood_in,&isgood[maxorder],nold*sizeof(int));
244	free(data);
245	free(isgood);
246	}
247
248	if (ar_coeff_in==NULL)
249	free(ar_coeff);
250
251	if (order)
252	*order = model_order;
253
254	return 0;
255	}
256
257
258	static void FILL_GAPS(int n, CTYPE data, int isgood, int order,
259	CTYPE *ar_coeff)
260	{
261	int iter, idx, i, shift, gap_idx, gap_idx_inv, *block;
262	int n_gap, n_data, n_blocks, maxblock;
263	CTYPE rhs, rhs2, data_gap, gamma2, *scratch;
264	RTYPE rnorm, tol;
265	void aarg[3], marg[3];
266	CTYPE *gamma;
267
268
269	// for (i=0; i<10; i++)
270	// printf("data[%d] = %f%+fi.\n",i,creal(data[i]),cimag(data[i]));
271
272	// Compute gamma containing the first row of T^H*T.
273	gamma = (CTYPE )malloc((order+1)sizeof(CTYPE));
274	rhs = (CTYPE )malloc(2(n-order)*sizeof(CTYPE));
275	for (shift = 0; shift<=order; shift++)
276	{
277	gamma[shift] = 0;
278	for (i=0; (i+shift)<=order; i++)
279	gamma[shift] += 2CONJ(ar_coeff[i+shift])ar_coeff[i];
280	// printf("gamma[%d] = %f%+fi\n",shift,creal(gamma[shift]),cimag(gamma[shift]));
281	}
282	// for (i=0; i<=order; i++)
283	// printf("ar_coeff[%d] = %f%+fi.\n",i,creal(ar_coeff[i]),cimag(ar_coeff[i]));
284
285
286	// Find gaps.
287	gap_idx = (int*)calloc(n,sizeof(int));
288	gap_idx_inv = (int*)calloc(n,sizeof(int));
289	block = (int*)calloc(n,sizeof(int));
290	n_data = 0;
291	n_gap = 0;
292	n_blocks = 0;
293	maxblock = 0;
294	for (i=0; i<n; i++)
295	{
296	if (isgood[i])
297	{
298	n_data++;
299	if (i>0 && !isgood[i-1])
300	n_blocks++;
301	}
302	else
303	{
304	// printf("gap at %4d\n",i);
305	gap_idx[i] = n_gap;
306	gap_idx_inv[n_gap] = i;
307	n_gap++;
308	block[n_blocks]++;
309	if (block[n_blocks] > maxblock)
310	maxblock = block[n_blocks];
311	}
312	}
313
314	/* Set up right-hand side for block Toeplitz system. */
315	for (shift=0; (shift+order)<n; shift++)
316	{
317	rhs[shift] = 0;
318	rhs[shift+(n-order)] = 0;
319	for (i=0; i<=order; i++)
320	{
321	if (isgood[i+shift])
322	{
323	rhs[shift] -= ar_coeff[order-i]*data[i+shift];
324	rhs[shift+(n-order)] -= CONJ(ar_coeff[i])*data[i+shift];
325	}
326	}
327	}
328
329	rhs2 = (CTYPE *)calloc(n_gap,sizeof(CTYPE));
330	for (shift=0; (shift+order)<n; shift++)
331	{
332	for (i=0; i<=order; i++)
333	{
334	if ( !isgood[i+shift] )
335	{
336	idx = gap_idx[i+shift];
337	// assert(idx>=0 && idx < n_gap);
338	rhs2[idx] += CONJ(ar_coeff[order-i])*rhs[shift];
339	rhs2[idx] += ar_coeff[i]*rhs[shift+(n-order)];
340
341	}
342	}
343	}
344
345	// printf("n_gap = %d\n",n_gap);
346	//for (i=0; i<n_gap; i++)
347	// printf("rhs2[%d] = %f%+fi.\n",i,creal(rhs2[i]),cimag(rhs2[i]));
348
349
350	data_gap = (CTYPE*)calloc(n_gap,sizeof(CTYPE));
351	if (n_gap == 1)
352	data_gap[0] = rhs2[0] / gamma[0];
353	else
354	{
355	// Prepare PCG inputs defining the matrix.
356	aarg[0] = (void *) &order;
357	aarg[1] = (void *) gap_idx_inv;
358	aarg[2] = (void *) gamma;
359
360	// Prepare PCG inputs defining the preconditioner.
361	gamma2 = (CTYPE *) calloc(maxblock,sizeof(CTYPE));
362	scratch = (CTYPE *) calloc(maxblock,sizeof(CTYPE));
363	for (i=0; i<=min(order,maxblock-1); i++)
364	gamma2[i] = gamma[i];
365	marg[0] = (void *) block;
366	marg[1] = (void *) gamma2;
367	marg[2] = (void *) scratch;
368
369	// Solve normal equations using PCG with a block Toeplitz preconditioner.
370	memset(data_gap, 0, n_gap*sizeof(CTYPE));
371	tol = sqrt(order)*PCG_TOL;
372	iter = PCG(n_gap, PCG_MAXIT+2, tol, TOEPLITZ_MULT,
373	PRECONDITION, rhs2, data_gap, &rnorm, aarg, marg);
374	if (rnorm>tol)
375	fprintf(stderr, "Warning: PCG did not converge. "
376	"After %d iterations (maxit=%d) the relative"
377	" residual was %e (tol=%e)\n",
378	iter, PCG_MAXIT+2, rnorm,tol);
379
380	/* printf("n_gap=%3d, maxblock=%3d, iter=%2d, rnorm=%e\n",
381	n_gap, maxblock, iter,rnorm); */
382	// for (i=0; i<n_gap; i++)
383	//printf("x[%d] = %f%+fi\n",i,creal(data_gap[i]),cimag(data_gap[i]));
384	free(scratch); free(gamma2);
385	}
386	// Fill gap the original time series with extrapolated data.
387	for (i=0; i<n_gap; i++)
388	data[gap_idx_inv[i]] = data_gap[i];
389
390	free(block); free(gap_idx); free(gap_idx_inv);
391	free(data_gap); free(gamma); free(rhs); free(rhs2);
392	}
393
394
395
396	static int GAPSTRUCTURE(int n, CTYPE data, int isgood, gapped_timeseries *ts,
397	int data_length, CTYPE *good_data, int minpercentage)
398	{
399	int i,j,first;
400
401	// Count number of good and bad points and
402	// number of good and bad intervals.
403	ts->n_data = 0; ts->m_data = 0;
404	ts->n_gap = 0; ts->m_gap = 0;
405	for (i=0; i<n-1; i++)
406	{
407	// assert(isgood[i]==1 \|\| isgood[i]==0);
408	if (isgood[i]==1)
409	{
410	ts->n_data++;
411	if (i==0)
412	{
413	ts->first_is_good = 1;
414	ts->m_data++;
415	}
416	if (isgood[i+1]==0)
417	ts->m_gap++;
418	}
419	else if (isgood[i]==0)
420	{
421	ts->n_gap++;
422	if (i==0)
423	{
424	ts->first_is_good = 0;
425	ts->m_gap++;
426	}
427	if (isgood[i+1]==1)
428	ts->m_data++;
429	}
430	else
431	{
432	fprintf(stderr, "ERROR: problem in routine gapstructure: value of window function must be either 0 or 1\n");
433	exit(1);
434	}
435	}
436	if (isgood[n-1]==1)
437	ts->n_data++;
438	else
439	ts->n_gap++;
440
441
442	#ifdef DEBUGOUT
443	printf("n_data = %d, m_data = %d\nn_gap = %d, m_gap = %d\n",
444	ts->n_data,ts->m_data,ts->n_gap,ts->m_gap);
445	#endif
446	/*
447	assert((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) \|\|
448	(isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) \|\|
449	(isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data));
450	*/
451
452	if (!((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) \|\|
453	(isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) \|\|
454	(isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data)))
455	{
456	fprintf(stderr, "ERROR: problem in routine gapstructure\n");
457	exit(1);
458	}
459
460
461	// Allocate space for data pointers and AR coefficients.
462	good_data = (CTYPE ) malloc(ts->m_datasizeof(CTYPE *));
463	data_length = (int ) malloc(ts->m_data*sizeof(int));
464	ts->gap_int = (interval ) malloc(ts->m_gapsizeof(interval));
465	ts->data_int = (interval ) malloc(ts->m_datasizeof(interval));
466
467
468	// Scan again and set up data structures describing gap structure.
469	j = 0;
470	first=0;
471	if (isgood[0]==1)
472	{
473	first=0;
474	ts->data_int[0].first = 0;
475	(*good_data)[0] = &data[0];
476	j++;
477	}
478	for (i=0; i<n-1; i++)
479	{
480	// first of an interval
481	if (isgood[i]==0 && isgood[i+1]==1)
482	{
483	first = i+1;
484	ts->data_int[j].first = first;
485	(*good_data)[j] = &(data[i+1]);
486	// printf("first[%d] = %d\n",j,i+1);
487	j++;
488	}
489	// end of an interval
490	if (isgood[i]==1 && isgood[i+1]==0)
491	{
492	(*data_length)[j-1] = i-first+1;
493	ts->data_int[j-1].last = i;
494	}
495	}
496	if (isgood[n-1]==1)
497	{
498	(*data_length)[j-1] = n-1-first+1;
499	ts->data_int[j-1].last = n-1;
500	}
501	if (ts->first_is_good)
502	{
503	for(j=0; j<ts->m_data;j++)
504	{
505	if (j<ts->m_gap)
506	ts->gap_int[j].first = ts->data_int[j].last+1;
507	if (j>0)
508	ts->gap_int[j-1].last = ts->data_int[j].first-1;
509	}
510	if (ts->m_data==ts->m_gap)
511	ts->gap_int[ts->m_gap-1].last = n-1;
512	}
513	else
514	{
515	ts->gap_int[0].first = 0;
516	for(j=0; j<ts->m_data;j++)
517	{
518	if (j<ts->m_gap)
519	ts->gap_int[j].last = ts->data_int[j].first-1;
520	if (j<ts->m_gap-1)
521	ts->gap_int[j+1].first = ts->data_int[j].last+1;
522	}
523	if (ts->m_gap==ts->m_data+1)
524	ts->gap_int[ts->m_gap-1].last = n-1;
525	}
526	return maxorder(ts, *data_length, minpercentage);
527	}
528
529
530
531	// Matrix-vector multiply with symmetric Toeplitz matrix.
532	// y = T(r) x
533
534	static void TOEPLITZ_MULT(int n, CTYPE x[n], CTYPE y[n], void **data)
535	{
536	int i,j,order,idx,ii;
537	int *gap_idx_inv;
538	CTYPE *r;
539
540	order = ((int ) data[0]);
541	gap_idx_inv = (int *) data[1];
542	r = (CTYPE *) data[2];
543
544	for (i=0; i<n; i++)
545	y[i] = r[0]*x[i];
546
547	for (i=0; i<n; i++)
548	{
549	ii = gap_idx_inv[i];
550	for (j=i+1; j<n; j++)
551	{
552	idx = ii-gap_idx_inv[j];
553	if (idx <= order && idx >= -order)
554	{
555	if (idx>=0)
556	{
557	y[i] += CONJ(r[idx])*x[j];
558	y[j] += r[idx]*x[i];
559	}
560	else
561	{
562	y[i] += r[-idx]*x[j];
563	y[j] += CONJ(r[-idx])*x[i];
564	}
565	}
566	}
567	}
568	}
569
570	// Solve M y = x;
571	static void PRECONDITION(int n, CTYPE b, CTYPE x, void **data)
572	{
573	int i,first,*block;
574	CTYPE gamma, scratch;
575	block = (int *) data[0];
576	gamma = (CTYPE *) data[1];
577	scratch = (CTYPE *) data[2];
578
579
580	for(i=0, first = 0; first<n; first+=block[i], i++)
581	{
582	memcpy(scratch,&b[first],block[i]*sizeof(CTYPE));
583	LEVINSON(block[i],gamma,scratch,&x[first]);
584	}
585	}
586
587
588	#undef GAPSTRUCTURE
589	#undef FAHLMAN_ULRYCH
590	#undef BURG
591	#undef FILL_GAPS
592	#undef TOEPLITZ_MULT
593	#undef PRECONDITION
594	#undef PCG
595	#undef LEVINSON
596	#undef EPSILON
597	#include "ctypes_undef.h"