add DAETS integrator

lib/OrdinaryDiffEqTaylorSeries/LICENSE.md

@@ -0,0 +1,24 @@
+The OrdinaryDiffEq.jl package is licensed under the MIT "Expat" License:
+
+> Copyright (c) 2016-2020: ChrisRackauckas, Yingbo Ma, Julia Computing Inc, and
+> other contributors:
+> 
+> https://github.com/SciML/OrdinaryDiffEq.jl/graphs/contributors
+> 
+> Permission is hereby granted, free of charge, to any person obtaining a copy
+> of this software and associated documentation files (the "Software"), to deal
+> in the Software without restriction, including without limitation the rights
+> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+> copies of the Software, and to permit persons to whom the Software is
+> furnished to do so, subject to the following conditions:
+> 
+> The above copyright notice and this permission notice shall be included in all
+> copies or substantial portions of the Software.
+> 
+> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+> SOFTWARE.

lib/OrdinaryDiffEqTaylorSeries/Project.toml

@@ -0,0 +1,59 @@
+name = "OrdinaryDiffEqTaylorSeries"
+uuid = "9c7f1690-dd92-42a3-8318-297ee24d8d39"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.1.0"
+
+[deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+Hungarian = "e91730f6-4275-51fb-a7a0-7064cfbd3b39"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
+
+[compat]
+DiffEqBase = "6.152.2"
+DiffEqDevTools = "2.44.4"
+FastBroadcast = "0.3.5"
+GLPK = "1.2.1"
+Hungarian = "0.7.0"
+JuMP = "1.24.0"
+LinearAlgebra = "<0.0.1, 1"
+MuladdMacro = "0.2.4"
+OrdinaryDiffEqCore = "1.1"
+PrecompileTools = "1.2.1"
+Preferences = "1.4.3"
+Random = "<0.0.1, 1"
+RecursiveArrayTools = "3.27.0"
+Reexport = "1.2.2"
+SafeTestsets = "0.1.0"
+SciMLBase = "2.72.2"
+Static = "1.1.1"
+SymbolicUtils = "3.15.0"
+Symbolics = "6.28.0"
+TaylorDiff = "0.3.1"
+Test = "<0.0.1, 1"
+TruncatedStacktraces = "1.4.0"
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]

lib/OrdinaryDiffEqTaylorSeries/src/DAETS_utils.jl

@@ -0,0 +1,166 @@
+using Symbolics
+using SymbolicUtils
+using Hungarian
+using JuMP
+using GLPK
+
+export compute_signature_matrix, compute_hvt, compute_offsets, compute_jacobian
+"""
+    signature_matrix(eqs, vars, t_ind) -> Matrix{Float64}
+
+Construct the signature matrix Σ for the system of equations `eqs` with respect to
+the variables `vars`. Each entry Σ[i,j] is the (highest) derivative order of `vars[j]`
+appearing in eqs[i], or -Inf if it does not appear.
+"""
+function signature_matrix(eqs::Vector, vars::Vector, t_ind)
+    # eqs[i] is something like f_i(x₁(t), x₂'(t), x₃''(t), ...)
+    n_eqs   = length(eqs)
+    n_vars  = length(vars)
+    Σ = fill(-Inf, n_eqs, n_vars) # Initialize with -Inf
+
+    for i in 1:n_eqs
+        # For each equation
+        for j in 1:n_vars
+            # Check for each variable in the equation what the highest derivative order is
+            order = max_derivative_order(eqs[i], vars[j], t_ind)
+            Σ[i,j] = order
+        end
+    end
+    return Σ
+end
+"""
+    max_derivative_order(ex, var, t_ind)
+
+Returns the highest derivative order of `var` that appears in the symbolic expression `ex`,
+using `t_ind` as the independent variable (e.g., time). If `var` does not appear, returns -Inf.
+"""
+function max_derivative_order(ex, var, t_ind)
+    # Base case: if ex is exactly var(t_ind), order is 0.
+    if isequal(ex, var(t_ind))
+        return 0
+    end
+
+    # If it's a number or an unrelated symbol we ignore
+    if ex isa Number || (ex isa Symbol && ex != var)
+        return -Inf
+    end
+
+    # Check if ex is a derivative expression.
+    if iscall(ex) && operation(ex) isa Symbolics.Differential
+        inner = arguments(ex)[1]  # The function being differentiated.
+        # Recursively check the order of the inner expression.
+        sub_order = max_derivative_order(inner, var, t_ind)
+        # If the inner expression is not related to var, return -Inf otherwise add 1 to derivative order.
+        return sub_order == -Inf ? -Inf : 1 + sub_order
+    end
+
+    # For composite expressions (e.g., sums, products), traverse their arguments.
+    if iscall(ex)
+        best = -Inf
+        for arg in arguments(ex)
+            # Recursively check the order of each component
+            best = max(best, max_derivative_order(arg, var, t_ind))
+        end
+        return best
+    end
+    return -Inf
+end
+
+"""
+    highest_value_transversal(Σ) -> (Vector{Tuple{Int, Int}}, Float64)
+
+Finds the highest value transversal (HVT) of the signature matrix `Σ` using the Hungarian algorithm.
+Returns the transversal as a vector of (row, column) indices and its value.
+"""
+function highest_value_transversal(Σ::Matrix{Float64})
+    n = size(Σ, 1)
+
+    # The Hungarian algorithm minimizes so multiply by -1 to max
+    cost_matrix = -Σ
+    assignment = hungarian(cost_matrix)[1]
+    # Extract transversal and its value.
+    transversal = [(i, assignment[i]) for i in 1:n]
+    value = sum(Σ[i, assignment[i]] for i in 1:n)
+
+    return transversal, value
+end
+
+
+"""
+    find_offsets(Σ::Matrix{Float64}, T::Vector{Tuple{Int, Int}}) -> (Vector{Int}, Vector{Int})
+
+Finds the canonical offsets `c` and `d` for the signature matrix `Σ` and the highest value transversal `T`.
+Returns the vectors `c` and `d`.
+"""
+function find_offsets(Σ::Matrix{Float64}, T::Vector{Tuple{Int, Int}})
+    n = size(Σ, 1)
+
+    # Create JuMP model with GLPK solver
+    model = Model(GLPK.Optimizer)
+
+    # Define variables for c and d
+    @variable(model, c[1:n] >= 0, Int)
+    @variable(model, d[1:n] >= 0, Int)
+
+    # Add constraints for all i, j: d[j] - c[i] >= Σ[i, j]
+    for i in 1:n
+        for j in 1:n
+            if Σ[i, j] != -Inf
+                @constraint(model, d[j] - c[i] >= Σ[i, j])
+            end
+        end
+    end
+
+    # Add constraints for equality over transversal
+    for (i, j) in T
+        @constraint(model, d[j] - c[i] == Σ[i, j])
+    end
+
+    # min sum c and d to find canonical offsets
+    @objective(model, Min, sum(c) + sum(d))
+    optimize!(model)
+    c_values = value.(c)
+    d_values = value.(d)
+    return c_values, d_values
+end
+
+
+"""
+    system_jacobian(eqs::Vector{SymbolicUtils.BasicSymbolic{Number}}, vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple{Number}, Number, Nothing}}}, t_ind::SymbolicUtils.BasicSymbolic{Number}, c::Vector{Int}, d::Vector{Int}, Σ::Matrix{Float64}) -> Matrix{SymbolicUtils.BasicSymbolic{Number}}
+
+Constructs the System Jacobian matrix J for the system of equations `eqs` with respect to the variables `vars`.
+The offsets `c` and `d` and the signature matrix `Σ` are used to determine the structure of the Jacobian.
+"""
+function system_jacobian(
+    eqs::Vector{SymbolicUtils.BasicSymbolic{Number}},
+    vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple{Number}, Number, Nothing}}},
+    t_ind::SymbolicUtils.BasicSymbolic{Number},
+    c::Vector{Int},
+    d::Vector{Int},
+    Σ::Matrix{Float64}
+)
+    n = length(eqs)
+    J = zeros(Symbolics.Num, n, n)
+
+    for i in 1:n
+        for j in 1:n
+            if d[j] - c[i] == Σ[i, j]
+                f_i = eqs[i]
+                x_j = vars[j]
+                σ_ij = Int(Σ[i, j])
+                # Compute the σ_ij-th derivative of x_j
+                x_j_deriv = x_j(t_ind)
+                for _ in 1:σ_ij
+                    x_j_deriv = Differential(t_ind)(x_j_deriv)
+                end
+                # Compute the partial derivative ∂f_i / ∂x_j^(σ_ij)
+                J[i, j] = expand_derivatives(Differential(x_j_deriv)(f_i))
+            else
+                # Set J[i, j] = 0 if d[j] - c[i] != Σ[i, j]
+                J[i, j] = 0
+            end
+        end
+    end
+
+    return J
+end

lib/OrdinaryDiffEqTaylorSeries/src/OrdinaryDiffEqTaylorSeries.jl

@@ -0,0 +1,63 @@
+module OrdinaryDiffEqTaylorSeries
+
+import OrdinaryDiffEqCore: alg_order, alg_stability_size, explicit_rk_docstring,
+                           OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqMutableCache,
+                           alg_cache,
+                           OrdinaryDiffEqConstantCache, @fold, trivial_limiter!,
+                           constvalue, @unpack, perform_step!, calculate_residuals, @cache,
+                           calculate_residuals!, _ode_interpolant, _ode_interpolant!,
+                           CompiledFloats, @OnDemandTableauExtract, initialize!,
+                           perform_step!, OrdinaryDiffEqAlgorithm,
+                           CompositeAlgorithm, _ode_addsteps!, copyat_or_push!,
+                           AutoAlgSwitch, get_fsalfirstlast,
+                           full_cache, DerivativeOrderNotPossibleError
+import Static: False
+import MuladdMacro: @muladd
+import FastBroadcast: @..
+import RecursiveArrayTools: recursivefill!, recursive_unitless_bottom_eltype
+import LinearAlgebra: norm
+using TruncatedStacktraces
+using TaylorDiff
+import DiffEqBase: @def
+import OrdinaryDiffEqCore
+
+using Reexport
+@reexport using DiffEqBase
+
+include("DAETS_utils.jl")
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("TaylorSeries_caches.jl")
+include("interp_func.jl")
+# include("interpolants.jl")
+include("TaylorSeries_perform_step.jl")
+
+import PrecompileTools
+import Preferences
+
+PrecompileTools.@compile_workload begin
+    lorenz = OrdinaryDiffEqCore.lorenz
+    lorenz_oop = OrdinaryDiffEqCore.lorenz_oop
+    solver_list = [ExplicitTaylor2()]
+    prob_list = []
+
+    if Preferences.@load_preference("PrecompileNoSpecialize", false)
+        push!(prob_list,
+            ODEProblem{true, SciMLBase.NoSpecialize}(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0)))
+        push!(prob_list,
+            ODEProblem{true, SciMLBase.NoSpecialize}(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0),
+                Float64[]))
+    end
+
+    for prob in prob_list, solver in solver_list
+        solve(prob, solver)(5.0)
+    end
+
+    prob_list = nothing
+    solver_list = nothing
+end
+
+export ExplicitTaylor2, ExplicitTaylor, DAETS
+
+end