htswanalysis/MACS/lib/gsl/gsl-1.11/bspline/bspline.c

   1 /* bspline/bspline.c
   2  *
   3  * Copyright (C) 2006, 2007, 2008 Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_errno.h>
  22 #include <gsl/gsl_bspline.h>
  23
  24 /*
  25  * This module contains routines related to calculating B-splines.
  26  * The algorithms used are described in
  27  *
  28  * [1] Carl de Boor, "A Practical Guide to Splines", Springer
  29  *     Verlag, 1978.
  30  */
  31
  32 static int bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
  33                             gsl_bspline_workspace *w);
  34
  35 static inline size_t bspline_find_interval(const double x, int *flag,
  36                                            gsl_bspline_workspace *w);
  37
  38 /*
  39 gsl_bspline_alloc()
  40   Allocate space for a bspline workspace. The size of the
  41 workspace is O(5k + nbreak)
  42
  43 Inputs: k      - spline order (cubic = 4)
  44         nbreak - number of breakpoints
  45
  46 Return: pointer to workspace
  47 */
  48
  49 gsl_bspline_workspace *
  50 gsl_bspline_alloc(const size_t k, const size_t nbreak)
  51 {
  52   if (k == 0)
  53     {
  54       GSL_ERROR_NULL("k must be at least 1", GSL_EINVAL);
  55     }
  56   else if (nbreak < 2)
  57     {
  58       GSL_ERROR_NULL("nbreak must be at least 2", GSL_EINVAL);
  59     }
  60   else
  61     {
  62       gsl_bspline_workspace *w;
  63
  64       w = (gsl_bspline_workspace *)
  65           calloc(1, sizeof(gsl_bspline_workspace));
  66       if (w == 0)
  67         {
  68           GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);
  69         }
  70
  71       w->k = k;
  72       w->km1 = k - 1;
  73       w->nbreak = nbreak;
  74       w->l = nbreak - 1;
  75       w->n = w->l + k - 1;
  76
  77       w->knots = gsl_vector_alloc(w->n + k);
  78       if (w->knots == 0)
  79         {
  80           gsl_bspline_free(w);
  81           GSL_ERROR_NULL("failed to allocate space for knots vector", GSL_ENOMEM);
  82         }
  83
  84       w->deltal = gsl_vector_alloc(k);
  85       w->deltar = gsl_vector_alloc(k);
  86       if (!w->deltal || !w->deltar)
  87         {
  88           gsl_bspline_free(w);
  89           GSL_ERROR_NULL("failed to allocate space for delta vectors", GSL_ENOMEM);
  90         }
  91
  92       w->B = gsl_vector_alloc(k);
  93       if (w->B == 0)
  94         {
  95           gsl_bspline_free(w);
  96           GSL_ERROR_NULL("failed to allocate space for temporary spline vector", GSL_ENOMEM);
  97         }
  98
  99       return (w);
 100     }
 101 } /* gsl_bspline_alloc() */
 102
 103 /* Return number of coefficients */
 104 size_t
 105 gsl_bspline_ncoeffs (gsl_bspline_workspace * w)
 106 {
 107   return w->n;
 108 }
 109
 110 /* Return order */
 111 size_t
 112 gsl_bspline_order (gsl_bspline_workspace * w)
 113 {
 114   return w->k;
 115 }
 116
 117 /* Return number of breakpoints */
 118 size_t
 119 gsl_bspline_nbreak (gsl_bspline_workspace * w)
 120 {
 121   return w->nbreak;
 122 }
 123
 124 double
 125 gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
 126 {
 127   size_t j = i + w->k - 1;
 128   return gsl_vector_get(w->knots, j);
 129 }
 130
 131 /*
 132 gsl_bspline_free()
 133   Free a bspline workspace
 134
 135 Inputs: w - workspace to free
 136
 137 Return: none
 138 */
 139
 140 void
 141 gsl_bspline_free(gsl_bspline_workspace *w)
 142 {
 143   if (w->knots)
 144     gsl_vector_free(w->knots);
 145
 146   if (w->deltal)
 147     gsl_vector_free(w->deltal);
 148
 149   if (w->deltar)
 150     gsl_vector_free(w->deltar);
 151
 152   if (w->B)
 153     gsl_vector_free(w->B);
 154
 155   free(w);
 156 } /* gsl_bspline_free() */
 157
 158 /*
 159 gsl_bspline_knots()
 160   Compute the knots from the given breakpoints:
 161
 162 knots(1:k) = breakpts(1)
 163 knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
 164 knots(n+1:n+k) = breakpts(l + 1)
 165
 166 where l is the number of polynomial pieces (l = nbreak - 1) and
 167 n = k + l - 1
 168 (using matlab syntax for the arrays)
 169
 170 The repeated knots at the beginning and end of the interval
 171 correspond to the continuity condition there. See pg. 119
 172 of [1].
 173
 174 Inputs: breakpts - breakpoints
 175         w        - bspline workspace
 176
 177 Return: success or error
 178 */
 179
 180 int
 181 gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w)
 182 {
 183   if (breakpts->size != w->nbreak)
 184     {
 185       GSL_ERROR("breakpts vector has wrong size", GSL_EBADLEN);
 186     }
 187   else
 188     {
 189       size_t i; /* looping */
 190
 191       for (i = 0; i < w->k; ++i)
 192         gsl_vector_set(w->knots, i, gsl_vector_get(breakpts, 0));
 193
 194       for (i = 1; i < w->l; ++i)
 195         {
 196           gsl_vector_set(w->knots, w->k - 1 + i,
 197                          gsl_vector_get(breakpts, i));
 198         }
 199
 200       for (i = w->n; i < w->n + w->k; ++i)
 201         gsl_vector_set(w->knots, i, gsl_vector_get(breakpts, w->l));
 202
 203       return GSL_SUCCESS;
 204     }
 205 } /* gsl_bspline_knots() */
 206
 207 /*
 208 gsl_bspline_knots_uniform()
 209   Construct uniformly spaced knots on the interval [a,b] using
 210 the previously specified number of breakpoints. 'a' is the position
 211 of the first breakpoint and 'b' is the position of the last
 212 breakpoint.
 213
 214 Inputs: a - left side of interval
 215         b - right side of interval
 216         w - bspline workspace
 217
 218 Return: success or error
 219
 220 Notes: 1) w->knots is modified to contain the uniformly spaced
 221           knots
 222
 223        2) The knots vector is set up as follows (using octave syntax):
 224
 225           knots(1:k) = a
 226           knots(k+1:k+l-1) = a + i*delta, i = 1 .. l - 1
 227           knots(n+1:n+k) = b
 228 */
 229
 230 int
 231 gsl_bspline_knots_uniform(const double a, const double b,
 232                           gsl_bspline_workspace *w)
 233 {
 234   size_t i;     /* looping */
 235   double delta; /* interval spacing */
 236   double x;
 237
 238   delta = (b - a) / (double) w->l;
 239
 240   for (i = 0; i < w->k; ++i)
 241     gsl_vector_set(w->knots, i, a);
 242
 243   x = a + delta;
 244   for (i = 0; i < w->l - 1; ++i)
 245     {
 246       gsl_vector_set(w->knots, w->k + i, x);
 247       x += delta;
 248     }
 249
 250   for (i = w->n; i < w->n + w->k; ++i)
 251     gsl_vector_set(w->knots, i, b);
 252
 253   return GSL_SUCCESS;
 254 } /* gsl_bspline_knots_uniform() */
 255
 256 /*
 257 gsl_bspline_eval()
 258   Evaluate the basis functions B_i(x) for all i. This is
 259 basically a wrapper function for bspline_eval_all()
 260 which formats the output in a nice way.
 261
 262 Inputs: x - point for evaluation
 263         B - (output) where to store B_i(x) values
 264             the length of this vector is
 265             n = nbreak + k - 2 = l + k - 1 = w->n
 266         w - bspline workspace
 267
 268 Return: success or error
 269
 270 Notes: The w->knots vector must be initialized prior to calling
 271        this function (see gsl_bspline_knots())
 272 */
 273
 274 int
 275 gsl_bspline_eval(const double x, gsl_vector *B,
 276                  gsl_bspline_workspace *w)
 277 {
 278   if (B->size != w->n)
 279     {
 280       GSL_ERROR("B vector length does not match n", GSL_EBADLEN);
 281     }
 282   else
 283     {
 284       size_t i;     /* looping */
 285       size_t idx = 0;
 286       size_t start; /* first non-zero spline */
 287
 288       /* find all non-zero B_i(x) values */
 289       bspline_eval_all(x, w->B, &idx, w);
 290
 291       /* store values in appropriate part of given vector */
 292
 293       start = idx - w->k + 1;
 294       for (i = 0; i < start; ++i)
 295         gsl_vector_set(B, i, 0.0);
 296
 297       for (i = start; i <= idx; ++i)
 298         gsl_vector_set(B, i, gsl_vector_get(w->B, i - start));
 299
 300       for (i = idx + 1; i < w->n; ++i)
 301         gsl_vector_set(B, i, 0.0);
 302
 303       return GSL_SUCCESS;
 304     }
 305 } /* gsl_bspline_eval() */
 306
 307 /****************************************
 308  *          INTERNAL ROUTINES           *
 309  ****************************************/
 310
 311 /*
 312 bspline_eval_all()
 313   Evaluate all non-zero B-splines B_i(x) using algorithm (8)
 314 of chapter X of [1]
 315
 316 The idea is something like this. Given x \in [t_i, t_{i+1}]
 317 with t_i < t_{i+1} and the t_i are knots, the values of the
 318 B-splines not automatically zero fit into a triangular
 319 array as follows:
 320                          0
 321             0
 322 0                        B_{i-2,3}
 323             B_{i-1,2}
 324 B_{i,1}                  B_{i-1,3}       ...
 325             B_{i,2}
 326 0                        B_{i,3}
 327             0
 328                          0
 329
 330 where B_{i,k} is the ith B-spline of order k. The boundary 0s
 331 indicate that those B-splines not in the table vanish at x.
 332
 333 To compute the non-zero B-splines of a given order k, we use
 334 Eqs. (4) and (5) of chapter X of [1]:
 335
 336 (4) B_{i,1}(x) = { 1, t_i <= x < t_{i+1}
 337                    0, else }
 338
 339 (5) B_{i,k}(x) =     (x - t_i)
 340                  ----------------- B_{i,k-1}(x)
 341                  (t_{i+k-1} - t_i)
 342
 343                       t_{i+k} - x
 344                  + ----------------- B_{i+1,k-1}(x)
 345                    t_{i+k} - t_{i+1}
 346
 347 So (4) gives us the first column of the table and we can use
 348 the recurrence relation (5) to get the rest of the columns.
 349
 350 Inputs: x   - point at which to evaluate splines
 351         B   - (output) where to store B-spline values (length k)
 352         idx - (output) B-spline function index of last output
 353               value (B_{idx}(x) is stored in the last slot of 'B')
 354         w   - bspline workspace
 355
 356 Return: success or error
 357
 358 Notes: 1) the w->knots vector must be initialized before calling
 359           this function
 360
 361        2) On output, B contains:
 362
 363        B = [B_{i-k+1,k}, B_{i-k+2,k}, ..., B_{i-1,k}, B_{i,k}]
 364
 365           where 'i' is stored in 'idx' on output
 366
 367        3) based on PPPACK bsplvb
 368 */
 369
 370 static int
 371 bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
 372                  gsl_bspline_workspace *w)
 373 {
 374   if (B->size != w->k)
 375     {
 376       GSL_ERROR("B vector not of length k", GSL_EBADLEN);
 377     }
 378   else
 379     {
 380       size_t i;  /* spline index */
 381       size_t j;  /* looping */
 382       size_t ii; /* looping */
 383       int flag = 0;  /* error flag */
 384       double saved;
 385       double term;
 386
 387       i = bspline_find_interval(x, &flag, w);
 388
 389       if (flag == -1)
 390         {
 391           GSL_ERROR("x outside of knot interval", GSL_EINVAL);
 392         }
 393       else if (flag == 1)
 394         {
 395           if (x <= gsl_vector_get(w->knots, i) + GSL_DBL_EPSILON)
 396             {
 397               --i;
 398             }
 399           else
 400             {
 401               GSL_ERROR("x outside of knot interval", GSL_EINVAL);
 402             }
 403         }
 404
 405       if (gsl_vector_get(w->knots, i) == gsl_vector_get(w->knots, i + 1))
 406         {
 407           GSL_ERROR("knot(i) = knot(i+1) will result in division by zero",
 408                     GSL_EINVAL);
 409         }
 410
 411       *idx = i;
 412
 413       gsl_vector_set(B, 0, 1.0);
 414
 415       for (j = 0; j < w->k - 1; ++j)
 416         {
 417           gsl_vector_set(w->deltar, j,
 418                          gsl_vector_get(w->knots, i + j + 1) - x);
 419           gsl_vector_set(w->deltal, j,
 420                          x - gsl_vector_get(w->knots, i - j));
 421
 422           saved = 0.0;
 423
 424           for (ii = 0; ii <= j; ++ii)
 425             {
 426               term = gsl_vector_get(B, ii) /
 427                      (gsl_vector_get(w->deltar, ii) +
 428                       gsl_vector_get(w->deltal, j - ii));
 429
 430               gsl_vector_set(B, ii,
 431                              saved +
 432                              gsl_vector_get(w->deltar, ii) * term);
 433
 434               saved = gsl_vector_get(w->deltal, j - ii) * term;
 435             }
 436
 437           gsl_vector_set(B, j + 1, saved);
 438         }
 439
 440       return GSL_SUCCESS;
 441     }
 442 } /* bspline_eval_all() */
 443
 444 /*
 445 bspline_find_interval()
 446   Find knot interval such that
 447
 448   t_i <= x < t_{i + 1}
 449
 450 where the t_i are knot values.
 451
 452 Inputs: x    - x value
 453         flag - (output) error flag
 454         w    - bspline workspace
 455
 456 Return: i (index in w->knots corresponding to left limit of interval)
 457
 458 Notes: The error conditions are reported as follows:
 459
 460 Condition                        Return value        Flag
 461 ---------                        ------------        ----
 462 x < t_0                               0               -1
 463 t_i <= x < t_{i+1}                    i                0
 464 t_i < x = t_{i+1} = t_{n+k-1}         i                0
 465 t_{n+k-1} < x                       l+k-1             +1
 466 */
 467
 468 static inline size_t
 469 bspline_find_interval(const double x, int *flag,
 470                       gsl_bspline_workspace *w)
 471 {
 472   size_t i;
 473
 474   if (x < gsl_vector_get(w->knots, 0))
 475     {
 476       *flag = -1;
 477       return 0;
 478     }
 479
 480   /* find i such that t_i <= x < t_{i+1} */
 481   for (i = w->k - 1; i < w->k + w->l - 1; ++i)
 482     {
 483       const double ti = gsl_vector_get(w->knots, i);
 484       const double tip1 = gsl_vector_get(w->knots, i + 1);
 485
 486       if (tip1 < ti)
 487         {
 488           GSL_ERROR("knots vector is not increasing", GSL_EINVAL);
 489         }
 490
 491       if (ti <= x && x < tip1)
 492         break;
 493
 494       if (ti < x && x == tip1 &&
 495           tip1 == gsl_vector_get(w->knots, w->k + w->l - 1))
 496         break;
 497     }
 498
 499   if (i == w->k + w->l - 1)
 500     *flag = 1;
 501   else
 502     *flag = 0;
 503
 504   return i;
 505 } /* bspline_find_interval() */