Merge branch 'master' into compathelper/new_version/2024-06-28-00-23-27-224-02392598976

Author: ChrisRackauckas

benchmarks/AutomaticDifferentiation/Project.toml

@@ -28,7 +28,7 @@ BenchmarkTools = "1"
 DataFrames = "1"
 DataFramesMeta = "0.15"
 DifferentiationInterface = "0.2"
-DifferentiationInterfaceTest = "0.2"
+DifferentiationInterfaceTest = "0.2, 0.5"
 Enzyme = "0.11, 0.12"
 FiniteDiff = "2"
 ForwardDiff = "0.10"

benchmarks/AutomaticDifferentiationSparse/Project.toml

@@ -15,3 +15,11 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 SparseDiffTools = "2"
+Symbolics = "5"
+DifferentiationInterface = "0.5"
+DataFrames = "1"
+SparseConnectivityTracer = "0.5"
+ForwardDiff = "0.10"
+BenchmarkTools = "1"
+SparseMatrixColorings = "0.3"
+Zygote = "0.6"

benchmarks/ModelingToolkit/Multibody_Robot.jmd

@@ -4,6 +4,7 @@ author: Fredrik Bagge Carlson
 ---
 
 This is a Multibody.jl model of a 6-axis industrial robot with 6 revolute joints. The first three joints are modeled as flexible. The robot is tracking a smooth, predefined reference trajectory.
+
 ```julia
 cd(@__DIR__)
 using Pkg
@@ -16,58 +17,108 @@ using JuliaSimCompiler
 using OrdinaryDiffEq
 using CairoMakie
 using Printf
+```
+
+## Setup
 
+```julia
 time_instantiate = @elapsed @named robot = Robot6DOF()
 robot = complete(robot)
 time_simplify = @elapsed ssys = structural_simplify(IRSystem(robot))
-time_prob = @elapsed prob = ODEProblem(ssys, [
-    robot.mechanics.r1.phi => deg2rad(-60)
-    robot.mechanics.r2.phi => deg2rad(20)
-    robot.mechanics.r3.phi => deg2rad(90)
-    robot.mechanics.r4.phi => deg2rad(0)
-    robot.mechanics.r5.phi => deg2rad(-110)
-    robot.mechanics.r6.phi => deg2rad(0)
-    robot.axis1.motor.Jmotor.phi => deg2rad(-60) * (-105) 
-    robot.axis2.motor.Jmotor.phi => deg2rad(20) * (210)
-    robot.axis3.motor.Jmotor.phi => deg2rad(90) * (60)
-], (0.0, 4.0))
-time_solve = @elapsed sol = solve(prob, Rodas5P(autodiff=false)); # With autodiff=true this takes over 150 seconds
-
-tv = 0:0.01:4
-time_extract_data = @elapsed data = sol(tv, idxs=[
-    robot.pathPlanning.controlBus.axisControlBus1.angle
-    robot.pathPlanning.controlBus.axisControlBus2.angle
-    robot.pathPlanning.controlBus.axisControlBus3.angle
-    robot.pathPlanning.controlBus.axisControlBus4.angle
-    robot.pathPlanning.controlBus.axisControlBus5.angle
-    robot.pathPlanning.controlBus.axisControlBus6.angle
-    robot.mechanics.r1.phi
-    robot.mechanics.r2.phi
-    robot.mechanics.r3.phi
-    robot.mechanics.r4.phi
-    robot.mechanics.r5.phi
-    robot.mechanics.r6.phi
-    robot.axis1.controller.feedback1.output.u
-    robot.axis2.controller.feedback1.output.u
-    robot.axis3.controller.feedback1.output.u
-    robot.axis4.controller.feedback1.output.u
-    robot.axis5.controller.feedback1.output.u
-    robot.axis6.controller.feedback1.output.u
-]);
-
+const C = JuliaSimCompiler.CTarget();
+const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 labels = ["Instantiate", "Simplify", "Problem creation", "Solve", "Extract data"]
-timings = [time_instantiate, time_simplify, time_prob, time_solve, time_extract_data]
+const u0 = [
+        robot.mechanics.r1.phi => deg2rad(-60)
+        robot.mechanics.r2.phi => deg2rad(20)
+        robot.mechanics.r3.phi => deg2rad(90)
+        robot.mechanics.r4.phi => deg2rad(0)
+        robot.mechanics.r5.phi => deg2rad(-110)
+        robot.mechanics.r6.phi => deg2rad(0)
+        robot.axis1.motor.Jmotor.phi => deg2rad(-60) * (-105) 
+        robot.axis2.motor.Jmotor.phi => deg2rad(20) * (210)
+        robot.axis3.motor.Jmotor.phi => deg2rad(90) * (60)
+    ]
+const tspan = (0.0, 4.0)
+```
+
+## Backend Benchmark Definition
+
+```julia
+function julia_multibody_timings(ssys, target)
+    time_prob = @elapsed prob = if target === JuliaSimCompiler.JuliaTarget()
+        ODEProblem(ssys, u0, tspan)
+    else
+        ODEProblem(ssys, target, u0, tspan)
+    end
+    time_solve = @elapsed sol = solve(prob, Rodas5P(autodiff=false)); # With autodiff=true this takes over 150 seconds
+    
+    tv = 0:0.01:4
+    time_extract_data = @elapsed data = sol(tv, idxs=[
+        robot.pathPlanning.controlBus.axisControlBus1.angle
+        robot.pathPlanning.controlBus.axisControlBus2.angle
+        robot.pathPlanning.controlBus.axisControlBus3.angle
+        robot.pathPlanning.controlBus.axisControlBus4.angle
+        robot.pathPlanning.controlBus.axisControlBus5.angle
+        robot.pathPlanning.controlBus.axisControlBus6.angle
+        robot.mechanics.r1.phi
+        robot.mechanics.r2.phi
+        robot.mechanics.r3.phi
+        robot.mechanics.r4.phi
+        robot.mechanics.r5.phi
+        robot.mechanics.r6.phi
+        robot.axis1.controller.feedback1.output.u
+        robot.axis2.controller.feedback1.output.u
+        robot.axis3.controller.feedback1.output.u
+        robot.axis4.controller.feedback1.output.u
+        robot.axis5.controller.feedback1.output.u
+        robot.axis6.controller.feedback1.output.u
+    ]);
+    timings = [time_instantiate, time_simplify, time_prob, time_solve, time_extract_data]
+end
+```
+
+```julia
+julia_timings = julia_multibody_timings(ssys, JuliaSimCompiler.JuliaTarget())
+
 f = Figure()
 xs = 1:length(labels)
-points = Makie.Point.(xs .- 0.25, timings .+ 1)
-anns = [@sprintf("%3.3g", t) for t in timings]
-barplot(f[1,1], timings, axis=(; title="Timings", xticks = (xs, labels), limits = (nothing, (0, maximum(timings)*1.2 + 2))))
+points = Makie.Point.(xs .- 0.25, julia_timings .+ 1)
+anns = [@sprintf("%3.3g", t) for t in julia_timings]
+barplot(f[1,1], julia_timings, axis=(; title="Julia Backend Timings", xticks = (xs, labels), limits = (nothing, (0, maximum(julia_timings)*1.2 + 2))))
+annotations!(f[1,1], anns, points)
+barplot(f[2,1], julia_timings, axis=(; title="Log Timings", xticks = (xs, labels), yscale=log10))
+f
+```
 
+```julia
+llvm_timings = julia_multibody_timings(ssys, LLVM)
+
+f = Figure()
+xs = 1:length(labels)
+points = Makie.Point.(xs .- 0.25, llvm_timings .+ 1)
+anns = [@sprintf("%3.3g", t) for t in llvm_timings]
+barplot(f[1,1], llvm_timings, axis=(; title="LLVM Backend Timings", xticks = (xs, labels), limits = (nothing, (0, maximum(llvm_timings)*1.2 + 2))))
+annotations!(f[1,1], anns, points)
+barplot(f[2,1], llvm_timings, axis=(; title="Log Timings", xticks = (xs, labels), yscale=log10))
+f
+```
+
+```julia
+c_timings = julia_multibody_timings(ssys, C)
+
+f = Figure()
+xs = 1:length(labels)
+points = Makie.Point.(xs .- 0.25, c_timings .+ 1)
+anns = [@sprintf("%3.3g", t) for t in c_timings]
+barplot(f[1,1], c_timings, axis=(; title="C Backend Timings", xticks = (xs, labels), limits = (nothing, (0, maximum(c_timings)*1.2 + 2))))
 annotations!(f[1,1], anns, points)
-barplot(f[2,1], timings, axis=(; title="Log Timings", xticks = (xs, labels), yscale=log10))
+barplot(f[2,1], c_timings, axis=(; title="Log Timings", xticks = (xs, labels), yscale=log10))
 f
 ```
 
+Since the C backend currently errors, it's omitted from the total timings.
+
 ## OpenModelica
 
 ```julia
@@ -91,26 +142,37 @@ OMJulia.quit(mod)
 ```julia
 colors = Makie.wong_colors()
 labels = ["Build", "Simulate", "Total"]
-julia_build = time_instantiate + time_simplify + time_prob
-julia_sim = time_solve + time_extract_data
+julia_build = sum(julia_timings[1:3])
+julia_sim = sum(julia_timings[4:5])
 julia_total = julia_build + julia_sim
+
+llvm_build = sum(llvm_timings[1:3])
+llvm_sim = sum(llvm_timings[4:5])
+llvm_total = llvm_build + llvm_sim
+
+# c_build = sum(c_timings[1:3])
+# c_sim = sum(c_timings[4:5])
+# c_total = c_build + c_sim
+
 dymola_build = sum([5.075, 3.912, 4.024])/3
 dymola_total = sum([5.267, 4.112, 4.255])/3
 dymola_sim = dymola_total - dymola_build
 
 data = [
     julia_build julia_sim julia_total
+    llvm_build llvm_sim llvm_total
+    # c_build c_sim c_total
     om_build om_sim om_total
     dymola_build dymola_sim dymola_total
 ]
 
-xs = repeat(1:length(labels), inner=3)
-group = repeat([1,2,3], outer=3)
+xs = repeat(1:length(labels), inner=4)
+group = repeat([1,2,3,4], outer=3)
 fig = Figure()
 barplot(fig[1,1], xs, vec(data), dodge=group, color=colors[group], axis=(; title="Timings", xticks = ([1,2,3], labels)))
 
 # Legend
-legendentries = ["Julia", "OpenModelica", "Dymola"]
+legendentries = ["Julia Backend", "LLVM Backend", #= "C Backend", =# "OpenModelica", "Dymola"]
 elements = [PolyElement(polycolor = colors[i]) for i in 1:length(legendentries)]
 title = "Contestants"
 

benchmarks/ModelingToolkit/Project.toml

@@ -26,11 +26,11 @@ JuliaSimCompiler = "0.1.9"
 JuliaSimCompilerRuntime = "1.0.0"
 LinearSolve = "2.30.0"
 Multibody = "0.1.0"
-ModelingToolkit = "9.19.0"
+ModelingToolkit = "9.20.0"
 ModelingToolkitStandardLibrary = "2"
 OMJulia = "0.3.1"
 OrdinaryDiffEq = "6.84.0"
 Polynomials = "4.0.8"
 PreferenceTools = "0.1.2"
 SciMLBenchmarks = "0.1.3"
-Symbolics = "5.30.4"
+Symbolics = "5.32"

docs/pages.jl

@@ -40,7 +40,7 @@ end
 # The result is in alphabetical order, change to the wanted order
 
 permute!(pages,
-    [1, 15, 12, 10, 21, 3, 19, 29, 5, 7, 17, 14, 27, 6, 8, 16, 20, 30, 18, 28, 11, 9, 23, 22, 26, 4, 13, 24, 25, 2, 31]
+    [1, 16, 13, 11, 22, 3, 4, 20, 30, 6, 8, 18, 15, 28, 7, 9, 17, 21, 31, 19, 29, 12, 10, 24, 23, 27, 5, 14, 25, 26, 2, 32]
 )
 
 names = [
@@ -50,6 +50,7 @@ names = [
     "Interval Rootfinding",
     "Nonlinear Solvers",
     "Automatic Differentiation",
+    "Sparse Automatic Differentiation",
     "Non-Stiff Ordinary Differential Equations (ODEs)",
     "Stiff Ordinary Differential Equations (ODEs)",
     "Biological Differential Equations",