@@ -8,138 +8,163 @@ author: Chris Rackauckas
8
8
using Distributed
9
9
addprocs(2)
10
10
11
- p1 = Vector{Any}(undef,3)
12
- p2 = Vector{Any}(undef,3)
13
- p3 = Vector{Any}(undef,3)
11
+ p1 = Vector{Any}(undef, 3)
12
+ p2 = Vector{Any}(undef, 3)
13
+ p3 = Vector{Any}(undef, 3)
14
14
15
15
@everywhere begin
16
- using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
17
- import SDEProblemLibrary: prob_sde_additive,
18
- prob_sde_linear, prob_sde_wave
16
+ using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots,
17
+ ParallelDataTransfer
18
+ import SDEProblemLibrary: prob_sde_additive,
19
+ prob_sde_linear, prob_sde_wave
19
20
end
20
21
21
22
using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
22
23
import SDEProblemLibrary: prob_sde_additive,
23
- prob_sde_linear, prob_sde_wave
24
+ prob_sde_linear, prob_sde_wave
24
25
25
- probs = Matrix{SDEProblem}(undef,3, 3)
26
+ probs = Matrix{SDEProblem}(undef, 3, 3)
26
27
## Problem 1
27
28
prob = prob_sde_linear
28
- probs[1,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
29
- probs[1,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
30
- probs[1,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
29
+ probs[1, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
30
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
31
+ probs[1, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
32
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
33
+ probs[1, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
34
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
31
35
## Problem 2
32
36
prob = prob_sde_wave
33
- probs[2,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
34
- probs[2,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
35
- probs[2,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
37
+ probs[2, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
38
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
39
+ probs[2, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
40
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
41
+ probs[2, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
42
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
36
43
## Problem 3
37
44
prob = prob_sde_additive
38
- probs[3,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
39
- probs[3,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
40
- probs[3,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
41
-
42
- fullMeans = Vector{Array}(undef,3)
43
- fullMedians = Vector{Array}(undef,3)
44
- fullElapsed = Vector{Array}(undef,3)
45
- fullTols = Vector{Array}(undef,3)
45
+ probs[3, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
46
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
47
+ probs[3, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
48
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
49
+ probs[3, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
50
+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
51
+
52
+ fullMeans = Vector{Array}(undef, 3)
53
+ fullMedians = Vector{Array}(undef, 3)
54
+ fullElapsed = Vector{Array}(undef, 3)
55
+ fullTols = Vector{Array}(undef, 3)
46
56
offset = 0
47
57
48
- Ns = [17,23,
49
- 17]
58
+ Ns = [17, 23,
59
+ 17]
50
60
```
51
61
52
62
Timings are only valid if no workers die. Workers die if you run out of memory.
53
63
54
64
```julia
55
- for k in 1:size(probs,1)
56
- global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
57
- println("Problem $k")
58
- ## Setup
59
- N = Ns[k]
60
-
61
- msims = Vector{Any}(undef,N)
62
- elapsed = Array{Float64}(undef,N,3)
63
- medians = Array{Float64}(undef,N,3)
64
- means = Array{Float64}(undef,N,3)
65
- tols = Array{Float64}(undef,N,3)
66
-
67
- #Compile
68
- prob = probs[k,1]
69
- ParallelDataTransfer.sendto(workers(), prob=prob)
70
- monte_prob = EnsembleProblem(prob)
71
- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
72
-
73
- println("RSwM1")
74
- for i=1+offset:N+offset
75
- tols[i-offset,1] = 2.0^(-i-1)
76
- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
77
- trajectories=1000,abstol=2.0^(-i-1),
78
- reltol=0,force_dtmin=true))
79
- elapsed[i-offset,1] = msims[i-offset].elapsedTime
80
- medians[i-offset,1] = msims[i-offset].error_medians[:final]
81
- means[i-offset,1] = msims[i-offset].error_means[:final]
82
- end
83
-
84
- println("RSwM2")
85
- prob = probs[k,2]
86
-
87
- ParallelDataTransfer.sendto(workers(), prob=prob)
88
- monte_prob = EnsembleProblem(prob)
89
- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
90
-
91
- for i=1+offset:N+offset
92
- tols[i-offset,2] = 2.0^(-i-1)
93
- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
94
- trajectories=1000,abstol=2.0^(-i-1),
95
- reltol=0,force_dtmin=true))
96
- elapsed[i-offset,2] = msims[i-offset].elapsedTime
97
- medians[i-offset,2] = msims[i-offset].error_medians[:final]
98
- means[i-offset,2] = msims[i-offset].error_means[:final]
99
- end
100
-
101
- println("RSwM3")
102
- prob = probs[k,3]
103
- ParallelDataTransfer.sendto(workers(), prob=prob)
104
- monte_prob = EnsembleProblem(prob)
105
- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
106
-
107
- for i=1+offset:N+offset
108
- tols[i-offset,3] = 2.0^(-i-1)
109
- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
110
- adaptive=true,trajectories=1000,abstol=2.0^(-i-1),
111
- reltol=0,force_dtmin=true))
112
- elapsed[i-offset,3] = msims[i-offset].elapsedTime
113
- medians[i-offset,3] = msims[i-offset].error_medians[:final]
114
- means[i-offset,3] = msims[i-offset].error_means[:final]
115
- end
116
-
117
- fullMeans[k] = means
118
- fullMedians[k] =medians
119
- fullElapsed[k] = elapsed
120
- fullTols[k] = tols
65
+ for k in 1:size(probs, 1)
66
+ global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
67
+ println("Problem $k")
68
+ ## Setup
69
+ N = Ns[k]
70
+
71
+ msims = Vector{Any}(undef, N)
72
+ elapsed = Array{Float64}(undef, N, 3)
73
+ medians = Array{Float64}(undef, N, 3)
74
+ means = Array{Float64}(undef, N, 3)
75
+ tols = Array{Float64}(undef, N, 3)
76
+
77
+ #Compile
78
+ prob = probs[k, 1]
79
+ ParallelDataTransfer.sendto(workers(), prob = prob)
80
+ monte_prob = EnsembleProblem(prob)
81
+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
82
+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
83
+
84
+ println("RSwM1")
85
+ for i in (1 + offset):(N + offset)
86
+ tols[i - offset, 1] = 2.0^(-i-1)
87
+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
88
+ trajectories = 1000, abstol = 2.0^(-i-1),
89
+ reltol = 0, force_dtmin = true))
90
+ elapsed[i - offset, 1] = msims[i - offset].elapsedTime
91
+ medians[i - offset, 1] = msims[i - offset].error_medians[:final]
92
+ means[i - offset, 1] = msims[i - offset].error_means[:final]
93
+ end
94
+
95
+ println("RSwM2")
96
+ prob = probs[k, 2]
97
+
98
+ ParallelDataTransfer.sendto(workers(), prob = prob)
99
+ monte_prob = EnsembleProblem(prob)
100
+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
101
+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
102
+
103
+ for i in (1 + offset):(N + offset)
104
+ tols[i - offset, 2] = 2.0^(-i-1)
105
+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
106
+ trajectories = 1000, abstol = 2.0^(-i-1),
107
+ reltol = 0, force_dtmin = true))
108
+ elapsed[i - offset, 2] = msims[i - offset].elapsedTime
109
+ medians[i - offset, 2] = msims[i - offset].error_medians[:final]
110
+ means[i - offset, 2] = msims[i - offset].error_means[:final]
111
+ end
112
+
113
+ println("RSwM3")
114
+ prob = probs[k, 3]
115
+ ParallelDataTransfer.sendto(workers(), prob = prob)
116
+ monte_prob = EnsembleProblem(prob)
117
+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
118
+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
119
+
120
+ for i in (1 + offset):(N + offset)
121
+ tols[i - offset, 3] = 2.0^(-i-1)
122
+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
123
+ adaptive = true, trajectories = 1000, abstol = 2.0^(-i-1),
124
+ reltol = 0, force_dtmin = true))
125
+ elapsed[i - offset, 3] = msims[i - offset].elapsedTime
126
+ medians[i - offset, 3] = msims[i - offset].error_medians[:final]
127
+ means[i - offset, 3] = msims[i - offset].error_means[:final]
128
+ end
129
+
130
+ fullMeans[k] = means
131
+ fullMedians[k] = medians
132
+ fullElapsed[k] = elapsed
133
+ fullTols[k] = tols
121
134
end
122
135
```
123
136
124
137
```julia
125
- gr(fmt= :svg)
138
+ gr(fmt = :svg)
126
139
lw=3
127
- leg=String["RSwM1","RSwM2","RSwM3"]
140
+ leg=String["RSwM1", "RSwM2", "RSwM3"]
128
141
129
142
titleFontSize = 16
130
143
guideFontSize = 14
131
- legendFontSize= 14
132
- tickFontSize = 12
133
-
134
- for k in 1:size(probs,1)
135
- global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
136
- p1[k] = Plots.plot(fullTols[k],fullMeans[k],xscale=:log10,yscale=:log10, xguide="Absolute Tolerance",yguide="Mean Final Error",title="Example $k" ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
137
- p2[k] = Plots.plot(fullTols[k],fullMedians[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Median Final Error",title="Example $k",linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
138
- p3[k] = Plots.plot(fullTols[k],fullElapsed[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Elapsed Time",title="Example $k" ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
144
+ legendFontSize = 14
145
+ tickFontSize = 12
146
+
147
+ for k in 1:size(probs, 1)
148
+ global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
149
+ p1[k] = Plots.plot(fullTols[k], fullMeans[k], xscale = :log10, yscale = :log10,
150
+ xguide = "Absolute Tolerance", yguide = "Mean Final Error",
151
+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
152
+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
153
+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
154
+ p2[k] = Plots.plot(fullTols[k], fullMedians[k], xscale = :log10, yscale = :log10,
155
+ xguide = "Absolute Tolerance", yguide = "Median Final Error",
156
+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
157
+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
158
+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
159
+ p3[k] = Plots.plot(fullTols[k], fullElapsed[k], xscale = :log10, yscale = :log10,
160
+ xguide = "Absolute Tolerance", yguide = "Elapsed Time",
161
+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
162
+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
163
+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
139
164
end
140
165
141
166
Plots.plot!(p1[1])
142
- Plots.plot(p1[1],p1[2],p1[3],layout= (3,1),size= (1000,800))
167
+ Plots.plot(p1[1], p1[2], p1[3], layout = (3, 1), size = (1000, 800))
143
168
```
144
169
145
170
```julia
@@ -148,17 +173,17 @@ Plots.plot(p1[1],p1[2],p1[3],layout=(3,1),size=(1000,800))
148
173
```
149
174
150
175
```julia
151
- plot(p3[1],p3[2],p3[3],layout= (3,1),size= (1000,800))
176
+ plot(p3[1], p3[2], p3[3], layout = (3, 1), size = (1000, 800))
152
177
#savefig("timevstol.png")
153
178
#savefig("timevstol.pdf")
154
179
```
155
180
156
181
```julia
157
- plot(p1[1],p3[1],p1[2],p3[2],p1[3],p3[3],layout= (3,2),size= (1000,800))
182
+ plot(p1[1], p3[1], p1[2], p3[2], p1[3], p3[3], layout = (3, 2), size = (1000, 800))
158
183
```
159
184
160
185
```julia
161
186
162
187
using SciMLBenchmarks
163
- SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
188
+ SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
164
189
```
0 commit comments