fill.c

/* movstat/fill.c
 * 
 * Copyright (C) 2018 Patrick Alken
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.h>

/*
gsl_movstat_fill()
  Fill window for sample 'idx' from x using given end conditions

Inputs: endtype     - how to handle end points
        x           - input vector, length n
        idx         - index of center sample in window in [0,n-1]
        H           - number of samples left of center to include
        J           - number of samples right of center to include
        window      - (output) window of samples centered on x_{idx},
                      W_{idx}^{H,J}, length H + J + 1

Return: size of filled window (<= H + J + 1)
*/

size_t
gsl_movstat_fill(const gsl_movstat_end_t endtype, const gsl_vector * x, const size_t idx,
                 const size_t H, const size_t J, double * window)
{
  if (idx >= x->size)
    {
      GSL_ERROR_VAL ("window center index must be between 0 and n - 1", GSL_EDOM, 0);
    }
  else
    {
      const int n = (int) x->size;
      const int iidx = (int) idx;
      const int iH = (int) H;
      const int iJ = (int) J;
      int idx1, idx2, j;
      size_t window_size;

      if (endtype == GSL_MOVSTAT_END_TRUNCATE)
        {
          idx1 = GSL_MAX(iidx - iH, 0);
          idx2 = GSL_MIN(iidx + iJ, n - 1);
        }
      else
        {
          idx1 = iidx - iH;
          idx2 = iidx + iJ;
        }

      window_size = (size_t) (idx2 - idx1 + 1);

      /* fill sliding window */
      for (j = idx1; j <= idx2; ++j)
        {
          int widx = j - idx1;

          if (j < 0)
            {
              /* initial condition */
              if (endtype == GSL_MOVSTAT_END_PADZERO)
                window[widx] = 0.0;
              else if (endtype == GSL_MOVSTAT_END_PADVALUE)
                window[widx] = gsl_vector_get(x, 0);
            }
          else if (j >= n)
            {
              if (endtype == GSL_MOVSTAT_END_PADZERO)
                window[widx] = 0.0;
              else if (endtype == GSL_MOVSTAT_END_PADVALUE)
                window[widx] = gsl_vector_get(x, n - 1);
            }
          else
            {
              window[widx] = gsl_vector_get(x, j);
            }
        }

      return window_size;
    }
}