recursive_mean

Author: ChrisRackauckas

src/convergence.jl

@@ -40,16 +40,16 @@ function test_convergence(dts::AbstractArray,prob::Union{AbstractRODEProblem,Abs
   # Now Calculate Weak Errors
   additional_errors = Dict()
   # Final
-  m_final = mean([s[end] for s in solutions],1)
-  m_final_analytic = mean([s.u_analytic[end] for s in solutions],1)
-  additional_errors[:weak_final] = mean.(vecvecapply((x)->abs.(x),m_final - m_final_analytic))
+  m_final = recursive_mean([s[end] for s in solutions],1)
+  m_final_analytic = recursive_mean([s.u_analytic[end] for s in solutions],1)
+  additional_errors[:weak_final] = recursive_mean.(vecvecapply((x)->abs.(x),m_final - m_final_analytic))
   if timeseries_errors
     l2_tmp = Vector{eltype(solutions[1][1])}(size(solutions,2))
     max_tmp = Vector{eltype(solutions[1][1])}(size(solutions,2))
     for i in 1:size(solutions,2)
       solcol = @view solutions[:,i]
-      m_errors = [mean([solcol[j][i] for j in 1:length(solcol)]) for i in 1:length(solcol[1])]
-      m_errors_analytic = [mean([solcol[j].u_analytic[i] for j in 1:length(solcol)]) for i in 1:length(solcol[1])]
+      m_errors = [recursive_mean([solcol[j][i] for j in 1:length(solcol)]) for i in 1:length(solcol[1])]
+      m_errors_analytic = [recursive_mean([solcol[j].u_analytic[i] for j in 1:length(solcol)]) for i in 1:length(solcol[1])]
       ts_weak_errors = [abs.(m_errors[i] - m_errors_analytic[i]) for i in 1:length(m_errors)]
       ts_l2_errors = [sqrt.(sum(abs2,err)/length(err)) for err in ts_weak_errors]
       l2_tmp[i] = sqrt(sum(abs2,ts_l2_errors)/length(ts_l2_errors))

src/test_solution.jl

@@ -2,16 +2,24 @@
 TestSolution
 
 """
-type TestSolution{hasinterp} <: DESolution
-  t
-  u
-  interp
-  dense
+type TestSolution{T,N,hasinterp,tType,uType,iType} <: AbstractTimeseriesSolution{T,N}
+  t::tType
+  u::uType
+  interp::iType
+  dense::Bool
 end
 (T::TestSolution)(t) = T.interp(t)
-TestSolution(t,u) = TestSolution{false}(t,u,nothing,false)
-TestSolution(t,u,interp) = TestSolution{true}(t,u,interp,true)
-TestSolution(interp::DESolution) = TestSolution{true}(nothing,nothing,interp,true)
+function TestSolution(t,u)
+  T = eltype(eltype(u))
+  N = length((size(u[1])..., length(u)))
+  TestSolution{T,N,false,typeof(t),typeof(u),Void}(t,u,nothing,false)
+end
+function TestSolution(t,u,interp)
+  T = eltype(eltype(u))
+  N = length((size(u[1])..., length(u)))
+  TestSolution{T,N,length(true),typeof(t),typeof(u),typeof(interp)}(t,u,interp,true)
+end
+TestSolution(interp::DESolution) = TestSolution{Void,0,true,Void,Void,typeof(interp)}(nothing,nothing,interp,true)
 
 """
 `appxtrue(sol::AbstractODESolution,sol2::TestSolution)`
@@ -25,17 +33,17 @@ function appxtrue{hasinterp}(sol::AbstractODESolution,sol2::TestSolution{hasinte
     sol2.u = sol2(sol.t)
     sol2.t = sol.t
   end
-  errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])))
+  errors = Dict(:final=>recursive_mean(abs.(sol[end]-sol2[end])))
   if sol2.dense
     timeseries_analytic = sol2(sol.t)
     errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol[:]-timeseries_analytic))
-    errors[:l2] = sqrt(mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic)))
+    errors[:l2] = sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic)))
     if !(typeof(sol) <: AbstractRODESolution) && sol.dense
       densetimes = collect(linspace(sol.t[1],sol.t[end],100))
       interp_u = sol(densetimes)
       interp_analytic = sol2(densetimes)
       interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),
-                           :L2=>sqrt(mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
+                           :L2=>sqrt(recursive_mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
       errors = merge(errors,interp_errors)
     end
   end
@@ -56,15 +64,15 @@ errors for sol. If sol2 has no interpolant, only the final error is
 calculated.
 """
 function appxtrue(sol::AbstractODESolution,sol2::AbstractODESolution)
-  errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])))
+  errors = Dict(:final=>recursive_mean(abs.(sol[end]-sol2[end])))
   if !(typeof(sol2) <: AbstractRODESolution) && sol2.dense
     timeseries_analytic = sol2(sol.t)
-    errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])),:l∞=>maximum(vecvecapply((x)->abs.(x),sol[:]-timeseries_analytic)),:l2=>sqrt(mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic))))
+    errors = Dict(:final=>recursive_mean(abs.(sol[end]-sol2[end])),:l∞=>maximum(vecvecapply((x)->abs.(x),sol[:]-timeseries_analytic)),:l2=>sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic))))
     if !(typeof(sol) <: AbstractRODESolution) && sol.dense
       densetimes = collect(linspace(sol.t[1],sol.t[end],100))
       interp_u = sol(densetimes)
       interp_analytic = sol2(densetimes)
-      interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),:L2=>sqrt(mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
+      interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),:L2=>sqrt(recursive_mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
       errors = merge(errors,interp_errors)
     end
   end