(file) Return to jtsslice.c CVS log (file) (dir) Up to [Development] / JSOC / proj / globalhs / apps

   1 tplarson 1.1 /* this JSOC module is a port of an SOI module written by Rasmus Larsen.
   2               * ported by Tim Larson.
   3               */
   4              
   5              /* 
   6               *  Detrending/gapfilling/differencing module
   7               *  ts_fiddle_rml upgrades ts_fiddle
   8               */
   9              
  10              #include <stdio.h>
  11              #include <stdlib.h>
  12              #include <strings.h>
  13              #include <time.h>
  14              #include <sys/time.h>
  15              #include <sys/resource.h>
  16              #include <math.h>
  17              #include <assert.h>
  18              //#include <complex.h>
  19              #include <fftw3.h>
  20              #include "fahlman_ulrych_C99.h"
  21              #include "detrend_C99.h"
  22 tplarson 1.1 //#include <xmem.h>
  23              
  24              #include "jsoc_main.h"
  25              #include "/home/jsoc/include/fitsio.h"
  26              
  27              #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
  28              #define kNOTSPECIFIED "not specified"
  29              
  30              char *module_name = "jtsslice";
  31              
  32              /* Command line arguments: */
  33              ModuleArgs_t module_args[] =
  34              {
  35                {ARG_STRING,   "in", "", "input data records"},
  36                {ARG_STRING,   "out", "", "output data series"},
  37                {ARG_STRING,   "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"},
  38                {ARG_STRING,   "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
  39                {ARG_STRING,   "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
  40                {ARG_STRING,   "srclink",  "SRCDATA", "name of link to source data"},
  41                {ARG_INT,      "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
  42                {ARG_INT,      "VERB", "1",  "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL},
  43 tplarson 1.1   {ARG_STRING,   "LOGFILE", "stdout", ""},
  44                {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},
  45                {ARG_STRING,   "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"},
  46              
  47                {ARG_INT,      "NDT", "-1", "", ""},
  48              //  {ARG_FLOAT,    "TSAMPLE", "60", "", ""},
  49                {ARG_INT,      "IFILL", "0", "", ""}, //changed from jtsfiddle
  50                {ARG_INT,      "IFIRST", "0", "", ""},
  51                {ARG_INT,      "FLIPM", "1", "", ""},
  52                {ARG_INT,      "MAX_AR_ORDER", "360", "", ""},
  53                {ARG_INT,      "FU_FIT_EDGES", "1", "", ""},
  54                {ARG_INT,      "FU_ITER", "3", "", ""},
  55                {ARG_INT,      "FU_MIN_PERCENT", "90", "", ""},
  56                {ARG_FLOAT,    "CDETREND", "0.0", "", ""}, //changed from jtsfiddle
  57                {ARG_INT,      "DETREND_LENGTH", "1600", "", ""},
  58                {ARG_INT,      "DETREND_SKIP", "1440", "", ""},
  59                {ARG_INT,      "TSOUT", "0", "", ""},
  60                {ARG_INT,      "FFTOUT", "0", "", ""},
  61                {ARG_INT,      "FFT1OUT", "0", "", ""},
  62                {ARG_INT,      "POWOUT", "0", "", ""},
  63                {ARG_INT,      "ARPOWOUT", "0", "", ""},
  64 tplarson 1.1   {ARG_INT,      "MAVGOUT", "0", "", ""},
  65                {ARG_INT,      "NAVG", "0", "", ""},
  66                {ARG_INT,      "NSKIP", "0", "", ""},
  67                {ARG_INT,      "NOFF", "0", "", ""},
  68                {ARG_INT,      "IMIN", "0", "", ""},
  69                {ARG_INT,      "IMAX", "-1", "", ""},
  70                {ARG_END},
  71              };
  72              
  73              #include "saveparm.c"
  74              #include "timing.c"
  75              
  76              #define DIE(code) { fprintf(stderr,"jtsfiddle died with error code %d\n",(code)); return 0;}
  77              
  78              /* global variables holding the values of the command line variables. */
  79              static char *logfile, *gapfile, *sectionfile;
  80              static int lmode, n_sampled_out, ifill, flipm;
  81              static int ar_order, max_ar_order, fu_iter, fu_fit_edges;
  82              static int fu_min_percent;
  83              static int detrend_length, detrend_skip;
  84              static float tsample, detrend_const;
  85 tplarson 1.1 static int ifirst, tsout,fftout, fft1out, powout, arpowout, mavgout;
  86              static int navg, nskip, noff, imin, imax;
  87              
  88              
  89              /* Prototypes for local functions. */
  90              static double splitting(int l, int m);
  91              
  92              /* Prototypes for external functions */
  93              extern void cffti_(int *, float *);
  94              extern void cfftb_(int *, float *, float *);
  95              static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
  96              		  float tsample_in, float *tsample_out,
  97              		  int *num_sections, int *sections);
  98              static void extract_data(int n_in, int *gaps, float *data_in, float *data_out);
  99              static int read_section_file(char *filename, int *num_section, int **section);
 100              
 101              /************ Here begins the main module **************/
 102              int DoIt(void)
 103              {
 104                int i;
 105                int start_index[2], counts[2], intervals[2];
 106 tplarson 1.1   int m, n_in, n_sampled_in;
 107                TIME epoch, step, start, stop;
 108                int gapsize, istart, istop, data_dim, detrend_order;
 109                int npow, mshift;
 110                float tsamplex, df1, c;
 111                int *agaps;
 112                int *gaps, *gaps2;
 113                float *msum, *pow;
 114                float *data, *out_data, *in_data;
 115                fftwf_plan plan;
 116                int num_sections, *sections;
 117                float *ar_coeff=NULL;
 118              
 119                int fstatus = 0;
 120                fitsfile *fitsptr;
 121                long fdims[1];
 122                long fpixel[]={1};
 123                int *anynul=0;
 124              
 125                int length[2], startind[2], endind[2];
 126                float *datarr;
 127 tplarson 1.1   int status=DRMS_SUCCESS;
 128                int newstat=0;
 129              
 130                DRMS_Segment_t *segin = NULL;
 131                DRMS_Segment_t *segout = NULL;
 132                DRMS_Array_t *inarr = NULL;
 133                DRMS_Array_t *outarr = NULL;
 134                DRMS_Record_t *inrec = NULL;
 135                DRMS_Record_t *outrec = NULL;
 136                DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
 137                DRMS_RecLifetime_t lifetime;
 138                long long histrecnum = -1;
 139                DRMS_Segment_t *gapseg = NULL;
 140                DRMS_Array_t *gaparr = NULL;
 141              
 142                HIterator_t outKey_list;
 143                DRMS_Keyword_t *outKey;
 144                TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
 145              
 146                int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
 147                int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
 148 tplarson 1.1 
 149                double tstart, tstartsave, tstop, tstopsave, tstep, tstartout, tstopout;
 150                double wt0, wt1, wt2, wt3, wt;
 151                double ut0, ut1, ut2, ut3, ut;
 152                double st0, st1, st2, st3, st;
 153                double ct0, ct1, ct2, ct3, ct;
 154              
 155                wt0=getwalltime();
 156                ct0=getcputime(&ut0, &st0);
 157              
 158                /* Read input parameters. */
 159                logfile     = cmdparams_save_str    (&cmdparams, "LOGFILE", &newstat);
 160                gapfile     = cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
 161                sectionfile = cmdparams_save_str    (&cmdparams, "SECTIONFILE", &newstat);
 162                n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);
 163              //  tsample     = cmdparams_save_float  (&cmdparams, "TSAMPLE", &newstat);
 164                ifill       = cmdparams_save_int    (&cmdparams, "IFILL", &newstat);
 165                max_ar_order= cmdparams_save_int    (&cmdparams, "MAX_AR_ORDER", &newstat);
 166                fu_fit_edges= cmdparams_save_int    (&cmdparams, "FU_FIT_EDGES", &newstat);
 167                fu_iter     = cmdparams_save_int    (&cmdparams, "FU_ITER", &newstat);
 168                fu_min_percent = cmdparams_save_int (&cmdparams, "FU_MIN_PERCENT", &newstat);
 169 tplarson 1.1   ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
 170                flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);
 171                detrend_const = cmdparams_save_float(&cmdparams, "CDETREND", &newstat);
 172                detrend_length = cmdparams_save_int (&cmdparams, "DETREND_LENGTH", &newstat);
 173                detrend_skip = cmdparams_save_int   (&cmdparams, "DETREND_SKIP", &newstat);
 174                tsout       = cmdparams_save_int    (&cmdparams, "TSOUT", &newstat);
 175                fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);
 176                fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);
 177                powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);
 178                arpowout    = cmdparams_save_int    (&cmdparams, "ARPOWOUT", &newstat);
 179                mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);
 180                navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);
 181                nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);
 182                noff        = cmdparams_save_int    (&cmdparams, "NOFF", &newstat);
 183                imin        = cmdparams_save_int    (&cmdparams, "IMIN", &newstat);
 184                imax        = cmdparams_save_int    (&cmdparams, "IMAX", &newstat);
 185              
 186                char *inRecQuery = cmdparams_save_str(&cmdparams, "in", &newstat);
 187                char *outSeries = cmdparams_save_str(&cmdparams, "out", &newstat);
 188                char *segnamein = cmdparams_save_str(&cmdparams, "segin", &newstat);
 189                char *segnameout = cmdparams_save_str(&cmdparams, "segout", &newstat);
 190 tplarson 1.1   int seginflag = strcmp(kNOTSPECIFIED, segnamein);
 191                int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
 192                char *histlinkname = cmdparams_save_str(&cmdparams, "histlink", &newstat);
 193                char *srclinkname = cmdparams_save_str(&cmdparams, "srclink", &newstat);
 194                int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
 195                int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
 196                if (permflag)
 197                  lifetime = DRMS_PERMANENT;
 198                else
 199                  lifetime = DRMS_TRANSIENT;
 200              
 201              #include "histinclude.c"
 202              
 203                /**** Sanity check of input parameters. ****/
 204              
 205                /* Default is to just output the timeseries. */
 206                if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0)) 
 207                  tsout=1;
 208                /* Only one type of output allowed at a time */
 209                if ((tsout+fftout+fft1out+powout+mavgout+arpowout) !=1) 
 210                {
 211 tplarson 1.1     fprintf(stderr, "ERROR: only one type of output allowed at a time\n");
 212                  return 0;
 213                }
 214                if (navg <= 0) 
 215                  navg=1;
 216                if (nskip <= 0) 
 217                  nskip=1;
 218                if ((navg != 1) && (nskip != 1))
 219                {
 220                  fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");
 221                  return 0;
 222                }
 223                if (noff < 0 || noff >= nskip) 
 224                {
 225                  fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
 226                  return 0;
 227                }
 228                if(arpowout && !ifill)
 229                {
 230                  fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");
 231                  return 0;
 232 tplarson 1.1   }
 233              
 234              
 235               // these must be present in the output dataseries and variable, not links or constants   
 236                char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX"};
 237              
 238                int itest;
 239                for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
 240                {
 241                  DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
 242                  if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)
 243                  {
 244                    fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
 245                    drms_close_record(tempoutrec, DRMS_FREE_RECORD);
 246                    return 0;
 247                  }
 248                }
 249              
 250                drms_close_record(tempoutrec, DRMS_FREE_RECORD);
 251              
 252              //  xmem_off();
 253 tplarson 1.1 
 254              
 255                if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
 256                {
 257                  DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status);
 258                  if (status != DRMS_SUCCESS || gaprecset == NULL) 
 259                  {
 260                    fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
 261                    return 0;
 262                  }
 263                  if (gaprecset->n == 0)
 264                  {
 265                    fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
 266                    return 0;
 267                  }
 268                  if (gaprecset->n > 1)
 269                  {
 270                    fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
 271                    return 0;
 272                  } 
 273                  gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);
 274 tplarson 1.1     if (gapseg != NULL)
 275                    gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status);
 276                  if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
 277                  {
 278                    fprintf(stderr, "ERROR: problem reading gaps segment: status = %d\n", status);
 279                    return 0;
 280                  }
 281                  agaps=(int *)(gaparr->data);
 282                  gapsize=gaparr->axis[0];
 283                }
 284                else
 285                {
 286                  fits_get_img_size(fitsptr, 1, fdims, &fstatus);
 287                  gapsize=fdims[0];
 288                  agaps=(int *)(malloc(gapsize*sizeof(int)));
 289                  fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus); 
 290                  fits_close_file(fitsptr, &fstatus);
 291                  if (fstatus) 
 292                  {
 293                    fprintf(stderr, "ERROR: ");
 294                    fits_report_error(stderr, fstatus);
 295 tplarson 1.1       return 0;
 296                  }
 297                }
 298              
 299                printf("gapfile read, gapsize = %d\n",gapsize);
 300              
 301                data_dim=gapsize;
 302                if (n_sampled_out>gapsize) 
 303                  data_dim=n_sampled_out;
 304              
 305                /* Read location of jumps. */
 306                if (!strcmp(sectionfile,kNOTSPECIFIED))
 307                {
 308                  /* No sections file specified. Assume that the whole 
 309                     time series in one section. */
 310                  num_sections=1;
 311                  sections = malloc(2*sizeof(int));
 312                  sections[0] = 0;
 313                  sections[1] = data_dim-1;
 314                }
 315                else
 316 tplarson 1.1   {
 317                  int secstat;
 318                  if (secstat=read_section_file(sectionfile, &num_sections, &sections))
 319                  {
 320                    DRMS_RecordSet_t *secrecset = drms_open_records(drms_env, sectionfile, &status);
 321                    if (status != DRMS_SUCCESS || secrecset == NULL) 
 322                    {
 323                      fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status);
 324                      return 0;
 325                    }
 326                    if (secrecset->n == 0)
 327                    {
 328                      fprintf(stderr,"ERROR: sectionfile recordset contains no records\n");
 329                      return 0;
 330                    }
 331                    if (secrecset->n > 1)
 332                    {
 333                      fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n");
 334                      return 0;
 335                    }
 336                    num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL);
 337 tplarson 1.1       if (num_sections < 1)
 338                    {
 339                      fprintf(stderr,"ERROR: Invalid NSECS keywords\n");
 340                      return 0;
 341                    }
 342                    sections = (int *)malloc(2*(num_sections)*sizeof(int));
 343                    char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL);
 344                    sscanf(strtok(sectkey," \n"),"%d",&(sections[0]));
 345                    sscanf(strtok(NULL," \n"),"%d",&(sections[1]));
 346                    int i;
 347                    for (i=2 ;i<2*(num_sections); i+=2)
 348                    {
 349                      sscanf(strtok(NULL," \n"), "%d", &(sections)[i]);
 350                      sscanf(strtok(NULL," \n"), "%d", &(sections)[i+1]);
 351              
 352                      if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1]))
 353                      {
 354                        fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n");
 355                        return 0;
 356                      }
 357                    }
 358 tplarson 1.1       free(sectkey);
 359                  }
 360                }
 361              
 362                printf("num_sections = %d\n",num_sections);
 363              
 364                /* allocate temporary storage */
 365                gaps=(int *)(malloc(gapsize*sizeof(int)));
 366                gaps2=(int *)(malloc(gapsize*sizeof(int)));
 367                data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
 368                ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float));
 369              
 370              
 371              
 372                DRMS_RecordSet_t *inRecSet = drms_open_records(drms_env, inRecQuery, &status);
 373              
 374                if (status != DRMS_SUCCESS || inRecSet == NULL) 
 375                {
 376                  fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
 377                  return 0;
 378                }
 379 tplarson 1.1 
 380                if (inRecSet->n == 0)
 381                {
 382                  printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]);
 383                  return 0;
 384                }
 385              
 386                if (inRecSet->n > 1)
 387                {
 388                  printf("ERROR: nRecs > 1 not yet supported.\n");
 389                  return 0;
 390                }
 391              
 392              
 393                inrec=inRecSet->records[0];
 394                char *inchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "T_STEP"};
 395              
 396                for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
 397                {
 398                  DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
 399                  if (!inkeytest)
 400 tplarson 1.1     {
 401                    fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
 402                    drms_close_records(inRecSet, DRMS_FREE_RECORD);
 403                    return 0;
 404                  }
 405                }
 406              
 407              
 408              // insert loop over input here
 409              int iRec=0;
 410                inrec = inRecSet->records[iRec];
 411              
 412                int lmin=drms_getkey_int(inrec, "LMIN", NULL);
 413                int lmax=drms_getkey_int(inrec, "LMAX", NULL);
 414                if (lmin != lmax)
 415                {
 416                  printf("ERROR: lmin != lmax not yet supported.\n");
 417                  return 0;
 418                }
 419                lmode=lmin;
 420              
 421 tplarson 1.1 
 422                tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL);
 423                tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL);
 424                tstep=drms_getkey_float(inrec, "T_STEP", NULL);
 425              
 426                n_in=(tstop-tstart)/tstep;
 427                if (n_in != gapsize)
 428                {
 429                  fprintf(stderr, "ERROR: gaps seem incompatible with time-series, gapsize=%d, n_in=%d\n", gapsize, n_in);
 430                  return 0;
 431                }
 432              
 433                tsample=tstep;
 434              
 435              
 436              
 437                // changed n_in to gapsize in following call
 438                if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
 439                {
 440                  fprintf(stderr, "ERROR: problem in extract_gaps\n");
 441                  return 0;
 442 tplarson 1.1   }
 443              
 444                /* Adjust detrend parameters according to the number of
 445                   available points. */
 446                if (n_sampled_in<detrend_length)
 447                {
 448                  detrend_length = n_sampled_in;
 449                  detrend_skip = detrend_length;
 450                }	
 451              
 452                if (n_sampled_out == -1) 
 453                  n_sampled_out = n_sampled_in;
 454                df1 = tsamplex*n_sampled_out;
 455              
 456                if (fftout || fft1out || powout || arpowout || mavgout)     
 457                  plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, 
 458              			    (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);
 459              
 460                /* Set sum to zero before starting */
 461                if (mavgout) 
 462                  msum = (float *)calloc(n_sampled_out/2,sizeof(float));
 463 tplarson 1.1 
 464              
 465                tstartout=tstart+ifirst*tsample;
 466                tstopout=tstartout+df1;
 467              
 468              
 469                /* Create output data set. */
 470                if (tsout || fftout) 
 471                {
 472                  length[0]=2*n_sampled_out;
 473                  length[1]=lmode+1;
 474                }
 475                else 
 476                {
 477                  /* Set first and last output index if not set as a parameter. */
 478                  if (imax < imin) 
 479                  {
 480                    imin=0;
 481                    imax=n_sampled_out/2-1;
 482                  }
 483                  npow=imax-imin+1;
 484 tplarson 1.1     pow=(float *)(malloc(n_sampled_out*sizeof(float)));
 485                  if (fft1out) 
 486                    length[0]=2*npow;
 487                  else 
 488                    length[0]=npow;
 489                  if (powout || arpowout || fft1out) 
 490                    length[1]=2*lmode+1;
 491                  else 
 492                    length[1]=1;
 493                }
 494              
 495                
 496                /* Loops through all the data and update the gaps to reflect missing 
 497                   data. */
 498                /* the above comment is from the original code.  it refers to an ifdef
 499                   that would read in all the data and run extract_data on it in case there 
 500                   were nan's not removed in the input window function.  the present version
 501                   assumes the window function is correct. */
 502              
 503                startind[0]=2*ifirst;
 504                endind[0]=2*(ifirst + n_sampled_out) - 1;
 505 tplarson 1.1   startind[1]=0;
 506                endind[1]=lmode;
 507              
 508                if (seginflag)
 509                  segin = drms_segment_lookup(inrec, segnamein);
 510                else
 511                  segin = drms_segment_lookupnum(inrec, 0);
 512                if (segin)
 513              //    inarr = drms_segment_read(segin, usetype, &status);
 514                  inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);
 515                if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL)
 516                {
 517                  fprintf(stderr, "ERROR: problem reading input segment: iRec = %d, status = %d\n", iRec, status);
 518                  return 0;
 519                }
 520                datarr=(float *)(inarr->data);
 521              
 522              
 523                outrec = drms_create_record(drms_env, 
 524                                            outSeries,
 525                                            lifetime, 
 526 tplarson 1.1                               &status);
 527              
 528                if (status != DRMS_SUCCESS || outrec==NULL)
 529                {
 530                   fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
 531                   return 0;
 532                }
 533              
 534                if (histlink)
 535                   drms_setlink_static(outrec, histlinkname,  histrecnum);
 536                if (srclink)
 537                   drms_setlink_static(outrec, srclinkname,  inrec->recnum);
 538              
 539                if (segoutflag)
 540                  segout = drms_segment_lookup(outrec, segnameout);
 541                else
 542                  segout = drms_segment_lookupnum(outrec, 0);
 543                outarr = drms_array_create(usetype, 2, length, NULL, &status);
 544                if (status != DRMS_SUCCESS || outarr == NULL || segout == NULL)
 545                {
 546                  fprintf(stderr,"ERROR: problem creating output array or segment: length = [%d, %d], status = %d", length[0], length[1], status);
 547 tplarson 1.1     return 0; 
 548                }
 549                out_data = outarr->data;
 550              
 551              // set ifirst=0 because data is read starting at ifirst by drms_segment_readslice
 552              // set gapoffset to ifirst to start reading gaps at same location
 553                int gapoffset = ifirst;
 554                ifirst=0;
 555                /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
 556                for (m=0; m<=lmode; m++) 
 557                {
 558              //    (*history)("### processing m = %d ###\n",m);
 559                  /* Read chunk of data corresponding to next m from the data file. */
 560              /*
 561                  start_index[1] = m;
 562                  if (!(tile_sds = vds_select(in_vds,0,0,1,1,start_index,counts,intervals))) 
 563                  {
 564                    fprintf(stderr,"vds_select failed");
 565                    DIE(9);
 566                  }
 567                  if (!sds_data(tile_sds)) 
 568 tplarson 1.1     {
 569                    fprintf(stderr,"*** lnu vds_select returned no data.\n");
 570                    sds_free(&tile_sds);
 571                    DIE(10);
 572                  }
 573                  in_data = (float *)(sds_data(tile_sds));
 574                  if (sds_dim0(tile_sds)<0) 
 575                  {
 576                    sds_free(&tile_sds);
 577                    DIE(soi_errno);
 578                  }
 579              */
 580              
 581              //    in_data=datarr + m*2*gapsize;
 582                  in_data=datarr + m*2*n_sampled_out;
 583              // changed this expression after using drms_segment_readslice above
 584              
 585              
 586                  /* Extracts data copy */
 587              //    extract_data(n_sampled_in, gaps, in_data, data);
 588                 extract_data(n_sampled_out, gaps+gapoffset, in_data, data);
 589 tplarson 1.1 // likewise
 590              
 591                  /* Detrend entire time series if required */
 592                  if (detrend_const != 0) 
 593                  {
 594                    detrend_order = floor(detrend_length/detrend_const)+2;
 595              
 596              //      for (i=0;i<num_sections; i++)
 597              //        printf("[%d:%d]\n",sections[2*i],sections[2*i+1]);
 598                    cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps+gapoffset, 
 599              			 detrend_order, detrend_length, detrend_skip, 
 600              			 num_sections,sections);
 601                  }
 602              
 603                  /* Fill gaps in entire time series if required */
 604                  memcpy(gaps2, gaps, n_sampled_in*sizeof(int));
 605                  if (ifill != 0 && max_ar_order > 0) 
 606                  {
 607              //      (*history)("Filling gaps, ");
 608                    cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2+gapoffset,
 609              		      fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, 
 610 tplarson 1.1 		      &ar_order, (_Complex float *)ar_coeff);
 611              /*
 612                    { 
 613              	FILE *fh;
 614              	unsigned char c;
 615              	fh = fopen("gaps.bin","w");
 616              	for(i=0; i<n_sampled_in; i++)
 617              	{
 618              	  c = (unsigned char) gaps2[i];
 619              	  fwrite(&c,1,1,fh);
 620              	}
 621              	fclose(fh);
 622                    }
 623              */
 624              //      (*history)("AR model order used = %d.\n",ar_order);
 625                  }
 626              
 627              
 628                  /* pad with zeros if the last point output (n_sampled_out+ifirst) 
 629                     is after the last data points read in (nread). */
 630                  if (arpowout)
 631 tplarson 1.1     {
 632                    memmove(data, ar_coeff, 2*(ar_order+1)*sizeof(float));
 633                    memset(&data[2*(ar_order+1)],0,2*(n_sampled_out-ar_order-1)*sizeof(float));
 634                  }
 635                  else
 636                  {
 637                    memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
 638                    if ((ifirst+n_sampled_out) >= n_sampled_in) 
 639              	memset(&data[2*n_sampled_in], 0, 
 640              	       2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));      
 641                  }
 642                  if (fftout || fft1out || powout || arpowout || mavgout) 
 643                  {
 644              //      (*history)("Computing FFT.\n");
 645                    fftwf_execute(plan);
 646                  }
 647                    
 648              
 649              //    (*history)("Preparing output.\n");
 650              
 651                  if (tsout || fftout) 
 652 tplarson 1.1       memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));
 653              
 654                  else if (fft1out) 
 655                  {
 656                    /* Do positive m */
 657                    memcpy(&out_data[2*(lmode+m)*npow], &data[2*imin], 2*npow*sizeof(float));
 658                    if (m > 0) 
 659                    {
 660              	/* Do negative m */
 661              	for (i=0; i<npow; i++) 
 662              	{
 663              	  /* Conjugate for negative m */
 664              	  out_data[2*((lmode-m)*npow+i)] = data[2*(n_sampled_out-imin-i)];
 665              	  out_data[2*((lmode-m)*npow+i)+1] = -data[2*(n_sampled_out-imin-i)+1];
 666              	}
 667                    }
 668                  }
 669                  else 
 670                  {
 671                    /* Compute power spectrum from complex FFT. */      
 672                    if (arpowout)
 673 tplarson 1.1       { 
 674              	for (i = 0; i < n_sampled_out; i++) 
 675              	  pow[i] = 1.0/(data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]);
 676                    }
 677                    else
 678                    {
 679              	for (i = 0; i < n_sampled_out; i++) 
 680              	  pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];
 681                    }
 682                    if (powout || arpowout) 
 683                    {
 684              	/* Do positive m */
 685              	memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));
 686              	if (m > 0) 
 687              	{
 688              	  /* Do negative m */
 689              	  if (imin)
 690              	    out_data[(lmode-m)*npow] = pow[n_sampled_out-imin];
 691              	  else
 692              	    out_data[(lmode-m)*npow] = pow[0];
 693              	  for (i=1; i<npow;i++) 
 694 tplarson 1.1 	    out_data[(lmode-m)*npow+i] = pow[n_sampled_out-imin-i];
 695              	}
 696                    }
 697                    else 
 698                    {
 699              	/* m-averaged output */
 700              	/* Sum all freqs, not only those in output. Should be fixed. */
 701              	/* Maybe should allow for wrapping around Nyquist. */
 702              	/* Do positive m */
 703              	mshift = floor(df1*splitting(lmode,m)+0.5f);
 704              	if (mshift >= 0) 
 705              	  for (i=0; i<n_sampled_out/2-mshift; i++)
 706              	    msum[i] += pow[mshift+i];
 707              	else
 708              	  for (i=0; i<n_sampled_out/2+mshift; i++)
 709              	    msum[i-mshift] += pow[i];
 710              	if (m > 0) {
 711              	  /* Do negative m */
 712              	  mshift = floor(df1*splitting(lmode,-m)+0.5f);
 713              	  if (mshift >=0)
 714              	    for (i=1; i<n_sampled_out/2-mshift;i++)
 715 tplarson 1.1 	      msum[i] += pow[n_sampled_out-mshift-i];
 716              	  else
 717              	    for (i=1; i<n_sampled_out/2+mshift; i++)
 718              	      msum[i-mshift] += pow[n_sampled_out-i];
 719              	}
 720                    }
 721                  }
 722              
 723                } /* End of m loop */
 724              
 725              /*
 726                hiter_new(&outKey_list, &outrec->keywords);
 727                while ( (outKey=(DRMS_Keyword_t *)hiter_getnext(&outKey_list)) )
 728                   {
 729                   char *keyName = outKey->info->name;
 730                   DRMS_Value_t inValue = {DRMS_TYPE_STRING, NULL};
 731                   inValue = drms_getkey_p(inrec, keyName, &status);
 732                   if (status) fprintf(stderr,"*** COPY drms_getkey_p %s status=%d\n",keyName,status);
 733                   drms_setkey_p(outrec, keyName, &inValue);
 734                   if ((inValue.type == DRMS_TYPE_STRING) && inValue.value.string_val)
 735                   free(inValue.value.string_val);
 736 tplarson 1.1      inValue.value.string_val = NULL;
 737                   }
 738              */
 739                drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
 740                drms_setkey_int(outrec, "QUALITY", 0);
 741                drms_setkey_time(outrec, "T_START", tstartout);
 742                drms_setkey_time(outrec, "T_STOP", tstopout);
 743                drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2);
 744                drms_setkey_float(outrec, "T_STEP", tsamplex);
 745                drms_setkey_int(outrec, "NDT", n_sampled_out);
 746              
 747              
 748               
 749                tnow = (double)time(NULL);
 750                tnow += UNIX_epoch;
 751                drms_setkey_time(outrec, "DATE", tnow);
 752              
 753                drms_free_array(inarr);
 754                if (inRecSet) 
 755                {
 756                  drms_close_records(inRecSet, DRMS_FREE_RECORD);
 757 tplarson 1.1   }
 758              
 759                if (mavgout) 
 760                {
 761                  c=1.0f/(2*lmode+1);
 762                  for (i=0; i<npow; i++) 
 763                    out_data[i] = msum[imin+i] * c;
 764                }
 765              
 766                if(ar_coeff)
 767                  free(ar_coeff);
 768              
 769                drms_segment_write(segout, outarr, 0);
 770                drms_close_record(outrec, DRMS_INSERT_RECORD);
 771              
 772                if (fftout || fft1out || powout || arpowout || mavgout)   
 773                  fftwf_destroy_plan(plan); 
 774              
 775                wt=getwalltime();
 776                ct=getcputime(&ut, &st);
 777                if (verbflag) 
 778 tplarson 1.1   {
 779              //    printf("number of records created = %d\n", nsuccess);
 780                  fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
 781                          wt-wt0, ct-ct0);
 782                  printf("module %s successful completion\n", cmdparams.argv[0]);
 783                }
 784              
 785              //  xmem_leakreport();
 786                return 0;
 787              }
 788              
 789              
 790              static double splitting(int l, int m)
 791              {
 792                double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0};
 793                double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8};
 794                double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5};
 795                /* f0 is the frequency used for generating the above a-coeff */
 796                double f0 = 1500.0;
 797                /* fcol is the frequency for which to optimize the collaps */
 798                double fcol = 800.0;
 799 tplarson 1.1 
 800                double ll,x,x3,x5,a1,a2,a3,a5,frac,lx;
 801                int ix;
 802              
 803                if (l == 0) 
 804                  ll = 1; 
 805                else 
 806                  ll = sqrt(l*(l+1.));
 807                x = m/ll;
 808                x3 = x*x*x;
 809                x5 = x3*x*x;
 810                lx = 5*log10(l*f0/fcol)-4;
 811                ix = floor(lx);
 812                frac = lx-ix;
 813                if (lx <= 0) {
 814                  ix = 0;
 815                  frac = 0.0;
 816                }
 817                if (lx >=8) {
 818                  ix = 7;
 819                  frac = 1.0;
 820 tplarson 1.1   }
 821                a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1];
 822                a2 = 0.0;
 823                a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1];
 824                a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1];
 825              
 826                return ll*(a1*x+a2*(1.5*x*x-0.5)+a3*(2.5*x3-1.5*x)+a5*(63./8.*x5-70./8.*x3+15./8.*x))/1e9;
 827              }
 828              
 829              
 830              
 831              /* Extract the data samples actually used by skipping or averaging
 832                 data. Replace missing data that are not marked as gaps with zero. */
 833              void extract_data(int n_in, int *gaps, float *data_in, float *data_out)
 834              {
 835                int i,j, nmiss = 0;
 836                assert(nskip==1 || navg==1);
 837                if ((navg == 1) && (nskip == 1)) 
 838                {
 839                  for (i=0; i<n_in; i++) 
 840                  {
 841 tplarson 1.1       if (gaps[i]==0 )
 842                    {
 843              	data_out[2*i] = 0.0f;
 844              	data_out[2*i+1] = 0.0f;
 845                    }
 846                    else
 847                    {
 848              	if (isnan(data_in[2*i]) || isnan(data_in[2*i+1]))
 849              	{
 850              	  data_out[2*i] = 0.0f;
 851              	  data_out[2*i+1] = 0.0f;
 852              	  gaps[i] = 0;
 853              	  nmiss++;
 854              	}
 855              	else
 856              	{
 857              	  data_out[2*i] = data_in[2*i];
 858              	  data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1];
 859              	}
 860                    }
 861                  }
 862 tplarson 1.1   }
 863                else if (nskip != 1) 
 864                {
 865                  for (i = 0; i<n_in; i++) 
 866                  {
 867                    if (gaps[i]==0 )
 868                    {
 869              	data_out[2*i] = 0.0f;
 870              	data_out[2*i+1] = 0.0f;
 871                    }
 872                    else
 873                    {
 874              	int ix = nskip*i+noff;
 875              	if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
 876              	{
 877              	  data_out[2*i] = 0.0f;
 878              	  data_out[2*i+1] = 0.0f;
 879              	  gaps[i] = 0;
 880              	  nmiss++;
 881              	}
 882              	else
 883 tplarson 1.1 	{
 884              	  data_out[2*i] = data_in[2*ix];
 885              	  data_out[2*i+1] =  flipm ? -data_in[2*ix+1] : data_in[2*ix+1];
 886              	}
 887                    }
 888                  }
 889                }
 890                else if (navg != 1) 
 891                {
 892                  for (i = 0; i<n_in; i++) 
 893                  {
 894                    if (gaps[i]==0 )
 895                    {
 896              	data_out[2*i] = 0.0f;
 897              	data_out[2*i+1] = 0.0f;
 898                    }
 899                    else
 900                    {
 901              	float avgr = 0.0f;
 902              	float avgi = 0.0f;
 903              	for (j=0; j<navg; j++) 
 904 tplarson 1.1 	{
 905              	  int ix = navg*i+j;
 906              	  if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
 907              	  {
 908              	    data_out[2*i] = 0.0f;
 909              	    data_out[2*i+1] = 0.0f;
 910              	    gaps[i] = 0;
 911              	    nmiss++;
 912              	    break;
 913              	  }
 914              	  else 
 915              	  {
 916              	    avgr += data_in[2*ix];
 917              	    avgi += data_in[2*ix+1];
 918              	  }
 919              	}
 920              	if (j == navg) 
 921              	{
 922              	  data_out[2*i] = avgr/navg;
 923              	  data_out[2*i+1] = avgi/navg;
 924              	}
 925 tplarson 1.1       }
 926                  }  
 927                }
 928              }
 929              
 930              /* Extract the windows function for the samples actually used after
 931                 subsampling or averaging. Modify the list of continous sections
 932                 accordingly, and make sure we do not average across section boundaries. */
 933              static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
 934              		  float tsample_in, float *tsample_out, 
 935              		  int *num_sections, int *sections)
 936              {
 937                assert(nskip==1 || navg==1);
 938                int i,j,sect, start,stop;
 939              
 940              
 941                if (*num_sections<1)
 942                {
 943                  fprintf(stderr,"Apparantly no sections of data available.");
 944                  return 1;
 945                }
 946 tplarson 1.1   /* Mask out data that in between sections if this hasn't already been
 947                   done in gapfile. */
 948                for (i=0; i<sections[0] && i<n_in; i++)
 949                  gaps_in[i] = 0;
 950                for (sect=1; sect<(*num_sections); sect++)
 951                {
 952                  for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++)
 953                    gaps_in[i] = 0;
 954                }
 955                for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++)
 956                  gaps_in[i] = 0;
 957                  
 958              
 959                if ((navg == 1) && (nskip == 1)) 
 960                {
 961                  *n_out = n_in;
 962                  *tsample_out = tsample_in;
 963                  for (i=0; i<*n_out; i++) 
 964                    gaps_out[i] = gaps_in[i];
 965                }
 966                else if (nskip != 1) 
 967 tplarson 1.1   {
 968                  *n_out = n_in/nskip;
 969                  *tsample_out = nskip*tsample_in;
 970                  for (i=0; i<*n_out; i++) 
 971                  {
 972                    gaps_out[i] = gaps_in[nskip*i+noff];
 973                  }
 974                  for (i=0; i<2*(*num_sections); i+=2)
 975                  {
 976                    start = (int)ceilf(((float)(sections[i]-noff))/nskip);
 977                    stop = (int)floorf(((float)(sections[i+1]-noff))/nskip);
 978                    if (start <= stop)
 979                    {
 980              	sections[i] = start;
 981              	sections[i+1] = stop;
 982                    }
 983                    else 
 984                    {
 985              	/* This section was skipped entirely. */
 986              	for (j=i; j< 2*(*num_sections-1); j+=2) 
 987              	{
 988 tplarson 1.1 	  sections[j] = sections[j+2];
 989              	  sections[j+1] = sections[j+3];
 990              	}
 991              	*num_sections -= 1;
 992                    }
 993                  }
 994                }
 995                else  if (navg != 1) 
 996                {
 997                  sect = 0;
 998                  *n_out = n_in/navg;
 999                  *tsample_out = tsample*navg;
1000                  for (i=0; i<*n_out; i++) 
1001                  {      
1002                    int igx = 1;
1003                    while (sect < *num_sections && 
1004              	     !(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1]))
1005              	sect++;
1006              
1007                    /* Don't allow averaging of data from different sections. */
1008                    if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1])
1009 tplarson 1.1       {
1010              	for (j=0; j<navg; j++) 
1011              	  igx *= gaps_in[navg*i+j];
1012              	gaps_out[i] = igx;
1013                    }
1014                    else
1015              	gaps_out[i] = 0;
1016                  }
1017                  for (i=0; i<2*(*num_sections); i+=2)
1018                  {
1019                    start = (int)ceilf(((float)sections[i])/navg);
1020                    stop = (int)floorf(((float)sections[i+1])/navg);
1021                    if (start <= stop)
1022                    {
1023              	sections[i] = start;
1024              	sections[i+1] = stop;
1025                    }
1026                    else 
1027                    {
1028              	/* This section was skipped entirely. */
1029              	for (j=i; j< 2*(*num_sections-1); j+=2) 
1030 tplarson 1.1 	{
1031              	  sections[j] = sections[j+2];
1032              	  sections[j+1] = sections[j+3];
1033              	}
1034              	*num_sections -= 1;
1035                    }
1036                  }
1037                }  
1038                return 0;
1039              }
1040              
1041              
1042              /*
1043                 Read a list of sections [start_sample:end_sample] that should
1044                 contain continuous data. When detrending, jumps are allow to
1045                 occur between sections. The sections file should be a text file 
1046                 of the form:
1047              
1048                 num
1049                 start_1  end_1
1050                 start_2  end_2
1051 tplarson 1.1    ...
1052                 start_num  end_num
1053              */
1054              static int read_section_file(char *filename, int *num_sections, int **sections)
1055              {
1056                FILE *fh;
1057                int i;
1058              
1059                fh = fopen(filename,"r");
1060                if (fh!=NULL)
1061                {
1062                  fscanf(fh,"%d",num_sections);
1063                  *sections = (int *)malloc(2*(*num_sections)*sizeof(int));
1064                  
1065                  for (i=0 ;i<2*(*num_sections); i+=2)
1066                  {
1067                    fscanf(fh,"%d",&(*sections)[i]);
1068                    fscanf(fh,"%d",&(*sections)[i+1]);
1069                  
1070                    if ((*sections)[i]>(*sections)[i+1] ||
1071              	  (i>0 && (*sections)[i]<=(*sections)[i-1]))
1072 tplarson 1.1       {
1073              	fprintf(stderr,"Invalid sections file, sections overlap.\n");
1074              	return 2;
1075                    }
1076                  }
1077                  fclose(fh);
1078                  return 0;
1079                }
1080                else 
1081                  return 1;
1082              }  
1083              

Karen Tian
Powered by
ViewCVS 0.9.4