Return to Question

Lately, we had this thread about interpolation where J. M. linked two interpolation methods. The background for my question is that I estimated a parameter in polar coordinates with dependence on the angle. Therefore, I have a set of values where I know that the values are periodic.

Lately, we had this thread about interpolation where J. M. linked two interpolation methods. The background for my question is that I estimated a parameter in polar coordinates with dependence on the angle. Therefore, I have a set of values where I know that the values are periodic.

pp = Dot @@@ Transpose[Map[Partition[Reverse[#], 2, 1, {-1, 1}] &, {h, del}]]/
  ListConvolve[{1, 1}, h, -1*{-1, 1}];

Options[IPCUMonotonicInterpolation] := {
  PeriodicInterpolation -> False
  }

steffenEnds[{{h1_, h2_}, {d1_, d2_}}] := 
 With[{p = d1 + h1 (d1 - d2)/(h1 + h2)}, (Sign[p] + Sign[d1]) Min[
    Abs[p]/2, Abs[d1]]]
    
IPCUMonotonicInterpolation[data_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
    IPCUMonotonicInterpolation[Transpose[{Range[Length[data]], data}], opts];

IPCUMonotonicInterpolation[data_?MatrixQ, OptionsPattern[]] := 
  Module[{dTrans = Transpose[data], del, h, m, pp, optPeriodic, overhangs},
  optPeriodic = OptionValue[PeriodicInterpolation];
  h = Differences[First[dTrans]]; 
  del = Differences[Last[dTrans]]/h;
  overhangs = If[optPeriodic === False, {1, -1}, {-1, 1}];
  (* Note that overhangs in Partition and ListConvolve are defined differently*)
  pp = Dot @@@ 
     Transpose[
      Map[Partition[Reverse[#], 2, 1, overhangs] &, {h, del}]]/
    ListConvolve[{1, 1}, h, -1*overhangs];
  If[optPeriodic === True,
   del = ArrayPad[del, 1, "Periodic"]
   ];
  m = ListConvolve[{1, 1}, 2 UnitStep[del] - 1] *
                MapThread[Min, {Partition[Abs[del], 2, 1], Abs[pp]/2}];
  Interpolation[
   {{#1}, ##2} & @@@ Transpose[Append[dTrans,
      If[optPeriodic === True,
       m,
       Flatten[{
         steffenEnds[#[[{1, 2}]] & /@ {h, del}],
         m,
         steffenEnds[#[[{-1, -2}]] & /@ {h, del}]
         }]
       ]
      ]],
   PeriodicInterpolation -> optPeriodic]
  ]

pp = Dot @@@ Transpose[ MapAt[Reverse, 
  Map[Partition[#, 2, 1, {-1, 1}] &, {h, del}], {1, All}]]/
  ListConvolve[{1, 1}, h, {1, -1}];

Options[IPCUMonotonicInterpolation] := {
  PeriodicInterpolation -> False
  }

steffenEnds[{{h1_, h2_}, {d1_, d2_}}] := 
 With[{p = d1 + h1 (d1 - d2)/(h1 + h2)}, (Sign[p] + Sign[d1]) Min[
    Abs[p]/2, Abs[d1]]]
    
IPCUMonotonicInterpolation[data_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
    IPCUMonotonicInterpolation[Transpose[{Range[Length[data]], data}], opts];

IPCUMonotonicInterpolation[data_?MatrixQ, OptionsPattern[]] := 
  Module[{dTrans = Transpose[data], del, h, m, pp, optPeriodic, overhangs},
  optPeriodic = OptionValue[PeriodicInterpolation];
  h = Differences[First[dTrans]]; 
  del = Differences[Last[dTrans]]/h;
  overhangs = If[optPeriodic === False, {1, -1}, {-1, 1}];
  (* Note that overhangs in Partition and ListConvolve are defined differently*)
  pp = Dot @@@ Transpose[MapAt[Reverse, 
    Map[Partition[#, 2, 1, overhangs] &, {h, del}], {1, All}]]/
    ListConvolve[{1, 1}, h, -1*overhangs];
  If[optPeriodic === True,
   del = ArrayPad[del, 1, "Periodic"]
   ];
  m = ListConvolve[{1, 1}, 2 UnitStep[del] - 1] *
                MapThread[Min, {Partition[Abs[del], 2, 1], Abs[pp]/2}];
  Interpolation[
   {{#1}, ##2} & @@@ Transpose[Append[dTrans,
      If[optPeriodic === True,
       m,
       Flatten[{
         steffenEnds[#[[{1, 2}]] & /@ {h, del}],
         m,
         steffenEnds[#[[{-1, -2}]] & /@ {h, del}]
         }]
       ]
      ]],
   PeriodicInterpolation -> optPeriodic]
  ]

Final solution

I'm accepting the Steffen interpolation in J.M.'s answer as solution for the following reason:

It shows nicely how one can supply derivative values and not only the interpolation values to Interpolation. It therefore does not implement a whole new interpolation but it only calculates adjusted derivatives and uses the internal Interpolation of Mathematica.

Note, that I made some changes in his function. First, the following part

pp = Apply[Dot, Transpose[MapAt[Map[Reverse, #] &,
  Map[Partition[#, 2, 1, {-1, 1}] &, {h, del}], 1]], 1]/
  ListConvolve[{1, 1}, h, {1, -1}];

can be (IMO) slightly simplified to

pp = Dot @@@ Transpose[Map[Partition[Reverse[#], 2, 1, {-1, 1}] &, {h, del}]]/
  ListConvolve[{1, 1}, h, -1*{-1, 1}];

which saves a whole pure function and a Map. Otherwise, I combined the periodic and the non-periodic interpolation and added a pattern so that data in the form {y1,y2,...} can be interpolated in the usual way. (Please change the function name. It's only called like that since I included it in a package)

Options[IPCUMonotonicInterpolation] := {
  PeriodicInterpolation -> False
  }

steffenEnds[{{h1_, h2_}, {d1_, d2_}}] := 
 With[{p = d1 + h1 (d1 - d2)/(h1 + h2)}, (Sign[p] + Sign[d1]) Min[
    Abs[p]/2, Abs[d1]]]
    
IPCUMonotonicInterpolation[data_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
    IPCUMonotonicInterpolation[Transpose[{Range[Length[data]], data}], opts];

IPCUMonotonicInterpolation[data_?MatrixQ, OptionsPattern[]] := 
  Module[{dTrans = Transpose[data], del, h, m, pp, optPeriodic, overhangs},
  optPeriodic = OptionValue[PeriodicInterpolation];
  h = Differences[First[dTrans]]; 
  del = Differences[Last[dTrans]]/h;
  overhangs = If[optPeriodic === False, {1, -1}, {-1, 1}];
  (* Note that overhangs in Partition and ListConvolve are defined differently*)
  pp = Dot @@@ 
     Transpose[
      Map[Partition[Reverse[#], 2, 1, overhangs] &, {h, del}]]/
    ListConvolve[{1, 1}, h, -1*overhangs];
  If[optPeriodic === True,
   del = ArrayPad[del, 1, "Periodic"]
   ];
  m = ListConvolve[{1, 1}, 2 UnitStep[del] - 1] *
                MapThread[Min, {Partition[Abs[del], 2, 1], Abs[pp]/2}];
  Interpolation[
   {{#1}, ##2} & @@@ Transpose[Append[dTrans,
      If[optPeriodic === True,
       m,
       Flatten[{
         steffenEnds[#[[{1, 2}]] & /@ {h, del}],
         m,
         steffenEnds[#[[{-1, -2}]] & /@ {h, del}]
         }]
       ]
      ]],
   PeriodicInterpolation -> optPeriodic]
  ]