simplfy cubic splines 1

misc/test_truncdiff.jl

@@ -0,0 +1,49 @@
+using Dolo: splines
+
+
+ranges = ((0.0, 1.0, 10), (0.0, 1.0, 20))
+methods(splines.interp)
+
+using StaticArrays
+
+x0 = SVector( 10,20)
+
+
+values = rand(typeof(x0), 10,20)
+
+values2 = rand(typeof(x0), 12,22)
+
+using ForwardDiff
+
+
+fun = u->splines.interp(ranges, values, u)
+
+ForwardDiff.jacobian(
+    fun,
+    x0
+)
+
+a = tuple((e[1] for e in ranges)...)
+b = tuple((e[2] for e in ranges)...)
+orders = tuple((e[3] for e in ranges)...)
+
+fun2 = u->splines.eval_UC_spline(a,b,orders, values2, u)
+fun2(x0)
+
+
+
+
+ForwardDiff.jacobian(
+    fun2,
+    x0
+)
+
+# import ForwardDiff: unsafe_trunc
+
+# function ForwardDiff.unsafe_trunc(TI, number::ForwardDiff.Dual{Tag, Tf}) where Tag where Tf<:Union{Float32, Float64}
+#     ForwardDiff.Dual{Tag}(
+#         unsafe_trunc(TI, ForwardDiff.value(number)),
+#         0*ForwardDiff.partials(number)
+#     )
+
+# end

src/splines/macros.jl

@@ -79,8 +79,9 @@ function create_Phi(d, extrap, diff; Tf=Float64)
         block = []
         rhs_1 = U("tp", i,1)
         rhs_2 = U("tp", i,2)
-        rhs_4 = U("tp", i,4)
         rhs_3 = U("tp", i,3)
+        rhs_4 = U("tp", i,4)
+
         if extrap == "none"
             for j=1:4
                 eq = :($(U("Phi_",i,j)) = ($(Ad[j,1])*$rhs_1 + $(Ad[j,2])*$rhs_2 + $(Ad[j,3])*$rhs_3 + $(Ad[j,4])*$rhs_4) )
@@ -221,6 +222,7 @@ function create_function(d,extrap="natural"; vectorize=true, Tf=Float64)
                 $(create_parameters(d; Tf=Tf)...)
                 N = size(S,1)
                 # @fastmath @inbounds @simd( # doesn't seem to make any difference
+
                 for n=1:N
                     $(create_local_parameters(d;Tf=Tf)...)
                     $(create_Phi(d,extrap,false;Tf=Tf)...)