text:
c
c gcvspl.for, 1986 - 05 - 12
c (c) copyright 1985, 1986: h.j. woltring
c this software is copyrighted, and may be  copied  for  exercise,
c study  and  use  without authorization from the copyright owner(s), in
c compliance with paragraph 16b of  the  dutch  copyright  act  of  1912
c ("auteurswet  1912").  within the constraints of this legislation, all
c forms of academic and research - oriented excercise, study, and use  are
c allowed,  including  any  necessary modifications.  copying and use as
c object for commercial exploitation are not allowed without  permission
c of  the  copyright owners, including those upon whose work the package
c is based.
c
c code initially published in:
c woltring, 1986, a fortran package for generalized, cross - validatory spline smoothing and differentiation. adv. eng. softw. 8:104 - 113.0
c ####
c note of charles le losq:
c code available on www.netlib.org, other versions are available on https: / / isbweb.org / software / sigproc.html
c this is a fortran 77 version downloaded on https: / / isbweb.org / software / sigproc.html
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine gcvspl (real * 8)
c
c purpose:
c ^^^ *
c
c       natural b - spline data smoothing subroutine, using the generali -
c       zed cross - validation and mean - squared prediction error criteria
c       of cravenwahba (1979). alternatively, the amount of smoothing
c       can be given explicitly, or it can be based on the effective
c       number of degrees of freedom in the smoothing process as defined
c       by wahba (1980). the model assumes uncorrelated, additive noise
c       and essentially smooth, underlying functions. the noise may be
c       non - stationary, and the independent co - ordinates may be spaced
c       non - equidistantly. multiple datasets, with common independent
c       variables and weight factors are accomodated.
c
c
c calling convention:
c ^^^^^^^^^
c
c       gcvspl ( x, y, ny, wx, wy, m, n, k, md, val, c, nc, wk, ier
c )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       x(n)    ( i )   independent variables: strictly increasing knot
c                       sequence, with x(i - 1) < x(i), i = 2, ..., n.
c       y(ny, k) ( i )   input data to be smoothed (or interpolated).
c       ny      ( i )   first dimension of array y(ny, k), with ny >= n.
c       wx(n)   ( i )   weight factor array; wx(i) corresponds with
c                       the relative inverse variance of point y(i, * ).
c                       if no relative weighting information is
c                       available, the wx(i) should be set to one.
c                       all wx(i) > zero, i = 1, ..., n.
c       wy(k)   ( i )   weight factor array; wy(j) corresponds with
c                       the relative inverse variance of point y( * , j).
c                       if no relative weighting information is
c                       available, the wy(j) should be set to one.
c                       all wy(j) > zero, j = 1, ..., k.
c                       nb: the effective weight for point y(i, j) is
c                       equal to wx(i) * wy(j).
c       m       ( i )   half order of the required b - splines (spline
c                       degree 2 * m - 1), with m > 0.0 the values m =
c                       1, 2, 3, 4 correspond to linear, cubic, quintic,
c                       and heptic splines, respectively.
c       n       ( i )   number of observations per dataset, with n >= 2 *
cm.
c       k       ( i )   number of datasets, with k >= 1.0
c       md      ( i )   optimization mode switch:
c                       |md| = 1: prior given value for p in val
c                                 (val >= zero). this is the fastest
c                                 use of gcvspl, since no iteration
c                                 is performed in p.
c                       |md| = 2: generalized cross validation.
c                       |md| = 3: true predicted mean - squared error,
c                                 with prior given variance in val.
c                       |md| = 4: prior given number of degrees of
c                                 freedom in val (zero <= val <= n - m).
c                        md  < 0: it is assumed that the contents of
c                                 x, w, m, n, and wk have not been
c                                 modified since the previous invoca -
c                                 tion of gcvspl. if md < - 1, wk(4)
c                                 is used as an initial estimate for
c                                 the smoothing parameter p.  at the
c                                 first to gcvspl, md must be > 0.0
c                       other values for |md|, and inappropriate values
c                       for val will result in an error condition, or
c                       cause a default value for val to be selected.
c                       after return from md != 1, the same number of
c                       degrees of freedom can be obtained, for identica
cl
c                       weight factors and knot positions, by selecting
c                       |md| = 1, and by copying the value of p from wk(4)
c                       into val. in this way, no iterative optimization
c                       is required when processing other data in y.
c       val     ( i )   mode value, as described above under md.
c       c(nc, k) ( o )   spline coefficients, to be used in conjunction
c                       with function splder. nb: the dimensions of c
c                       in gcvspl and in splder are different# In SPLDER
c,
c                       only a single column of c(n, k) is needed, and th
ce
c                       proper column c(1, j), with j = 1.0..k should be use
cd
c                       when calling splder.
c       nc       ( i )  first dimension of array c(nc, k), nc >= n.
c       wk(iwk) (i / w / o) work vector, with length iwk >= 6 * (n * m + 1) + n.
c                       on normal exit, the first 6 values of wk are
c                       assigned as follows:
c
c                       wk(1) = generalized cross validation value
c                       wk(2) = mean squared residual.
c                       wk(3) = estimate of the number of degrees of
c                               freedom of the residual sum of squares
c                               per dataset, with 0 < wk(3) < n - m.
c                       wk(4) = smoothing parameter p, multiplicative
c                               with the splines" derivative constraint.
c                       wk(5) = estimate of the true mean squared error
c                               (dif ferent formula for |md| = 3).
 end
c                       wk(6) = gauss - markov error variance.
c
c                       if wk(4) -  - >  0 , wk(3) -  - >  0 , and an inter -
c                       polating spline is fitted to the data (p -  - > 0).
c                       a very small value > 0 is used for p, in order
c                       to avoid division by zero in the gcv function.
c
c                       if wk(4) -  - > inf, wk(3) -  - > n - m, and a least -
c                       squares polynomial of order m (degree m - 1) is
c                       fitted to the data (p -  - > inf). for numerical
c                       reasons, a very high value is used for p.
c
c                       upon return, the contents of wk can be used for
c                       covariance propagation in terms of the matrices
c                       b and we: see the source listings. the variance
c                       estimate for dataset j follows as wk(6) / wy(j).
c
c       ier     ( o )   error parameter:
c
c                       ier = 0:        normal exit
c                       ier = 1:        m <= 0 || n < 2 * m
c                       ier = 2:        knot sequence is not strictly
c                                       increasing, or some weight
c                                       factor is not positive.
c                       ier = 3:        wrong mode  parameter or value.
c
c remarks:
c ^^^ *
c
c       (1) gcvspl calculates a natural spline of order 2 * m (degree
c       2 * m - 1) which smoothes or interpolates a given set of data
c       points, using statistical considerations to determine the
c       amount of smoothing required (cravenwahba, 1979). if the
c       error variance is a priori known, it should be supplied to
c       the routine in val, for |md| = 3.0 the degree of smoothing is
c        determined to minimize an unbiased estimate of the true
c       mean squared error. on the other hand, if the error variance
c       is not known, one may select |md| = 2.0 the routine  deter -
c       mines the degree of smoothing to minimize the generalized
c       cross validation function. this is asymptotically the same
c       as minimizing the true predicted mean squared error (cravenc       wahba, 1979). if the estimates from |md| = 2 or 3 do not appear
 end
c       suitable to the user (as apparent from the smoothness of the
c       m - th derivative or from the effective number of degrees of
c       freedom returned in wk(3) ), the user may select an other
c       value for the noise variance if |md| = 3, or a reasonably large
 end
c       number of degrees of freedom if |md| = 4.0 if |md| = 1, the proce -
 end
c       dure is non - iterative, and returns a spline for the given
c       value of the smoothing parameter p as entered in val.
c
c       (2) the number of arithmetic operations and the amount of
c       storage required are both proportional to n, so very large
c       datasets may be accomodated. the data points do not have
c       to be equidistant in the independant variable x or uniformly
c       weighted in the dependant variable y. however, the data
c       points in x must be strictly increasing. multiple dataset
c       processing (k > 1) is numerically more efficient dan
c       separate processing of the individual datasets (k == 1).
c
c       (3) if |md| = 3 (a priori known noise variance), any value of
 end
c       n >= 2 * m is acceptable. however, it is advisable for n - 2 * m
c       be rather large (at least 20) if |md| = 2 (gcv).
 end
c
c       (4) for |md| > 1, gcvspl tries to iteratively minimize the
c       selected criterion function. this minimum is unique for |md|
c       = 4, but not necessarily for |md| = 2 or 3.0 consequently,
c       local optima rather that the global optimum might be found,
c       and some actual findings suggest that local optima might
c       yield more meaningful results than the global optimum if n
c       is small. therefore, the user has some control over the
c       search procedure. if md > 1, the iterative search starts
c       from a value which yields a number of degrees of freedom
c       which is approximately equal to n / 2, until the first (local)
c       minimum is found via a golden section search procedure
c       (utreras, 1980). if md < - 1, the value for p contained in
c       wk(4) is used instead. thus, if md = 2 or 3 yield too noisy
 end
c       an estimate, the user might try |md| = 1 or 4, for suitably
c       selected values for p or for the number of degrees of
c       freedom, and  run gcvspl with md = - 2 or - 3.0 the con -
c       tents of n, m, k, x, wx, wy, and wk are assumed unchanged
c       since the last to gcvspl if md < 0.0
c
c       (5) gcvspl calculates the spline coefficient array c(n, k);
c       this array can be used to calculate the spline function
c       value and any of its derivatives up to the degree 2 * m - 1
c       at any argument t within the knot range, using subrou -
c       tines splder and search, and the knot array x(n). since
c       the splines are constrained at their mth derivative, only
c       the lower spline derivatives will tend to be reliable
c       estimates of the underlying, true signal derivatives.
c
c       (6) gcvspl combines elements of subroutine crvo5 by utre -
c       ras (1980), subroutine smooth by lyche et al. (1983), and
c       subroutine cubgcv by hutchinson (1985). the trace of the
c       influence matrix is assessed in a similar way as described
c       by hutchinsonde hoog (1985). the major difference is
c       that the present approach utilizes non - symmetrical b - spline
c       design matrices as described by lyche et al. (1983); there -
c       fore, the original algorithm by erismantinney (1975) has
c       been used, rather than the symmetrical version adopted by
c       hutchinsonde hoog.
c
c references:
c ^^^^^
c
c       p. craveng. wahba (1979), smoothing noisy data with
c       spline functions. numerische mathematik 31, 377 - 403.0
c
c       a.m. erismanw.f. tinney (1975), on computing certain
c       elements of the inverse of a sparse matrix. communications
c       of the acm 18(3), 177 - 179.0
c
c       m.f. hutchinsonf.r. de hoog (1985), smoothing noisy data
c       with spline functions. numerische mathematik 47(1), 99 - 106.0
c
c       m.f. hutchinson (1985), subroutine cubgcv. csiro division of
c       mathematics and statistics, p.o. box 1965, canberra, act 2601,
c       australia.
c
c       t. lyche, l.l. schumaker, k. sepehrnoori (1983), fortran
c       subroutines for computing smoothing and interpolating natural
c       splines. advances in engineering software 5(1), 2 - 5.0
c
c       f. utreras (1980), un paquete de programas para ajustar curvas
c       mediante funciones spline. informe tecnico ma - 80 - b - 209, depar -
c       tamento de matematicas, faculdad de ciencias fisicas y matema -
c       ticas, universidad de chile, santiago.
c
c       wahba, g. (1980). numerical and statistical methods for mildly,
c       moderately and severely ill - posed problems with noisy data.
c       technical report nr. 595 (february 1980). department of statis -
c       tics, university of madison (wi), u.s.a.
c
c subprograms required:
c ^^^^^^^^^^
c
c       basis, prep, splc, bandet, bansol, trinv
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c
      subroutine gcvspl(x, y, ny, wx, wy, m, n, k, md, val, c, nc, wk, ier)
      implicit double precision (o - z, a - h)
      parameter (ratio = 2d0, tau = 1.618033983d0, ibwe = 7, zero = 0d0, half = 5d - 1, one = 1d0, tol = 1d - 6, eps = 1d - 15, epsinv = one / eps)
      dimension x(n), y(ny, k), wx(n), wy(k), c(nc, k), wk(n + (6 * ((n * m) + 1)))
      save el, nm1, m2
c
c^ *  parameter check and work array initialization
c
c^ *  check on mode parameter
      data m2 / 0 /
      data nm1 / 0 /
      data el / 0d0 /
      ier = 0
      if (((((iabs(md) > 4) || (md == 0)) || ((iabs(md) == 1) && (val < zero))) || ((iabs(md) == 3) && (val < zero))) || ((iabs(md) == 4) && ((val < zero) || (val > (n - m)))))
cwrong mode value
      ier = 3
      return
c^ *  check on m and n
      end
      if (md > 0)
      m2 = 2 * m
      nm1 = n - 1
      else
      if ((m2 != (2 * m)) || (nm1 != (n - 1)))
cm or n modified since previous
      ier = 3
      return
      end
      end
      if ((m <= 0) || (n < m2))
cm or n invalid
      ier = 1
      return
c^ *  check on knot sequence and weights
      end
      if (wx(1) <= zero) ier = 2
 end
      for 10 i = 2: n
      if ((wx(i) <= zero) || (x(i - 1) >= x(i))) ier = 2
 end
      if (ier != 0) return
   10 continue
      for 15 j = 1: k
      if (wy(j) <= zero) ier = 2
 end
      if (ier != 0) return
c
c^ *  work array parameters (address information for covariance
c^ *  propagation by means of the matrices stat, b, and we). nb:
c^ *  bwe cannot be used since it is modified by function trinv.
c
   15 continue
      nm2p1 = n * (m2 + 1)
c     istat = 1            #Statistics array STAT(6)
c     ibwe  = istat + 6      #Smoothing matrix BWE( - M:M  , N)
      nm2m1 = n * (m2 - 1)
cdesign matrix    b  (1 - m:m - 1, n)
      ib = ibwe + nm2p1
c     iwk   = iwe   + nm2p1      #Total work array length N + 6 * (N * M + 1)
c
c^ *  compute the design matrices b and we, the ratio
c^ *  of their l1 - norms, and check for iterative mode.
c
cdesign matrix    we ( - m:m  , n)
      iwe = ib + nm2m1
      if (md > 0)
      basis(m, n, x, wk(ib), r1, wk(ibwe))
      prep(m, n, x, wx, wk(iwe), el)
cl1 - norms ratio (saved upon return)
      el = el / r1
      end
c^ *     prior given value for p
      if (iabs(md) != 1) #FIXME goto 20
      r1 = val
c
c^ *  iterate to minimize the gcv function (|md| = 2),
c^ *  the mse function (|md| = 3), or to obtain the prior
c^ *  given number of degrees of freedom (|md| = 4).
c
      #FIXME goto 100
   20 if (md < ( - 1))
cuser - determined starting value
      r1 = wk(4)
      else
cdefault (dof ~ 0.5)
      r1 = one / el
      end
      r2 = r1 * ratio
      gf2 = splc(m, n, k, y, ny, wx, wy, md, val, r2, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
   40 gf1 = splc(m, n, k, y, ny, wx, wy, md, val, r1, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      if (gf1 > gf2) #FIXME goto 50
cinterpolation
      if (wk(4) <= zero) #FIXME goto 100
      r2 = r1
      gf2 = gf1
      r1 = r1 / ratio
      #FIXME goto 40
   50 r3 = r2 * ratio
   60 gf3 = splc(m, n, k, y, ny, wx, wy, md, val, r3, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      if (gf3 > gf2) #FIXME goto 70
cleast - squares polynomial
      if (wk(4) >= epsinv) #FIXME goto 100
      r2 = r3
      gf2 = gf3
      r3 = r3 * ratio
      #FIXME goto 60
   70 r2 = r3
      gf2 = gf3
      α = (r2 - r1) / tau
      r4 = r1 + α
      r3 = r2 - α
      gf3 = splc(m, n, k, y, ny, wx, wy, md, val, r3, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      gf4 = splc(m, n, k, y, ny, wx, wy, md, val, r4, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
   80 if (gf3 <= gf4)
      r2 = r4
      gf2 = gf4
      err = (r2 - r1) / (r1 + r2)
      if ((((err * err) + one) == one) || (err <= tol)) #FIXME goto 90
      r4 = r3
      gf4 = gf3
      α = α / tau
      r3 = r2 - α
      gf3 = splc(m, n, k, y, ny, wx, wy, md, val, r3, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      else
      r1 = r3
      gf1 = gf3
      err = (r2 - r1) / (r1 + r2)
      if ((((err * err) + one) == one) || (err <= tol)) #FIXME goto 90
      r3 = r4
      gf3 = gf4
      α = α / tau
      r4 = r1 + α
      gf4 = splc(m, n, k, y, ny, wx, wy, md, val, r4, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      end
      #FIXME goto 80
c
c^ *  calculate final spline coefficients
c
   90 r1 = half * (r1 + r2)
c
c^ *  ready
c
  100 gf1 = splc(m, n, k, y, ny, wx, wy, md, val, r1, eps, c, nc, wk, wk(ib), wk(iwe), el, wk(ibwe))
      return
c basis.for, 1985 - 06 - 03
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine basis (real * 8)
c
c purpose:
c ^^^ *
c
c       subroutine to assess a b - spline tableau, stored in vectorized
c       form.
c
c calling convention:
c ^^^^^^^^^
c
c       basis ( m, n, x, b, bl, q )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       m               ( i )   half order of the spline (degree 2 * m - 1),
c                               m > 0.0
c       n               ( i )   number of knots, n > = 2 * m.
c       x(n)            ( i )   knot sequence, x(i - 1) < x(i), i = 2, n.
c       b(1 - m:m - 1, n)    ( o )   output tableau. element b(j, i) of array
c                               b corresponds with element b(i, i + j) of
c                               the tableau matrix b.
c       bl              ( o )   l1 - norm of b.
c       q(1 - m:m)        ( w )   internal work array.
c
c remark:
c ^^^
c
c       this subroutine is an adaptation of subroutine basis from the
c       paper by lyche et al. (1983). no checking is performed on the
c       validity of m and n. if the knot sequence is not strictly in -
c       creasing, division by zero may occur.
c
c reference:
c ^^^^ *
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      subroutine basis(m, n, x, b, bl, q)
      implicit double precision (o - z, a - h)
      parameter (zero = 0d0, one = 1d0)
c
      dimension x(n), b(1 - m:m - 1, n), q(1 - m:m)
c^ *         linear spline
      if (m == 1)
      for 3 i = 1: n
      b(0, i) = one
    3 continue
      bl = one
      return
c
c^ *  general splines
c
      end
      mm1 = m - 1
      mp1 = m + 1
      m2 = 2 * m
c^ *     1st row
      for 15 l = 1: n
      for 5 j = - mm1: m
      q(j) = zero
    5 continue
      q(mm1) = one
c^ *     successive rows
      if ((l != 1) && (l != n)) q(mm1) = one / (x(l + 1) - x(l - 1))
 end
      arg = x(l)
      for 13 i = 3: m2
      ir = mp1 - i
      v = q(ir)
c^ *               left - hand b - splines
      if (l < i)
      for 6 j = l + 1: i
      u = v
      v = q(ir + 1)
      q(ir) = u + ((x(j) - arg) * v)
      ir = ir + 1
    6 continue
      end
      j1 = max0((l - i) + 1, 1)
      j2 = min0(l - 1, n - i)
c^ *               ordinary b - splines
      if (j1 <= j2)
      if (i < m2)
      for 8 j = j1: j2
      y = x(i + j)
      u = v
      v = q(ir + 1)
      q(ir) = u + (((v - u) * (y - arg)) / (y - x(j)))
      ir = ir + 1
    8 continue
      else
      for 10 j = j1: j2
      u = v
      v = q(ir + 1)
      q(ir) = ((arg - x(j)) * u) + ((x(i + j) - arg) * v)
      ir = ir + 1
   10 continue
      end
      end
      nmip1 = (n - i) + 1
c^ *           right - hand b - splines
      if (nmip1 < l)
      for 12 j = nmip1: l - 1
      u = v
      v = q(ir + 1)
      q(ir) = ((arg - x(j)) * u) + v
      ir = ir + 1
   12 continue
      end
   13 continue
      for 14 j = - mm1: mm1
      b(j, l) = q(j)
   14 continue
c
c^ *  zero unused parts of b
c
   15 continue
      for 17 i = 1: mm1
      for 16 k = i: mm1
      b( - k, i) = zero
      b(k, (n + 1) - i) = zero
   16 continue
c
c^ *  assess l1 - norm of b
c
   17 continue
      bl = 0d0
      for 19 i = 1: n
      for 18 k = - mm1: mm1
      bl = bl + abs(b(k, i))
   18 continue
   19 continue
c
c^ *  ready
c
      bl = bl / n
      return
c prep.for, 1985 - 07 - 04
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine prep (real * 8)
c
c purpose:
c ^^^ *
c
c       to compute the matrix we of weighted divided difference coeffi -
c       cients needed to set up a linear system of equations for sol -
c       ving b - spline smoothing problems, and its l1 - norm el. the matrix
c       we is stored in vectorized form.
c
c calling convention:
c ^^^^^^^^^
c
c       prep ( m, n, x, w, we, el )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       m               ( i )   half order of the b - spline (degree
c                               2 * m - 1), with m > 0.0
c       n               ( i )   number of knots, with n > = 2 * m.
c       x(n)            ( i )   strictly increasing knot array, with
c                               x(i - 1) < x(i), i = 2, n.
c       w(n)            ( i )   weight matrix (diagonal), with
c                               w(i) > 0.0, i = 1, n.
c       we( - m:m, n)      ( o )   array containing the weighted divided
c                               difference terms in vectorized format.
c                               element we(j, i) of array e corresponds
c                               with element e(i, i + j) of the matrix
c                               w^ - 1 * e.
c       el              ( o )   l1 - norm of we.
c
c remark:
c ^^^
c
c       this subroutine is an adaptation of subroutine prep from the pap
cer
c       by lyche et al. (1983). no checking is performed on the validity
c       of m and n. division by zero may occur if the knot sequence is
c       not strictly increasing.
c
c reference:
c ^^^^ *
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      subroutine prep(m, n, x, w, we, el)
      implicit double precision (o - z, a - h)
      parameter (zero = 0d0, one = 1d0)
c
c^ *  calculate the factor f1
c
cwe( - m:m, n)
      dimension x(n), w(n), we(((2 * m) + 1) * n)
      m2 = 2 * m
      mp1 = m + 1
      m2m1 = m2 - 1
      m2p1 = m2 + 1
      nm = n - m
      f1 = - one
      if (m != 1)
      for 5 i = 2: m
      f1 = - (f1 * i)
    5 continue
      for 6 i = mp1: m2m1
      f1 = f1 * i
    6 continue
c
c^ *  columnwise evaluation of the unweighted design matrix e
c
      end
      i1 = 1
      i2 = m
      jm = mp1
      for 17 j = 1: n
      inc = m2p1
      if (j > nm)
      f1 = - f1
      f = f1
      else
      if (j < mp1)
      inc = 1
      f = f1
      else
      f = f1 * (x(j + m) - x(j - m))
      end
      end
      if (j > mp1) i1 = i1 + 1
 end
      if (i2 < n) i2 = i2 + 1
 end
c^ *     loop for divided difference coefficients
      jj = jm
      ff = f
      y = x(i1)
      i1p1 = i1 + 1
      for 11 i = i1p1: i2
      ff = ff / (y - x(i))
   11 continue
      we(jj) = ff
      jj = jj + m2
      i2m1 = i2 - 1
      if (i1p1 <= i2m1)
      for 14 l = i1p1: i2m1
      ff = f
      y = x(l)
      for 12 i = i1: l - 1
      ff = ff / (y - x(i))
   12 continue
      for 13 i = l + 1: i2
      ff = ff / (y - x(i))
   13 continue
      we(jj) = ff
      jj = jj + m2
   14 continue
      end
      ff = f
      y = x(i2)
      for 16 i = i1: i2m1
      ff = ff / (y - x(i))
   16 continue
      we(jj) = ff
      jj = jj + m2
      jm = jm + inc
c
c^ *  zero the upper left and lower right corners of e
c
   17 continue
      kl = 1
      n2m = (m2p1 * n) + 1
      for 19 i = 1: m
      ku = (kl + m) - i
      for 18 k = kl: ku
      we(k) = zero
      we(n2m - k) = zero
   18 continue
      kl = kl + m2p1
c
c^ *  weighted matrix we = w^ - 1 * e and its l1 - norm
c
   19 continue
   20 jj = 0
      el = 0d0
      for 22 i = 1: n
      wi = w(i)
      for 21 j = 1: m2p1
      jj = jj + 1
      we(jj) = we(jj) / wi
      el = el + abs(we(jj))
   21 continue
   22 continue
c
c^ *  ready
c
      el = el / n
      return
c splc.for, 1985 - 12 - 12
c
c author: h.j. woltring
c
c organizations: university of nijmegen, and
c                philips medical systems, eindhoven
c                (the netherlands)
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c function splc (real * 8)
c
c purpose:
c ^^^ *
c
c       to assess the coefficients of a b - spline and various statistical
c       parameters, for a given value of the regularization parameter p.
c
c calling convention:
c ^^^^^^^^^
c
c       fv = splc ( m, n, k, y, ny, wx, wy, mode, val, p, eps, c, nc,
c       1           stat, b, we, el, bwe)
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       splc            ( o )   gcv function value if |mode| == 2,
c                               mse value if |mode| == 3, and absolute
c                               difference with the prior given number o
cf
c                               degrees of freedom if |mode| == 4.0
c       m               ( i )   half order of the b - spline (degree 2 * m - 1
c),
c                               with m > 0.0
c       n               ( i )   number of observations, with n > = 2 * m.
c       k               ( i )   number of datasets, with k > = 1.0
c       y(ny, k)         ( i )   observed measurements.
c       ny              ( i )   first dimension of y(ny, k), with ny >= n
c.
c       wx(n)           ( i )   weight factors, corresponding to the
c                               relative inverse variance of each measur
ce -
c                               ment, with wx(i) > 0.0.0
c       wy(k)           ( i )   weight factors, corresponding to the
c                               relative inverse variance of each datase
ct,
c                               with wy(j) > 0.0.0
c       mode            ( i )   mode switch, as described in gcvspl.
c       val             ( i )   prior variance if |mode| == 3, and
c                               prior number of degrees of freedom if
c                               |mode| == 4.0 for other values of mode,
c                               val is not used.
c       p               ( i )   smoothing parameter, with p > = 0.0.0 if
c                               p == 0.0, an interpolating spline is
c                               calculated.
c       eps             ( i )   relative rounding tolerance * 10.0.0 eps is
c                               the smallest positive number such that
c                               eps / 10.0 + 1.0 != 1.0.0
c       c(nc, k)         ( o )   calculated spline coefficient arrays. nb
c:
c                               the dimensions of in gcvspl and in splde
cr
c                               are different# In SPLDER, only a single
c                               column of c(n, k) is needed, and the prop
cer
c                               column c(1, j), with j = 1.0..k, should be u
csed
c                               when calling splder.
c       nc              ( i )   first dimension of c(nc, k), with nc >= n
c.
c       stat(6)         ( o )   statistics array. see the description in
c                               subroutine gcvspl.
c       b (1 - m:m - 1, n)   ( i )   b - spline tableau as evaluated by subrout
cine
c                               basis.
c       we( - m:m  , n)   ( i )   weighted b - spline tableau (w^ - 1 * e) as
c                               evaluated by subroutine prep.
c       el              ( i )   l1 - norm of the matrix we as evaluated by
c                               subroutine prep.
c       bwe( - m:m, n)     ( o )   central 2 * m + 1 bands of the inverted
c                               matrix ( b  +  p * w^ - 1 * e )^ - 1
c
c remarks:
c ^^^ *
c
c       this subroutine combines elements of subroutine splc0 from the
c       paper by lyche et al. (1983), and of subroutine spfit1 by
c       hutchinson (1985).
c
c references:
c ^^^^^
c
c       m.f. hutchinson (1985), subroutine cubgcv. csiro division of
c       mathematics and statistics, p.o. box 1965, canberra, act 2601,
c       australia.
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      function splc(m, n, k, y, ny, wx, wy, mode, val, p, eps, c, nc, stat, b, we, el, bwe)
      implicit double precision (o - z, a - h)
      parameter (zero = 0d0, one = 1d0, two = 2d0)
c
c^ *  check on p - value
c
      dimension y(ny, k), wx(n), wy(k), c(nc, k), stat(6), b(1 - m:m - 1, n), we( - m:m, n), bwe( - m:m, n)
      dp = p
      stat(4) = p
c^ *  pseudo - interpolation if p is too small
      pel = p * el
      if (pel < eps)
      dp = eps / el
      stat(4) = zero
c^ *  pseudo least - squares polynomial if p is too large
      end
      if ((pel * eps) > one)
      dp = one / (el * eps)
      stat(4) = dp
c
c^ *  calculate  bwe  =  b  +  p * w^ - 1 * e
c
      end
      for 40 i = 1: n
      km = - min0(m, i - 1)
      kp = min0(m, n - i)
      for 30 l = km: kp
      if (iabs(l) == m)
      bwe(l, i) = dp * we(l, i)
      else
      bwe(l, i) = b(l, i) + (dp * we(l, i))
      end
   30 continue
c
c^ *  solve bwe * c = y, and assess trace [ b * bwe^ - 1 ]
c
   40 continue
      bandet(bwe, m, n)
      bansol(bwe, y, ny, c, nc, m, n, k)
ctrace * p = res. d.o.
      stat(3) = trinv(we, bwe, m, n) * dp
c
c^ *  compute mean - squared weighted residual
c
      trn = stat(3) / n
      esn = zero
      for 70 j = 1: k
      for 60 i = 1: n
      dt = - y(i, j)
      km = - min0(m - 1, i - 1)
      kp = min0(m - 1, n - i)
      for 50 l = km: kp
      dt = dt + (b(l, i) * c(i + l, j))
   50 continue
      esn = esn + (((dt * dt) * wx(i)) * wy(j))
   60 continue
   70 continue
c
c^ *  calculate statistics and function value
c
      esn = esn / (n * k)
cestimated variance
      stat(6) = esn / trn
cgcv function value
      stat(1) = stat(6) / trn
c     stat(3) = trace [p * b * bwe^ - 1] #Estimated residuals" d.o.f.
c     stat(4) = p                     #Normalized smoothing factor
cmean squared residual
      stat(2) = esn
c^ *     unknown variance: gcv
      if (iabs(mode) != 3)
      stat(5) = stat(6) - esn
      if (iabs(mode) == 1) splc = zero
 end
      if (iabs(mode) == 2) splc = stat(1)
 end
      if (iabs(mode) == 4) splc = dabs(stat(3) - val)
 end
c^ *     known variance: estimated mean squared error
      else
      stat(5) = esn - (val * ((two * trn) - one))
      splc = stat(5)
c
      end
      return
c bandet.for, 1985 - 06 - 03
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine bandet (real * 8)
c
c purpose:
c ^^^ *
c
c       this subroutine computes the lu decomposition of an n * n matrix
c       e. it is assumed that e has m bands above and m bands below the
c       diagonal. the decomposition is returned in e. it is assumed that
c       e can be decomposed without pivoting. the matrix e is stored in
c       vectorized form in the array e( - m:m, n), where element e(j, i) of
c       the array e corresponds with element e(i, i + j) of the matrix e.
c
c calling convention:
c ^^^^^^^^^
c
c       bandet ( e, m, n )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       e( - m:m, n)       (i / o)   matrix to be decomposed.
c       m, n            ( i )   matrix dimensioning parameters,
c                               m > = 0, n > = 2 * m.
c
c remark:
c ^^^
c
c       no checking on the validity of the input data is performed.
c       if (m <= 0), no action is taken.
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      subroutine bandet(e, m, n)
      implicit double precision (o - z, a - h)
c
      dimension e( - m:m, n)
      if (m <= 0) return
      for 40 i = 1: n
      di = e(0, i)
      mi = min0(m, i - 1)
      if (mi >= 1)
      for 10 k = 1: mi
      di = di - (e( - k, i) * e(k, i - k))
   10 continue
      e(0, i) = di
      end
      lm = min0(m, n - i)
      if (lm >= 1)
      for 30 l = 1: lm
      dl = e( - l, i + l)
      km = min0(m - l, i - 1)
      if (km >= 1)
      du = e(l, i)
      for 20 k = 1: km
      du = du - (e( - k, i) * e(l + k, i - k))
      dl = dl - (e(( - l) - k, l + i) * e(k, i - k))
   20 continue
      e(l, i) = du
      end
      e( - l, i + l) = dl / di
   30 continue
      end
c
c^ *  ready
c
   40 continue
      return
c bansol.for, 1985 - 12 - 12
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine bansol (real * 8)
c
c purpose:
c ^^^ *
c
c       this subroutine solves systems of linear equations given an lu
c       decomposition of the design matrix. such a decomposition is pro -
c       vided by subroutine bandet, in vectorized form. it is assumed
c       that the design matrix is not singular.
c
c calling convention:
c ^^^^^^^^^
c
c       bansol ( e, y, ny, c, nc, m, n, k )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       e( - m:m, n)       ( i )   input design matrix, in lu - decomposed,
c                               vectorized form. element e(j, i) of the
c                               array e corresponds with element
c                               e(i, i + j) of the n * n design matrix e.
c       y(ny, k)         ( i )   right hand side vectors.
c       c(nc, k)         ( o )   solution vectors.
c       ny, nc, m, n, k ( i )   dimensioning parameters, with m > = 0,
c                               n > 2 * m, and k > = 1.0
c
c remark:
c ^^^
c
c       this subroutine is an adaptation of subroutine bansol from the
c       paper by lyche et al. (1983). no checking is performed on the
c       validity of the input parameters and data. division by zero may
c       occur if the system is singular.
c
c reference:
c ^^^^ *
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      subroutine bansol(e, y, ny, c, nc, m, n, k)
      implicit double precision (o - z, a - h)
c
c^ *  check on special cases: m = 0, m = 1, m>1
c
      dimension e( - m:m, n), y(ny, k), c(nc, k)
      nm1 = n - 1
c
c^ *  m = 0: diagonal system
c
      if (m - 1) 10, 40, 80
   10 for 30 i = 1: n
      for 20 j = 1: k
      c(i, j) = y(i, j) / e(0, i)
   20 continue
   30 continue
c
c^ *  m = 1: tridiagonal system
c
      return
   40 for 70 j = 1: k
      c(1, j) = y(1, j)
cforward sweep
      for 50 i = 2: n
      c(i, j) = y(i, j) - (e( - 1, i) * c(i - 1, j))
   50 continue
      c(n, j) = c(n, j) / e(0, n)
cbackward sweep
      for 60 i = nm1, 1: - 1
      c(i, j) = (c(i, j) - (e(1, i) * c(i + 1, j))) / e(0, i)
   60 continue
   70 continue
c
c^ *  m > 1: general system
c
      return
   80 for 130 j = 1: k
      c(1, j) = y(1, j)
cforward sweep
      for 100 i = 2: n
      mi = min0(m, i - 1)
      d = y(i, j)
      for 90 l = 1: mi
      d = d - (e( - l, i) * c(i - l, j))
   90 continue
      c(i, j) = d
  100 continue
      c(n, j) = c(n, j) / e(0, n)
cbackward sweep
      for 120 i = nm1, 1: - 1
      mi = min0(m, n - i)
      d = c(i, j)
      for 110 l = 1: mi
      d = d - (e(l, i) * c(i + l, j))
  110 continue
      c(i, j) = d / e(0, i)
  120 continue
  130 continue
c
      return
c trinv.for, 1985 - 06 - 03
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c function trinv (real * 8)
c
c purpose:
c ^^^ *
c
c       to calculate trace [ b * e^ - 1 ], where b and e are n * n
c       matrices with bandwidth 2 * m + 1, and where e is a regular matrix
c       in lu - decomposed form. b and e are stored in vectorized form,
c       compatible with subroutines bandet and bansol.
c
c calling convention:
c ^^^^^^^^^
c
c       trace = trinv ( b, e, m, n )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       b( - m:m, n)       ( i ) input array for matrix b. element b(j, i)
c                             corresponds with element b(i, i + j) of the
c                             matrix b.
c       e( - m:m, n)       (i / o) input array for matrix e. element e(j, i)
c                             corresponds with element e(i, i + j) of the
c                             matrix e. this matrix is stored in lu -
c                             decomposed form, with l unit lower tri -
c                             angular, and u upper triangular. the unit
c                             diagonal of l is not stored. upon return,
c                             the array e holds the central 2 * m + 1 bands
c                             of the inverse e^ - 1, in similar ordering.
c       m, n            ( i ) array and matrix dimensioning parameters
c                             (m > 0, n >= 2 * m + 1).
c       trinv           ( o ) output function value trace [ b * e^ - 1 ]
c
c reference:
c ^^^^ *
c
c       a.m. erismanw.f. tinney, on computing certain elements of the
c       inverse of a sparse matrix. communications of the acm 18(1975),
c       nr. 3, pp. 177 - 179.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      double precision function trinv(b, e, m, n)
      implicit double precision (o - z, a - h)
      parameter (zero = 0d0, one = 1d0)
c
c^ *  assess central 2 * m + 1 bands of e^ - 1 and store in array e
c
      dimension b( - m:m, n), e( - m:m, n)
cnth pivot
      e(0, n) = one / e(0, n)
      for 40 i = n - 1, 1: - 1
      mi = min0(m, n - i)
c^ *     save ith column of l and ith row of u, and normalize u row
cith pivot
      dd = one / e(0, i)
      for 10 k = 1: mi
cith row of u (normalized)
      e(k, n) = e(k, i) * dd
cith column of l
      e( - k, 1) = e( - k, k + i)
   10 continue
c^ *     invert around ith pivot
      dd = dd + dd
      for 30 j = mi, 1: - 1
      du = zero
      dl = zero
      for 20 k = 1: mi
      du = du - (e(k, n) * e(j - k, i + k))
      dl = dl - (e( - k, 1) * e(k - j, i + j))
   20 continue
      e(j, i) = du
      e( - j, j + i) = dl
      dd = dd - ((e(j, n) * dl) + (e( - j, 1) * du))
   30 continue
      e(0, i) = 5d - 1 * dd
c
c^ *  assess trace [ b * e^ - 1 ] and clear working storage
c
   40 continue
      dd = zero
      for 60 i = 1: n
      mn = - min0(m, i - 1)
      mp = min0(m, n - i)
      for 50 k = mn: mp
      dd = dd + (b(k, i) * e( - k, k + i))
   50 continue
   60 continue
      trinv = dd
      for 70 k = 1: m
      e(k, n) = zero
      e( - k, 1) = zero
c
c^ *  ready
c
   70 continue
      return
c splder.for, 1985 - 06 - 11
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c function splder (real * 8)
c
c purpose:
c ^^^ *
c
c       to produce the value of the function (ider == 0) or of the
c       iderth derivative (ider > 0) of a 2m - th order b - spline at
c       the point t. the spline is described in terms of the half
c       order m, the knot sequence x(n), n >= 2 * m, and the spline
c       coefficients c(n).
c
c calling convention:
c ^^^^^^^^^
c
c       svider = splder ( ider, m, n, t, x, c, l, q )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       splder  ( o )   function or derivative value.
c       ider    ( i )   derivative order required, with 0 <= ider
c                       and ider <= 2 * m. if ider == 0, the function
c                       value is returned; otherwise, the ider - th
c                       derivative of the spline is returned.
c       m       ( i )   half order of the spline, with m > 0.0
c       n       ( i )   number of knots and spline coefficients,
c                       with n >= 2 * m.
c       t       ( i )   argument at which the spline or its deri -
c                       vative is to be evaluated, with x(1) <= t
c                       and t <= x(n).
c       x(n)    ( i )   strictly increasing knot sequence array,
c                       x(i - 1) < x(i), i = 2, ..., n.
c       c(n)    ( i )   spline coefficients, as evaluated by
c                       subroutine gvcspl.
c       l       (i / o)   l contains an integer such that:
c                       x(l) <= t and t < x(l + 1) if t is within
c                       the range x(1) <= t and t < x(n). if
c                       t < x(1), l is set to 0, and if t >= x(n),
c                       l is set to n. the search for l is facili -
c                       tated if l has approximately the right
c                       value on entry.
c       q(2 * m)  ( w )   internal work array.
c
c remark:
c ^^^
c
c       this subroutine is an adaptation of subroutine splder of
c       the paper by lyche et al. (1983). no checking is performed
c       on the validity of the input parameters.
c
c reference:
c ^^^^ *
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      double precision function splder(ider, m, n, t, x, c, l, q)
      implicit double precision (o - z, a - h)
      parameter (zero = 0d0, one = 1d0)
cf2py intent(out) :: splder
cf2py intent(hide) :: q
c
c^ *  derivatives of ider >= 2 * m are alway zero
c
      dimension x(n), c(n), q(2 * m)
      m2 = 2 * m
      k = m2 - ider
      if (k < 1)
      splder = zero
      return
c
c^ *  search for the interval value l
c
      end
c
c^ *  initialize parameters and the 1st row of the b - spline
c^ *  coefficients tableau
c
      search(n, x, t, l)
      tt = t
      mp1 = m + 1
      npm = n + m
      m2m1 = m2 - 1
      k1 = k - 1
      nk = n - k
      lk = l - k
      lk1 = lk + 1
      lm = l - m
      jl = l + 1
      ju = l + m2
      ii = n - m2
      ml = - l
      for 2 j = jl: ju
      if ((j >= mp1) && (j <= npm))
      q(j + ml) = c(j - m)
      else
      q(j + ml) = zero
      end
c
c^ *  the following loop computes differences of the b - spline
c^ *  coefficients. if the value of the spline is required,
c^ *  differencing is not necessary.
c
    2 continue
      if (ider > 0)
      jl = jl - m2
      ml = ml + m2
      for 6 i = 1: ider
      jl = jl + 1
      ii = ii + 1
      j1 = max0(1, jl)
      j2 = min0(l, ii)
      mi = m2 - i
      j = j2 + 1
      if (j1 <= j2)
      for 3 jin = j1: j2
      j = j - 1
      jm = ml + j
      q(jm) = (q(jm) - q(jm - 1)) / (x(j + mi) - x(j))
    3 continue
      end
      if (jl >= 1) #FIXME goto 6
      i1 = i + 1
      j = ml + 1
      if (i1 <= ml)
      for 5 jin = i1: ml
      j = j - 1
      q(j) = - q(j - 1)
    5 continue
      end
    6 continue
      for 7 j = 1: k
      q(j) = q(j + ider)
    7 continue
c
c^ *  compute lower half of the evaluation tableau
c
      end
ctableau ready if ider == 2 * m - 1
      if (k1 >= 1)
      for 14 i = 1: k1
      nki = nk + i
      ir = k
      jj = l
      ki = k - i
c^ *        right - hand b - splines
      nki1 = nki + 1
      if (l >= nki1)
      for 9 j = nki1: l
      q(ir) = q(ir - 1) + ((tt - x(jj)) * q(ir))
      jj = jj - 1
      ir = ir - 1
    9 continue
c^ *        middle b - splines
      end
      lk1i = lk1 + i
      j1 = max0(1, lk1i)
      j2 = min0(l, nki)
      if (j1 <= j2)
      for 11 j = j1: j2
      xjki = x(jj + ki)
      z = q(ir)
      q(ir) = z + (((xjki - tt) * (q(ir - 1) - z)) / (xjki - x(jj)))
      ir = ir - 1
      jj = jj - 1
   11 continue
c^ *        left - hand b - splines
      end
      if (lk1i <= 0)
      jj = ki
      lk1i1 = 1 - lk1i
      for 13 j = 1: lk1i1
      q(ir) = q(ir) + ((x(jj) - tt) * q(ir - 1))
      jj = jj - 1
      ir = ir - 1
   13 continue
      end
   14 continue
c
c^ *  compute the return value
c
      end
c^ *  multiply with factorial if ider > 0
      z = q(k)
      if (ider > 0)
      for 16 j = k: m2m1
      z = z * j
   16 continue
      end
c
c^ *  ready
c
      splder = z
      return
c search.for, 1985 - 06 - 03
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
c subroutine search (real * 8)
c
c purpose:
c ^^^ *
c
c       given a strictly increasing knot sequence x(1) < ... < x(n),
c       where n > = 1, and a real number t, this subroutine finds the
c       value l such that x(l) < = t < x(l + 1).  if t < x(1), l = 0;
 end
c       if x(n) < = t, l = n.
 end
c
c calling convention:
c ^^^^^^^^^
c
c       search ( n, x, t, l )
c
c meaning of parameters:
c ^^^^^^^^^^ *
c
c       n       ( i )   knot array dimensioning parameter.
c       x(n)    ( i )   stricly increasing knot array.
c       t       ( i )   input argument whose knot interval is to
c                       be found.
c       l       (i / o)   knot interval parameter. the search procedure
c                       is facilitated if l has approximately the
c                       right value on entry.
c
c remark:
c ^^^
c
c       this subroutine is an adaptation of subroutine search from
c       the paper by lyche et al. (1983). no checking is performed
c       on the input parameters and data; the algorithm may fail if
c       the input sequence is not strictly increasing.
c
c reference:
c ^^^^ *
c
c       t. lyche, l.l. schumaker, k. sepehrnoori, fortran subroutines
c       for computing smoothing and interpolating natural splines.
c       advances in engineering software 5(1983)1, pp. 2 - 5.0
c
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *
c
      end
c
      subroutine search(n, x, t, l)
      implicit double precision (o - z, a - h)
c
      dimension x(n)
c^ *     out of range to the left
      if (t < x(1))
      l = 0
      return
      end
c^ *     out of range to the right
      if (t >= x(n))
      l = n
      return
c^ *  validate input value of l
      end
      l = max0(l, 1)
c
c^ *  often l will be in an interval adjoining the interval found
c^ *  in a previous to search
c
      if (l >= n) l = n - 1
 end
      if (t >= x(l)) #FIXME goto 5
      l = l - 1
c
c^ *  perform bisection
c
      if (t >= x(l)) return
      il = 1
    3 iu = l
    4 l = (il + iu) / 2
      if ((iu - il) <= 1) return
      if (t < x(l)) #FIXME goto 3
      il = l
      #FIXME goto 4
    5 if (t < x(l + 1)) return
      l = l + 1
      if (t < x(l + 1)) return
      il = l + 1
      iu = n
c
      #FIXME goto 4
      end

max_stars_repo_name: charlesll/gcvspl

max_issues_repo_name: charlesll/gcvspl

max_forks_repo_name: charlesll/gcvspl


text:
using JLD
using PyPlot
function plot_tangent_sens_cum_avg()
	n_files = 10000
	dJds_avg = zeros(2)
	dJds_var = zeros(2)
	dJds_arr = zeros(2,n_files)
	nSteps = 2000
	#vsh_arr = zeros(d, nSteps, n_files)
	for n = 1:10000
		filename = string("../data/",
						  "rijke_tangent_sensitivities/",
						  "corrected/",
						  "vsh_and_dJds_",
						  string(n),
						  ".jld")

		X = load(filename)
		dJds = X["dJds"]
		#vsh = X["vsh"]
		dJds_arr[:,n] = dJds
		#vsh_arr[:,:,n] = vsh
	end
	#=
	for n = 4001:n_files
		filename = string("../data/",
						  "rijke_tangent_sensitivities/",
						  "withf/",
						  "vsh_and_dJds_",
						  string(n + 1000),
						  ".jld")

		X = load(filename)
		dJds = X["dJds"]
		#vsh = X["vsh"]
		dJds_arr[:,n] = dJds
		#vsh_arr[:,:,n] = vsh
	end
	=#
	#indices = -10.0 .<= dJds_arr[1,:] .<= 10.0
	#n_files = sum(indices)
	#dJds_arr = dJds_arr[:,indices]
	#X = load("../data/sensitivities/sensitivity_t2000_dt01_n500.jld")
	#dJds = X["dJdbeta"]
	#indices = -10.0 .<= dJds[1,:] .<= 10.0
	#n_files += sum(indices)
	#dJds = dJds[:,indices]
	#dJds_arr = [dJds_arr dJds]
	n_files = size(dJds_arr)[2]
	dJds_avg = sum(dJds_arr, dims=2)/n_files
	dJds_var = sum((dJds_arr .- dJds_avg).^2.0)/n_files
	fig, ax = subplots(1,1)
	fig1, ax1 = subplots(1,1)

	ax.xaxis.set_tick_params(labelsize=28)
	ax1.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax1.yaxis.set_tick_params(labelsize=28)
	T = 20.0
	ax.plot(T*(1:n_files),
			cumsum(dJds_arr[1,:])./(1:n_files),
			"r.",  label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \beta}$",ms=4.0)
	#ax.set_ylim([-4.,-2.75])
	ax1.plot(T*(1:n_files),
			 cumsum(dJds_arr[2,:])./(1:n_files),
			 "b.", label=L"$\dfrac{d \langle J_{\rm ray}\rangle}{d \beta}$",ms=4.0)
	ax.plot(T*(1:n_files), ones(n_files)*dJds_avg[1],"r",lw=2.0)
	ax1.plot(T*(1:n_files), ones(n_files)*dJds_avg[2],"b",lw=2.0)
	ax.set_xlabel("Time",fontsize=28)
	ax.set_ylabel(L"$\dfrac{d\langle J_{\rm ac}\rangle}{d\beta}$",fontsize=28)
	ax.grid(true)
	ax1.set_xlabel("Time",fontsize=28)
	ax1.set_ylabel(L"$\dfrac{d\langle J_{\rm ray}\rangle}{d\beta}$",fontsize=28)
	ax1.set_ylim([-1.6,-0.9])
	ax.set_ylim([-4.6,-3.6])
	ax1.grid(true)
	lgnd = fig.legend(loc=(0.9,0.7),fontsize=28)
	lgnd.legendHandles[1]._legmarker.set_markersize(10)
	lgnd = fig1.legend(loc=(0.9,0.2),fontsize=28)
	lgnd.legendHandles[1]._legmarker.set_markersize(10)

	return dJds_arr
end
function plot_adjoint_sens_cum_avg()
	n_files = 1000
	dJds_avg = zeros(2)
	dJds_var = zeros(2)
	dJds_arr = zeros(2,n_files)
	nSteps = 2000
	#vsh_arr = zeros(d, nSteps, n_files)
	for n = 1:n_files
		filename = string("../data/",
						  "rijke_adjoint_sensitivities/",
						  "corrected/",
						  "vsh_and_dJds_",
						  string(n),
						  ".jld")

		X = load(filename)
		dJds = X["dJds"]
		#vsh = X["vsh"]
		#dJds *= nSteps
		dJds_arr[:,n] = dJds
		#vsh_arr[:,:,n] = vsh
	end
	indices = -100.0 .<= dJds_arr[1,:] .<= 100.0
	n_files = sum(indices)
	dJds_arr = dJds_arr[:, indices]
	dJds_avg = sum(dJds_arr, dims=2)/n_files
	dJds_var = sum((dJds_arr .- dJds_avg).^2.0, dims=2)/n_files
	fig, ax = subplots(1,1)
	fig1, ax1 = subplots(1,1)

	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax1.xaxis.set_tick_params(labelsize=28)
	ax1.yaxis.set_tick_params(labelsize=28)
	T = 20.
	ax.plot(T*(1:n_files),
			cumsum(dJds_arr[1,:])./(1:n_files), "r.", label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \beta}$",ms=4.0)
	ax1.plot(T*(1:n_files),
			 cumsum(dJds_arr[2,:])./(1:n_files), "g.", label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \tau}$",ms=4.0)
	ax.plot(T*(1:n_files), ones(n_files)*dJds_avg[1],"r",lw=2.0)
	ax1.plot(T*(1:n_files), ones(n_files)*dJds_avg[2],"g",lw=2.0)

	#ax.set_ylim([-4.9,-3.1])
	#ax1.set_ylim([100,200])
	ax.set_xlabel("Time",fontsize=28)
	ax.set_ylabel(L"$\dfrac{d\langle J_{\rm ac}\rangle}{d\beta}$",fontsize=28)
	ax.grid(true)
	ax1.set_xlabel("Time",fontsize=28)
	ax1.set_ylabel(L"$\dfrac{d\langle J_{\rm ac}\rangle}{d\tau}$",fontsize=28)
	#ax1.set_ylabel("Shadowing sensitivities",fontsize=28)
	ax1.grid(true)
	#lgnd = fig.legend(loc=(0.5,0.6),fontsize=28)
	#lgnd.legendHandles[1]._legmarker.set_markersize(10)
	#lgnd = fig1.legend(loc=(0.5,0.6),fontsize=28)
	#lgnd.legendHandles[1]._legmarker.set_markersize(10)
	#plt.xticks(LinRange(0,200000,5))
	#plt.locator_params(axis='x', nbins=6)

	return dJds_arr
end

#=
fig1, ax1 = subplots(1,1)
ax1.xaxis.set_tick_params(labelsize=18)
ax1.yaxis.set_tick_params(labelsize=18)

ax1.errorbar(n_arr*0.01, dJds[1,:], sqrt.(var_dJds[1,:]), marker="o", mfc="blue", mec="blue", label=L"$\dfrac{d \langle J_{\rm ac}\rangle}{d \beta}$")
ax1.set_xlabel("time segment length",fontsize=18)
ax1.set_ylabel("Shadowing sensitivities",fontsize=18)

fig1.legend(fontsize=18)
ax1.grid(true)

fig2, ax2 = subplots(1,1)
ax2.xaxis.set_tick_params(labelsize=18)
ax2.yaxis.set_tick_params(labelsize=18)

ax2.errorbar(n_arr*0.01, dJds[2,:], sqrt.(var_dJds[2,:]), marker="o", mfc="red", mec="red", color="red", label=L"$\dfrac{d \langle J_{\rm ray}\rangle}{d \beta}$")

ax2.set_xlabel("time segment length",fontsize=18)
ax2.set_ylabel("Shadowing sensitivities",fontsize=18)

fig2.legend(fontsize=18)
ax2.grid(true)
=#

ext: jl

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
abstract type AbstractTrainingData end
abstract type AbstractTestingData end

"""
Represents training dataset for GP regression problem.
"""
struct trainingData <: AbstractTrainingData
    x::Array{Float64, 2} #matrix of horizontal location vectors stacked vertically
    Fx::Array{Float64, 2} #matrix of covariates
    y::Array{Float64, 1} #array of labels, normalized
    ymax::Float64
    d::Int64 #dimension of location vectors in x
    p::Int64 #dimension of covariate vectors in Fx
    n:: Int64 # number of incorporated points
    function trainingData(x, Fx, y)
        @assert Base.size(x, 1) == Base.size(Fx, 1)
        @assert Base.size(Fx, 1) == length(y)
        ymax = maximum(y)
        y ./= ymax
        return new(x, Fx, y, ymax, Base.size(x, 2), Base.size(Fx, 2), Base.size(x, 1))
    end
end

"""
Represents mutable training dataset which can be extended with each newly incorporated observation
First inner constructor initializes and returns empty extensible container for training data
Second inner constructor initializes filled extensible container using x, Fx, y
"""
mutable struct extensible_trainingData<:AbstractTrainingData
    x::Array{Float64, 2} #locations
    Fx::Array{Float64, 2} #covariates
    y::Array{Float64, 1} #labels
    d::Int64 #dimension of data points in container
    p::Int64 #dimension of covariate vectors associated with each data point
    n::Int64 #number of incorporated points so far
    capacity::Int64 #amount of preallocated space
    function extensible_trainingData(d::Int64, p::Int64, capacity=300)::extensible_trainingData
        x = Array{Float64}(undef, capacity, d)
        Fx = Array{Float64}(undef, capacity, p)
        y = Array{Float64}(undef, capacity)
        n = 0 #number incorporated points
        new(x, Fx, y, d, p, n, capacity)
    end
    function extensible_trainingData(x, Fx, y, capacity=300)
        @assert Base.size(x, 1) == Base.size(Fx, 1)
        @assert Base.size(Fx, 1) == length(y)
        d = size(x, 2)
        n = size(x, 1)
        p = size(Fx, 2)
        x_full = Array{Float64}(undef, capacity, d)
        Fx_full = Array{Float64}(undef, capacity, p)
        y_full = Array{Float64}(undef, capacity)
        x_full[1:n, :] = x
        Fx_full[1:n, :] = Fx
        y_full[1:n] = y
        new(x_full, Fx_full, y_full, d, p, n, capacity)
    end
end

getLabel(td::trainingData) = td.y
getPosition(td::trainingData) = td.x
getCovariates(td::trainingData) = td.Fx

getLabel(td::extensible_trainingData) = td.y[1:td.n]
getPosition(td::extensible_trainingData) = td.x[1:td.n, :]
getCovariates(td::extensible_trainingData) = td.Fx[1:td.n, :]

getDimension(td::AbstractTrainingData) = td.d
getCovDimension(td::AbstractTrainingData) = td.p
getNumPts(td::AbstractTrainingData) = td.n

unpack(t::trainingData) = (t.x, t.Fx, t.y, t.d, t.n, t.p)
unpack(t::extensible_trainingData) = (t.x[1:t.n, :], t.Fx[1:t.n, :], t.y[1:t.n], t.d, t.n, t.p)

function print(data::AbstractTrainingData)
    println("\n=============== TRAINING DATA =============== ")
    println("Number of training points:     $(getNumPts(data))")
    println("Dimension of training points:  $(getDimension(data))")
    println("Dimension of covariates:       $(getCovDimension(data))")
    flush(stdout)

    # println("\nTraining points position: ")
    # flush(stdout)
    # display(getPosition(data)')
    # println("Covariates: ")
    # flush(stdout)
    # display(getCovariates(data)')
end

"""
Represents a set of testing data. Currently supports single-point prediction.
"""
mutable struct testingData<:AbstractTestingData
    x0::Array{T} where T<:Real
    Fx0::Array{T} where T<:Real
    y0_true::Union{Array{T}, Nothing} where T<:Real
    d::Int64
    p::Int64
    k::Int64
    testingData(x0::Array{T, 2}, Fx0::Array{T, 2}; y0_true=nothing) where T<:Real = (@assert size(x0, 1)==size(Fx0, 1);  new(x0, Fx0, y0_true, size(x0, 2), size(Fx0, 2), size(x0, 1)))
    testingData() = new()
end
getPosition(td::testingData) = td.x0
getCovariates(td::testingData) = td.Fx0
getDimension(td::testingData) = td.d
getCovDimension(td::testingData) = td.p
getNumPts(td::testingData) = td.k

function print(data::AbstractTestingData)
    println("\n=============== TESTING DATA =============== ")
    println("Number of testing points:     $(getNumPts(data))")
    println("Dimension of testing points:  $(getDimension(data))")
    println("Dimension of covariates:      $(getCovDimension(data))")
    flush(stdout)

    # println("\nTesting points position: ")
    # flush(stdout)
    # display(getPosition(data)')
    # println("Covariates: ")
    # flush(stdout)
    # display(getCovariates(data)')
end
#unpack(t::testingData) = (t.x0, t.Fx0, t.k) #never used, because we will always supply this data when calling pdf, cdf, etc.

"""
Updates extensible training dataset with new data points (locations, covariates, labels)
"""
function update!(e::extensible_trainingData, x0, Fx0, y0)
    @assert typeof(x0)<:Array{T, 2} where T<:Real
    @assert typeof(Fx0)<:Array{T, 2} where T<:Real
    @assert typeof(y0)<:Array{T, 1} where T<:Real
    k = size(x0, 1)
    if e.n + size(x0, 1) > e.capacity
        throw(BoundsError)
    end
    e.x[e.n + 1 : e.n + k, :] = x0
    e.Fx[e.n + 1 : e.n + k, :] = Fx0
    e.y[e.n + 1 : e.n + k] = y0
    e.n += k
    return nothing
end

"""
Replaces prediction location and covariates in testingData with new
new locations and covariates
"""
function update!(e:: AbstractTestingData, x0, Fx0)
    try
        @assert typeof(x0)<:Array{T, 2} where T<:Real
    catch e
        @warn "update x0 value is not 2-D array"
        #@warn "x0", x0
    end
    @assert typeof(Fx0)<:Array{T, 2} where T<:Real
    @assert size(x0, 1) == size(Fx0, 1)
    if (bool = [x0[i] == NaN for i = 1:length(x0)]; reduce((x, y)-> x || y, bool)) == true
        @warn "entry of x0 is NaN"
        @info "x0", x0
        @info "Fx0", Fx0
    end
    #@info("x0", x0)
    e.x0 = x0
    e.Fx0 = Fx0
    e.d = size(x0, 2)
    e.p = size(Fx0, 2)
    e.k = size(x0, 1)
    return nothing
end

function update_needed(e::AbstractTestingData, x0, Fx0)
    try
        a = getPosition(e) != x0
        b = getCovariates(e) != Fx0
        return ( a || b )
    catch UndefRefError #test buffer has not yet been initialized
        return true
    end
end

"""
Get capacity of extensible training object or number of data points in vanilla training object
"""
function getCapacity(e::AbstractTrainingData)
    typeof(e)<:extensible_trainingData ? e.capacity : e.n
end

function checkCompatibility(x::AbstractTrainingData, y::AbstractTestingData)
    #println("getDimension(x): ", getDimension(x))
    #println("getDimension(y): ", getDimension(y))
    #println("getCovDimension(x): ", getCovDimension(x))
    #println("getCovDimension(y): ", getCovDimension(y))

    return getDimension(x) == getDimension(y) && getCovDimension(x) == getCovDimension(y)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp

"""
    gen_grad_func(vi::VarInfo, sampler::Sampler, model)

Generate a function that takes a vector of reals `θ` and compute the logpdf and
gradient at `θ` for the model specified by `(vi, sampler, model)`.
"""
function gen_grad_func(vi::VarInfo, sampler::Sampler, model)
    return θ::AbstractVector{<:Real}->gradient(θ, vi, model, sampler)
end

"""
    gen_lj_func(vi::VarInfo, sampler::Sampler, model)

Generate a function that takes `θ` and returns logpdf at `θ` for the model specified by
`(vi, sampler, model)`.
"""
function gen_lj_func(vi::VarInfo, sampler::Sampler, model)
    return function(θ::AbstractVector{<:Real})
        vi[sampler] = θ
        return runmodel!(model, vi, sampler).logp
    end
end

"""
  gen_rev_func(vi::VarInfo, sampler::Sampler)

Generate a function on `(θ, logp)` that sets the variables referenced by `sampler` to `θ`
and the current `vi.logp` to `logp`.
"""
function gen_rev_func(vi::VarInfo, sampler::Sampler)
    return function(θ::AbstractVector{<:Real}, logp::Real)
        vi[sampler] = θ
        setlogp!(vi, logp)
    end
end

"""
    gen_log_func(sampler::Sampler)

Generate a function that takes no argument and performs logging for the number of leapfrog
steps used in `sampler`.
"""
function gen_log_func(sampler::Sampler)
    return function()
        sampler.info[:lf_num] += 1
        sampler.info[:total_lf_num] += 1
    end
end

function sample_momentum(vi::VarInfo, sampler::Sampler)
    d = length(getranges(vi, sampler))
    stds = sampler.info[:wum][:stds]
    return _sample_momentum(d, stds)
end

_sample_momentum(d::Int, stds::Vector) = randn(d) ./ stds

function runmodel!(model::Function, vi::VarInfo, spl::Union{Nothing,Sampler})
    setlogp!(vi, zero(Real))
    if spl != nothing && :total_eval_num ∈ keys(spl.info)
        spl.info[:total_eval_num] += 1
    end
    Base.invokelatest(model, vi, spl)
    return vi
end

function leapfrog(
    θ::AbstractVector{<:Real},
    p::AbstractVector{<:Real},
    τ::Int,
    ϵ::Real,
    model::Function,
    vi::VarInfo,
    sampler::Sampler,
)
    lp_grad = gen_grad_func(vi, sampler, model)
    rev = gen_rev_func(vi, sampler)
    logger = gen_log_func(sampler)
    return _leapfrog(θ, p, τ, ϵ, lp_grad; rev_func=rev, log_func=logger)
end

function _leapfrog(
    θ::AbstractVector{<:Real},
    p::AbstractVector{<:Real},
    τ::Int,
    ϵ::Real,
    lp_grad_func::Function;
    rev_func=nothing,
    log_func=nothing,
)
    _, grad = lp_grad_func(θ)
    verifygrad(grad) || (return θ, p, 0)

    p, θ, τ_valid = deepcopy(p), deepcopy(θ), 0

    p .-= ϵ .* grad ./ 2
    for t in 1:τ

        log_func != nothing && log_func()

        θ .+= ϵ .* p
        logp, grad = lp_grad_func(θ)

        # If gradients explode, tidy up and return.
        if ~verifygrad(grad)
            θ .-= ϵ .* p
            rev_func != nothing && rev_func(θ, logp)
            break
        end

        p .-= ϵ .* grad
        τ_valid += 1
    end

    # Undo half a step in the momenta.
    p .+= ϵ .* grad ./ 2

    return θ, p, τ_valid
end

# Compute the Hamiltonian
function find_H(p::AbstractVector{<:Real}, model::Function, vi::VarInfo, sampler::Sampler)
    logpdf_func_float = gen_lj_func(vi, sampler, model)
    return _find_H(vi[sampler], p, logpdf_func_float, sampler.info[:wum][:stds])
end

function _find_H(
    θ::AbstractVector{<:Real},
    p::AbstractVector{<:Real},
    logpdf_func_float::Function,
    stds::AbstractVector{<:Real},
)
    return _find_H(θ, p, logpdf_func_float(θ), stds)
end

function _find_H(
    θ::AbstractVector{<:Real},
    p::AbstractVector{<:Real},
    logp::Real,
    stds::AbstractVector{<:Real},
)
    H = sum(abs2, p .* stds) / 2 - logp
    return isnan(H) ? Inf : H
end

# TODO: remove used Turing-wrapper functions

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/base_hmc.hpp
function find_good_eps(model::Function, vi::VarInfo, sampler::Sampler{T}) where T
    println("[Turing] looking for good initial eps...")
    ϵ = 0.1

    p = sample_momentum(vi, sampler)
    H0 = find_H(p, model, vi, sampler)

    θ = vi[sampler]
    θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, sampler)
    h = τ == 0 ? Inf : find_H(p_prime, model, vi, sampler)

    delta_H = H0 - h
    direction = delta_H > log(0.8) ? 1 : -1

    iter_num = 1

    # Heuristically find optimal ϵ
    while (iter_num <= 12)

        p = sample_momentum(vi, sampler)
        H0 = find_H(p, model, vi, sampler)

        θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, sampler)
        h = τ == 0 ? Inf : find_H(p_prime, model, vi, sampler)
        @debug "direction = $direction, h = $h"

        delta_H = H0 - h

        if ((direction == 1) && !(delta_H > log(0.8)))
            break
        elseif ((direction == -1) && !(delta_H < log(0.8)))
            break
        else
            ϵ = direction == 1 ? 2.0 * ϵ : 0.5 * ϵ
        end

        iter_num += 1
    end

    while h == Inf  # revert if the last change is too big
        ϵ = ϵ / 2               # safe is more important than large
        θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, sampler)
        h = τ == 0 ? Inf : find_H(p_prime, model, vi, sampler)
    end
    println("\r[$T] found initial ϵ: ", ϵ)
    return ϵ
end

"""
    mh_accept(H::Real, H_new::Real)
    mh_accept(H::Real, H_new::Real, log_proposal_ratio::Real)

Peform MH accept criteria. Returns a boolean for whether or not accept and the
acceptance ratio in log space.
"""
mh_accept(H::Real, H_new::Real) = log(rand()) + H_new < min(H_new, H), min(0, -(H_new - H))
function mh_accept(H::Real, H_new::Real, log_proposal_ratio::Real)
    return log(rand()) + H_new < H + log_proposal_ratio, min(0, -(H_new - H))
end


text:
include("../src/TableSchema.jl")
using TableSchema

import HTTP: request
import DelimitedFiles: readdlm

# Either import the functions, or use TableSchema.read(<Table>) in the code
import TableSchema: read, is_valid, validate

REMOTE_URL = "https://raw.githubusercontent.com/frictionlessdata/tableschema-jl/master/data/data_simple.csv"

println( "Fetching remote data ..." )
t = Table()
req = request("GET", REMOTE_URL)
data = readdlm(req.body, ',')
tr = read(t, data=data, cast=false) # Array{Any,2}

column1 = tr[:,1]
println( "The length is ", length(column1) ) # 3
println( "Fun cities are ", join([ row for row in column1 ], ",") ) # london,paris,rome

println( "Fetching remote schema ..." )
s = Schema("https://raw.githubusercontent.com/frictionlessdata/tableschema-jl/master/data/schema_valid_simple.json")

if is_valid(s); println("A valid Schema is ready"); end

t.schema = s
if validate(t); println("The table is valid according to the Schema"); end

for err in t.errors
    println(string(err.field.name, " has error: ", err.message))
end


text:
@testset "OrderBook" begin

    @testset "Cancel orders" begin
        max_price = UInt128(110)
        min_price = UInt128(90)
        orderbook = OrderBook(max_price, min_price)

        trades, order_id = limitorder!(orderbook, BUY, UInt128(100), UInt128(100), UInt128(1))

        cancel!(orderbook, order_id)
        @test best_bid(orderbook) == 90
        @test best_ask(orderbook) == 110

        trades, order_id = limitorder!(orderbook, SELL, UInt128(100),
                                       UInt128(100), UInt128(1))

        cancel!(orderbook, order_id)
        @test best_bid(orderbook) == 90
        @test best_ask(orderbook) == 110

    end

    @testset "Update orders" begin

        max_price = UInt128(110)
        min_price = UInt128(90)
        orderbook = OrderBook(max_price, min_price)

        trades, order_id = limitorder!(orderbook, BUY, UInt128(100), UInt128(100),
                                       UInt128(1))

        update!(orderbook, order_id, UInt128(50))
        @test best_bid_size(orderbook) == 50

        trades, order_id = limitorder!(orderbook, SELL, UInt128(101), UInt128(100),
                                       UInt128(1))

        update!(orderbook, order_id, UInt128(20))
        @test best_ask_size(orderbook) == 20
    end

    @testset "Passive Orders" begin
        max_price = UInt128(100)
        min_price = UInt128(1)
        orderbook = OrderBook(max_price, min_price)
        order_id_expected = UInt128(0)
        @testset "Insert Buy orders" begin
            for p in 1:20
                price = UInt128(p)
                size = UInt128(10)
                trader_id = UInt128(1)
                side = BUY # Buy
                trades, order_id = limitorder!(orderbook, side,
                                               price, size, trader_id)

                # Check order_id
                order_id_expected = order_id_expected + 1
                @test order_id == order_id_expected

                # Check bid and bid size
                bid = best_bid(orderbook)
                @test bid == price
                bid_size = best_bid_size(orderbook)
                @test bid_size == size
            end
        end
        @testset "Insert Sell orders" begin
            for p in 100:-1:80
                price = UInt128(p)
                size = UInt128(10)
                trader_id = UInt128(1)
                side = SELL # Buy
                trades, order_id = limitorder!(orderbook, side,
                                               price, size, trader_id)

                # Check order_id
                order_id_expected = order_id_expected + 1
                @test order_id == order_id_expected

                # Check bid and bid size
                ask = best_ask(orderbook)
                @test ask == price
                ask_size = best_ask_size(orderbook)
                @test ask_size == size
            end
        end
    end

    @testset "Aggressive Orders" begin

        @testset "Buy Order" begin

            max_price = UInt128(15)
            min_price = UInt128(7)
            orderbook = OrderBook(max_price, min_price)
            order_id_expected = UInt128(0)

            # Place passive sell order
            limitorder!(orderbook, SELL, UInt128(10), UInt128(100), UInt128(1))

            # Match sell order
            trades, order_id = limitorder!(orderbook, BUY, UInt128(10),
                                          UInt128(10), UInt128(2))

            @test amount(trades) == 10*10

            trades, order_id = limitorder!(orderbook, BUY, UInt128(11),
                                           UInt128(40), UInt128(2))

            @test amount(trades) == 10*40

            # Not there exists 50 left of the sell order
            ask = best_ask(orderbook)
            ask_volume = best_ask_size(orderbook)
            @test ask == 10
            @test ask_volume == 50

            # We insert another sell order at price 11
            limitorder!(orderbook, SELL, UInt128(11), UInt128(100), UInt128(1))

            # If we place a buy order at 11 for 151 we should get 50 for 10 and
            # 100 for 11
            trades, order_id = limitorder!(orderbook, BUY, UInt128(11),
                                           UInt128(151), UInt128(2))

            @test amount(trades) == 10*50 + 11*100
        end

        @testset "Sell Order" begin
            max_price = UInt128(15)
            min_price = UInt128(7)
            orderbook = OrderBook(max_price, min_price)
            order_id_expected = UInt128(0)

            # Place passive buy order
            limitorder!(orderbook, BUY, UInt128(10), UInt128(100), UInt128(1))
            # Match buy order
            trades, order_id = limitorder!(orderbook, SELL, UInt128(10),
                                           UInt128(10), UInt128(2))

            @test amount(trades) == 10*10

            trades, order_id = limitorder!(orderbook, SELL, UInt128(9),
                                           UInt128(40), UInt128(2))

            @test amount(trades) == 10*40

            # Not there exists 50 left of the buy order
            bid = best_bid(orderbook)
            bid_volume = best_bid_size(orderbook)
            @test bid == 10
            @test bid_volume == 50

            # We insert another buy order at price 9
            limitorder!(orderbook, BUY, UInt128(9), UInt128(100), UInt128(1))

            # If we place a sell order at 9 for 151 we should get 50 for 10 and
            # 100 for 9
            trades, order_id = limitorder!(orderbook, SELL, UInt128(9),
                                           UInt128(151), UInt128(2))

            @test amount(trades) == 10*50 + 9*100
        end

        @testset "Aggressive limit order becomes bid and ask" begin
            max_price = UInt128(15)
            min_price = UInt128(7)
            orderbook = OrderBook(max_price, min_price)
            # Place passive sell order to me matched with
            limitorder!(orderbook, SELL, UInt128(10), UInt128(10), UInt128(1))
            # Matches with this order and adds as a passive buy order
            trades, order_id = limitorder!(orderbook, BUY, UInt128(11), UInt128(20),
                                           UInt128(2))
            @test amount(trades) == 10*10
            @test best_bid(orderbook) == 11
            @test best_bid_size(orderbook) == 10
            @test best_ask(orderbook) == max_price

            # Places aggressive sell order that becomes passive
            trades, order_id = limitorder!(orderbook, SELL, UInt128(9), UInt128(20),
                                           UInt128(1))
            @test amount(trades) == 11*10
            @test best_ask(orderbook) == 9
            @test best_ask_size(orderbook) == 10
            @test best_bid(orderbook) == min_price
        end
    end
end


text:
# Linear interpolation routines

# ValsVector type

struct ValsVector{T<:Real}
    data::Vector{T}
end

Base.length(vals::ValsVector) = length(vals.data)

@inline function (vals::ValsVector{TV})(i::Integer, w::TW) where {TV<:Real,TW<:Real}
    T = promote_type(TV, TW)
    n = length(vals)
    1 <= i <= n || throw(DomainError())
    i == n && return convert(T, vals.data[n])
    return vals.data[i] * (1 - w) + vals.data[i+1] * w
end

# search_index_weight method

@inline function search_index_weight(a::AbstractVector, x::Real)
    T = Float64
    n = length(a)
    i = searchsortedlast(a, x)
    i <= 0 && return 1, zero(T)
    i >= n && return n, zero(T)

    w::T = (x - a[i]) / (a[i+1] - a[i])
    return i, w
end

@inline function _search_index_weight(start, step, stop, x::Real)
    T = Float64
    x <= start && return 1, zero(T)
    x >= stop && (x = convert(typeof(x), stop))

    w::T = (x - start) / step
    i = floor(Int, w)
    w -= i
    i += 1
    return i, w
end

search_index_weight(a::Range, x::Real) =
    _search_index_weight(first(a), step(a), last(a), x)

# MyLinInterp type

struct MyLinInterp{TG<:AbstractVector,TV<:Real}
    grid::TG
    vals::ValsVector{TV}
end

MyLinInterp(grid::AbstractVector{T1}, vals::Vector{T2}) where {T1<:Real,T2<:Real} =
    MyLinInterp(grid, ValsVector(vals))

@inline function (f::MyLinInterp)(x::Real)
    t = search_index_weight(f.grid, x)
    return f.vals(t[1], t[2])
end


text:
# TODO: Custom "range" of Frequencies to form a spectrum for, e.g. lightning sferic

"""
    Frequency

Electromagnetic wave defined by frequency `f`, but also carrying angular wave frequency `ω`,
wavenumber `k`, and wavelength `λ`, all in SI units.
"""
struct Frequency
    f::Float64
    ω::Float64
    k::Float64
    λ::Float64
end

"""
    Frequency(f)

Return `Frequency` given electromagnetic wave frequency `f` in Hertz.
"""
function Frequency(f)
    ω = 2π*f
    k = ω/C0
    λ = C0/f

    Frequency(f, ω, k, λ)
end
Frequency(f::Frequency) = f

########

"""
    Emitter

Abstract type for types that energize the waveguide with electromagnetic energy.
"""
abstract type Emitter end

function transmitterchecks(alt)
    alt < 0 && error("Transmitter altitude cannot be negative.")

    return true
end

"""
    Transmitter{A<:Antenna} <: Emitter

Typical ground-based Transmitter.

# Fields

- `name::String`: transmitter name.
- `latitude::Float64`: transmitter geographic latitude in degrees.
- `longitude::Float64`: transmitter geographic longitude in degrees.
- `antenna::Antenna`: transmitter antenna.
- `frequency::Frequency`: transmit frequency.
- `power::Float64`: transmit power in Watts.
"""
struct Transmitter{A<:Antenna} <: Emitter
    name::String
    latitude::Float64
    longitude::Float64
    antenna::A
    frequency::Frequency
    power::Float64

    Transmitter{A}(n, la, lo, an, fr, po) where {A<:Antenna} = new(n, la, lo, an, fr, po)
end

altitude(t::Transmitter) = 0.0
frequency(t::Transmitter) = t.frequency
azimuth(t::Transmitter) = azimuth(t.antenna)
inclination(t::Transmitter) = inclination(t.antenna)
power(t::Transmitter) = t.power

"""
    Transmitter(name, lat, lon, frequency)

Return a `Transmitter` with a `VerticalDipole` antenna and transmit power of 1000 W.
"""
Transmitter(name, lat, lon, frequency) =
    Transmitter{VerticalDipole}(name, lat, lon, VerticalDipole(), Frequency(frequency), 1000)

"""
    Transmitter(frequency)

Return a `Transmitter` with zeroed geographic position, 1000 W transmit power, and
`VerticalDipole` antenna.
"""
Transmitter(frequency) = Transmitter("", 0, 0, frequency)

"""
    Transmitter(antenna, frequency, power)

Return a `Transmitter` with zeroed geographic position.
"""
Transmitter(antenna::A, frequency, power) where {A<:Antenna} =
    Transmitter{A}("", 0, 0, antenna, Frequency(frequency), power)

"""
    AirborneTransmitter{A<:Antenna} <: Emitter

Similar to `Transmitter` except it also contains an `altitude` field in meters.
"""
struct AirborneTransmitter{A<:Antenna} <: Emitter
    name::String
    latitude::Float64
    longitude::Float64
    altitude::Float64
    antenna::A
    frequency::Frequency
    power::Float64

    AirborneTransmitter{A}(n, la, lo, al, an, fr, po) where {A<:Antenna} =
        transmitterchecks(al) && new(n, la, lo, al, an, fr, po)
end

altitude(t::AirborneTransmitter) = t.altitude
frequency(t::AirborneTransmitter) = t.frequency
azimuth(t::AirborneTransmitter) = azimuth(t.antenna)
inclination(t::AirborneTransmitter) = inclination(t.antenna)
power(t::AirborneTransmitter) = t.power


text:
using SafeTestsets
@safetestset "Benchmark Tests" begin include("safe_test.jl") end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IdentityUserAssignedIdentities <: SwaggerModel
    principalId::Any # spec type: Union{ Nothing, String } # spec name: principalId
    clientId::Any # spec type: Union{ Nothing, String } # spec name: clientId

    function IdentityUserAssignedIdentities(;principalId=nothing, clientId=nothing)
        o = new()
        validate_property(IdentityUserAssignedIdentities, Symbol("principalId"), principalId)
        setfield!(o, Symbol("principalId"), principalId)
        validate_property(IdentityUserAssignedIdentities, Symbol("clientId"), clientId)
        setfield!(o, Symbol("clientId"), clientId)
        o
    end
end # type IdentityUserAssignedIdentities

const _property_map_IdentityUserAssignedIdentities = Dict{Symbol,Symbol}(Symbol("principalId")=>Symbol("principalId"), Symbol("clientId")=>Symbol("clientId"))
const _property_types_IdentityUserAssignedIdentities = Dict{Symbol,String}(Symbol("principalId")=>"String", Symbol("clientId")=>"String")
Base.propertynames(::Type{ IdentityUserAssignedIdentities }) = collect(keys(_property_map_IdentityUserAssignedIdentities))
Swagger.property_type(::Type{ IdentityUserAssignedIdentities }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IdentityUserAssignedIdentities[name]))}
Swagger.field_name(::Type{ IdentityUserAssignedIdentities }, property_name::Symbol) =  _property_map_IdentityUserAssignedIdentities[property_name]

function check_required(o::IdentityUserAssignedIdentities)
    true
end

function validate_property(::Type{ IdentityUserAssignedIdentities }, name::Symbol, val)
end


text:
include("./superres.jl")

using Plots
pyplot();

let ntests=15000,seed=0,Ωₘ=2000
    let ℓ=3,R=4,clusterFun=Ζ, amplitudeFun = imag1s
        let pertFun=𝔉ᵣ(Ωₘ,seed)
            let ϵR = 10.0.^LinRange(-14,-5,150), ΩR = NrangeLog(2,3,150), ΔR = 10.0.^LinRange(-4,-2,150)
                let res = mapreduce(transpose,vcat,
                        [runSingle(ℓ,R,clusterFun, pertFun,amplitudeFun, ϵR,ΩR,ΔR) for _ in 1:ntests])
                    let clr = map(cc,res[:,2]), sh = map(ss,res[:,2]),idxsSucc = findall(res[:,2].==1),idxsFail=findall(res[:,2].==0)
                        #p = scatter(log10.(res[:,1]),log10.(inv.(res[:,3])),
                        #    mc=clr,msc=:white,ms=7,shape=sh,markerstrokewidth=0.4,markeralpha=0.5,
                        #    xlabel="log SRF",ylabel="log 1/ϵ",
                        #    color_palette=:magma,lab="")
                        p=scatter(log10.(res[idxsSucc,1]),log10.(inv.(res[idxsSucc,3])),
                            mc=:blue,msc=:white,ms=7,shape=:dtriangle,markerstrokewidth=0.1,markeralpha=0.1,
                            xlabel="log SRF",ylabel="log 1/ϵ",lab="Success",color_palette=:magma)
                        scatter!(log10.(res[idxsFail,1]),log10.(inv.(res[idxsFail,3])),
                            mc=:red,msc=:white,ms=7,shape=:circle,markerstrokewidth=0.1,markeralpha=1.0,
                            xlabel="log SRF",ylabel="log 1/ϵ",lab="Fail",color_palette=:magma)
                        #res[:,2]
                        plot!(t->(2*ℓ-1)*t+4,w=2,l=:dash,lc=:black,
                            lab="\$ SRF^{2p-1}\$")
                        ylims!(5,13)
                        xlims!(0,2)
                        title!("p=$(ℓ),d=$(R), \$ {\\bf S_1} \$")
                        savefig(p,"./figures/fig5-2.pdf")
                    end
                end
            end
        end
    end
end

max_stars_repo_name: batenkov-team/iai-code

max_issues_repo_name: batenkov-team/iai-code

max_forks_repo_name: batenkov-team/iai-code


text:
function parseroman(rnum::AbstractString)
    romandigits = Dict('I' => 1, 'V' => 5, 'X' => 10, 'L' => 50,
                       'C' => 100, 'D' => 500, 'M' => 1000)
    mval = accm = 0
    for d in reverse(uppercase(rnum))
        val = try
            romandigits[d]
        catch
            throw(DomainError())
        end
        if val > mval maxval = val end
        if val < mval
            accm -= val
        else
            accm += val
        end
    end
    return accm
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct SolidifyResult
    sys::System
    eqs::Vector{EquilibResult}
end

function Base.display(res::SolidifyResult)
    N = length(res.eqs)
    N == 0 && return
    Tl = res.eqs[1].T
    Tu = res.eqs[end].T
    println("Solidification:")
    println("\t$Tl ≤ T ≤ $Tu")
    println("\t$N calculations")
end

function Base.print(res::SolidifyResult)
    println("Solidification Result:\n")
    for eq in res.eqs
        @printf "T = %.2f\n" eq.T
        print(eq)
    end
end

function solidify(sys::System, X::Vector{<:Real} = equiatom(sys),
                  Trange::AbstractRange = defaulttemprange(sys))
    eqs = Vector{Calmato.EquilibResult}()
    for T in Trange
        eq = equilib(sys, X, T)
        push!(eqs, eq)
    end
    return SolidifyResult(sys, eqs)
end

function defaulttemprange(sys::System)
    Tl, Tu = sys.temp
    N::Int = div(Tu - Tl, 100)
    return LinRange(Tl, Tl + 100N, N)
end

@recipe function plot(res::SolidifyResult)
    sys = res.sys
    eqs = res.eqs
    K = sys.nphas
    N = length(eqs)

    names = reshape([ sys.phass[k].name for k in 1:K ], 1,K)
    Ts = [ eqs[n].T for n in 1:N ]
    xss = [ [ eqs[n].Y[k] / sum(eqs[n].Y) for n in 1:N ] for k in 1:K ]

    xguide --> "T (K)"
    yguide --> "mole fraction"
    label --> names

    (Ts, xss)
end


text:
using GridPotentialFlow
using Plots
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

R = Lx/8;
b_left = PotentialFlowBody(Circle(R,Δx),Γ=1.0)
T = RigidTransform((-Lx/4,0.0),0.0)
T(b_left);
b_right = PotentialFlowBody(Circle(R,Δx))
T = RigidTransform((Lx/4,0.0),0.0)
T(b_right);

model = VortexModel(g,bodies=[b_left,b_right])
s = streamfunction(model)
plot(s,g)
plot!(b_left,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
plot!(b_right,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)

Δx = 0.004
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

a = Lx/4 # semi-major axis
AR_list = [1.0,2.0,3.0]
ellipses = [Ellipse(a,a/AR,Δx) for AR in AR_list]
bodies = PotentialFlowBody[]
push!(bodies,[PotentialFlowBody(SplinedBody(hcat(e.x,e.y),3*cellsize(g))) for e in ellipses]...)
push!(bodies,PotentialFlowBody(Plate(2*a,3*Δx)));
colors = [:blue,:red,:green,:black]
plot()
for i in 1:length(bodies)
    plot!(bodies[i],fillrange=0,fillalpha=0.0,linecolor=colors[i],linewidth=2)
end
plot!(xlim=xlim,ylim=ylim,xlabel="x",ylabel="y")

f₀_list = []
for i in 1:length(bodies)
    Δs = dlength(bodies[i]);
    model = VortexModel(g,bodies=[bodies[i]]);
    ones = ScalarData(length(bodies[i]))
    ones .= 1.0
    f₀ = model.system.ibp.Sfact\ones
    push!(f₀_list,f₀)
end

append!(f₀_list[end],reverse!(f₀_list[end]));
Δs_plate = dlengthmid(bodies[end]);
append!(Δs_plate,reverse!(Δs_plate));
plot(ylim=(-2,0),xlabel="p/N",title="f₀ for ellipses of different aspect ratios")
for i in 1:length(bodies)-1
    plot!((1:length(f₀_list[i]))/length(f₀_list[i]),f₀_list[i]./dlength(bodies[i]),linecolor=colors[i],label="AR=$(AR_list[i])")
end
plot!((1:length(f₀_list[end]))/length(f₀_list[end]),0.5*f₀_list[end]./Δs_plate,linecolor=colors[end],label="AR=inf")

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
const CTXS = Dict{Symbol, UnionAll}()

"""
    mock(f::Function[, ctx::Symbol], args...; filters::Vector{<:Function}=Function[])

Run `f` with specified functions mocked out.

!!! note
    Mocking functions with keyword arguments is only partially supported.
    See the "Keyword Arguments" section below for more details.

## Examples

Mocking a single function:

```julia
f(x) = x + 1
mock(+) do plus
    @assert plus isa Mock
    @assert f(0) != 1  # The call to + is mocked.
    @assert called_once_with(plus, 0, 1)
end
```

Mocking a function with a custom [`Mock`](@ref):

```julia
mock((+) => Mock(1)) do plus
    @assert 1 + 1 == 1
    @assert called_once_with(plus, 1, 1)
end
```

Mocking methods that match a given signature:

```julia
mock((+, Float64, Float64) => Mock((a, b) -> 2a + 2b)) do plus
    @assert 1 + 1 == 2
    @assert 2.0 + 2.0 == 8
    @assert called_once_with(plus, 2.0, 2.0)
end
```

Mocking with something other than a `Mock`:

```julia
mock((+) => (a, b) -> 2a + 2b) do _plus
    @assert 1 + 2 == 6
end
```

## Using Filters

Oftentimes, you mock a function with a very specific idea of where you want that mocking to happen.
It can be confusing when a call you didn't anticipate gets mocked somewhere deep in the call stack, botching everything.
To avoid this, you can use filter functions like so:

```julia
f(x, y) = x + y
g(x, y) = f(x, y)
mock((+) => Mock((a, b) -> 2a + 2b); filters=[max_depth(2)]) do plus
    @assert f(1, 2) == 6  # The call depth of + here is 2.
    @assert g(3, 4) == 7  # Here, it's 3.
    @assert called_once_with(plus, 1, 2)
end
```

Filter functions take a single argument of type [`Metadata`](@ref).
If any filter rejects, then mocking is not performed.
See [Filter Functions](@ref) for a list of included filters, as well as building blocks for you to create your own.

## Reusing `Context`s

Under the hood, this function creates a new [Cassette `Context`](https://jrevels.github.io/Cassette.jl/stable/api.html#Cassette.Context) on every call by default.
This provides a nice clean mocking environment, but it can be slow to create and call new types and methods over and over.
If you find yourself repeatedly mocking the same set of functions, you can specify a context name to reuse that context like so:

```julia
ctx = gensym()
mock(g -> @assert(!called(g)), ctx, get)
# This one is faster, especially when there's a lot going on in your mock blocks.
mock(g -> @assert(!called(g)), ctx, get)
```

## Keyword Arguments

Mocking of functions with keyword arguments is fully supported when the following conditions are met:

- Filter functions are not used
- The context in use has no previously-mocked functions that are now unmocked

if a filter function rejects and calls a mocked function (instead of its mock), that call will have no keywords.

If you reuse a context that has previously mocked some function, unmocked calls to that function will have no keywords.
For example:

```julia
kwfunc(; kwargs...) = nothing
calls_kwfunc(; kwargs...) = kwfunc(; kwargs...)
ctx = gensym()
mock(ctx, calls_kwfunc) do c
    calls_kwfunc(; x=1, y=2)
    @assert called_once_with(c; x=1, y=2)
end
mock(ctx, kwfunc) do k
    calls_kwfunc(; x=1, y=2)               # This will issue a warning.
    @assert called_once_with(k; x=1, y=2)  # This will fail!
end
```

In short, avoid using filters and reusing contexts when mocking functions that accept keywords.
"""
mock(f::Function, args...; filters::Vector{<:Function}=Function[]) =
    mock(f, gensym(), args...; filters=filters)

function mock(f::Function, ctx::Symbol, args...; filters::Vector{<:Function}=Function[])
    mocks = map(sig2mock, args)  # ((f, sig) => mock).
    isempty(mocks) && throw(ArgumentError("At least one function must be mocked"))

    # Create the new context type if it doesn't already exist.
    ctx_is_new = !haskey(CTXS, ctx)
    ctx_is_new && make_context(ctx)
    Ctx = CTXS[ctx]

    # Implement the overdubs, but only if they aren't already implemented.
    has_new_od = false
    foreach(map(first, mocks)) do k
        fun = k[1]
        sig = k[2:end]
        if ctx_is_new || !overdub_exists(Ctx, fun, sig)
            make_overdub(Ctx, fun, sig)
            has_new_od = true
        end
    end

    # Only use `invokelatest` if the Context/overdub implementations are new.
    meta = Metadata(Dict(mocks), filters)
    c = ctx_is_new ? invokelatest(Ctx; metadata=meta) : Ctx(; metadata=meta)
    od_args = [c, f, map(last, mocks)...]
    return has_new_od ? invokelatest(overdub, od_args...) : overdub(od_args...)
end

# Output (f, sig) => mock.
sig2mock(p::Pair{<:Tuple}) = p
sig2mock(p::Pair) = (p.first, Vararg{Any}) => p.second
sig2mock(t::Tuple) = t => Mock()
sig2mock(f) = (f, Vararg{Any}) => Mock()

# Create a new context type.
make_context(name::Symbol) = @eval begin
    Ctx = CTXS[$(QuoteNode(name))] = @context $(gensym())

    # TODO: Maybe these should be inlined, but it slows down compilation a lot.

    @noinline function Cassette.prehook(ctx::Ctx{Metadata{true}}, f, args...)
        @nospecialize f args
        update!(ctx.metadata, prehook, f, args...)
    end

    @noinline function Cassette.posthook(ctx::Ctx{Metadata{true}}, v, f, args...)
        @nospecialize v f args
        update!(ctx.metadata, posthook, f, args...)
    end
end

# Has a given function and signature already been overdubbed?
function overdub_exists(::Type{Ctx}, ::F, sig::Tuple) where {Ctx <: Context, F}
    squashed = foldl(sig; init=[]) do acc, T
        if T isa DataType && T.name.name === :Vararg
            append!(acc, repeat([T.parameters[1]], T.parameters[2]))
        else
            push!(acc, T)
        end
    end
    return any(m -> m.sig === Tuple{typeof(overdub), Ctx, F, squashed...}, methods(overdub))
end

# Implement `overdub` for a given Context, function, and signature.
function make_overdub(::Type{Ctx}, f::F, sig::Tuple) where {Ctx <: Context, F}
    sig_exs = Expr[]
    sig_names = []

    foreach(sig) do T
        T_uw = unwrap_unionall(T)
        name = gensym()

        if VERSION >= v"1.7" && T_uw isa TypeofVararg
            # T === T_uw.
            if isdefined(T, :T) && isdefined(T, :N)  # Vararg{T, N}.
                foreach(1:T.N) do _i
                    push!(sig_exs, :($name::$(T.T)))
                    push!(sig_names, name)
                    name = gensym()
                end
            else  # Vararg{T, N} where {T, N} or Vararg{T, N} where N.
                T = Vararg{Any}
                push!(sig_exs, :($name::$(T.T)...))
                push!(sig_names, :($name...))
            end
        elseif T_uw.name.name === :Vararg
            if T isa UnionAll && T.body isa UnionAll  # Vararg{T, N} where {T, N}.
                T = Vararg{Any}
            end
            if T isa UnionAll  # Vararg{T, N} where N.
                push!(sig_exs, :($name::$(T_uw.parameters[1])...))
                push!(sig_names, :($name...))
            else  # Vararg{T, N}.
                foreach(1:T.parameters[2]) do _i
                    push!(sig_exs, :($name::$(T.parameters[1])))
                    push!(sig_names, name)
                    name = gensym()
                end
            end
        else
            push!(sig_exs, :($name::$T))
            push!(sig_names, name)
        end
    end

    @eval begin
        @inline function Cassette.overdub(ctx::$Ctx, f::$F, $(sig_exs...); kwargs...)
            method = (f, $(sig...))
            if should_mock(ctx.metadata, method)
                ctx.metadata.mocks[method]($(sig_names...); kwargs...)
            else
                isempty(kwargs) || @warn "Discarding keyword arguments" f kwargs
                recurse(ctx, f, $(sig_names...))
            end
        end

        # https://github.com/jrevels/Cassette.jl/issues/48#issuecomment-440605481
        @inline Cassette.overdub(ctx::$Ctx, ::kwftype($F), kwargs, f::$F, $(sig_exs...)) =
            overdub(ctx, f, $(sig_names...); kwargs...)
    end
end


text:
@testset "Sets constant" begin
    @test MOIU.getconstant(MOI.EqualTo(3)) == 3
    @test MOIU.getconstant(MOI.GreaterThan(6)) == 6
    @test MOIU.getconstant(MOI.LessThan(2)) == 2
end

@testset "Sets shifts" begin
    @test MOIU.shift_constant(MOI.EqualTo(3), 1) == MOI.EqualTo(4)
    @test MOIU.shift_constant(MOI.GreaterThan(6), -1) == MOI.GreaterThan(5)
    @test MOIU.shift_constant(MOI.LessThan(2), 2) == MOI.LessThan(4)
    @test MOIU.shift_constant(MOI.Interval(-2, 3), 1) == MOI.Interval(-1, 4)
end

@testset "Set dimension" begin
    @test MOI.dimension(MOI.EqualTo(3.0)) === 1
    @test MOI.dimension(MOI.Reals(8)) === 8
    @test MOI.dimension(MOI.DualExponentialCone()) === 3
    @test MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(4)) === 10
    @test MOI.dimension(MOI.PositiveSemidefiniteConeSquare(5)) === 25
    @test MOI.dimension(MOI.RootDetConeTriangle(6)) === 22
    @test MOI.dimension(MOI.LogDetConeTriangle(6)) === 22
    @test MOI.dimension(MOI.RootDetConeSquare(4)) === 17
    @test MOI.dimension(MOI.LogDetConeSquare(4)) === 17
    @test MOI.dimension(MOI.SOS2(collect(1 : 6))) === 6
end


text:
# by Lazaro Alonso
# used like a set that contains beta skeletons in : # hide
# Geometrical and spectral study of β-skeleton graphs # hide
# L. Alonso, J. A. Méndez-Bermúdez, and Ernesto Estrada # hide
# Phys. Rev. E 100, 062309 – Published 19 December 2019 # hide

using CairoMakie, Random
using GR: delaunay
using Makie.GeometryBasics: Polygon
CairoMakie.activate!() # hide

let
    Random.seed!(123)
    function pointsDelaunay(nodes::Int64)
        points = rand(nodes, 2)
        n, tri = delaunay(points[:,1], points[:,2])
        return (n, tri, points)
    end
    n, tri, points = pointsDelaunay(2000)
    # vector of polygons
    polys = [Polygon([Point2f0(points[tri[k,:],:][i,:]) for i in 1:3])
        for k in 1:size(tri)[1]]

    #poly(points, tri, strokewidth = 1) # also works, but just one color
    poly(polys, color=rand(size(tri)[1]), strokewidth = 0.5,
        colormap = :linear_worb_100_25_c53_n256,
        figure = (;resolution=(800,800), fontsize = 22),
        axis = (;aspect = DataAspect(), title = "Delaunay triangulation"))
    limits!(0,1,0,1)
    current_figure()
    save(joinpath(@__DIR__, "output", "delaunayTriangulation.png"), current_figure(), px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "GR", "Makie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie