//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** @author John Miller, Michael Cotterell, Hao Peng, Dong-Yu Yu * @version 1.6 * @date Thu Sep 22 21:45:58 EDT 2016 * @see LICENSE (MIT style license file). * * @see en.wikipedia.org/wiki/B-spline * @see cran.r-project.org/web/packages/crs/vignettes/spline_primer.pdf * @see http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node17.html * @see open.uct.ac.za/bitstream/item/16664/thesis_sci_2015_essomba_rene_franck.pdf?sequence=1 */ package scalation.calculus import scala.math.sqrt import scalation.linalgebra.{MatrixD, VectoD, VectorD} import scalation.math.ExtremeD.EPSILON import scalation.math.double_exp import scalation.plot.{Plot, PlotM} import scalation.stat.Statistic import scalation.util.{Error, timed} //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `B_Spline` class provides B-Spline basis functions for various orders 'm', * where the order is one more than the degree. A spline function is a piecewise * polynomial function where the pieces are stitched together at knots with the * goal of maintaining continuity and differentiability. B-Spline basis functions * form a popular form of basis functions used in Functional Data Analysis. * @see http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node17.html *----------------------------------------------------------------------------- * @param ττ the time-points of the original knots in the time dimension * @param mMax the maximum order, allowing splines orders from 1 to mMax * @param clamp flag for augmenting ττ */ class B_Spline (ττ: VectoD, mMax: Int = 4, clamp: Boolean = true) extends BasisFunction with Error { protected val DEBUG = false // debug flag protected val l = ττ.dim - 1 // the number of intervals protected val τ = if (clamp) B_Spline.clamp (mMax, ττ, true) // augment to accommodate knots else ττ protected val head = τ(0) protected val tail = τ.last if (mMax < 1 || mMax > 10) flaw ("constructor", "B_Spline order restricted to 1 thru 10") if (DEBUG) { println (s"B_Spline (ττ = $ττ, mMax = $mMax)") println ("intervals l = " + l) println ("augmented τ = " + τ) } // if //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Range of "usable" splines when using the `bs` function. This is needed, * since extra splines may be generated by the general B-spline recurrence. * @param m the order of the spline */ def range (m: Int = mMax) = 0 to l + m - 2 //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The number of usable spline basis functions for a specified order, * given the configured knot vector. * @param m the order of the spline */ def size (m: Int = mMax) = l + m - 1 //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Order one 'm = 1' B-Spline basis functions (flat functions). * @param j indicates which spline function 0 to l-1 * @param t the time parameter */ private def bb1 (j: Int, t: Double): Double = { if (j == size()-1 && t =~ tail) 1.0 else if ((τ(j) < t || τ(j) =~ t) && (t < τ(j+1))) 1.0 else 0.0 } // bb1 //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Order 'm' B-Spline basis functions (general recurrence). * @param m the order of the spline function (degree = order - 1) * @param j indicates which spline function * @param t the time parameter */ def bb (m: Int)(j: Int)(t: Double): Double = { if (mMax < m) flaw ("bb", s"mMax = $mMax can't be less than m = $m") if (m == 1) return bb1 (j, t) val f1 = bb (m-1)(j)(t) val f2 = bb (m-1)(j+1)(t) val a = if (f1 =~ 0) 0 else (t - τ(j)) * f1 / (τ(j+m-1) - τ(j)) val b = if (f2 =~ 0) 0 else (τ(j+m) - t) * f2 / (τ(j+m) - τ(j+1)) a + b } // bb //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Adjusted order 'm' B-Spline basis functions (general recurrence). These * are adjusted so that the first "usable" spline is at `j = 0`. The valid * range of usable splines is defined in `range`. * @param m the order of the spline function (degree = order - 1) * @param j indicates which spline function * @param t the time parameter */ def bfr (m: Int)(j: Int)(t: Double): Double = bb (m)(j+mMax-m)(t) //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Evaluate the 'j+1'-th order `m` B-spline basis function * at `t` non-recursively, using an efficient dynamic programming approach. * ==Example== * {{{ * val m = 4 // order m (degree m-1) * val t = VectorD (1, 2, 3, 4) // original time points * val N = B_Spline (t, k) // B_Spline instance * val j = 0 // spline j in N.range(m) * val z = N.bf (m)(j)(t(2)) // evaluate at t(2) * println ("order k basis function j at t(2) = z") * }}} * ==Implementation Details== * Let `N(m)(i)` denote the `i`-th order `k` basis function evaluated at * `t`. Then each `N(m)(i)` depends on the evaluation of `N(m-1)(i)` and * `N(m-1)(i+1)`. Here is an example, given a set of order `m=4` B-spline * basis functions and knot vector `τ` of length `n`: * {{{ * N(1)(0), N(1)(1), N(1)(2), N(1)(3) * |/ |/ |/ * N(2)(0), N(2)(1), N(2)(2) * |/ |/ * N(3)(0), N(3)(1) * |/ * N(4)(0) * }}} * This algorithm applies the procedure described above using O(k) * storage and O(k*k) floating point operations. * @param t point to evaluate * @author Michael Cotterell, Hao Peng * @see Carl de Boor (1978). A Practical Guide to Splines. Springer-Verlag. ISBN 3-540-90356-9. */ def bf (m: Int)(j: Int)(t: Double): Double = { val N = Array.ofDim [Double] (m) // all evaluations; overwrite as needed for (k <- 0 until m) { // build for each order 1 <= (k+1) <= m val nf = m - k // number of basis functions for order (k+1) if (j == 0 && t =~ head) N(0) = 1.0 // trivial case: near first knot else if (j == size(m)-1 && t =~ tail) N(0) = 1.0 // trivial case: near last knot else if (k == 0) { // evaluate order 1 for (i <- 0 until nf) { N(i) = if ((τ(i+j) < t || τ(i+j) =~ t) && (t < τ(i+j+1))) 1.0 else 0.0 } // for } else if (k > 0) { // evaluate orders (k+1) <= k for (i <- 0 until nf) { val bs1 = N(i) val bs2 = N(i+1) // using "== 0" instead of "=~ 0" should be safe val bf1 = if (bs1 == 0) 0 else (t - τ(i+j)) * bs1 / (τ(i+j+k) - τ(i+j)) val bf2 = if (bs2 == 0) 0 else (τ(i+j+k+1) - t) * bs2 / (τ(i+j+k+1) - τ(i+j+1)) N(i) = bf1 + bf2 } // for } // if } // for N(0) } // bf //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Evaluate each of the `τ.dim-1`-many order `m` B-spline basis functions * at `t` non-recursively, using an efficient dynamic programming approach. * ==Example== * {{{ * val m = 4 // order m (degree m-1) * val t = VectorD (1, 2, 3, 4) // original time points * val N = B_Spline (t, k) // B_Spline instance * val z = N.abf_ (m)(t) // vector of evaluations * println ("order k basis functions at t = z") * }}} * ==Implementation Details== * Let `N(m)(i)` denote the `i`-th order `k` basis function evaluated at * `t`. Then each `N(m)(i)` depends on the evaluation of `N(m-1)(i)` and * `N(m-1)(i+1)`. Here is an example, given a set of order `m=4` B-spline * basis functions and knot vector `τ` of length `n`: * {{{ * N(1)(0), N(1)(1), N(1)(2), N(1)(3), ..., N(1)(n-1) * |/ |/ |/ |/ * N(2)(0), N(2)(1), N(2)(2), ..., N(2)(n-2) * |/ |/ |/ * N(3)(0), N(3)(1), ..., N(3)(n-3) * |/ |/ * N(4)(0), ..., N(4)(n-4) * }}} * This algorithm applies the procedure described above using O(n) * storage and O(k*n) floating point operations. * @param t point to evaluate * @author Michael Cotterell * @see Carl de Boor (1978). A Practical Guide to Splines. Springer-Verlag. ISBN 3-540-90356-9. */ override def abf_ (m: Int = mMax)(t: Double): VectorD = { val N = Array.ofDim [Double] (τ.dim-1) // all evaluations; overwrite as needed for (k <- 0 until m) { // build for each order 1 <= (k+1) <= m val nf = τ.dim - k - 1 // number of basis functions for order (k+1) if (t =~ head) N( 0) = 1.0 // trivial case: near first knot else if (t =~ tail) N(nf-1) = 1.0 // trivial case: near last knot else if (k == 0) { // evaluate order 1 for (i <- 0 until nf-1) { N(i) = if ((τ(i) < t || τ(i) =~ t) && (t < τ(i+1))) 1.0 else 0.0 } // for } else if (k > 0) { // evaluate orders (k+1) <= k for (i <- 0 until nf) { val bs1 = N(i) val bs2 = N(i+1) // using "== 0" instead of "=~ 0" should be safe val bf1 = if (bs1 == 0) 0 else (t - τ(i)) * bs1 / (τ(i+k) - τ(i)) val bf2 = if (bs2 == 0) 0 else (τ(i+k+1) - t) * bs2 / (τ(i+k+1) - τ(i+1)) N(i) = bf1 + bf2 } // for } // if } // for new VectorD (τ.dim - m, N) // do not allocate new storage } // abf_ //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Retrieve the order of the this B_Spline. */ def getOrder = mMax } // B_Spline class //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Companion object for `B_Spline` class provides functions for clamping the * the ends of a spline and running timing benchmarks. */ object B_Spline { //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return a "clamped" version of the input vector, augmented to ensure * that each end has `m`-many repeated points. If `isInclusive` is true, * then the first and last values of the vector are repeated; otherwise, * the values `t(0)-sqrt(EPSILON)` and `t(t.dim-1)+sqrt(EPSILON)` are * repeated for the beginning and end, respectively. * @param m intended B-spline order (degree = m - 1) * @param t non-decreasing vector of time points * @param isInclusive repeat end points (default = true) * @author Michael Cotterell */ def clamp (m: Int, t: VectoD, isInclusive: Boolean = true): VectorD = { val e = sqrt (EPSILON) val head = t(0) val tail = t.last if (isInclusive) VectorD.fill (m-1)(head-e) ++ t ++ VectorD.fill (m-1)(tail+e) else VectorD.fill (m)(head-e) ++ t ++ VectorD.fill (m)(tail+e) } // clamp //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Gather timing statistics for repeated executions of a block, prints the * statistics, and returns the value produced by the block. * @tparam R result type of block * @param reps number of replications * @param useSeconds record time as seconds instead of milliseconds (default = false) * @param title title for statistic * @param block block of code to execute */ def benchmark [R] (reps: Int = 100, useSeconds: Boolean = false, title: String = "benchmark") (block: => R): R = { val (r, stat) = benchmarked (reps, useSeconds, title) (block) println(Statistic.labels) println(stat) r } // benchmark //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Gather timing statistics for repeated executions of a block and returns * a tuple containing the return value produced by the block and the * `Statistic` instance used to gather the statistics. * @tparam R result type of block * @param reps number of replications * @param useSeconds record time as seconds instead of milliseconds (default = false) * @param title title for statistic * @param block block of code to execute */ def benchmarked [R] (reps: Int = 100, useSeconds: Boolean = false, title: String = "benchmark") (block: => R): (R, Statistic) = { val stat = new Statistic (title) for (i <- 0 until reps-1) { val (r, ms) = timed (block) stat.tally (if (useSeconds) ms / 1000 else ms) } // for val (r, ms) = timed(block) stat.tally (if (useSeconds) ms / 1000 else ms) (r, stat) } // benchmarked } // B_Spline object //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `B_SplineTest` object is used to test the `B_Spline` class. * It tests the B-Spline functions using the general recurrence. * > runMain scalation.calculus.B_SplineTest */ object B_SplineTest extends App { val mM = 4 // maximum order to test val τ = VectorD (0.0, 20.0, 40.0, 60.0, 80.0, 100.0) // knot time-points val bs = new B_Spline (τ, mM) // B-Spline generator val n = 100 // number of time-points for plotting val t = VectorD.range (0, n) // time-points for plotting for (m <- 1 to mM) { //--------------------------------------------------------------------- // order m B-Splines (polynomial functions) val y1 = new VectorD (n) val y2 = new VectorD (n) for (i <- 0 until n) { y1(i) = bs.bf (m)(1)(t(i)) // first "interesting" B-Spline y2(i) = bs.bf (m)(2)(t(i)) // next B-Spline } // for new Plot (t, y1, y2, "B-Spline order " + m) } // for } // B_SplineTest object //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `B_SplineTest2` object is used to test the `B_Spline` class. * It tests the B-Spline functions using the general recurrence and plots * several basis functions using `PlotM`. * > runMain scalation.calculus.B_SplineTest2 */ object B_SplineTest2 extends App { val mM = 5 // maximum order to test val τ = VectorD (0.0, 50.0, 100.0) // knot time-points val bs = new B_Spline (τ, mM) // B-Spline generator val n = 100 // number of time-points for plotting val t = VectorD.range (0, n)/(n-1)*n // time-points for plotting val k = 0 // index for initial B-Spline for (m <- 2 to 4) { //--------------------------------------------------------------------- // order m B-Splines (polynomial functions) val y = new MatrixD (τ.dim + mM, n) // matrix to hold initial B-Splines for (i <- 0 until n; j <- 0 to bs.range(m).last) { y(j, i) = bs (m)(j)(t(i)) } // for new PlotM (t, y, null, "B-Spline order " + m) } // for } // B_SplineTest2 object //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `B_SplineTest3` object is used to test the `B_Spline` class. * It tests the B-Spline functions using the general recurrence. * > runMain scalation.calculus.B_SplineTest3 */ object B_SplineTest3 extends App { import scalation.calculus.BasisFunction.plot val mM = 4 // maximum order to test val n = 100 // number of time points val τ = VectorD (0.0, 0.5 * (n-1), (n-1)) / (n) // knot vector (unaugmented) val bs = new B_Spline (τ, mM) // B-Spline generator val t = VectorD.range (0, n) / n // time-points for plotting for (m <- 1 to mM) plot (bs, m, t) // plot } // B_SplineTest3 object //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `B_SplineTest4` object is used to test the `B_Spline` class. * It tests the B-Spline functions using the general recurrence. * > runMain scalation.calculus.B_SplineTest4 */ object B_SplineTest4 extends App { val mM = 4 // maximum order to test val t = VectorD.range (1, 11) / 10 // time-points for plotting val bs = new B_Spline (t, mM) // B-Spline generator for (j <- bs.range(mM)) { println (s"j = $j") val m1 = VectorD (for (i <- t.indices) yield bs.bf (mM)(j)(t(i))) val m2 = VectorD (for (i <- t.indices) yield bs.bfr(mM)(j)(t(i))) println (s"m1 = $m1") println (s"m2 = $m2") val diff = m1 - m2 val e2 = diff dot diff println (s"e2 = $e2") } // for } // B_SplineTest4 object