%BPATH1 Brownian path simulation
function [t,W]=BPATH1(T,N)
randn('state',100)
dt = T/N;
dW = zeros(10,N);
W = zeros(10,N);

dW(:,1) = sqrt(dt)*randn(1,10);
W(:,1) = dW(:,1);

for j = 2:N
    dW(:,j) = sqrt(dt)*randn(1,10);
    W(:,j) = W(:,j-1) + dW(:,j);
end
t = [0:dt:T-dt];
end

% In Another File...

f = @() BPATH1(1,10000);
timeit(f,2)
[t,W] = BPATH1(1,10000);
plot(t,W,'r-')
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

using Random
function bpath(T,N)
    Random.seed!(100)
    dt = T/N
    dW = zeros(10,N)
    W = zeros(10,N)

    dW[:,1] = sqrt(dt)*randn(10)
    W[:,1] = dW[:,1]

    for j = 2:N
        dW[:,j] .= sqrt(dt)*randn(10)
        W[:,j] .= W[:,j-1] .+ dW[:,j]
    end
    [0:dt:T-dt],W
end # Translation took < 1 minute
@time t,W = bpath(1,100000)
@time t,W = bpath(1,100000)
# 10x speedup over MATLAB
using Plots; gr(fmt = :png)
plot(t,W',color=:red,xlabel="t",ylabel="W(t)")

f(x,y) = x + y
@code_llvm f(1,1)
@code_llvm f(1.0,1.0)

using DifferentialEquations
function f(du,u,p,t)
  x, y = u
  α,β,γ,δ = p
  du[1] = α*x - β*x*y
  du[2] = -γ*y + δ*x*y
end
p = (1.5,1.0,3.0,1.0); u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob1 = ODEProblem(f,u0,tspan,p)
sol = solve(prob1)
using Plots; plot(sol)

f = @ode_def LotkaVolterra begin
  dx  = α*x  - β*x*y
  dy = -γ*y + δ*x*y
end α β γ δ

f = @ode_def LotkaVolterra begin
  d🐁  = α*🐁 - β*🐁*🐈
  d🐈 = -γ*🐈 + δ*🐁*🐈
end α β γ δ

sir_model = @reaction_network SIR begin
    c1, s + i --> 2i
    c2, i --> r
end c1 c2

using StaticArrays
static_u0 = @SVector [1.0,1.0]
function f2(u,p,t)
  x, y = u
  α,β,γ,δ = p
  @SVector [α*x - β*x*y, -γ*y + δ*x*y]
end
prob2 = ODEProblem(f2,static_u0,tspan,p)
sol = solve(prob2);

using BenchmarkTools
@benchmark solve(prob1)

@benchmark solve(prob2)

@benchmark solve(prob2,Vern9())

using LSODA
@benchmark solve(prob1,lsoda()) # Default method of SciPy, deSolve, etc.

using Measurements
u0 = [1.0 ± 0.0 ,1.0 ± 0.0]
p = (1.5 ± 0.0,1.0 ± 0.1,3.0 ± 0.2,1.0 ± 0.1)
prob_error = ODEProblem(f,u0,tspan,p)
sol = solve(prob_error,Tsit5(), saveat=0.2)
plot(sol)

using ForwardDiff: Dual
p1dual = Dual{Float64}(1.5, (1.0, 0.0))
p2dual = Dual{Float64}(1.0, (0.0, 1.0))
pdual = (p1dual, p2dual, 3.0, 1.0)
u0 = [Dual{Float64}(1.0, (0.0, 0.0)),Dual{Float64}(1.0, (0.0, 0.0))]
prob_dual = ODEProblem(f,u0,tspan,pdual)
sol_dual = solve(prob_dual,Tsit5(), saveat=0.2)

timepoints = [i for i in sol_dual.t]
sensitivity_forward_diff = [i[1].partials.values[1] for i in sol_dual.u]
Plots.plot(timepoints,sensitivity_forward_diff,title="dx/da",lw=3,xaxis="Time")

p = big.((1.5,1.0,3.0,1.0)); u0 = big.([1.0,1.0])
tspan = big.((0.0,10.0))
prob1 = ODEProblem(f,u0,tspan,p)
sol = solve(prob1,Vern9(),abstol=1e-25,reltol=1e-25)
sol[10]

using Unitful
p = (1.5,1.0,3.0,1.0)./u"s"; u0 = [1.0u"kg",1.0u"kg"]
tspan = (0.0u"s",10.0u"s")
prob_units = ODEProblem(f,u0,tspan,p)
sol = solve(prob_units,Tsit5())

#= 
f = @ode_def LotkaVolterra begin
  d🐁  = α*🐁 - β*🐁*🐈
  d🐈 = -γ*🐈 + δ*🐁*🐈
end α β γ δ
=#
p = (1.5,1.0/u"kg",3.0,1.0/u"kg")./u"s"; u0 = [1.0u"kg",1.0u"kg"]
tspan = (0.0u"s",10.0u"s")
prob_units = ODEProblem(f,u0,tspan,p)
sol = solve(prob_units,Tsit5())
sol.t[10],sol[10]

sol(5.5u"s")

using DifferentialEquations
function g(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.2u[2]
end
p = (1.5,1.0,3.0,1.0); u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob1 = SDEProblem(f,g,u0,tspan,p)
sol = solve(prob1)
using Plots; plot(sol)

function f(du,u,h,p,t)
  x, y = u
  α,β,γ,δ = p
  du[1] = α*h(p,t-1)[1] - β*x*y
  du[2] = -γ*y + δ*x*y
end
p = (1.5,1.0,3.0,1.0); u0 = [1.0;1.0]
tspan = (0.0,10.0)
_h(p,t) = ones(2)
prob1 = DDEProblem(f,u0,_h,tspan,p)
sol = solve(prob1)
using Plots; plot(sol)