Skip to main content

splines - How to make the result of Piecewise be a closed interval?


Description


This question comes from two questions. Namely Q1 and Q2


The defintion of B-Spline basis function as shown below:


Let $\vec{U}=\{u_0,u_1,\ldots,u_m\}$ a nondecreasing sequence of real numbers,i.e, $u_i\leq u_{i+1}\quad i=0,1,2\ldots m-1$


$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u

  NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),u_] /;
   i <= Length[knots] - 2 :=

   Piecewise[
    {{1, knots[[i + 1]] <= u < knots[[i + 2]]}},0]

  coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
  coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

  NBSpline[i_Integer, p_Integer, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), 
    u_] /;p > 0 && i + p <= Length[knots] - 2 :=
   Module[{init, res},
    init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];

    res = First@
    Nest[
     Dot @@@
      (Thread@
        {Partition[#, 2, 1],
         With[{m = p + 2 - Length@#},
          Table[
           {coeff[u, k + 1, k + m + 1, knots],
            coeff[u, k + m + 2, k + 2, knots]}, {k, i, i + Length@# - 2}]]}) &,
    init, p]

  ]

Compare to built-in function BSplineBasis


 sortResult[x_ /; x == 0] := 0;
 sortResult[res_] := MapAt[SortBy[#, Last] &, res, 1]

 pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
 knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};

Comparsion 1



 NBSpline[#, 0, knots, u] & /@ Range[0, 8] // Simplify


enter image description here



 PiecewiseExpand@BSplineBasis[{0, knots}, #, u] & /@ Range[0, 8]


enter image description here




Comparsion 2


sortResult /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7] // Simplify)


enter image description here



sortResult /@ (PiecewiseExpand@BSplineBasis[{1, knots}, #, u] & /@ 
Range[0, 7] // Simplify)



enter image description here




By the comparsion, I found that the result of BSplineBasis is always a closed interval, whereas the result my function NBSpline is a half closed- half open interval.


My trail:


postProcess[res_ /; res == 0] := 0;

postProcess[res_] :=
 Module[{interval, pos, expr},
  interval =

   First@Simplify@res;
  pos =
   Position[
   SortBy[interval, Last], Less];
  expr =
   ReplacePart[
    SortBy[interval, Last], pos[[-1]] -> LessEqual];
  Piecewise[expr, 0]
]


Test


 postProcess /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7])


enter image description here



Question:


Is there any good/simple method(strategy) to deal with this problem? I think my method is complex.



enter image description here





Answer



a fix, specially treating the end interval:


 NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),
      u_] /; i <= Length[knots] - 2 := (
   If[ knots[[i + 2]] != knots[[-1]],
       Piecewise[{{1, knots[[i + 1]] <= u < knots[[i + 2]]}}, 0],
       Piecewise[{{1, knots[[i + 1]] <= u <= knots[[i + 2]]}}, 0] ]);
 coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
 coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

 NBSpline[i_Integer, p_Integer, 
      knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), u_] /; 
           p > 0 && i + p <= Length[knots] - 2 := 
 Module[{init, res}, 
      init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];
      res = First@
      Nest[Dot @@@ (Thread@{Partition[#, 2, 1], 
          With[{m = p + 2 - Length@#}, 
              Table[{coeff[u, k + 1, k + m + 1, knots], 
              coeff[u, k + m + 2, k + 2, knots]}, {k, i, 

                  i + Length@# - 2}]]}) &, init, p]]
  pts = {{1/10, 2/10}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
  knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};
  BSplinePlot1[pts : {{_, _} ..}, knots_, opts : OptionsPattern[Plot]] :=
  Module[{p = Length@First@Split[knots] - 1, a, b},
        {a, b} = {First[knots], Last[knots]};
    f[u_] = 
        Evaluate@
           Simplify@
               Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 

                      pts] /. x_Piecewise :> Piecewise[x[[1]], Null];
    ParametricPlot[f[u], {u, a, b}, opts, PlotPoints -> 1000]]
    BSplinePlot1[pts, knots, ImageSize -> 600]

this bit: /. x_Piecewise :> Piecewise[x[[1]], Null] fixes the catchall behavior of Piecewise so you will get an error if you try to evaluate out of bounds. (not necessary for this example )


here is f[u] , note that we have closed just the end interval..


enter image description here


note as a practical matter @belisarius approach works as well, so long as you apply at the last stage:


 f[u_] = Evaluate@
          Simplify@

             Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 
                pts] /. Less :> LessEqual /. x_Piecewise :> Piecewise[x[[1]], Null]

enter image description here


since the function is continuous it doesn't really matter if you close all the internal intervals.


Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more