Monotonicity preserving C1 interpolation (continued)

[modelica-trac-importer]

Reported by beutlich on 2 Nov 2015 15:56 UTC In #1717 I added the univariate Fritsch-Butland-spline interpolation to the MSL trunk.

Now I noticed that GSL 2.0 also added support for monotonicity preserving C1 interpolation based on "Steffen, M.; A Simple Method for Monotonic Interpolation in One Dimension; Astronomy and Astrophysics, Vol. 239, NO. NOV(II), P. 443, 1990".

What should be done for MSL 3.2.2?

1. Do nothing and keep Fritsch-Butland as it is now?
2. Replace Fritsch-Butland by Steffen.
3. Keep Fritsch-Butland and also add Steffen.

@haumer @otter Any feedback is appreciated.

---

Migrated-From: https://trac.modelica.org/Modelica/ticket/1814

[modelica-trac-importer]

Modified by beutlich on 2 Nov 2015 15:58 UTC

[modelica-trac-importer]

Comment by beutlich on 23 Nov 2015 21:27 UTC I started a poll about the potential change alternatives to gain more feedback from the MSL shapers and users. Feel free to contribute (by the end of next week).

http://doodle.com/poll/txgxv2i58rs84hh8

[modelica-trac-importer]

Comment by hansolsson on 24 Nov 2015 16:15 UTC I believe the most important question is to understand in what way the different schemes differ.

This is both necessary to make an informed decision in this poll and to present it in a reasonable way for the users.

Looking at http://adsabs.harvard.edu/full/1990A%26A...239..443S it seems that Steffens' variant looks better than Akima - while avoiding some problems with Fritsch-Butland (but I couldn't find any direct comparison between them). However, some of the issues seem related to derivatives at knot-points - which may be an advanced option that we do not use.

So, can anyone answer the following: Is there a problem with Fritsch-Butland for this use-case? If there is a problem with Fritsch-Butland, what is the reason for keeping it (if we only want to keep it for backward compatibility reasons we might use annotation choices to hide it).

[modelica-trac-importer]

Comment by hansolsson on 26 Nov 2015 11:32 UTC Ok, I have now read the paper by Steffen in more detail. I still feel that it seems as if we hand over the decision to users since we haven't understood it ourselves.

Let's see if I understand it correctly - that in addition to Akima-interpolation the following interpolation alternatives exist:

1 Akima 2 Fritsch&Carlsson 1980 3 Fritsch&Butland 1984/Brodlie 1980 4 Steffen 1990

As I understand it 1 is already present, 2 is not considered. The question is if we should add 3 and/or 4. Is that correct?

Akima can produce local extrema inside intervals. Akima method is not always forward stable (i.e. infinitesimal changes in y-values can cause large changes in results) However, Akima interpolation is already used and should thus be kept (even if problematic).

Fritsch&Carlsson potentially has the issues that it 2-pass, and non-local; thus we should not add it. Fritsch&Butland does not match 2nd order polynomials Akima completely matches 2nd order polynomials Steffen completely matches 2nd order polynomials - if the extrema is on a knot-point

The difference between F&B and Steffen seems to be the difference between equation 30 and 11; and F&B seems more elegant.

But is it important to match 2nd order polynomials? If so is there any reason to assume that extrema are on knot-points? Are there any other considerations that we need to consider? E.g. are any of the interpolations linear (and is that important)? (I.e such that interpolation([x,y1+a*y2], t)=interpolation([x,y1],t)+a*interpolation([x,y2],t) ).

[modelica-trac-importer]

Comment by beutlich on 26 Nov 2015 11:53 UTC Replying to [comment:4 hansolsson]:

The question is if we should add 3 and/or 4. Is that correct?

Exactly.

I studied why the GSL implemented Steffen and not F&B. They discussed it: https://www.sourceware.org/ml/gsl-discuss/2014-q1/msg00036.html It looks like GSL preferred Steffen over F&B due to the three reasons Steffen stated in his paper.

I have the patch for Steffen interpolation ready. The differences (both in sense of source and the resulting interpolant) are marginal only. And for that reason I cannot see why Steffen claims his three aspects.

So the difference in slope approximation simply is:

--- fb.c    Do Nov 26 12:51:08 2015
+++ s.c Do Nov 26 13:03:38 2015
@@ -1,13 +1,12 @@
-static CubicHermite1D* fritschButlandSpline1DInit(const double* table,
-                                                  size_t nRow, size_t nCol,
-                                                  const int* cols,
-                                                  size_t nCols) {
+static CubicHermite1D* steffenSpline1DInit(const double* table,
+                                           size_t nRow, size_t nCol,
+                                           const int* cols,
+                                           size_t nCols) {
   /* Reference:

-     Frederick N. Fritsch and Judy Butland. A method for constructing local
-     monotone piecewise cubic interpolants. SIAM Journal on Scientific and
-     Statistical Computing, 5(2), 300-304, June 1984.
-     (http://dx.doi.org/10.1137/0905021)
+     Matthias Steffen. A simple method for monotonic interpolation in one
+     dimension. Astronomy and Astrophysic, 239, 443-450, August 1990.
+     (https://trac.modelica.org/Modelica/attachment/ticket/1814/1990A%2BA___239__443S.pdf)
   */

     CubicHermite1D* spline = NULL;
@@ -55,7 +54,15 @@
             }
             else {
                 const double dx_ = TABLE_COL0(i + 2) - TABLE_COL0(i + 1);
-                c2 = 3*(dx + dx_)/((dx + 2*dx_)/d[i] + (dx_ + 2*dx)/d[i + 1]);
+                double abs_c2_d2, abs_di, abs_di1;
+                c2 = (d[i]*dx_ + d[i + 1]*dx)/(dx + dx_);
+                abs_c2_d2 = 0.5*fabs(c2);
+                abs_di = fabs(d[i]);
+                abs_di1 = fabs(d[i + 1]);
+                if (abs_c2_d2 > abs_di || abs_c2_d2 > abs_di1) {
+                    const double two_a = d[i] > 0 ? 2 : -2;
+                    c2 = two_a*(abs_di < abs_di1 ? abs_di : abs_di1);
+                }
             }
             c[1] = (3*d[i] - 2*c[2] - c2)/dx;
             c[0] = (c[2] + c2 - 2*d[i])/(dx*dx);

[modelica-trac-importer]

Comment by hansolsson on 26 Nov 2015 14:41 UTC Replying to [comment:5 beutlich]:

Replying to [comment:4 hansolsson]:

The question is if we should add 3 and/or 4. Is that correct? Exactly.

I studied why the GSL implemented Steffen and not F&B. They discussed it: https://www.sourceware.org/ml/gsl-discuss/2014-q1/msg00036.html It looks like GSL preferred Steffen over F&B due to the three reasons Steffen stated in his paper.

Now I am getting confused again:

I thought the three reasons was why Steffen is preferred over F&C, but compared to F&B there is only the "doesn't match 2nd order polynomial where extrema is at knotpoint" argument.

[modelica-trac-importer]

Comment by beutlich on 26 Nov 2015 15:08 UTC OK, I always thought that F&C (1980) gives the idea, but lacks a real implementation, which than was added by F&B (1984).

[modelica-trac-importer]

Comment by hansolsson on 26 Nov 2015 15:42 UTC Replying to [comment:7 beutlich]:

OK, I always thought that F&C (1980) gives the idea, but lacks a real implementation, which than was added by F&B (1984).

And I think that F&C (1980) had one idea, which was found wanting and a better solution proposed in F&B (1984).

Steffen (1990) independently found another solution without knowing about these; the reviewers found them for him - and he added references, but in such a way that it could be mis-interpreted as above. (I might be too jaded.)

Obviously this needs to be resolved.

[modelica-trac-importer]

Comment by beutlich on 26 Nov 2015 15:57 UTC I doubt that we'll find a resolution. I also tried to reproduce Steffen's Akima figures with his datasets (reconstructed from the figures) but failed. I neither observed the numerical instabilities nor the oscillations he mentioned.

I also doubt that we'll see some other user's input since most simply do not care in such a detail. See #1153 for an example.

For me there is no major difference between F&B and Steffen and I also cannot judge if one is superior than the other. You said that F&B does not match 2nd order polynomials whereas Steffen does. This might be the only reason to go for one, the other or even both. I mean I can simply add both and let the user have the choice (that we cannot or do not want to take).

[modelica-trac-importer]

Comment by hansolsson on 26 Nov 2015 16:14 UTC Replying to [comment:9 beutlich]:

I doubt that we'll find a resolution. I also tried to reproduce Steffen's Akima figures with his datasets (reconstructed from the figures) but failed. I neither observed the numerical instabilities nor the oscillations he mentioned.

Even more confusing.

I also doubt that we'll see some other user's input since most simply do not care in such a detail. See #1153 for an example.

Well, there are two kinds of user input:

1. I prefer x over y due to z.
2. This is too complicated, can you give me a simpler variant?

We get less of the 2nd kind, but I still see some.

For me there is no major difference between F&B and Steffen and I also cannot judge if one is superior than the other. You said that F&B does not match 2nd order polynomials whereas Steffen does. This might be the only reason to go for one, the other or even both. I mean I can simply add both and let the user have the choice (that we cannot or do not want to take).

Well, even that claim is less complete: Steffen claims his method matches 2nd order polynomials IF the extrema is at a knot-point; otherwise not.

Letting users choose between two different interpolation formulas where we don't know the real difference and we doubt that users care is not good library design to me.

[modelica-trac-importer]

Comment by beutlich on 26 Nov 2015 20:03 UTC Replying to [comment:10 hansolsson]:

Even more confusing.

Steffen mentions Akima (1970) that also is used as reference for the MSL tables (in function akimaSpline1DInit). Somewhere I read that there is an improved Akima variant. That might be one reason. For sake of completeness, this is the (reconstructed) data set of figure 3: [0,8.5;1,8;2,7.5;3,7;4,6.5;4.25,3.25;4.5,0.5;7,0.5;8,0.5;9,0.5;10,0.5]

Letting users choose between two different interpolation formulas where we don't know the real difference and we doubt that users care is not good library design to me.

You are right, to be honest.

I read Steffen's conclusion again. The paragraph after formula (30) shows the differences he observed with respect to F&B (1984). He states that the approximated slopes of F&B are always lower in comparison to his algorithm. And the difference can be up to 25%.

And on page 2 Steffen states the his algorithm gives exact results if the data points are derived from a 2nd order polynomial (which is not true for F&B/Brodlie).

Thus, if we won't let the user decide between F&B or Steffen we could reason that we chose Steffen over F&B for exactly this point.

[modelica-trac-importer]

Comment by beutlich on 30 Nov 2015 12:37 UTC In http://doodle.com/poll/txgxv2i58rs84hh8 there are now three votes to have both F&B and Steffen. Hans, would you agree on replacing F&B by Steffen?

[modelica-trac-importer]

Comment by otter on 10 Dec 2015 10:54 UTC Interpolation has always the difficulty that it is not possible to know which is suitable/best for the problem at hand. The only way is to concretly try it out with the concrete task and then compare what looks best for this problem at hand. Therefore, if there are already two implementations, I would include them both in MSL.

[modelica-trac-importer]

Modified by beutlich on 10 Dec 2015 11:12 UTC

[modelica-trac-importer]

Comment by beutlich on 16 Dec 2015 13:07 UTC Option 3. resolved by 7631ba7a0c5dab286bb065d198bd36fd5bfb1274.

CombiTable2D gives a not-implemented error for smoothness=MonotoneContinuousDerivative1 or smoothness=MonotoneContinuousDerivative2.

[modelica-trac-importer]

Changelog modified by beutlich on 16 Dec 2015 13:07 UTC Added univariate Steffen-spline interpolation (for CombiTimeTable, CombiTable1D and CombiTable1Ds)