text:
#=
            # P is the precision used for this value
type ArfFloat{P}  <: Real
  exponentOf2::Int # fmpz
  nwords_sign::UInt # mp_size_t
  significand1::UInt # significand_struct
  significand2::UInt
end
=#

precision(x::ArfFloat{P}) where {P} = P
precision(::Type{ArfFloat{P}}) where {P} = P
precision(::Type{ArfFloat}) = ArbFloatPrecision[1]
setprecision(::Type{ArfFloat}, x::Int) = setprecision(ArbFloat, x)

# a type specific hash function helps the type to 'just work'
const hash_arffloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4
const hash_0_arffloat_lo = hash(zero(UInt), hash_arffloat_lo)
hash(z::ArfFloat{P}, h::UInt) where {P} =
    hash(reinterpret(UInt,z.significand1)$z.exponentOf2,
         (h $ hash(z.significand2$(~reinterpret(UInt,P)), hash_arffloat_lo) $ hash_0_arffloat_lo))

weakcopy(x::ArfFloat{P}) where {P} = WeakRef(x)

function copy(x::ArfFloat{P}) where {P}
    z = initializer(ArfFloat{P})
    ccall(@libarb(arf_set), Cvoid, (Ref{ArfFloat{P}}, Ref{ArfFloat{P}}), z, x)
    return z
end

# initialize and zero a variable of type ArfFloat
function release(x::ArfFloat{P}) where {P}
    ccall(@libarb(arf_clear), Cvoid, (Ref{ArfFloat{P}}, ), x)
    return nothing
end

function initializer(::Type{ArfFloat{P}}) where {P}
    z = ArfFloat{P}(zero(Int), zero(UInt64), zero(Int64), zero(Int64))
    ccall(@libarb(arf_init), Cvoid, (Ref{ArfFloat{P}}, ), z)
    finalizer(release, z)
    return z
end

# empty constructor
@inline function ArfFloat()
     P = precision(ArfFloat)
     return initializer(ArfFloat{P})
end

function deepcopy(x::ArfFloat{P}) where {P}
    z = initializer(ArfFloat{P})
    ccall(@libarb(arf_set), Cvoid, (Ref{ArfFloat{P}}, Ref{ArfFloat{P}}), z, x)
    return z
end

function zero(::Type{ArfFloat{P}}) where {P}
    z = initializer( ArbFloat{P} )
    ccall(@libarb(arf_zero), Cvoid, (Ref{ArfFloat{P}}, ), z)
    return z
end

function one(::Type{ArfFloat{P}}) where {P}
    z = initializer( ArbFloat{P} )
    ccall(@libarb(arf_one), Cvoid, (Ref{ArfFloat{P}}, ), z)
    return z
end

function isnan(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_nan), Cint, (Ref{ArfFloat},), x)
end
function isinf(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_inf), Cint, (Ref{ArfFloat},), x)
end
function isposinf(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_pos_inf), Cint, (Ref{ArfFloat},), x)
end
function isneginf(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_neg_inf), Cint, (Ref{ArfFloat},), x)
end
function iszero(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_zero), Cint, (Ref{ArfFloat},), x)
end
function isone(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_one), Cint, (Ref{ArfFloat},), x)
end
function isfinite(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_finite), Cint, (Ref{ArfFloat},), x)
end
function isnormal(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_normal), Cint, (Ref{ArfFloat},), x)
end
function isspecial(x::T) where {T <: ArfFloat}
    zero(Cint) != ccall(@libarb(arf_is_special), Cint, (Ref{ArfFloat},), x)
end

function notnan(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_nan), Cint, (Ref{ArfFloat},), x)
end
function notinf(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_inf), Cint, (Ref{ArfFloat},), x)
end
function notposinf(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_pos_inf), Cint, (Ref{ArfFloat},), x)
end
function notneginf(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_neg_inf), Cint, (Ref{ArfFloat},), x)
end
function notzero(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_zero), Cint, (Ref{ArfFloat},), x)
end
function notone(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_one), Cint, (Ref{ArfFloat},), x)
end
function notfinite(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_finite), Cint, (Ref{ArfFloat},), x)
end
function notnormal(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_normal), Cint, (Ref{ArfFloat},), x)
end
function notspecial(x::T) where {T <: ArfFloat}
    zero(Cint) == ccall(@libarb(arf_is_special), Cint, (Ref{ArfFloat},), x)
end

function convert(::Type{BigFloat}, x::ArfFloat{P}) where {P}
    z = BigFloat(0)
    r = ccall(@libarb(arf_get_mpfr), Int, (Ref{BigFloat}, Ref{ArfFloat{P}}, Cint), z, x, 0)
    return z
end

function convert(::Type{ArfFloat{P}}, x::BigFloat) where {P}
    z = initializer(ArfFloat{P})
    ccall(@libarb(arf_set_mpfr), Cvoid, (Ref{ArfFloat{P}}, Ref{BigFloat}), z, x)
    z
end
convert(::Type{ArfFloat}, x::BigFloat) = convert(ArfFloat{precision(ArfFloat)}, x)

convert(::Type{T}, x::BigInt) where {T <: ArfFloat} = convert(T, convert(BigFloat, x))
convert(::Type{ArfFloat{P}}, x::BigInt) where {P} = convert(ArfFloat{P}, convert(BigFloat, x))

function convert(::Type{ArfFloat{P}}, x::Int64) where {P}
    z = initializer(ArfFloat{P})
    ccall(@libarb(arf_set_si), Cvoid, (Ref{ArfFloat{P}}, Ref{Int64}), z, x)
    z
end
convert(::Type{ArfFloat}, x::Int64) = convert(ArfFloat{precision(ArfFloat)}, x)

function convert(::Type{ArfFloat{P}}, x::Float64) where {P}
    z = initializer(ArfFloat{P})
    ccall(@libarb(arf_set_d), Cvoid, (Ref{ArfFloat{P}}, Ref{Float64}), z, x)
    z
end
function convert(::Type{ArfFloat}, x::Float64)
     P = precision(ArfFloat)
     return convert(ArfFloat{P}, x)
end

midpoint(x::ArfFloat{P}) where {P} = x

radius(x::ArfFloat{P}) where {P} = ArfFloat{P}(0)

#=
#define ARF_RND_DOWN FMPR_RND_DOWN
#define ARF_RND_UP FMPR_RND_UP
#define ARF_RND_FLOOR FMPR_RND_FLOOR
#define ARF_RND_CEIL FMPR_RND_CEIL
#define ARF_RND_NEAR FMPR_RND_NEAR
=#

function round(x::T, prec::Int64) where {T <: ArfFloat}
    P = precision(T)
    z = ArfFloat{P}()
    ccall(@libarb(arf_set_round), Int, (Ref{ArfFloat}, Ref{ArfFloat}, Int64, Int), z, x, prec, 2)
    return z
end

#=
function frexp{P}(x::ArfFloat{P})
   significand = ArfFloat{P}()
   exponentOf2 = zero(Int64)
   ccall(@libarb(arf_frexp), Void, (Ref{ArfFloat{P}}, Int64, Ref{ArfFloat{P}}), &significand, exponentOf2, x)
   significand, exponentOf2
end
=#

function min(x::T, y::T) where {T <: ArfFloat}
    c = ccall(@libarb(arf_cmp), Cint, (Ref{ArfFloat}, Ref{ArfFloat}), x, y)
    return (c < 0) ? x : y
end

function max(x::T, y::T) where {T <: ArfFloat}
    c = ccall(@libarb(arf_cmp), Cint, (Ref{ArfFloat}, Ref{ArfFloat}), x, y)
    return (c > 0) ? x : y
end


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 786a23ad-ee13-4299-b58f-b4bc2cff50f7
using Plots

# ╔═╡ f84ba617-f84e-4072-94a1-9d85fd6e2ef0
html"
Laplace.jl

Solving the Laplace Equation using Jacobi, Gauss-Seidel, and Simultaneous Over Relaxation methods on a square grid

First, set N, the number of grid points on a side:
"

# ╔═╡ 7cea9966-828e-48b1-aa2d-7f3fcb21cadc
@bind N html"<input type='range' min='3' max='30' step='1' value='20'>"

# ╔═╡ 8b743ed0-ac8e-4135-8461-ad0ce14da8db
N

# ╔═╡ 34e94d32-6b28-4cf7-bc22-20f9cfd0d346
html"
Setting up the grid
"

# ╔═╡ cd6f8f6c-1e4c-4eb4-b7ae-be1c62e5f3ad
begin
L = 1 				# Size of system (length of side normalized to be 1)
h = L/(N-1) 		# grid spacing
x = y = -L/2:h:L/2 	# x and y locations
end

# ╔═╡ 084cd4d2-2e7f-4700-862e-c88f311dc633
begin
ϕ = zeros(N, N) 	# Matrix potential
ϕ₀ = 1 				# Potential at y = L
end

# ╔═╡ 09a56e33-ef32-45de-822e-299a8b209180
html"
Resetting functions

The Initial Condition is defined as:"

# ╔═╡ db7f9c67-0c89-4b63-86de-aa7a7da608c0
function reset!(ψ::Array, ψ₀::Real = ϕ₀, n::Integer = N)
	for i ∈ 1:n, j ∈ 1:n
		ψ[i, j] = ψ₀ * ( 4/(π*sinh(π)) ) * sin(π*x[i]/L) * sinh(π*y[j]/L)
	end
end

# ╔═╡ 98412958-f1d0-4ee3-8439-815c982d9e5d
html"Boundary Conditions: ϕ = ϕ₀ along the rightmost column, 0 along the other sides:"

# ╔═╡ 10913cfe-fbb1-432e-bc34-a22dfe3602fe
function setboundaryconditions!(ψ::Array, ψ₀::Real = ϕ₀)
	ψ[1, :] .= ψ[end, :] .= ψ[:, 1] .= 0
	ψ[:, end] .= fill(ψ₀, size(ψ, 1))
end

# ╔═╡ 0955e4b8-19e7-49e5-bdcf-187996e337c5
html"
Step functions

Jacobi Method"

# ╔═╡ 31a399e8-fb24-4210-b0a7-58a9a6192abe
function jacobistep!(ψ::Array, changesum::Real, n::Integer = N)
	ψₙ = copy(ψ)
	for i ∈ 2:(n-1), j ∈ 2:(n-1)
		ψₙ[i, j] = .25*(ψ[i+1, j] + ψ[i-1, j] + ψ[i, j+1] + ψ[i, j-1])
		changesum += abs(1 - ψ[i, j]/ψₙ[i, j])
	end
	ψ[:, :] = ψₙ[:, :]
end

# ╔═╡ 35eefeef-7d13-4e4e-9656-0e72214ef7ef
html"Gauss-Seidel Method"

# ╔═╡ f214b84c-4fdb-4567-904e-4df485dca12b
function gs_step!(ψ::Array, changesum::Real, n::Integer = N)
	for i ∈ 2:(n-1), j ∈ 2:(n-1)
		newval = .25*(ψ[i+1, j] + ψ[i-1, j] + ψ[i, j+1] + ψ[i, j-1])
		changesum += abs(1 - ψ[i, j]/newval)
		ψ[i, j] = newval
	end
end

# ╔═╡ b01bef7f-c9c3-4e38-9aef-c7e727739bb8
html"Simultaneous Over-Relaxation Method"

# ╔═╡ 927e449c-699e-40d6-ad93-d4da9a77a392
Ω₀ = 2 / (1+sin(π/N)) 	# optimizes convergence time for SOR method

# ╔═╡ 1eb90fd0-6be6-4cbe-b513-ed8eb9bc81a1
function sorstep!(ψ::Array, changesum::Real, n::Integer = N, Ω::Real = Ω₀)
	for i ∈ 2:(n-1), j ∈ 2:(n-1)
		newval = .25Ω * (ψ[i+1, j] + ψ[i-1, j] + ψ[i, j+1] + ψ[i, j-1]) +
			(1 - Ω) * ψ[i, j]
		changesum += abs(1 - ψ[i, j]/newval)
		ψ[i, j] = newval
	end
end

# ╔═╡ 6b360f13-ec44-4a9e-a064-38eaadb13c0e
html"
Other preparations

Maximum number of steps (so we don't just loop forever):"

# ╔═╡ 667382b6-c6aa-4655-bec5-9c00f3b7b463
imax = N^2

# ╔═╡ 10bdbf0f-32c7-4df3-869a-56d266f56fac
html"Iteration ceases once the average correction is < 10-4:"

# ╔═╡ ce3419a8-3984-488e-85e9-0d0786bbae2d
changedesired = 10^-4

# ╔═╡ 1d13f124-1d16-4374-81c5-4ef8a0c8cfa8
html"
Plotting
"

# ╔═╡ 08463463-a613-45dc-8bf6-8ba2569d9764
let
reset!(ϕ)
setboundaryconditions!(ϕ)
for _ ∈ 1:imax
	changesum = 0 		# Total change in system this iteration
	jacobistep!(ϕ, changesum) 	# Then take the average change per (N-2)^2 points:
	changesum/(N-2)^2 < changedesired && continue
end
surface(x, y, rotr90(ϕ), title="Jacobi Method", xticks=[], yticks=[],
		xlabel="x", ylabel="y")
end

# ╔═╡ e425201e-acb5-4975-bc70-47d80fc5686e
let
reset!(ϕ)
setboundaryconditions!(ϕ)
for _ ∈ 1:imax
	changesum = 0
	gs_step!(ϕ, changesum)
	changesum/(N-2)^2 < changedesired && continue
end
surface(x, y, rotr90(ϕ), title="Gauss-Seidel Method", xticks=[], yticks=[],
		xlabel="x", ylabel="y")
end

# ╔═╡ ee21c706-b9f3-43ba-8c28-ea152ef4a9bb
let
reset!(ϕ)
setboundaryconditions!(ϕ)
for _ ∈ 1:imax
	changesum = 0
	sorstep!(ϕ, changesum)
	changesum/(N-2)^2 < changedesired && continue
end
surface(x, y, rotr90(ϕ), title="Simultaneous Over-Relaxation Method",
		xticks=[], yticks=[], xlabel="x", ylabel="y")
end

# ╔═╡ Cell order:
# ╟─f84ba617-f84e-4072-94a1-9d85fd6e2ef0
# ╟─7cea9966-828e-48b1-aa2d-7f3fcb21cadc
# ╠═8b743ed0-ac8e-4135-8461-ad0ce14da8db
# ╟─34e94d32-6b28-4cf7-bc22-20f9cfd0d346
# ╠═cd6f8f6c-1e4c-4eb4-b7ae-be1c62e5f3ad
# ╠═084cd4d2-2e7f-4700-862e-c88f311dc633
# ╟─09a56e33-ef32-45de-822e-299a8b209180
# ╠═db7f9c67-0c89-4b63-86de-aa7a7da608c0
# ╟─98412958-f1d0-4ee3-8439-815c982d9e5d
# ╠═10913cfe-fbb1-432e-bc34-a22dfe3602fe
# ╟─0955e4b8-19e7-49e5-bdcf-187996e337c5
# ╠═31a399e8-fb24-4210-b0a7-58a9a6192abe
# ╟─35eefeef-7d13-4e4e-9656-0e72214ef7ef
# ╠═f214b84c-4fdb-4567-904e-4df485dca12b
# ╟─b01bef7f-c9c3-4e38-9aef-c7e727739bb8
# ╠═927e449c-699e-40d6-ad93-d4da9a77a392
# ╠═1eb90fd0-6be6-4cbe-b513-ed8eb9bc81a1
# ╟─6b360f13-ec44-4a9e-a064-38eaadb13c0e
# ╠═667382b6-c6aa-4655-bec5-9c00f3b7b463
# ╟─10bdbf0f-32c7-4df3-869a-56d266f56fac
# ╠═ce3419a8-3984-488e-85e9-0d0786bbae2d
# ╟─1d13f124-1d16-4374-81c5-4ef8a0c8cfa8
# ╠═786a23ad-ee13-4299-b58f-b4bc2cff50f7
# ╠═08463463-a613-45dc-8bf6-8ba2569d9764
# ╠═e425201e-acb5-4975-bc70-47d80fc5686e
# ╠═ee21c706-b9f3-43ba-8c28-ea152ef4a9bb

ext: jl

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
using GLAbstraction, Colors, GeometryTypes, GLVisualize, Reactive

if !isdefined(:runtests) # We use the examples to generate docs and tests
    window = glscreen()
    timesignal = bounce(linspace(0.0, 1.0, 360))
end

description = """
Coloured 3D line example.
"""

n = 400 # The number of points per line
nloops = 20 # The number of loops
# The scalar parameter for each line
TL = linspace(-2f0 * pi, 2f0 * pi, n)
# We create a list of positions and connections, each describing a line.
# We will collapse them in one array before plotting.
xyz = Point3f0[]
intensities = Float32[]

# Create each line one after the other in a loop
for i = 1:nloops
    append!(xyz, [Point3f0(sin(t), cos((2 + .02 * i) * t), cos((3 + .03 * i) * t)) for t in TL])
    append!(intensities, fill(i, length(TL)))
end

# map comes from Reactive.jl and allows you to map any Signal to another.
# In this case we create a rotation matrix from the timesignal signal.

rotation = map(timesignal) do t
    rotationmatrix_z(Float32(t*2pi)) # -> 4x4 Float32 rotation matrix
end

# creates a color map from which we can sample for each line
# and add some transparency
cmap = map(x-> RGBA{Float32}(x, 0.4), colormap("Blues", nloops))

lines3d = visualize(
    xyz, :lines,
    intensity = intensities,
    color_map = cmap,
    color_norm = Vec2f0(0, nloops), # normalize intensities. Lookup in cmap will be between 0-1
    model = rotation
)

_view(lines3d, window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
export fixed

include("systems/mixed_system.jl")
include("systems/affine_chart_system.jl")
include("systems/composition_system.jl")
include("systems/fixed_parameter_system.jl")
include("systems/randomized_system.jl")
include("systems/sliced_system.jl")
include("systems/start_pair_system.jl")

"""
    fixed(F::System; compile::Union{Bool,Symbol} = $(COMPILE_DEFAULT[]))

Constructs either a [`CompiledSystem`](@ref) (if `compile = :all`), an
[`InterpretedSystem`](@ref) (if `compile = :none`) or a [`MixedSystem`](@ref) (`compile = :mixed`).
"""
function fixed(F::System; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[], kwargs...)
    if compile == true || compile == :all
        CompiledSystem(F; kwargs...)
    elseif compile == false || compile == :none
        InterpretedSystem(F; kwargs...)
    elseif compile == :mixed
        MixedSystem(F; kwargs...)
    else
        error("Unknown argument $compile for keyword `compile`.")
    end
end
fixed(F::AbstractSystem; kwargs...) = F

set_solution!(x, ::AbstractSystem, y) = (x .= y; x)


text:
using SugarKelp, Test, Interpolations, DataFrames

@info "These tests check if the code has been broken but does not thoroughly check if it produces correct results (for now at least). They are also so monalithic you can't actually tell which bit you've broken."

const arr_lat = [1:3;]
const arr_lon = [1:3;]
const arr_dep = [1:5;]
const arr_t = [1:200;]

const temp = permutedims(repeat([9:(20 - 9) / 199:20;], 1, 5, 4, 3), (4, 3, 2, 1))
const no3 = repeat([5.0], 3, 4, 5, 200)
par = repeat([2.0], 3, 4, 92)
const par_t = [1:2:184;]
const par_fill = -1.0
par[2,3,4] = -1.0
kd = repeat([2.0], 3, 4, 110)
const kd_t = [1:2:220;]
const kd_fill = -1.0
kd[2,3,4] = -1
const beta = repeat([0.1], 3, 4, 5, 200)
const u = repeat([0.15], 3, 4, 5, 200)

t_i = 1.0
nd = 184
a_0 = 0.1;n_0 = 0.022;c_0 = 0.3;

# haveing to use 2013 parameters as cant re define constants but need to be able to run with resp_model 2
@testset "Grid with kd" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (kd, kd_t, kd_fill), nothing, "../src/parameters/2013.jl")
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.323811260597751
    @test all_results[2,2,3,4,123] ≈ 0.02159999999999999
    @test all_results[3,2,3,4,123] ≈ 0.010000000000000002
    @test all_results[4,2,3,4,123] ≈ 0.0036635181806166307
end

@testset "Grid with beta" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl")
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.33853815235193496
    @test all_results[2,2,3,4,123] ≈ 0.02159999999999999
    @test all_results[3,2,3,4,123] ≈ 0.010000000000000002
    @test all_results[4,2,3,4,123] ≈ 0.003952014050125887
end

@testset "Grid with beta, specifying resp model to be 1" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 1)
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.33853815235193496
    @test all_results[2,2,3,4,123] ≈ 0.02159999999999999
    @test all_results[3,2,3,4,123] ≈ 0.010000000000000002
    @test all_results[4,2,3,4,123] ≈ 0.003952014050125887
end

@testset "Grid with beta, resp model 2" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 2)
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.6949592729765186
    @test all_results[2,2,3,4,123] ≈ 0.0215977671314431
    @test all_results[3,2,3,4,123] ≈ 0.01000005496349408
    @test all_results[4,2,3,4,123] ≈ 0.01014604344770969
end

@testset "Grid with beta, resp model 2, specifying dt=1" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 1, 1)
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.33853815235193496
    @test all_results[2,2,3,4,123] ≈ 0.02159999999999999
    @test all_results[3,2,3,4,123] ≈ 0.010000000000000002
    @test all_results[4,2,3,4,123] ≈ 0.003952014050125887
end
@testset "Grid with beta, resp model 2, specifying dt=2" begin
    # This configuration is unstable so area goes below zero and it stops, that would be a useful behaviour to test for but its not done here
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 1, 2)
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), round(nd / 2 + 1))
    @test all_results[1,2,3,4,5] ≈ 0.19669013321363188
    @test all_results[2,2,3,4,5] ≈ 0.01614271696738769
    @test all_results[3,2,3,4,5] ≈ 0.02894612990034886
    @test all_results[4,2,3,4,5] ≈ 0.0011375879549344894
end
@testset "Grid with beta, resp model 2, specifying dt=1, specifying to show progress" begin
    all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 1, 1, true)
    @test size(all_results) == (4, length(arr_lon), length(arr_lat), length(arr_dep), nd + 1)
    @test all_results[1,2,3,4,123] ≈ 0.33853815235193496
    @test all_results[2,2,3,4,123] ≈ 0.02159999999999999
    @test all_results[3,2,3,4,123] ≈ 0.010000000000000002
    @test all_results[4,2,3,4,123] ≈ 0.003952014050125887
end
# This is a test for progress supression but there is no way to actually test it
# all_results = SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_t, no3, temp, u, (par, par_t, par_fill), (nothing, nothing, nothing), beta, "../src/parameters/2013.jl", 1, 1, false)

u_itp = Interpolations.LinearInterpolation(arr_t, u[1,1,1,:])
temp_itp = Interpolations.LinearInterpolation(arr_t, temp[1,1,1,:])
no3_itp = Interpolations.LinearInterpolation(arr_t, no3[1,1,1,:])
par_itp = Interpolations.LinearInterpolation(par_t, par[1,1,:], extrapolation_bc=Flat())

@testset "Single point, dataframe output" begin
    solution, a_result = SugarKelp.solve(t_i, nd, u_itp, temp_itp, par_itp, no3_itp, 55, a_0, n_0, c_0, "../src/parameters/2013.jl")
    @test typeof(a_result) == DataFrames.DataFrame
    @test size(a_result) == (185, 5)
    @test a_result.area[17] ≈ 0.36857308566331204
    @test a_result.nitrogen[17] ≈ 0.01579264171780362
    @test a_result.carbon[17] ≈ 0.027214793958727508
end

@testset "Depth interpolator utility" begin
    tnew = SugarKelp.interp_deps(temp, arr_dep, [1:0.5:5;], -1)
    @test size(tnew) == (3, 4, 9, 200)
    @test tnew[2,3,4,5] ≈ 9.22110552763819
end

@testset "Indexer" begin @test SugarKelp.get_ind(2, arr_lat, 0) == 2; @test SugarKelp.get_ind(2.34, arr_lat, 1) == 2 end


text:
function main_stats(args)
    println("="^70)
    println("Début de l'action stats\n")
    inst = Instance(args[:infile])

    if lg3()
        println("="^70)
        println("Regéneration de l'instance au format ampl")
        println(to_s_long(inst, format="ampl"))
        println("="^70)
        println("Regéneration de l'instance au format alp")
        println(to_s_long(inst, format="alp"))
        println("="^70)
        println("Regéneration de l'instance au format alpx")
        println(to_s_long(inst, format="alpx"))
    end

    println("="^70)
    println(to_s_stats(inst))
    # println("...")
    #
    # println("="^70)
    # println("Vérification de l'inégalité triangulaire de l'instance")
    # print_triangular_inequality_viols(inst; doprint=true)

    println("Fin de l'action stats")
end

main_stats(Args.args)


text:
using Maker

empty!(Maker.TARGETS)
Maker.rm(".maker-cache.jld")

include("example.jl")
include("file.jl")
include("variable.jl")
include("test-targets.jl")
include("utils.jl")
include("glob.jl")
include("function-dependency.jl")

rm(".maker-cache.jld")


text:
using CSV
using MortalityTables

select_filepath = joinpath(@__DIR__, "validate_sel.csv")
select_csv = CSV.File(select_filepath)

# test the select values
for row in select_csv
    table_id, age, dur, val = row
    # load the table
    local table = MortalityTables.table(table_id)
    if !(table.select[age][age + dur - 1] == val)
        error("The values don't match")
    end
end

ultimate_filepath = joinpath(@__DIR__, "validate_ult.csv")
ultimate_csv = CSV.File(ultimate_filepath)

for row in ultimate_csv
    table_id, age, val = row
    # load the table
    local table = MortalityTables.table(table_id)
    if !(table.ultimate[age] == val)
        error("The values don't match")
    end
end

max_stars_repo_name: actuarialopensource/validations

max_issues_repo_name: actuarialopensource/validations

max_forks_repo_name: actuarialopensource/validations


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = (1.0, 1.0, 1.0)
equations = LinearScalarAdvectionEquation3D(advectionvelocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

initial_condition = initial_condition_convergence_test
boundary_condition = BoundaryConditionDirichlet(initial_condition)

boundary_conditions = Dict(
  :all => boundary_condition
)

# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping.
# The mapping will be interpolated at tree level, and then refined without changing
# the geometry interpolant. The original mapping applied to this unstructured mesh
# causes some Jacobians to be negative, which makes the mesh invalid.
function mapping(xi, eta, zeta)
  # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
  # xi = 1.5 * xi_ + 1.5
  # eta = 1.5 * eta_ + 1.5
  # zeta = 1.5 * zeta_ + 1.5

  y = eta + 1/6 * (cos(1.5 * pi * (2 * xi - 3)/3) *
                   cos(0.5 * pi * (2 * eta - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  x = xi + 1/6 * (cos(0.5 * pi * (2 * xi - 3)/3) *
                  cos(2 * pi * (2 * y - 3)/3) *
                  cos(0.5 * pi * (2 * zeta - 3)/3))

  z = zeta + 1/6 * (cos(0.5 * pi * (2 * x - 3)/3) *
                    cos(pi * (2 * y - 3)/3) *
                    cos(0.5 * pi * (2 * zeta - 3)/3))

  return SVector(x, y, z)
end

# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3
mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp",
                              mesh_file)

mesh = P4estMesh{3}(mesh_file, polydeg=3,
                    mapping=mapping,
                    initial_refinement_level=2)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 0.1
ode = semidiscretize(semi, (0.0, 0.1));

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

# The SaveRestartCallback allows to save a file from which a Trixi simulation can be restarted
save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=1.2)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
using JLD2
# using PyCallJLD
# using PyCall
using ScikitLearn: fit!, score, predict, @sk_import
@sk_import linear_model:RidgeClassifier

classifier = RidgeClassifier()

save_object("asdf.jld2", classifier)

my_new = load_object("asdf.jld2")


text:
include("$(pwd())/startup.jl")

file = "/home/wiebke/Schreibtisch/2207_ppb.csv"

dataset, header = readdlm(file, header = true)
datetime = unix2datetime.(dataset[:,1])
for i in findall( x -> occursin("C10H16", x), header)
	semilogy(datetime, dataset[:,i[2]], label = header[i[2]])
end
legend()


text:
# TODO: find/replace all [PkgName] with CamelCase

# [ADD BACKEND WEBSITE]

function _initialize_backend(::[PkgName]Backend; kw...)
    @eval begin
        import [PkgName]
        export [PkgName]
        # todo: other initialization that needs to be eval-ed
    end
    # todo: other initialization
end

# ---------------------------------------------------------------------------

# Create the window/figure for this backend.
function _create_backend_figure(plt::Plot{[PkgName]Backend})
    nothing
end

# # this is called early in the pipeline, use it to make the plot current or something
# function _prepare_plot_object(plt::Plot{[PkgName]Backend})
# end

# Set up the subplot within the backend object.
function _initialize_subplot(plt::Plot{[PkgName]Backend}, sp::Subplot{[PkgName]Backend})
end

# ---------------------------------------------------------------------------

# Add one series to the underlying backend object.
function _series_added(plt::Plot{[PkgName]Backend}, series::Series)
end

# When series data is added/changed, this callback can do dynamic updates to the backend object.
# note: if the backend rebuilds the plot from scratch on display, then you might not do anything here.
function _series_updated(plt::Plot{[PkgName]Backend}, series::Series)
end

# ---------------------------------------------------------------------------

# called just before updating layout bounding boxes... in case you need to prep
# for the calcs
function _before_layout_calcs(plt::Plot{[PkgName]Backend})
end

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot{[PkgName]Backend})
    sp.minpad = (20mm, 5mm, 2mm, 10mm)
end

# ----------------------------------------------------------------

# Override this to update plot items (title, xlabel, etc), and add annotations (d[:annotations])
function _update_plot_object(plt::Plot{[PkgName]Backend})
end

# ----------------------------------------------------------------

# Write a png to io.  You could define methods for:
    # "application/eps"         => "eps",
    # "image/eps"               => "eps",
    # "application/pdf"         => "pdf",
    # "image/png"               => "png",
    # "application/postscript"  => "ps",
    # "image/svg+xml"           => "svg"
function _show(io::IO, ::MIME"image/png", plt::Plot{[PkgName]Backend})
end

# Display/show the plot (open a GUI window, or browser page, for example).
function _display(plt::Plot{[PkgName]Backend})
end


text:
module TestCommutativeDiagrams
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra.FinSets
using Catlab.CategoricalAlgebra.CommutativeDiagrams

A, B, C, D = Ob(FreeCategory, :A, :B, :C, :D)

# Commutative squares
#####################

l, t, b, r = Hom(:lef, A,B), Hom(:top, A, C), Hom(:bot, B, D), Hom(:rht, C,D)
sq1 = SquareDiagram(t, b, l, r)

@test_throws AssertionError composeH(sq1, sq1)

l, t, b, r = Hom(:lef, A,B), Hom(:top, A, A), Hom(:bot, B, B), Hom(:rht, A,B)
rr = Hom(:rr, A,B)
sq2 = SquareDiagram(t, b, l, r)
sq3 = composeH(sq2, SquareDiagram(t, b, r, rr))
@test left(sq3)   == l
@test top(sq3)    == compose(t,t)
@test bottom(sq3) == compose(b,b)
@test right(sq3)  == rr

@test_throws AssertionError composeV(sq2, sq2)

ll = Hom(:ll, B, A)
rr = Hom(:rr, B, A)
sq4 = composeV(sq2, SquareDiagram(b, t, ll, rr))
@test left(sq4)    == compose(l, ll)
@test right(sq4)  == compose(r, rr)
@test top(sq4)   == t
@test bottom(sq4) == t

@test hom(sq4) == [top(sq4), bottom(sq4), left(sq4), right(sq4)]

# Double category of squares
#---------------------------

t  = FinFunction([1,2,4], 5)
l  = FinFunction([1,2,3], 3)
b  = FinFunction([2,3,1], 3)
f  = FinFunction([2,3,2,1,3], 3)
α = SquareDiagram(t, b, l, f)

t₁ = FinFunction([1,1,1,1,1], 1)
r  = FinFunction([1], 1)
b₁ = FinFunction([1,1,1], 1)
β = SquareDiagram(t₁, b₁, f, r)

γ = composeH(α, β)
@test top(γ) == t⋅t₁
@test bottom(γ) == b⋅b₁
@test left(γ) == left(α)
@test right(γ) == right(β)

t = FinFunction([1,2,4], 5)
b = FinFunction([2,1,3], 3)
l = FinFunction([2,1,3], 3)
r = FinFunction([1,2,2,3,3], 3)
α = SquareDiagram(t,b,l,r)
β = SquareDiagram(FinFunction([1,2,1,2,1], 2),
                  FinFunction([1,1,1],1),
                  FinFunction([1,2,2,3,3], 3),
                  FinFunction([1,1],1)
)
γ = composeH(α, β)
@test collect(top(γ)) == [1,2,2]
@test collect(bottom(γ)) == [1,1,1]

α = SquareDiagram(l,r, t,b)
β = SquareDiagram(FinFunction([1,2,2,3,3], 3),
                  FinFunction([1,1],1),
                  FinFunction([1,2,1,2,1], 2),
                  FinFunction([1,1,1],1)
)
γ = composeV(α, β)
@test_throws AssertionError composeH(α, β)
@test collect(left(γ)) == [1,2,2]
@test collect(right(γ)) == [1,1,1]

@test collect(idH(FinSet(3))) == collect(1:3)
@test collect(idV(FinSet(5))) == collect(1:5)
@test collect(idH(FinSet(0))) == []

@test collect(composeH(FinFunction([2,1], 3), FinFunction([2, 1, 3]))) == [1,2]
@test collect(composeV(FinFunction([2,1], 3), FinFunction([2, 1, 3]))) == [1,2]

end


text:
module UnitfulLsqFit

using Unitful:
    AbstractQuantity, Unit, Units, NoDims, dimension, unit, ustrip, @u_str, @unit, register
import LsqFit: curve_fit, LsqFitResult
import Base: *

function __init__()
    return register(UnitfulLsqFit)
end

# One {{{
@unit one "" One 1 false
register(UnitfulLsqFit)
*(q::AbstractQuantity, ::Units{(Unit{:One,NoDims}(0, 1),),NoDims,nothing}) = q
function *(
    a::AbstractQuantity, b::T
) where {
    T<:AbstractQuantity{<:Number,NoDims,<:Units{(Unit{:One,NoDims}(0, 1),),NoDims,nothing}},
}
    return a * b.val
end
function *(
    b::T, a::AbstractQuantity
) where {
    T<:AbstractQuantity{<:Number,NoDims,<:Units{(Unit{:One,NoDims}(0, 1),),NoDims,nothing}},
}
    return b.val * a
end
# One }}}

function curve_fit(
    model,
    xdata::AbstractArray{<:AbstractQuantity},
    ydata::AbstractArray,
    p0::AbstractArray{<:AbstractQuantity};
    normalize=false,
    kwargs...,
)
    # unit check
    dimension(model(first(xdata), p0)...) === dimension(first(ydata)) ||
        error("Model and ydata dimensions incompatible")
    xunit = unit(first(xdata))
    yunit = unit(first(ydata))
    punits = unit.(p0)

    X = ustrip.(xunit, xdata)
    Y = ustrip.(yunit, ydata)
    P = ustrip.(punits, p0)

    auxmodel(x, p) = ustrip.(yunit, model(x .* xunit, p .* punits))
    auxfit = curve_fit(auxmodel, X, Y, P; kwargs...)

    return LsqFitResult(
        auxfit.param .* punits,
        auxfit.resid * yunit,
        auxfit.jacobian * yunit ./ reshape(punits, 1, :),
        auxfit.converged,
        auxfit.wt,
    )
end

# Orthogonals {{{
function curve_fit( # 1 0
    model,
    xdata::AbstractArray{<:AbstractQuantity},
    ydata::AbstractArray,
    p0::AbstractArray;
    kwargs...,
)
    return curve_fit(model, xdata, ydata * u"one", p0 * u"one"; kwargs...)
end

function curve_fit( # 0 1
    model,
    xdata::AbstractArray,
    ydata::AbstractArray,
    p0::AbstractArray{<:AbstractQuantity};
    kwargs...,
)
    return curve_fit(model, xdata * u"one", ydata * u"one", p0; kwargs...)
end

end


text:
using ArgParse
using CUDAapi
using Dates
using FileIO
using JLD2
using LinearAlgebra
using Logging
using MPI
using Printf
using Requires
using CLIMAParameters

using ..Atmos
using ..Callbacks
using ..ColumnwiseLUSolver
using ..ConfigTypes
using ..Courant
using ..Diagnostics
using ..DGmethods
using ..DGmethods: vars_state, vars_aux, update_aux!
using ..DGmethods.NumericalFluxes
using ..GenericCallbacks
using ..HydrostaticBoussinesq
using ..Mesh.Grids
using ..Mesh.Topologies
using ..MoistThermodynamics
using ..MPIStateArrays
using ..ODESolvers
using ..TicToc
using ..VariableTemplates
using ..VTK

# Note that the initial values specified here are overwritten by the
# command line argument defaults in `parse_commandline()`.
Base.@kwdef mutable struct CLIMA_Settings
    disable_gpu::Bool = false
    show_updates::Bool = true
    update_interval::Int = 60
    enable_diagnostics::Bool = false
    diagnostics_interval::Int = 10000
    enable_vtk::Bool = false
    vtk_interval::Int = 10000
    enable_wall_clock::Bool = false
    wall_clock_interval::Int = 10000
    monitor_courant_numbers::Bool = false
    monitor_courant_interval::Int = 10
    checkpoint_walltime::Int = -1
    checkpoint_keep_one::Bool = true
    checkpoint_at_end::Bool = false
    checkpoint_dir::String = "checkpoint"
    restart_from_num::Int = -1
    log_level::String = "INFO"
    output_dir::String = "output"
    integration_testing::Bool = false
    array_type
end

const Settings = CLIMA_Settings(array_type = Array)

array_type() = Settings.array_type

const cuarray_pkgid =
    Base.PkgId(Base.UUID("3a865a2d-5b23-5a0f-bc46-62713ec82fae"), "CuArrays")

@init begin
    if get(ENV, "CLIMA_GPU", "") != "false" && CUDAapi.has_cuda_gpu()
        CuArrays = Base.require(cuarray_pkgid)
    end
end

function _init_array(::Type{Array})
    return nothing
end

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CuArrays, .CuArrays.CUDAdrv, .CuArrays.CUDAnative
    function _init_array(::Type{CuArray})
        comm = MPI.COMM_WORLD
        # allocate GPUs among MPI ranks
        local_comm = MPI.Comm_split_type(
            comm,
            MPI.MPI_COMM_TYPE_SHARED,
            MPI.Comm_rank(comm),
        )
        # we intentionally oversubscribe GPUs for testing: may want to disable this for production
        CUDAnative.device!(MPI.Comm_rank(local_comm) % length(devices()))
        CuArrays.allowscalar(false)
        return nothing
    end
end

function gpu_allowscalar(val)
    if haskey(Base.loaded_modules, CLIMA.cuarray_pkgid)
        Base.loaded_modules[CLIMA.cuarray_pkgid].allowscalar(val)
    end
    return nothing
end

"""
    parse_commandline()
"""
function parse_commandline()
    exc_handler =
        isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler
    s = ArgParseSettings(exc_handler = exc_handler)

    @add_arg_table! s begin
        "--disable-gpu"
        help = "do not use the GPU"
        action = :store_true
        "--no-show-updates"
        help = "do not show simulation updates"
        action = :store_true
        "--update-interval"
        help = "interval in seconds for showing simulation updates"
        arg_type = Int
        default = 60
        "--enable-diagnostics"
        help = "enable the collection of diagnostics to <output-dir>"
        action = :store_true
        "--diagnostics-interval"
        help = "override the interval for gathering diagnostics (in simulation steps)"
        arg_type = Int
        default = 10000
        "--enable-vtk"
        help = "output VTK to <output-dir> every <vtk-interval> simulation steps"
        action = :store_true
        "--vtk-interval"
        help = "interval in simulation steps for VTK output"
        arg_type = Int
        default = 10000
        "--enable-wall-clock"
        help = "output wall-clock time per time-step every <wall-clock-interval> simulation steps"
        action = :store_true
        "--wall-clock-interval"
        help = "interval in simulation steps for wall-clock time per time-step output"
        arg_type = Int
        default = 10000
        "--monitor-courant-numbers"
        help = "output acoustic, advective, and diffusive Courant numbers"
        action = :store_true
        "--monitor-courant-interval"
        help = "interval in Courant number calculations"
        arg_type = Int
        default = 10
        "--checkpoint-walltime"
        help = "interval in seconds at which to create a checkpoint"
        arg_type = Int
        default = -1
        "--checkpoint-keep-all"
        help = "keep all checkpoints (instead of just the most recent)"
        action = :store_true
        "--checkpoint-at-end"
        help = "create a checkpoint at the end of the simulation"
        action = :store_true
        "--checkpoint-dir"
        help = "the directory in which to store checkpoints"
        arg_type = String
        default = "checkpoint"
        "--restart-from-num"
        help = "checkpoint number from which to restart (in <checkpoint-dir>)"
        arg_type = Int
        default = -1
        "--log-level"
        help = "set the log level to one of debug/info/warn/error"
        arg_type = String
        default = "info"
        "--output-dir"
        help = "directory for output data"
        arg_type = String
        default = "output"
        "--integration-testing"
        help = "enable integration testing"
        action = :store_true
    end

    return parse_args(s)
end

"""
    CLIMA.init(; disable_gpu=false)

Initialize MPI, allocate GPUs among MPI ranks if using GPUs, parse command
line arguments for CLIMA, and return a Dict of any additional command line
arguments.
"""
function init(; disable_gpu = false)
    # initialize MPI
    if !MPI.Initialized()
        MPI.Init()
    end

    # set up timing mechanism
    tictoc()

    # parse command line arguments
    try
        parsed_args = parse_commandline()
        Settings.disable_gpu = disable_gpu || parsed_args["disable-gpu"]
        Settings.show_updates = !parsed_args["no-show-updates"]
        Settings.update_interval = parsed_args["update-interval"]
        Settings.enable_diagnostics = parsed_args["enable-diagnostics"]
        Settings.diagnostics_interval = parsed_args["diagnostics-interval"]
        Settings.enable_vtk = parsed_args["enable-vtk"]
        Settings.vtk_interval = parsed_args["vtk-interval"]
        Settings.enable_wall_clock = parsed_args["enable-wall-clock"]
        Settings.wall_clock_interval = parsed_args["wall-clock-interval"]
        Settings.output_dir = parsed_args["output-dir"]
        Settings.monitor_courant_numbers =
            parsed_args["monitor-courant-numbers"]
        Settings.monitor_courant_interval =
            parsed_args["monitor-courant-interval"]
        Settings.checkpoint_walltime = parsed_args["checkpoint-walltime"]
        Settings.checkpoint_keep_one = !parsed_args["checkpoint-keep-all"]
        Settings.checkpoint_at_end = parsed_args["checkpoint-at-end"]
        Settings.checkpoint_dir = parsed_args["checkpoint-dir"]
        Settings.restart_from_num = parsed_args["restart-from-num"]
        Settings.integration_testing = parsed_args["integration-testing"]
        Settings.log_level = uppercase(parsed_args["log-level"])
    catch
        Settings.disable_gpu = disable_gpu
    end

    # set up the array type appropriately depending on whether we're using GPUs
    if !Settings.disable_gpu &&
       get(ENV, "CLIMA_GPU", "") != "false" &&
       CUDAapi.has_cuda_gpu()
        atyp = CuArrays.CuArray
    else
        atyp = Array
    end
    _init_array(atyp)
    Settings.array_type = atyp

    # create the output directory if needed
    if Settings.enable_diagnostics || Settings.enable_vtk
        mkpath(Settings.output_dir)
    end
    if Settings.checkpoint_walltime > 0 || Settings.checkpoint_at_end
        mkpath(Settings.checkpoint_dir)
    end

    # set up logging
    loglevel = Settings.log_level == "DEBUG" ? Logging.Debug :
        Settings.log_level == "WARN" ? Logging.Warn :
        Settings.log_level == "ERROR" ? Logging.Error : Logging.Info
    # TODO: write a better MPI logging back-end and also integrate Dlog for large scale
    logger_stream = MPI.Comm_rank(MPI.COMM_WORLD) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    return nothing
end

include("driver_configs.jl")
include("solver_configs.jl")
include("diagnostics_configs.jl")

"""
    CLIMA.invoke!(solver_config::SolverConfiguration;
                  diagnostics_config       = nothing,
                  user_callbacks           = (),
                  check_euclidean_distance = false,
                  adjustfinalstep          = false,
                  user_info_callback       = (init)->nothing)

Run the simulation defined by the `solver_config`.

Keyword Arguments:

The `user_callbacks` are passed to the ODE solver as callback functions; see
[`ODESolvers.solve!]@ref().

If `check_euclidean_distance` is `true, then the Euclidean distance between the
final solution and initial condition function evaluated with
`solver_config.timeend` is reported.

The value of 'adjustfinalstep` is passed to the ODE solver; see
[`ODESolvers.solve!]@ref().

The function `user_info_callback` is called after the default info callback
(which is called every `Settings.update_interval` seconds of wallclock time).
The single input argument `init` is `true` when the callback is called
called for initialization before time stepping begins and `false` when called
during the actual ODE solve; see [`GenericCallbacks`](@ref) and
[`ODESolvers.solve!]@ref().
"""
function invoke!(
    solver_config::SolverConfiguration;
    diagnostics_config = nothing,
    user_callbacks = (),
    check_euclidean_distance = false,
    adjustfinalstep = false,
    user_info_callback = (init) -> nothing,
)
    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balancelaw
    Q = solver_config.Q
    FT = eltype(Q)
    timeend = solver_config.timeend
    init_on_cpu = solver_config.init_on_cpu
    init_args = solver_config.init_args
    solver = solver_config.solver

    # set up callbacks
    callbacks = ()
    if Settings.show_updates
        # set up the information callback
        upd_starttime = Ref(now())
        cbinfo = GenericCallbacks.EveryXWallTimeSeconds(
            Settings.update_interval,
            mpicomm,
        ) do (init = false)
            if init
                upd_starttime[] = now()
            else
                runtime = Dates.format(
                    convert(Dates.DateTime, Dates.now() - upd_starttime[]),
                    Dates.dateformat"HH:MM:SS",
                )
                energy = norm(Q)
                @info @sprintf(
                    """
         Update
             simtime = %8.2f / %8.2f
             runtime = %s
             norm(Q) = %.16e""",
                    ODESolvers.gettime(solver),
                    timeend,
                    runtime,
                    energy
                )
                isnan(energy) && error("norm(Q) is NaN")
            end
            user_info_callback(init)
        end
        callbacks = (callbacks..., cbinfo)
    end

    dgn_starttime = ""
    if Settings.enable_diagnostics && diagnostics_config !== nothing
        dgn_starttime = replace(string(now()), ":" => ".")
        Diagnostics.init(mpicomm, dg, Q, dgn_starttime, Settings.output_dir)

        # set up a callback for each diagnostics group
        dgncbs = ()
        for dgngrp in diagnostics_config.groups
            if Settings.diagnostics_interval > 0
                interval = Settings.diagnostics_interval
            else
                interval = dgngrp.interval
            end
            fn = GenericCallbacks.EveryXSimulationSteps(
                interval,
            ) do (init = false)
                currtime = ODESolvers.gettime(solver)
                @info @sprintf(
                    """Diagnostics: %s
                    %s at %s""",
                    dgngrp.name,
                    init ? "initializing" : "collecting",
                    string(currtime)
                )
                dgngrp(currtime, init = init)
                nothing
            end
            dgncbs = (dgncbs..., fn)
        end
        callbacks = (callbacks..., dgncbs...)
    end

    if Settings.enable_vtk
        # set up VTK output callback
        vtknum = Ref(1)
        cbvtk =
            GenericCallbacks.EveryXSimulationSteps(Settings.vtk_interval) do (
                init = false
            )
                save_vtk(solver_config, vtknum[])
                vtknum[] += 1
                nothing
            end
        callbacks = (callbacks..., cbvtk)
    end

    if Settings.enable_wall_clock
        cbwall = Callbacks.wall_clock_time_per_time_step(
            Settings.wall_clock_interval,
            Settings.array_type,
            mpicomm,
        )
        callbacks = (callbacks..., cbwall)
    end

    if Settings.monitor_courant_numbers
        # set up the callback for Courant number calculations
        cbcfl =
            GenericCallbacks.EveryXSimulationSteps(Settings.monitor_courant_interval) do (
                init = false
            )
                show_courant(solver_config)
                nothing
            end
        callbacks = (callbacks..., cbcfl)
    end

    if Settings.checkpoint_walltime > 0
        cpnum = Ref(1)
        cbcheckpoint = GenericCallbacks.EveryXWallTimeSeconds(
            Settings.checkpoint_walltime,
            mpicomm,
        ) do (init = false)
            write_checkpoint(solver_config, cpnum[])
            if Settings.checkpoint_keep_one
                rm_checkpoint(solver_config, cpnum[] - 1)
            end
            cpnum[] += 1
            nothing
        end
        callbacks = (callbacks..., cbcheckpoint)
    end

    callbacks = (user_callbacks..., callbacks...)

    # initial condition norm
    eng0 = norm(Q)
    @info @sprintf(
        """
Starting %s
    dt              = %.5e
    timeend         = %8.2f
    number of steps = %d
    norm(Q)         = %.16e""",
        solver_config.name,
        solver_config.dt,
        solver_config.timeend,
        solver_config.numberofsteps,
        eng0
    )

    # run the simulation
    @tic solve!
    solve!(
        Q,
        solver;
        timeend = timeend,
        callbacks = callbacks,
        adjustfinalstep = adjustfinalstep,
    )
    @toc solve!

    # write checkpoint
    if Settings.checkpoint_at_end
        write_checkpoint(solver_config, solver_config.numberofsteps)
    end

    # fini diagnostics groups
    if Settings.enable_diagnostics
        currtime = ODESolvers.gettime(solver)
        for dgngrp in diagnostics_config.groups
            dgngrp(currtime, fini = true)
        end
    end

    engf = norm(Q)
    @info @sprintf(
        """Finished
        norm(Q)            = %.16e
        norm(Q) / norm(Q₀) = %.16e
        norm(Q) - norm(Q₀) = %.16e""",
        engf,
        engf / eng0,
        engf - eng0
    )

    if check_euclidean_distance
        Qe =
            init_ode_state(dg, timeend, init_args...; init_on_cpu = init_on_cpu)
        engfe = norm(Qe)
        errf = euclidean_distance(Q, Qe)
        @info @sprintf(
            """Euclidean distance
            norm(Q - Qe)            = %.16e
            norm(Q - Qe) / norm(Qe) = %.16e""",
            errf,
            errf / engfe
        )
    end

    return engf / eng0
end

function save_vtk(solver_config::SolverConfiguration, num::Int)
    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balancelaw
    Q = solver_config.Q
    FT = eltype(Q)

    vprefix = @sprintf(
        "%s_mpirank%04d_num%04d",
        solver_config.name,
        MPI.Comm_rank(mpicomm),
        num
    )
    outprefix = joinpath(Settings.output_dir, vprefix)

    statenames = flattenednames(vars_state(bl, FT))
    auxnames = flattenednames(vars_aux(bl, FT))

    writevtk(outprefix, Q, dg, statenames, dg.auxstate, auxnames)

    # Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        # name of the pvtu file
        pprefix = @sprintf("%s_num%04d", solver_config.name, num)
        pvtuprefix = joinpath(Settings.output_dir, pprefix)

        # name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_num%04d", solver_config.name, i - 1, num)
        end
        writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...))
    end
end

function show_courant(solver_config::SolverConfiguration)
    Δt = solver_config.dt
    c_v = DGmethods.courant(
        nondiffusive_courant,
        solver_config,
        direction = VerticalDirection(),
    )
    c_h = DGmethods.courant(
        nondiffusive_courant,
        solver_config,
        direction = HorizontalDirection(),
    )
    ca_v = DGmethods.courant(
        advective_courant,
        solver_config,
        direction = VerticalDirection(),
    )
    ca_h = DGmethods.courant(
        advective_courant,
        solver_config,
        direction = HorizontalDirection(),
    )
    cd_v = DGmethods.courant(
        diffusive_courant,
        solver_config,
        direction = VerticalDirection(),
    )
    cd_h = DGmethods.courant(
        diffusive_courant,
        solver_config,
        direction = HorizontalDirection(),
    )
    simtime = ODESolvers.gettime(solver_config.solver)
    @info @sprintf """
Courant numbers at simtime: %8.2f, Δt = %8.2f s
    Acoustic (vertical) Courant number    = %.2g
    Acoustic (horizontal) Courant number  = %.2g
    Advection (vertical) Courant number   = %.2g
    Advection (horizontal) Courant number = %.2g
    Diffusion (vertical) Courant number   = %.2g
    Diffusion (horizontal) Courant number = %.2g
    """ simtime Δt c_v c_h ca_v ca_h cd_v cd_h

    return nothing
end

function write_checkpoint(solver_config::SolverConfiguration, num::Int)
    nm = replace(solver_config.name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(solver_config.mpicomm),
        num,
    )
    cfull = joinpath(Settings.checkpoint_dir, cname)
    @info @sprintf(
        """
Checkpoint
    saving to %s""",
        cfull
    )

    dg = solver_config.dg
    Q = solver_config.Q
    if Array ∈ typeof(Q).parameters
        h_Q = Q.realdata
        h_aux = dg.auxstate.realdata
    else
        h_Q = Array(Q.realdata)
        h_aux = Array(dg.auxstate.realdata)
    end
    t = ODESolvers.gettime(solver_config.solver)
    @save cfull h_Q h_aux t

    return nothing
end

function rm_checkpoint(solver_config::SolverConfiguration, num::Int)
    nm = replace(solver_config.name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(solver_config.mpicomm),
        num,
    )
    rm(joinpath(Settings.checkpoint_dir, cname), force = true)

    return nothing
end

function read_checkpoint(name::String, mpicomm::MPI.Comm, num::Int)
    nm = replace(name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(mpicomm),
        num,
    )
    cfull = joinpath(Settings.checkpoint_dir, cname)
    if !isfile(cfull)
        error("Cannot restore from checkpoint in %s, file not found")
    end

    @load cfull h_Q h_aux t

    ArrayType = array_type()
    return (ArrayType(h_Q), ArrayType(h_aux), t)
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# imports
using FMI
using FMIFlux
using Flux
using DifferentialEquations: Tsit5
import Plots

pathFMU = joinpath(dirname(@__FILE__), "../model/SpringPendulumExtForce1D.fmu")
println("FMU path: ", pathFMU)

tStart = 0.0
tStep = 0.01
tStop = 5.0
tSave = tStart:tStep:tStop

referenceFMU = fmiLoad(pathFMU)
fmiInstantiate!(referenceFMU; loggingOn=false)
fmiInfo(referenceFMU)

fmiSetupExperiment(referenceFMU, tStart, tStop)
fmiSetReal(referenceFMU, "mass_s0", 1.3)   # increase amplitude, invert phase
fmiEnterInitializationMode(referenceFMU)
fmiExitInitializationMode(referenceFMU)

x₀ = fmiGetContinuousStates(referenceFMU)

vrs = ["mass.s", "mass.v", "mass.a"]
_, referenceSimData = fmiSimulate(referenceFMU, tStart, tStop; recordValues=vrs, setup=false, reset=false, saveat=tSave)
fmiPlot(referenceFMU, vrs, referenceSimData)

posReference = collect(data[1] for data in referenceSimData.saveval)
velReference = collect(data[2] for data in referenceSimData.saveval)
accReference = collect(data[3] for data in referenceSimData.saveval)

fmiReset(referenceFMU)
defaultFMU = referenceFMU

fmiSetupExperiment(defaultFMU, tStart, tStop)
fmiEnterInitializationMode(defaultFMU)
fmiExitInitializationMode(defaultFMU)

x₀ = fmiGetContinuousStates(defaultFMU)

_, defaultSimData = fmiSimulate(defaultFMU, tStart, tStop; recordValues=vrs, setup=false, reset=false, saveat=tSave)
fmiPlot(defaultFMU, vrs, defaultSimData)

posDefault = collect(data[1] for data in defaultSimData.saveval)
velDefault = collect(data[2] for data in defaultSimData.saveval)
accDefault = collect(data[3] for data in defaultSimData.saveval)

function extForce(t)
    return [0.0]
end

# loss function for training
function lossSum()
    solution = csNeuralFMU(extForce, tStep)

    accNet = collect(data[1] for data in solution)

    Flux.Losses.mse(accReference, accNet)
end

# callback function for training
global counter = 0
function callb()
    global counter += 1

    if counter % 20 == 1
        avgLoss = lossSum()
        @info "Loss [$counter]: $(round(avgLoss, digits=5))"
    end
end

# NeuralFMU setup
numInputs = length(defaultFMU.modelDescription.inputValueReferences)
numOutputs = length(defaultFMU.modelDescription.outputValueReferences)

net = Chain(inputs -> fmiInputDoStepCSOutput(defaultFMU, tStep, inputs),
            Dense(numOutputs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numOutputs))

csNeuralFMU = CS_NeuralFMU(defaultFMU, net, (tStart, tStop); saveat=tSave);

solutionBefore = csNeuralFMU(extForce, tStep)
Plots.plot(tSave, collect(data[1] for data in solutionBefore), label="acc CS-NeuralFMU", linewidth=2)

# train
paramsNet = Flux.params(csNeuralFMU)

optim = ADAM()
Flux.train!(lossSum, paramsNet, Iterators.repeated((), 300), optim; cb=callb)

# plot results mass.a
solutionAfter = csNeuralFMU(extForce, tStep)

fig = Plots.plot(xlabel="t [s]", ylabel="mass acceleration [m/s^2]", linewidth=2,
                 xtickfontsize=12, ytickfontsize=12,
                 xguidefontsize=12, yguidefontsize=12,
                 legendfontsize=8, legend=:topright)

accNeuralFMU = collect(data[1] for data in solutionAfter)

Plots.plot!(fig, tSave, accDefault, label="defaultFMU", linewidth=2)
Plots.plot!(fig, tSave, accReference, label="referenceFMU", linewidth=2)
Plots.plot!(fig, tSave, accNeuralFMU, label="CS-NeuralFMU (300 eps.)", linewidth=2)
fig

fmiUnload(defaultFMU)


text:
export DP

# Stick breaking construction.

@pp function pickstick(sticks, j)
    flip(sticks(j)) ? j : pickstick(sticks, j+1)
end

@pp function makesticks(alpha)
    local sticks = mem(index -> ~Beta(1.0, alpha))
    () -> pickstick(sticks,1)
end

@pp function DP(thunk, alpha)
    local atoms = mem(index -> thunk())
    local sticks = makesticks(alpha)
    () -> atoms(sticks())
end


text:
using Test, LinearMaps, LinearAlgebra

@testset "function maps" begin
    N = 100
    function myft(v::AbstractVector)
        # not so fast fourier transform
        N = length(v)
        w = zeros(complex(eltype(v)), N)
        for k = 1:N
            kappa = (2*(k-1)/N)*pi
            for n = 1:N
                w[k] += v[n]*exp(kappa*(n-1)*im)
            end
        end
        return w
    end
    MyFT = @inferred LinearMap{ComplexF64}(myft, N) / sqrt(N)
    U = Matrix(MyFT) # will be a unitary matrix
    @test @inferred U'U ≈ Matrix{eltype(U)}(I, N, N)

    CS = @inferred LinearMap(cumsum, 2)
    @test size(CS) == (2, 2)
    @test @inferred !issymmetric(CS)
    @test @inferred !ishermitian(CS)
    @test @inferred !isposdef(CS)
    @test @inferred !(LinearMaps.ismutating(CS))
    @test @inferred Matrix(CS) == [1. 0.; 1. 1.]
    @test @inferred Array(CS) == [1. 0.; 1. 1.]
    CS = @inferred LinearMap(cumsum, 10; ismutating=false)
    v = rand(10)
    cv = cumsum(v)
    @test CS * v == cv
    @test *(CS, v) == cv
    @test_throws ErrorException CS' * v
    CS = @inferred LinearMap(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false)
    cv = cumsum(v)
    @test @inferred CS * v == cv
    @test @inferred *(CS, v) == cv
    @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))

    CS! = @inferred LinearMap(cumsum!, 10; ismutating=true)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException transpose(CS!) * v

    CS! = @inferred LinearMap{ComplexF64}(cumsum!, 10; ismutating=true)
    v = rand(ComplexF64, 10)
    cv = cumsum(v)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException adjoint(CS!) * v
    CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), adjoint(CS), v) == reverse!(cumsum(reverse(v)))

    # Test fallback methods:
    L = @inferred LinearMap(x -> x, x -> x, 10)
    v = randn(10)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v
    L = @inferred LinearMap{ComplexF64}(x -> x, x -> x, 10)
    v = rand(ComplexF64, 10)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v
end


text:
# Define a custom Normal distribution, for illustrative puspose.
struct CustomNormal{T<:Real} <: Distributions.ContinuousUnivariateDistribution
    m::T
end
CustomNormal() = CustomNormal(0)

Distributions.rand(rng::AbstractRNG, d::CustomNormal) = d.m + randn(rng)


text:
using Flux
using Flux: Chain, gradient, logitcrossentropy, params
using Flux.Data: DataLoader
using Flux.Optimise: Momentum
using JLD

"""
    generate_permutations(model, dataset, [optimizer])

Generate a JLD permutations file which shifts high-entropy weight at the end of
the array. To do so, a model and dataset are trained on a few eopchs to generate
the entropy pattern of the weights.
"""
function generate_permutations(model::Chain, dataset::DataLoader; optimizer=Momentum(0.01, 0.5))
    loss(x, y) = logitcrossentropy(model(x), y)

    num_epochs = 5
    weights = Vector{Vector{Float32}}(undef, num_epochs+1)
    weights[1] = flatten_model(model)

    for i = 2:num_epochs+1
        println("epochs [$(i-1)/$(num_epochs)]")
        local train_loss = 0.0f0
        for (i, batch) in enumerate(dataset)
            grads = gradient(params(model)) do
                train_loss += loss(batch...)
            end
            Flux.update!(optimizer, params(model), grads)
        end
        train_loss /= length(dataset)
        @show train_loss

        weights[i] = flatten_model(model)
    end

    weights_diff = sum([abs.(weights[i-1] - weights[i]) for i in 2:num_epochs+1]) / num_epochs

    perms = sortperm(weights_diff)

    JLD.save("permutations.jld", "permutations", perms)
end


text:
using Test
@testset "RandomNumbers" begin
    for (i, testfile) in enumerate((
        "generic.jl",
        "wrapped_rng.jl",
        "randfloat.jl",
        "PCG/runtests.jl",
        "MersenneTwisters/runtests.jl",
        "Xorshifts/runtests.jl"))
        @eval module $(Symbol("T$i"))
            include("common.jl")
            include($testfile)
        end
    end
end


text:
function normalise01(x::AbstractArray)
    μ = mean(x)
    return (x .- μ) ./ (stdm(x, μ; corrected=false) .+ eps(eltype(x)))
end

Zygote.@adjoint function stdm(x, μ; corrected::Bool=true)
    σ = stdm(x, μ; corrected=corrected)
    function stdm_sensitivity(Δ)
        xμ = x .- μ
        ΔdivNσ = Δ ./ (σ .* (length(x) - corrected))
        return (xμ .* ΔdivNσ, zero(μ))
    end
    return σ, stdm_sensitivity
end


text:
using OrdinaryDiffEq, DiffEqBase, Test, DiffEqProblemLibrary, DiffEqDevTools,
      Random

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_bigfloatlinear,
                                               prob_ode_bigfloat2Dlinear,
                                               prob_ode_2Dlinear

## Convergence Testing
println("Convergence Test on Linear")

testTol = 1
prob = prob_ode_2Dlinear
println("Feagin RKs")
dts = (1//2) .^ (4:-1:2)
sol = solve(prob,Feagin10(),dt=dts[1])
prob = remake(prob_ode_bigfloat2Dlinear,tspan=(big(0)//1,big(1)//1))
sol = solve(prob,Feagin10(),dt=dts[1])

prob = remake(prob_ode_bigfloat2Dlinear,tspan=(big(0.0),big(1.0)))
dts = (1//2) .^ (4:-1:2)
sim = test_convergence(dts,prob,Feagin10())
@test abs(sim.𝒪est[:final]-8) < testTol #Lowered due to low test dt

sim = test_convergence(dts,prob,Feagin12())
@test abs(sim.𝒪est[:final]-12) < testTol

sim = test_convergence(dts,prob,Feagin14())
@test abs(sim.𝒪est[:final]-15) < testTol #Upped to 15 for test

prob = prob_ode_bigfloatlinear

dts = (1//2) .^ (6:-1:3)
sim = test_convergence(dts,prob,Feagin10())
@test abs(sim.𝒪est[:final]-10) < testTol

dts = (1//2) .^ (4:-1:2)
sim = test_convergence(dts,prob,Feagin12())
@test abs(sim.𝒪est[:final]-12) < testTol

sim = test_convergence(dts,prob,Feagin14())
@test abs(sim.𝒪est[:final]-15) < testTol #Upped to 15 for test

prob = prob_ode_bigfloat2Dlinear

#compile
sol =solve(prob,Feagin10(),dt=dts[1])
sol =solve(prob,Feagin12(),dt=dts[1])
sol =solve(prob,Feagin14(),dt=dts[1])

#test
@time sol =solve(prob,Feagin10(),dt=dts[1])
@time sol =solve(prob,Feagin12(),dt=dts[1])
@time sol =solve(prob,Feagin14(),dt=dts[1])


text:
# Example for 2D FWI of Camembert model using SPG
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, opesciSLIM.TimeModeling, NLopt, opesciSLIM.SLIM_optim, HDF5

# Load camembert model
~isfile("camembert.h5") && run(`wget ftp://ftp.slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/camembert.h5`)
n,d,o,m = read(h5open("camembert.h5","r"), "n", "d", "o", "m")

# Set up model structure
model = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m)

# Bound constraints and initial model
m0 = ones(Float32,model.n)*0.25f0
vmin = ones(Float32,model.n) + 0.3f0
vmax = ones(Float32,model.n) + 4f0
mmin = vec((1f0./vmax).^2)
mmax = vec((1f0./vmin).^2)

# Setup info and model structure
nsrc = 10	# number of sources
model0 = Model(model.n,model.d,model.o,m0)

## Set up receiver geometry
nxrec = 100
xrec = ones(Float32,100)*50f0
yrec = 0f0
zrec = linspace(50f0,1950f0,nxrec)

# receiver sampling and recording time
timeR = 3000f0	# receiver recording time [ms]
dtR = 4f0	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(linspace(1950f0, 1950f0, nsrc))
ysrc = convertToCell(linspace(0f0, 0f0, nsrc))
zsrc = convertToCell(linspace(50f0, 1950f0, nsrc))

# source sampling and number of time steps
timeS = 3000f0
dtS = 4f0	# receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.004f0
wavelet = source(timeS,dtS,f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(model.n),nsrc,ntComp)

###################################################################################################

# Combined operator Pr*F*Ps'
F = joModeling(info,model,srcGeometry,recGeometry)
q = joData(srcGeometry,wavelet)

# Nonlinear modeling
dobs = F*q

# Optimization library
opt = "minConf"	# or "NLopt"

###################################################################################################

if opt=="minConf"

	# Optimization parameters
	iterations = 20
	batchsize = 10

	function objective_function(x)
		model0.m = reshape(x,model0.n);

		# select batch
		idx = randperm(dobs.nsrc)[1:batchsize]

		# fwi function value and gradient
		fval,grad = fwi_objective(model0,q[idx],dobs[idx])
		grad = reshape(grad,model0.n); grad[1:20,:] = 0.; grad[end-20:end,:] = 0.;

 	   return fval, vec(grad)
	end

	ProjBound(x) = vec(reshape(median([vec(mmin) vec(x) vec(mmax)],2), model.n))

	# FWI with SPG
	#options = spg_options(verbose=3, maxIter=iterations, memory=3)
	#x, fsave, funEvals= minConf_SPG(objective_function, vec(model0.m), ProjBound, options)

	# FWI with PQN
	options = pqn_options(verbose=3, maxIter=iterations, corrections=15)
	x, fsave, funEvals= minConf_PQN(objective_function, vec(model0.m), ProjBound, options)

####################################################################################################

elseif opt=="NLopt"

	# NLopt objective function
	batchsize = 10;	# number of sources per iteration
	count = 0
	function f!(x,grad)
		# calcualte gradient and function value
		model0.m = convert(Array{Float32,2},reshape(x,model0.n))
		idx = randperm(dobs.nsrc)[1:batchsize]
		fval,gradient = fwi_objective(model0,subsample(q,idx),subsample(dobs,idx))

		# reset water column to zero
		gradient = reshape(gradient,model0.n); gradient[1:20,:] = 0f0; gradient[180:end,:] = 0f0
		grad[1:end] = vec(gradient)

		global count; count += 1
		return convert(Float64,fval)
	end

	# Optimization parameters
	opt = Opt(:LD_LBFGS, info.n)
	lower_bounds!(opt, mmin)
	upper_bounds!(opt, mmax)
	maxeval!(opt,30)
	maxtime!(opt,100)
	min_objective!(opt,f!)
	(minf, minx, ret) = optimize(opt, vec(model0.m))

else
	throw("wrong optimization keyword")
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
begin
    begin
        g = Core.Box()
        begin
              #= /home/redbq/Desktop/stable-dev/julia-community/GeneralizedGenerated.jl/src/GeneralizedGenerated.jl:41 =#
            begin
                  #= /home/redbq/Desktop/stable-dev/julia-community/GeneralizedGenerated.jl/test/runtests.jl:168 =#
                g.contents = (x, r = 0)->begin
                    begin
                        begin
                            if (Main).:(===)(x, 0)
                                r
                            else
                                g.contents = g.contents
                                g.contents((Main).:-(x, 1), (Main).:+(r, x))
                            end
                        end
                    end
                end
                  #= /home/redbq/Desktop/stable-dev/julia-community/GeneralizedGenerated.jl/test/runtests.jl:172 =#
                println(g.contents(10))
            end
        end
    end
end


text:
module Blox

using ArrayViews

export view

include("blockvector.jl")

end # module