text:
#
# Minimize the acceleration of a composite Bézier curve on the Sphere S2 with approximation
#
# This example appeared in Sec. 5.2, second example, of
#
# > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
# > by minimizing the acceleration of a Bézier curve_.
# > Frontiers in Applied Mathematics and Statistics, 2018.
# > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
# > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
#
using Manopt, Manifolds, Colors, ColorSchemes
experiment_name = "Bezier_Approximation"
results_folder = joinpath(@__DIR__, "Minimize_Acceleration_Bezier")
!isdir(results_folder) && mkdir(results_folder)
asy_export = true #export data and results to asyExport
λ = 10.0

curve_samples = [range(0, 3; length=101)...] # sample curve for the gradient
curve_samples_plot = [range(0, 3; length=201)...] # sample curve for asy exports
camera_position = (-1.0, -0.7, 0.3)
curveColor = RGBA{Float64}(colorant"#000000")
sColor = RGBA{Float64}(colorant"#BBBBBB")
dColor = RGBA{Float64}(colorant"#EE7733") # data Color: Tol Vibrant Orange
pColor = RGBA{Float64}(colorant"#0077BB") # control point data color: Tol Virbant Blue
ξColor = RGBA{Float64}(colorant"#33BBEE") # tangent vector: Tol Vibrant blue
bColor = RGBA{Float64}(colorant"#009988") # inner control points: Tol Vibrant teal
#
# Data
#
M = Sphere(2)
B = artificial_S2_composite_bezier_curve()
cP = de_casteljau(M, B, curve_samples_plot)
dataP = get_bezier_junctions(M, B)
# export original data.
if asy_export
    asymptote_export_S2_signals(
        joinpath(results_folder, experiment_name * "-orig.asy");
        curves=[cP],
        points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)],
        tangent_vectors=[[
            Tuple(a) for a in
            zip(get_bezier_junctions(M, B, true), get_bezier_junction_tangent_vectors(M, B))
        ]],
        colors=Dict(
            :curves => [curveColor], :points => [dColor, bColor], :tvectors => [ξColor]
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 1.5],
        dot_size=4.0,
    )
    render_asymptote(joinpath(results_folder, experiment_name * "-orig.asy"); render=4)
end
pB = get_bezier_points(M, B, :differentiable)
N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
function F(pB)
    return cost_L2_acceleration_bezier(
        M, pB, get_bezier_degrees(M, B), curve_samples, λ, dataP
    )
end
∇F(pB) = ∇L2_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples, λ, dataP)
x0 = pB
pB_opt = gradient_descent(
    N,
    F,
    ∇F,
    x0;
    stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001), # use Armijo lineSearch
    stopping_criterion=StopWhenAny(
        StopWhenChangeLess(10.0^(-5)),
        StopWhenGradientNormLess(10.0^-7),
        StopAfterIteration(300),
    ),
    debug=[
        :Stop,
        :Iteration,
        " | ",
        :Cost,
        " | ",
        DebugGradientNorm(),
        " | ",
        DebugStepsize(),
        " | ",
        :Change,
        "\n",
        25,
    ],
)
B_opt = get_bezier_segments(M, pB_opt, get_bezier_degrees(M, B), :differentiable)
res_curve = de_casteljau(M, B_opt, curve_samples_plot)
#res_curve ./= norm.(res_curve)
if asy_export
    asymptote_export_S2_signals(
        joinpath(results_folder, experiment_name * "-result.asy");
        curves=[res_curve, cP],
        points=[get_bezier_junctions(M, B_opt), get_bezier_inner_points(M, B_opt)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B_opt, true),
                get_bezier_junction_tangent_vectors(M, B_opt),
            )
        ]],
        colors=Dict(
            :curves => [curveColor, pColor],
            :points => [dColor, bColor],
            :tvectors => [ξColor],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 0.75, 1.5],
        dot_size=4.0,
    )
    render_asymptote(joinpath(results_folder, experiment_name * "-result.asy"); render=4)
end


text:
function BinomPutDelta(K, T, S, sig, r, N)

  PutArray, AssetArray = (GenAmericanBinomialPut(K, T, S, sig, r, N))

  delta = (PutArray[1,2] - PutArray[2,2])/(AssetArray[1,2]-AssetArray[2,2])
  delta

end

ext: jl


text:
"""
Based on MATLAB's meshgrid

X,Y = meshgrid(x,y) returns 2-D grid coordinates based on the coordinates contained in vectors x and y.
X is a matrix where each row is a copy of x, and Y is a matrix where each column is a copy of y.
The grid represented by the coordinates X and Y has length(y) rows and length(x) columns.

X,Y,Z = meshgrid(x,y,z) returns 3-D grid coordinates based on the coordinates contained in vectors x, y and z.
"""

function meshgrid(x::AbstractVector{T}, y::AbstractVector{T}) where T
    m, n = length(y), length(x)
    x = reshape(x, 1, n)
    y = reshape(y, m, 1)
    (repeat(x, m, 1), repeat(y, 1, n))
end

function meshgrid(x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}) where T
    l, m, n = length(z), length(y), length(x)
    x = reshape(x, 1, 1, n)
    y = reshape(y, 1, m, 1)
    z = reshape(z, l, 1, 1)
    (repeat(x, l, m, 1), repeat(y, l, 1, n), repeat(z, 1, m, n))
end

max_stars_repo_name: timkimd/plnde

max_issues_repo_name: timkimd/plnde

max_forks_repo_name: timkimd/plnde


text:
const _VI = MOI.VariableIndex
const _SV = MOI.SingleVariable
const _VOV = MOI.VectorOfVariables
const _SAT = MOI.ScalarAffineTerm{Float64}
const _SAF = MOI.ScalarAffineFunction{Float64}
const _ET = MOI.EqualTo{Float64}
const _GT = MOI.GreaterThan{Float64}
const _LT = MOI.LessThan{Float64}
const _IN = MOI.Interval{Float64}
const _CI = MOI.ConstraintIndex
const _ZO = MOI.ZeroOne
const _GI = MOI.Integer

const _CVI = Cerberus.VariableIndex
const _CCI = Cerberus.ConstraintIndex
const _CSAF = Cerberus.ScalarAffineFunction

function _build_polyhedron()
    return Cerberus.Polyhedron(
        [
            Cerberus.AffineConstraint(
                _CSAF([1.0, 2.1, 3.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
                _ET(3.0),
            ),
            Cerberus.AffineConstraint(
                _CSAF([-3.5, 1.2], [_CVI(1), _CVI(2)], 0.0),
                _LT(4.0),
            ),
        ],
        [_IN(0.5, 1.0), _IN(-1.3, 2.3), _IN(0.0, 1.0)],
    )
end

function _build_dmip_formulation()
    return Cerberus.DMIPFormulation(
        _build_polyhedron(),
        [_ZO(), nothing, _ZO()],
        _CSAF([1.0, -1.0], [_CVI(1), _CVI(2)], 0.0),
    )
end

function _build_gi_polyhedron()
    return Cerberus.Polyhedron(
        [
            Cerberus.AffineConstraint(
                _CSAF([1.3, 3.7, 2.4], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
                _LT(5.5),
            ),
        ],
        [_IN(0.0, 4.5), _IN(0.0, 1.0), _IN(0.0, 3.0)],
    )
end

function _build_gi_dmip_formulation()
    return Cerberus.DMIPFormulation(
        _build_gi_polyhedron(),
        [nothing, _ZO(), _GI()],
        _CSAF(),
    )
end

# Indices correspond to what Gurobi.jl, not Cerberus, uses
const DMIP_BASIS = Cerberus.Basis(
    [MOI.NONBASIC_AT_LOWER, MOI.BASIC, MOI.NONBASIC_AT_LOWER],
    [MOI.BASIC],
    MOI.BasisStatusCode[],
    [MOI.NONBASIC],
    MOI.BasisStatusCode[],
    MOI.BasisStatusCode[],
)

function _test_is_equal_to_dmip_basis(basis::Cerberus.Basis)
    @test basis.base_var_constrs == DMIP_BASIS.base_var_constrs
    @test basis.base_lt_constrs == DMIP_BASIS.base_lt_constrs
    @test basis.base_gt_constrs == DMIP_BASIS.base_gt_constrs
    @test basis.base_et_constrs == DMIP_BASIS.base_et_constrs
    @test basis.branch_lt_constrs == DMIP_BASIS.branch_lt_constrs
    @test basis.branch_gt_constrs == DMIP_BASIS.branch_gt_constrs
end

function _is_equal(u::_CSAF, v::_CSAF)
    u_terms = [Cerberus.index(u.indices[i]) => u.coeffs[i] for
     i in 1:length(u.indices)]
    v_terms = [Cerberus.index(v.indices[i]) => v.coeffs[i] for
     i in 1:length(v.indices)]
    u_t = sort(u_terms, lt = (x, y) -> x[1] < y[1])
    v_t = sort(v_terms, lt = (x, y) -> x[1] < y[1])
    if !(u.constant ≈ v.constant)
        return false
    end
    u_idx = [t[1] for t in u_t]
    v_idx = [t[1] for t in v_t]
    if u_idx != v_idx
        return false
    end
    u_coeffs = [t[2] for t in u_t]
    v_coeffs = [t[2] for t in v_t]
    @assert length(u_coeffs) == length(v_coeffs)
    for i in 1:length(u_coeffs)
        if !(u_coeffs[i] ≈ v_coeffs[i])
            return false
        end
    end
    return true
end

function _is_equal(u::_SAF, v::_SAF)
    return _is_equal(convert(Cerberus.CSAF, u), convert(Cerberus.CSAF, v))
end

function _build_formulation_with_single_disjunction(
    method::DisjunctiveConstraints.AbstractDisjunctiveFormulation = DisjunctiveConstraints.NaiveBigM(
        DisjunctiveConstraints.IntervalArithmetic(),
    ),
)
    # min  y + 1.2
    # s.t. (-2 <= x <= -1 & y = -x - 1) or
    #          (-1 <= x <= +1 & y = 0) or
    #          (+1 <= x <= +2 & y = x - 1)
    #      x + y >= -0.5
    #      x in [-1.5, 3.0]
    #      y in [0.0, 0.5]
    form = Cerberus.DMIPFormulation(
        Cerberus.Polyhedron(
            [
                Cerberus.AffineConstraint(
                    _CSAF([1.0, 1.0], [_CVI(1), _CVI(2)], 0.0),
                    _GT(-0.5),
                ),
            ],
            [_IN(-1.5, 3.0), _IN(0.0, 0.5)],
        ),
        [nothing, nothing],
        _CSAF([1.0], [_CVI(2)], 1.2),
    )
    disjunction = Cerberus.Disjunction(
        [
            _CSAF([1.0], [_CVI(1)], 0.0),
            _CSAF([1.0, 1.0], [_CVI(1), _CVI(2)], 0.0),
            _CSAF([1.0], [_CVI(2)], 0.0),
            _CSAF([-1.0, 1.0], [_CVI(1), _CVI(2)], 0.0),
        ],
        DisjunctiveConstraints.DisjunctiveSet(
            [
                -2.0 -1.0 +1.0
                -1.0 -Inf -Inf
                -Inf 0.0 -Inf
                -Inf -Inf -1.0
            ],
            [
                -1.0 +1.0 +2.0
                -1.0 Inf Inf
                Inf 0.0 Inf
                Inf Inf -1.0
            ],
        ),
    )
    formulater = Cerberus.DisjunctiveFormulater(disjunction, method)
    Cerberus.attach_formulater!(form, formulater)
    return form
end

function _test_roundtrip_model(
    model::MOI.ModelLike,
    expected_bounds::Vector{Tuple{Float64,Float64}},
    expected_lt_acs::Vector{Tuple{_SAF,_LT}},
    expected_gt_acs::Vector{Tuple{_SAF,_GT}},
    expected_et_acs::Vector{Tuple{_SAF,_ET}},
)
    n = length(expected_bounds)
    @test MOI.get(model, MOI.NumberOfVariables()) == n
    for i in 1:n
        @test MOIU.get_bounds(model, Float64, _VI(i)) == expected_bounds[i]
    end

    expected_loc = Set([])
    if !isempty(expected_bounds)
        push!(expected_loc, (_SV, _IN))
    end
    if !isempty(expected_lt_acs)
        push!(expected_loc, (_SAF, _LT))
    end
    if !isempty(expected_gt_acs)
        push!(expected_loc, (_SAF, _GT))
    end
    if !isempty(expected_et_acs)
        push!(expected_loc, (_SAF, _ET))
    end
    @test Set(MOI.get(model, MOI.ListOfConstraints())) == expected_loc

    lt_constr_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_LT}())
    @test length(lt_constr_cis) == length(expected_lt_acs)
    for (i, ci) in enumerate(lt_constr_cis)
        f_actual = MOI.get(model, MOI.ConstraintFunction(), ci)
        s_actual = MOI.get(model, MOI.ConstraintSet(), ci)
        f_expected, s_expected = expected_lt_acs[i]
        @test _is_equal(f_actual, f_expected)
        @test s_actual == s_expected
    end

    gt_constr_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_GT}())
    @test length(gt_constr_cis) == length(expected_gt_acs)
    for (i, ci) in enumerate(gt_constr_cis)
        f_actual = MOI.get(model, MOI.ConstraintFunction(), ci)
        s_actual = MOI.get(model, MOI.ConstraintSet(), ci)
        f_expected, s_expected = expected_gt_acs[i]
        @test _is_equal(f_actual, f_expected)
        @test s_actual == s_expected
    end

    et_constr_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_ET}())
    @test length(et_constr_cis) == length(expected_et_acs)
    for (i, ci) in
        enumerate(MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_ET}()))
        f_actual = MOI.get(model, MOI.ConstraintFunction(), ci)
        s_actual = MOI.get(model, MOI.ConstraintSet(), ci)
        f_expected, s_expected = expected_et_acs[i]
        @test _is_equal(f_actual, f_expected)
        @test s_actual == s_expected
    end
end

function _CurrentState(; primal_bound = Inf)
    state = Cerberus.CurrentState(primal_bound = primal_bound)
    state.gurobi_env = GRB_ENV
    return state
end


text:
using Statistics

"""
`CALCULATE_AFFINE` - Compute left-hand affine transform from two Nxd point sets

```
tform = calculate_affine(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being homogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being homogeneous coords of a point, of the
    reference point set.
* tform: affine transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_affine(moving_pts, fixed_pts)
    return moving_pts \ fixed_pts
end

"""
`CALCULATE_RIGID` - Compute right-hand rigid transform from two Nxd point sets

```
tform = calculate_rigid(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being homogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being homogeneous coords of a point, of the
    reference point set.
* tform: rigid transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_rigid(moving_pts, fixed_pts)
    moving_pts = moving_pts[:, 1:end-1]
    fixed_pts = fixed_pts[:, 1:end-1]
    n, dim = size(moving_pts)
    moving_bar = mean(moving_pts, dims = 1)
    fixed_bar = mean(fixed_pts, dims = 1)
    moving_centered = moving_pts .- moving_bar
    fixed_centered = fixed_pts .- fixed_bar
    C = transpose(fixed_centered) * moving_centered / n
    U, s, V = svd(C)
  # Rotation R in least squares sense:
  # moving_pts - moving_bar = (fixed_pts - fixed_bar)*R
    R = transpose(U * diagm(0 => vcat(ones(dim - 1), det(U * transpose(V)))) * transpose(V))
    t = fixed_bar - moving_bar * R
    return [R zeros(dim); t 1]
end

"""
`CALCULATE_TRANSLATION` - Compute left-hand translation from two Nxd point sets

```
tform = calculate_translation(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being nonhomogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being nonhomogeneous coords of a point, of the
    reference point set.
* tform: translation transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_translation(moving_pts, fixed_pts)
    n, dim = size(moving_pts)
    moving_bar = mean(moving_pts, dims = 1)
    fixed_bar = mean(fixed_pts, dims = 1)
    t = fixed_bar - moving_bar
    return [Matrix(1.0I, 2, 2) zeros(dim); t 1]
end

function calculate_affine(matches::Matches)
    return calculate_affine(matches.src_points, matches.dst_points)
end

function calculate_rigid(matches::Matches)
    return calculate_rigid(matches.src_points, matches.dst_points)
end

function calculate_translation(matches::Matches)
    return calculate_translation(matches.src_points, matches.dst_points)
end

function calculate_affine(mesh::Mesh)
    return calculate_affine(mesh.src_nodes, mesh.dst_nodes)
end

function calculate_rigid(mesh::Mesh)
    return calculate_rigid(mesh.src_nodes, mesh.dst_nodes)
end

function calculate_translation(mesh::Mesh)
    return calculate_translation(mesh.src_nodes, mesh.dst_nodes)
end

"""
Matrix for rotation in degrees (clockwise rotation with left-handed Cartesian)
"""
function make_rotation_matrix(theta, image_size = nothing)
    angle = deg2rad(theta)
    tform = [cos(angle) sin(angle) 0; -sin(angle) cos(angle) 0; 0 0 1]
    tform[abs(tform).<eps] = 0
    if image_size != nothing
        rotated_bb = snap_bb(tform_bb(sz_to_bb([image_size[1], image_size[2]]), tform))
        rotation_offset = [rotated_bb.i, rotated_bb.j]
        translation_matrix = make_translation_matrix(-rotation_offset)
        tform = tform * translation_matrix
    end
    return tform
end

function make_translation_matrix(offset)
    return [1 0 0; 0 1 0; offset[1] offset[2] 1]
end

function make_scale_matrix(scale)
    return [scale 0 0; 0 scale 0; 0 0 1]
end


text:
using PyPlot, PyCall
using NetCDF
pyimport_conda("cartopy.crs", "cartopy", "conda-forge")
ssl = pyimport_conda("ssl","ssl")
ssl._create_default_https_context = ssl._create_unverified_context

"""
    plotmap(fname, var; ...)

Easy plotting of gridded datasets on a global map.
Plots variable `var` from a netcdf file `fname`.
Optional arguments are available to control the details of the plot.

# Arguments
* data::Array{Float32,2} : data to plot. If a 3D array is passed, only the first frame is plotted: `data[:,:,1]`
* lon::Array{Float64,1} : longitudes
* lat::Array{Float64,1} : latitudes

Optional:
* title :       title string ("")
* cstep :       divisions of the colorbar axis ([])
* cmap :        colormap ("RdBu_r")
* proj :        projection. One of ["platecarree", "robinson", "mollweide"]. Defaults to "platecarree".
* cpad :        padding (shift) of the colorbar (0.08)
* sub :         matplotlib subplot option (e.g. "221" for the first panel of 2x2 subplots) ("111")
* clabel :      label of the colorbar (defaults to the units string read from the netcdf file)
* cdir :        direction of the colorbar. One of ["horizontal", "vertical"]. ("horizontal")
* cscale :      scaling of the colorbar (0.65)
* cfs :         colorbar ticks font size (12)
* lfs :         colorbar label font size (12)
* tfs :         title font size (14)
* borders :     national borders (false)
* tpad :        padding (shift) of the title string
* tweight :     weight of title font. One of [ 'normal" "bold" "heavy" "light" "ultrabold" "ultralight"]. ("normal")
* grid :        grid spacing (defaults to [60,30]). set to empty [] to remove gridlines.
* region::NTuple{4,Int64}:       region to plot in format (lon1, lon2, lat1, lat2). Defaults to global.
* style :       one of ["pcolormesh" "contourf"]. Defaults to "pcolormesh".
* levels :      contour plot levels. Can be an array or a number of levels (auto)
* extend :      plot levels outside `levels` range. One of ["neither", "both", "min", "max"]. Default: "neither".

"""
function plotmap(fname::String, var::String; lon="lon", lat="lat", lonb="lon_bnds", latb="lat_bnds", titles="", title="", cstep=[], cmap="RdBu_r", proj="", cpad=0.08, tpad=24, sub=111, clabel="NONE", cdir="horizontal", cscale=0.65, tfs=14, cfs=12, lfs=12, tweight="normal", grid=[60,30], region=(), style="pcolormesh", levels=0, extend="neither", borders=false)

# pcolormesh needs cell boundaries
if style=="pcolormesh"
    try
        lonb=ncread(fname, lonb);
        lonv=vcat(lonb[1,:],lonb[2,end])
    catch
        lonv=ncread(fname, lon);
    end
    try
        latb=ncread(fname, latb);
        latv=vcat(latb[1,:],latb[2,end])
    catch
        latv=ncread(fname, lat);
    end
else
    lonv=ncread(fname, lon);
    latv=ncread(fname, lat);
end

data=ncread(fname, var);
units=ncgetatt(fname, var, "units")
if clabel=="NONE" clabel=units end

plotmap(lonv, latv, data; title=title, titles=titles, cstep=cstep, cmap=cmap, proj=proj, cpad=cpad, tpad=tpad, sub=sub, clabel=clabel, cdir=cdir, cscale=cscale, tfs=tfs, cfs=cfs, lfs=lfs, tweight=tweight, grid=grid, region=region, style=style, levels=levels, extend=extend, borders=borders)

end

"""
    plotmap(lon, lat, data; ...)

Easy plotting of gridded datasets on a global map.
Plots data in 2D array `data` with longitudes `lon` and latitudes `lat`.
Optional arguments are available to control the details of the plot.

# Arguments
* data::Array{Float32,2} : data to plot. If a 3D array is passed, only the first frame is plotted: `data[:,:,1]`
* lon::Array{Float64,1} : longitudes
* lat::Array{Float64,1} : latitudes

Optional:
* title :       title string ("")
* cstep :       divisions of the colorbar axis ([])
* cmap :        colormap ("RdBu_r")
* proj :        projection. One of ["platecarree", "robinson", "mollweide"]. Defaults to "platecarree".
* cpad :        padding (shift) of the colorbar (0.08)
* sub :         matplotlib subplot option (e.g. "221" for the first panel of 2x2 subplots) ("111")
* clabel :      label of the colorbar (defaults to the units string read from the netcdf file)
* cdir :        direction of the colorbar. One of ["horizontal", "vertical"]. ("horizontal")
* cscale :      scaling of the colorbar (0.65)
* cfs :         colorbar ticks font size (12)
* lfs :         colorbar label font size (12)
* tfs :         title font size (14)
* borders :     national borders (false)
* tpad :        padding (shift) of the title string
* tweight :     weight of title font. One of [ 'normal" "bold" "heavy" "light" "ultrabold" "ultralight"]. ("normal")
* grid :        grid spacing (defaults to [60,30]). set to empty [] to remove gridlines.
* region::NTuple{4,Int64}:       region to plot in format (lon1, lon2, lat1, lat2). Defaults to global.
* style :       one of ["pcolormesh" "contourf"]. Defaults to "pcolormesh".
* levels :      contour plot levels. Can be an array or a number of levels (auto)
* extend :      plot levels outside `levels` range. One of ["neither", "both", "min", "max"]. Default: "neither".

"""
function plotmap(lon, lat, data; title="", titles="", cstep=[], cmap="RdBu_r", proj="", cpad=0.08, tpad=24, sub=111, clabel="", cdir="horizontal", cscale=0.65, tfs=14, cfs=12, lfs=12, tweight="normal", grid=[60,30], region=(), style="pcolormesh", levels=0, extend="neither", borders=false)

#cf = pyimport("cartopy.feature")

dd = size(data)

if length(dd)==3 data=data[:,:,1] end
if style=="pcolormesh"
    if length(lon) in dd
        #println("pcolormesh needs cell boundaries, reconstructing lon")
        lonb=zeros(2,length(lon))
        lonb[1,2:end]=0.5*(lon[2:end]+lon[1:(end-1)])
        lonb[1,1]=lon[1]-(lon[2]-lon[1])*0.5
        lonb[2,end]=lon[end]+(lon[end]-lon[end-1])*0.5
        lon=vcat(lonb[1,:],lonb[2,end])
    end
    if length(lat) in dd
        #println("pcolormesh needs cell boundaries, reconstructing lat")
        latb=zeros(2,length(lat))
        latb[1,2:end]=0.5*(lat[2:end]+lat[1:(end-1)])
        latb[1,1]=lat[1]-(lat[2]-lat[1])*0.5
        latb[2,end]=lat[end]+(lat[end]-lat[end-1])*0.5
        if latb[1,1]>89; latb[1,1]=90 ; end
        if latb[1,1]<-89; latb[1,1]=-90 ; end
        if latb[2,end]>89; latb[2,end]=90 ; end
        if latb[2,end]<-89; latb[2,end]=-90 ; end
        lat=vcat(latb[1,:],latb[2,end])
    end
    if length(lon)==(dd[1]+1) data=data' end
else
    if length(lon)==dd[1] data=data' end
end

ccrs = pyimport("cartopy.crs")
cutil = pyimport("cartopy.util")
cfeat = pyimport("cartopy.feature")

if proj=="robinson"
    proj=ccrs.Robinson()
    dlabels=false
elseif proj == "mollweide"
    proj=ccrs.Mollweide()
    dlabels=false
elseif proj == "polar_north"
    proj=ccrs.NorthPolarStereo()
    region = (-180, 180, 45, 90)
    dlabels=false
elseif proj == "polar_south"
    proj=ccrs.SouthPolarStereo()
    region = (-180, 180, -90, -45)
    dlabels=false
else
    proj=ccrs.PlateCarree()
    dlabels=true
end

ax = subplot(sub, projection=proj)
if length(region)>0 ax.set_extent(region, crs=ccrs.PlateCarree()) end
ax.coastlines()
if borders ax.add_feature(cfeat.BORDERS) end
xlocvec=vcat(-vcat(grid[1]:grid[1]:180)[end:-1:1], vcat(0:grid[1]:180))
ylocvec=vcat(-vcat(grid[2]:grid[2]:90)[end:-1:1], vcat(0:grid[2]:90))

if dlabels
    ax.gridlines(linewidth=1, color="gray", alpha=0.5, linestyle="--", draw_labels=true, xlocs=xlocvec, ylocs=ylocvec)
else
    ax.gridlines(linewidth=1, color="gray", alpha=0.5, linestyle="--", xlocs=xlocvec, ylocs=ylocvec)
end

if style=="contourf"
    data_cyc, lon_cyc = cutil.add_cyclic_point(data, coord=lon)
    if levels==0
        contourf(lon_cyc, lat, data_cyc, transform=ccrs.PlateCarree(), cmap=cmap, extend=extend)
    else
        contourf(lon_cyc, lat, data_cyc, transform=ccrs.PlateCarree(), cmap=cmap, levels=levels, extend=extend)
    end
else
    pcolormesh(lon, lat, data, transform=ccrs.PlateCarree(), cmap=cmap)
end

if length(cstep)>0 clim(cstep[1],cstep[end]); end
if length(title)>1 PyPlot.title(title, pad=tpad, fontsize=tfs, weight=tweight) end
if length(titles)>1 PyPlot.title(titles, pad=tpad, fontsize=tfs, weight=tweight) end
cbar=colorbar(orientation=cdir, extend="both", pad=cpad, label=clabel, shrink=cscale)
cbar.set_label(label=clabel,size=lfs)
cbar.ax.tick_params(labelsize=cfs)
if length(cstep)>0 cbar.set_ticks(cstep) end
tight_layout()

end


text:
#
# rules-tests.jl -
#
# Test algebraic rules and simplifications.
#
module TestingLazyAlgebraRules

using LazyAlgebra
using LazyAlgebra.Foundations
import LazyAlgebra: ⋅,
    Adjoint, Inverse, InverseAdjoint, Jacobian,
    Scaled, Sum, Composition
using Test

function to_string(A::Mapping)
    io = IOBuffer()
    show(io, A)
    String(take!(io))
end

@testset "Algebraic rules" begin
    include("common.jl")
    identical = LazyAlgebra.identical
    Composition = LazyAlgebra.Composition
    dims = (3,4,5)

    M = SymbolicMapping(:M)
    Q = SymbolicMapping(:Q)
    R = SymbolicMapping(:R)
    A = SymbolicLinearMapping(:A)
    B = SymbolicLinearMapping(:B)
    C = SymbolicLinearMapping(:C)
    D = SymbolicLinearMapping(:D)

    # Test properties.
    @test M !== R
    @test A !== B
    @test identical(A, B) == false
    let E = A
        @test identical(A, E) == true
    end
    @test is_linear(M) == false
    @test is_linear(A) == true
    @test is_linear(A + B) == true
    @test is_linear(A') == true
    @test is_linear(inv(A)) == true
    @test is_linear(inv(M)) == false
    @test is_linear(A + B + C) == true
    @test is_linear(A' + B + C) == true
    @test is_linear(A + B' + C) == true
    @test is_linear(A + B + C') == true
    @test is_linear(M + B + C) == false
    @test is_linear(A + M + C) == false
    @test is_linear(A + B + M) == false
    @test is_linear(A*B*C) == true
    @test is_linear(A'*B*C) == true
    @test is_linear(A*B'*C) == true
    @test is_linear(A*B*C') == true
    @test is_linear(M*B*C) == false
    @test is_linear(A*M*C) == false
    @test is_linear(A*B*M) == false
    @test is_endomorphism(M) == false
    @test is_endomorphism(A) == false
    @test is_selfadjoint(M) == false
    @test is_selfadjoint(A) == false
    @test is_selfadjoint(A'*A) == true # FIXME: was broken
    @test is_selfadjoint(A*A') == true # FIXME: was broken
    @test_broken is_selfadjoint(B'*A*A'*B) == true

    # Test identity.
    @test Id === LazyAlgebra.Id
    @test Id' === Id
    @test inv(Id) === Id
    @test 1Id === Id
    @test Id*M*Id === M
    @test Id*Id === Id
    @test Id*Id*Id === Id
    @test Id\Id === Id
    @test Id/Id === Id
    @test Id + Id === 2Id
    @test Id - Id === 0Id
    @test Id + 2Id === 3Id
    @test 2Id + 3Id === 5Id
    @test -4Id + (Id + 2Id) === -Id
    @test Id + Id - 2Id === 0Id
    @test 2Id - (Id + Id) === 0Id
    @test inv(3Id) === (1/3)*Id
    @test SelfAdjointType(Id) === SelfAdjoint()
    @test MorphismType(Id) === Endomorphism()
    @test DiagonalType(Id) === DiagonalMapping()
    for T in (Float32, Float64)
        atol, rtol = zero(T), sqrt(eps(T))
        x = randn(T, dims)
        y = randn(T, dims)
        @test Id*x === x # test same object
        for P in (Direct, Adjoint, Inverse, InverseAdjoint)
            @test apply(P,Id,x) === x # test same object
            test_api(P, Id, x, y)
        end
    end

    # Neutral elements.
    @test isone(Id)
    @test iszero(Id) == false
    @test iszero(A - A)
    @test iszero(-M + M)
    @test one(A) === Id
    @test one(M) === Id
    @test zero(A) === 0*A
    @test zero(M) === 0*M

    # Basic methods for sums and compositions of mappings.
    let A1 = A + B + M, A2 = C*B*M*Q
        @test eltype(A1) <: Mapping
        @test eltype(A2) <: Mapping
        @test Tuple(A1) === terms(A1)
        @test Tuple(A2) === terms(A2)
        @test length(A1) == 3
        @test length(A2) == 4
        @test firstindex(A1) == 1
        @test firstindex(A2) == 1
        @test lastindex(A1) == length(A1)
        @test lastindex(A2) == length(A2)
        @test first(A1) === A1[1]
        @test first(A2) === A2[1]
        @test last(A1) === A1[end]
        @test last(A2) === A2[end]
    end

    @test B + A === A + B
    @test B + A + Q === Q + A + B
    @test B + (M + A + Q) === A + B + M + Q
    @test (-1)*A === -A
    @test -(-A) === A
    @test 2\A === (1/2)*A
    @test A + A === 2A
    @test A + B - A === B
    @test 2R + A - Q + B - 2A + Q + 1R === 3R + B - A
    # FIXME: for the following to work, we must impose the value of the
    #        first multiplier of a sum
    #@test A + B*(M - Q) + A - 2B*(Q - M) === 2*A + 3*B*(M - Q)
    @test A + B*(M - Q) + 3A - 3B*(M - Q) === 4*A - 2*B*(M - Q)
    @test A*2M === 2*(A*M)
    @test A*2M === 2*(A*M)
    @test 3A*2M === 6A*M
    @test 3R*2M !== 6R*M

    # Test adjoint and Jacobian.
    x = nothing
    @test (A*A')' === A*A'
    for X in (A*A', A'*A, B'*A'*A*B, B'*A*A'*B)
        @test X' === X
    end
    @test A' === Adjoint(A)
    @test A' isa Adjoint
    @test A'' === (A')' === A
    @test adjoint(A) === A'
    @test jacobian(A,x) === A
    @test ∇(A,x) === A
    @test (3A)' === 3*(A')
    @test (A + 2B)' - A' === 2*B'
    @test_throws ArgumentError Jacobian(A,x)
    @test_throws ArgumentError M'
    @test_throws ArgumentError adjoint(M)
    @test_throws ArgumentError Adjoint(M)
    @test jacobian(M,x) isa Jacobian
    @test ∇(M,x) === jacobian(M,x)
    @test ∇(3M,x) === 3*∇(M,x)
    @test ∇(M,x) + ∇(2M,x) === 3∇(M,x)

    # Inverse.
    @test inv(M) === Id/M
    @test inv(inv(M)) === M
    @test Id/A === inv(A)
    @test Id\A === A
    @test M/Q === M*inv(Q)
    @test M\Q === inv(M)*Q
    @test inv(2M) === inv(M)*(2\Id)
    @test inv(2M) === inv(M)*((1/2)*Id)
    @test inv(2A) === 2\inv(A)
    @test inv(2A) === (1/2)*inv(A)
    @test inv(A*B) === inv(B)*inv(A)
    @test A/3B === 3\A/B
    @test 2/A === 2*inv(A)
    @test 4A/4B === A/B
    @test 4A\4B === A\B
    @test inv(A*M*B*Q) === inv(Q)*inv(B)*inv(M)*inv(A)
    @test inv(M)*M === Id
    @test M*inv(M) === Id
    let D = M*Q*(A - B)
        @test inv(D)*D === Id
        @test D*inv(D) === Id
    end
    let D = A + 2B - C
        @test inv(D)*D === Id
        @test D*inv(D) === Id
    end
    @test inv(M*Q*(A - B)) === inv(A - B)*inv(Q)*inv(M)
    @test inv(A*3M) === inv(M)*(3\inv(A))
    @test inv(A*3B) === 3\inv(A*B) === 3\inv(B)*inv(A)
    @test inv(A*B*3M) === inv(M)*(3\inv(B))*inv(A)
    @test inv(A*B*3C) === 3\inv(A*B*C) === 3\inv(C)*inv(B)*inv(A)

    # Inverse-adjoint.
    @test inv(A)' === inv(A')
    @test inv(A')*A' === Id
    @test A'*inv(A') === Id
    let E = inv(A'), F = inv(A)'
        @test inv(E) === A'
        @test inv(F) === A'
    end

    # Test aliases for composition.
    @test isa(M*R, Composition)
    @test M⋅R === M*R
    @test M∘R === M*R

    # Test associativity of sum and composition.
    @test (M + R) + Q === M + R + Q
    @test M + (R + Q) === (M + R) + Q
    @test (M*R)*Q === M*R*Q
    @test M*(R*Q) === (M*R)*Q

    # Test adjoint of sums and compositions.
    @test (A*B)' === (B')*(A') === B'*A'
    @test (A'*B)' === B'*A
    @test (A*B')' === B*A'
    @test (A*B*C)' === C'*B'*A'
    @test (A'*B*C)' === C'*B'*A
    @test (A*B'*C)' === C'*B*A'
    @test (A*B*C')' === C*B'*A'
    @test (A + B)' === A' + B'
    @test (A' + B)' === A + B'
    @test (A + B')' === A' + B
    @test (A' + B + C)' === A + B' + C'
    @test (A + B' + C)' === A' + B + C'
    @test (A + B + C')' === A' + B' + C

    # Test inverse of sums and compositions.

    # Test unary plus and negation.
    @test +M === M
    @test -(-M) === M
    @test -M === (-1)*M
    @test 2A + A === 3A
    @test 10M - 3M === 7M
    @test M + A + M === A + 2M

    # Test forbidden calls to constructors because they should yield an
    # instance of a different type if simplification rules were applied.
    @test_throws ArgumentError Adjoint(Id)
    @test_throws ArgumentError Adjoint(A')
    @test_throws ArgumentError Adjoint(inv(A))
    @test_throws ArgumentError Adjoint(inv(A'))
    @test_throws ArgumentError Adjoint(3A)
    @test_throws ArgumentError Adjoint(A + B)
    @test_throws ArgumentError Adjoint(A*B)

    @test_throws ArgumentError Inverse(Id)
    @test_throws ArgumentError Inverse(A')
    @test_throws ArgumentError Inverse(inv(A))
    @test_throws ArgumentError Inverse(inv(A'))
    @test_throws ArgumentError Inverse(3A)
    @test_throws ArgumentError Inverse(A*B)

    @test_throws ArgumentError InverseAdjoint(Id)
    @test_throws ArgumentError InverseAdjoint(A')
    @test_throws ArgumentError InverseAdjoint(inv(A))
    @test_throws ArgumentError InverseAdjoint(inv(A'))
    @test_throws ArgumentError InverseAdjoint(3A)
    @test_throws ArgumentError InverseAdjoint(A*B)

    @test_throws ArgumentError Jacobian(3M,x)

    @test_throws ArgumentError Scaled(2,3M)

    # Check that sums and compositions must have at least 2 terms
    @test_throws ArgumentError Sum()
    @test_throws ArgumentError Sum(A)
    @test Sum(A,B) isa Sum
    @test_throws ArgumentError Composition()
    @test_throws ArgumentError Composition(A)
    @test Composition(A,B) isa Composition

    # Test the `show` method.
    @test to_string(A) == "A"
    @test to_string(A') == "A'"
    @test to_string(A + A) == "2⋅A"
    @test to_string(A' + A') == "2⋅A'"
    @test to_string(Id/A) == "inv(A)"
    @test to_string(Id/(A + B)) == "inv($(to_string(A + B)))"
    @test to_string(M) == "M"
    @test to_string(M + M) == "2⋅M"
    @test to_string(∇(M,x)) == "∇(M,x)"
    @test to_string(∇(M,x) + ∇(M,x)) == "2⋅∇(M,x)"
end # testset
nothing

end # module


text:
mult(x, y) = x * y
typeof(mult) #=> Function

plustwo = function (x)
              x + 2
          end
plustwo(3) #=> 5

function quadratic(f)
       # calculations (for full code see end of chapter)
    return (-b + sqrt(b^2 -4a*c + 0im)) / 2a, (-b - sqrt(b^2 -4a*c + 0im)) / 2a
end
r1, r2 = quadratic(x -> 2x^2 + 30x + 9)

function derivative(f)
    return function(x)
        # definiton of anonymous function calculating derivative
    end
end

# nested:
function a(x)
    z = x * 2
    function b(z)
        z += 1
    end
    b(z)
end

d = 5
a(d) #=> 11

# closure:
function counter()
    n = 0
    () -> n +=1, () -> n = 0
end
(addOne, reset) = counter()
addOne()
addOne()
addOne()
reset()

# recursive:
sum(n) =  n > 1 ? sum(n-1) + n : n
sum(100) #=> 5050
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
fib(1000) #=>

max_stars_repo_path: Module 1/Chapter03/Old/first_class.jl

max_issues_repo_path: Module 1/Chapter03/Old/first_class.jl

max_forks_repo_path: Module 1/Chapter03/Old/first_class.jl


text:
Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y
Base.:/(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x / y

struct ValidationError <: Exception
    message::String
end

"Throw exception on invalid unit types, otherwise return argument."
function screen_unit(result)
    result isa Unitful.Unitlike || throw(ValidationError("Unit must be a subtype of Unitful.Unitlike, not $(typeof(result))."))
    result isa Unitful.ScalarUnits || throw(ValidationError("Non-scalar units such as $result are not supported. Use a scalar unit instead."))
    result == u"°" && throw(ValidationError("Degrees are not supported. Use radians instead."))
    result
end

"""Test unit equivalence.

Example of implemented behavior:
```julia
using ModelingToolkit, Unitful
MT = ModelingToolkit
@parameters γ P [unit = u"MW"] E [unit = u"kJ"] τ [unit = u"ms"]
@test MT.equivalent(u"MW" ,u"kJ/ms") # Understands prefixes
@test !MT.equivalent(u"m", u"cm") # Units must be same magnitude
@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E/τ)^γ)) # Handles symbolic exponents
```
"""
equivalent(x,y) = isequal(1*x,1*y)
unitless = Unitful.unit(1)

#For dispatching get_unit
Literal = Union{Sym,Symbolics.ArrayOp,Symbolics.Arr,Symbolics.CallWithMetadata}
Conditional = Union{typeof(ifelse),typeof(IfElse.ifelse)}
Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}

"Find the unit of a symbolic item."
get_unit(x::Real) = unitless
get_unit(x::Unitful.Quantity) = screen_unit(Unitful.unit(x))
get_unit(x::AbstractArray) = map(get_unit,x)
get_unit(x::Num) = get_unit(value(x))
get_unit(x::Literal) = screen_unit(getmetadata(x,VariableUnit, unitless))
get_unit(op::Differential, args) = get_unit(args[1]) / get_unit(op.x)
get_unit(op::Difference, args) =   get_unit(args[1]) / get_unit(op.t)
get_unit(op::typeof(getindex),args) = get_unit(args[1])
get_unit(x::SciMLBase.NullParameters) = unitless

function get_unit(op,args) # Fallback
    result = op(1 .* get_unit.(args)...)
    try
        unit(result)
    catch
        throw(ValidationError("Unable to get unit for operation $op with arguments $args."))
    end
end

function get_unit(op::Integral,args)
    unit = 1
    if op.domain.variables isa Vector
        for u in op.domain.variables
            unit *= get_unit(u)
        end
    else
        unit *= get_unit(op.domain.variables)
    end
    return get_unit(args[1]) * unit
end

function get_unit(x::Pow)
    pargs = arguments(x)
    base,expon = get_unit.(pargs)
    @assert expon isa Unitful.DimensionlessUnits
    if base == unitless
        unitless
    else
        pargs[2] isa Number ? base^pargs[2] : (1*base)^pargs[2]
    end
end

function get_unit(x::Add)
    terms = get_unit.(arguments(x))
    firstunit = terms[1]
    for other in terms[2:end]
        termlist = join(map(repr, terms), ", ")
        equivalent(other, firstunit) || throw(ValidationError(", in sum $x, units [$termlist] do not match."))
    end
    return firstunit
end

function get_unit(op::Conditional, args)
    terms = get_unit.(args)
    terms[1] == unitless || throw(ValidationError(", in $x, [$(terms[1])] is not dimensionless."))
    equivalent(terms[2], terms[3]) || throw(ValidationError(", in $x, units [$(terms[2])] and [$(terms[3])] do not match."))
    return terms[2]
end

function get_unit(op::typeof(Symbolics._mapreduce),args)
    if args[2] == +
        get_unit(args[3])
    else
        throw(ValidationError("Unsupported array operation $op"))
    end
end

function get_unit(op::Comparison, args)
    terms = get_unit.(args)
    equivalent(terms[1], terms[2]) || throw(ValidationError(", in comparison $op, units [$(terms[1])] and [$(terms[2])] do not match."))
    return unitless
end

function get_unit(x::Symbolic)
    if SymbolicUtils.istree(x)
        op = operation(x)
        if op isa Sym || (op isa Term && operation(op) isa Term) # Dependent variables, not function calls
            return screen_unit(getmetadata(x, VariableUnit, unitless)) # Like x(t) or x[i]
        elseif op isa Term && !(operation(op) isa Term)
            gp = getmetadata(x,Symbolics.GetindexParent,nothing) # Like x[1](t)
            return screen_unit(getmetadata(gp, VariableUnit, unitless))
        end  # Actual function calls:
        args = arguments(x)
        return get_unit(op, args)
    else # This function should only be reached by Terms, for which `istree` is true
        throw(ArgumentError("Unsupported value $x."))
    end
end

"Get unit of term, returning nothing & showing warning instead of throwing errors."
function safe_get_unit(term, info)
    side = nothing
    try
        side = get_unit(term)
    catch err
        if err isa Unitful.DimensionError
            @warn("$info: $(err.x) and $(err.y) are not dimensionally compatible.")
        elseif err isa ValidationError
            @warn(info*err.message)
        elseif err isa MethodError
            @warn("$info: no method matching $(err.f) for arguments $(typeof.(err.args)).")
        else
            rethrow()
        end
    end
    side
end

function _validate(terms::Vector, labels::Vector{String}; info::String = "")
    valid = true
    first_unit = nothing
    first_label = nothing
    for (term,label) in zip(terms,labels)
        equnit = safe_get_unit(term, info*label)
        if equnit === nothing
            valid = false
        elseif !isequal(term,0)
            if first_unit === nothing
                first_unit = equnit
                first_label = label
            elseif !equivalent(first_unit, equnit)
                valid = false
                @warn("$info: units [$(first_unit)] for $(first_label) and [$(equnit)] for $(label) do not match.")
            end
        end
    end
    valid
end

function validate(jump::Union{ModelingToolkit.VariableRateJump, ModelingToolkit.ConstantRateJump}, t::Symbolic; info::String = "")
    newinfo = replace(info,"eq."=>"jump")
    _validate([jump.rate, 1/t], ["rate", "1/t"], info = newinfo) && # Assuming the rate is per time units
    validate(jump.affect!,info = newinfo)
end

function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
    left_symbols = [x[1] for x in jump.reactant_stoch] #vector of pairs of symbol,int -> vector symbols
    net_symbols = [x[1] for x in jump.net_stoch]
    all_symbols = vcat(left_symbols,net_symbols)
    allgood = _validate(all_symbols, string.(all_symbols), info = info)
    n = sum(x->x[2],jump.reactant_stoch,init = 0)
    base_unitful = all_symbols[1] #all same, get first
    allgood && _validate([jump.scaled_rates, 1/(t*base_unitful^n)], ["scaled_rates", "1/(t*reactants^$n))"], info = info)
end

function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Symbolic)
    labels = ["in Mass Action Jumps,", "in Constant Rate Jumps,", "in Variable Rate Jumps,"]
    all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
end

validate(eq::ModelingToolkit.Equation; info::String = "") = _validate([eq.lhs, eq.rhs], ["left", "right"], info = info)
validate(eq::ModelingToolkit.Equation, term::Union{Symbolic,Unitful.Quantity,Num}; info::String = "") = _validate([eq.lhs, eq.rhs, term], ["left","right","noise"], info = info)
validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "") = _validate(vcat([eq.lhs, eq.rhs], terms), vcat(["left", "right"], "noise  #".*string.(1:length(terms))), info = info)

"Returns true iff units of equations are valid."
validate(eqs::Vector; info::String = "") = all([validate(eqs[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, noise::Vector; info::String = "") = all([validate(eqs[idx], noise[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, noise::Matrix; info::String = "") = all([validate(eqs[idx], noise[idx, :], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(eqs::Vector, term::Symbolic; info::String = "") = all([validate(eqs[idx], term, info = info*" in eq. #$idx") for idx in 1:length(eqs)])
validate(term::Symbolics.SymbolicUtils.Symbolic) = safe_get_unit(term,"") !== nothing

"Throws error if units of equations are invalid."
check_units(eqs...) = validate(eqs...) || throw(ValidationError("Some equations had invalid units. See warnings for details."))
all_dimensionless(states) = all(x->safe_get_unit(x,"") in (unitless,nothing),states)


text:
"""
    most_common_colors(img::String)
img: path to local image file
n: how many common colors to return
"""
function most_common_colors(img::String, n::Int)

    colorthief = pyimport("colorthief")
    color_thief = colorthief.ColorThief(img)
    cols = color_thief.get_palette(color_count=n)

	function getone(x)
		RGB(cols[x][1]/255, cols[x][2]/255,  cols[x][3]/255)
	end

	output = []
	for i in 1:n-1
		push!(output,getone(i))
	end

    return output

end


text:
# Defines equations of state which are combinations of others.

using JLD  # for accessing data files in .jld files (HDF5-based)
using ProgressMeter  # for showing progress bars during calculations

export PressurePiecewiseEOS, StitchedEOS

# This EOS is used for signalling we're outside the domain
struct OutOfDomainEOS <: EOS end
(::OutOfDomainEOS)(args...) = throw(DomainError())

# Piecewise one-dimensional EOS

"Equation of state which is stored piecewise in some coordinate"
abstract type PiecewiseEOS <: EOS end

""" Equation of state which is stored piecewise in the pressure coordinate

    * `eoses`: Vector of `EOS`es. Any EOS can be chosen, though it's assumed you
      will build a `PressurePiecewiseEOS` from one-dimensional `EOS`es.
    * `edges`: Vector of pressure values that bracket each individual EOS. For
      example, the array [0, 1, 2] would serve to define an EOS with a domain
      [0, 2], piecewise from [0, 1] and (1, 2].

    Calling the `PressurePiecewiseEOS` will evaluate the correct EOS at a given
    pressure. """
struct PressurePiecewiseEOS <: PiecewiseEOS
    eoses::Vector{EOS}
    edges::Vector{Float64}

    function PressurePiecewiseEOS(eoses, edges)
        @assert(length(edges) == length(eoses) + 1,
                "EOS/edge mismatch: need one more edge than number of EOSes")
        new(eoses, edges)
    end
end

"Find the appropriate individual EOS in a 1D piecewise `eos` at point `x`"
function extracteos(eos::PiecewiseEOS, x::Number)
    if isoutside(x, first(eos.edges), last(eos.edges))
        if x ≈ first(eos.edges) || x ≈ last(eos.edges)
            # might be caused by floating point problems
            warn("""Attempted to evaluate a piecewise EOS just outside its domain.
                    This may be a floating point precision error.
                    Consider widening the domain slightly.""")
            end
        return OutOfDomainEOS()
    end

    for (layer_num, layer_end) in enumerate(eos.edges[2:end])
        if x <= layer_end
            # it's within the half-open interval (x, y]
            return eos.eoses[layer_num]
        end
    end
end
extracteos(eos::EOS, args...) = eos  # fallback for when it's not piecewise

# calling a PressurePiecewiseEOS just evaluates the appropriate EOS
(eos::PressurePiecewiseEOS)(P) = extracteos(eos, P)(P)
(eos::PressurePiecewiseEOS)(P, T) = extracteos(eos, P)(P, T)

""" Save piecewise EOSes to `eos-piecewise.jld`:

    * The water EOS used by Sara Seager in her 2007 paper, consisting of an ice
    VII BME at low pressures, an intermediate density functional theory
    calculation, and the TFD at high pressures.
    * A copy of the above which is bounded to the region of (P,T) space defined
    by ice X, continuing to higher pressures.
    * The MgSiO3 EOS from that paper, consisting of a 4th-order BME at low
    pressures and the TFD at high pressures.
    * The Fe EOS from that paper, consisting of a Vinet fit at low pressures and
    the TFD at high pressures."""
function save_piecewise_eoses!()
    funcs = load("$(config.datadir)/eos-functional.jld")
    tables = load("$(config.datadir)/eos-tabular.jld")

    phaseregions = load_phase_boundaries()["regions"]

    "Seager's H2O piecewise EOS"
    h2o = let
        eoses = [funcs["seager"]["h2o"],
                 tables["seager_dft"],
                 funcs["seager"]["h2o_tfd"]]
        transitionP = [0, 44.3e9, 7686e9, Inf]
        raw = PressurePiecewiseEOS(eoses, transitionP)
        P = logspace(5, 19)
        ρ = map(raw, P)
        LineEOS(P, ρ)
    end

    "Seager's MgSiO3 piecewise EOS"
    mgsio3 = let
        eoses = [funcs["seager"]["mgsio3_pv"],
                 funcs["seager"]["mgsio3_tfd"]]
        transitionP = [0, 1.35e13, Inf]
        raw = PressurePiecewiseEOS(eoses, transitionP)
        P = logspace(5, 19)
        ρ = map(raw, P)
        LineEOS(P, ρ)
    end

    "Seager's Fe piecewise EOS"
    fe = let
        eoses = [funcs["seager"]["fe_eps"],
                 funcs["seager"]["fe_tfd"]]
        transitionP = [0, 2.09e14, Inf]
        raw = PressurePiecewiseEOS(eoses, transitionP)
        P = logspace(5, 19)
        ρ = map(raw, P)
        LineEOS(P, ρ)
    end

    iceXbound = phaseregions["X"]
    beyondXbound = BoundingBox(4e10, Inf, -Inf, Inf)
    iceX = BoundedEOS(h2o, iceXbound)
    iceX_beyond = BoundedEOS(h2o, beyondXbound)

    jldopen("$(config.datadir)/eos-piecewise.jld", "w") do file
        addrequire(file, WaterData)
        write(file, "h2o", h2o)
        write(file, "mgsio3", mgsio3)
        write(file, "fe", fe)
        write(file, "iceX", iceX)
        write(file, "iceX_beyond", iceX_beyond)
    end
end

# Stitched (piecewise 2D) EOS

"Equation of state which consists of several other EOS stitched together"
struct StitchedEOS <: EOS
    eoses::Vector{EOS}
end
StitchedEOS(a::EOS, b...) = StitchedEOS([a, b...])

"Bounding box for this EOS"
function BoundingBox(s::StitchedEOS)
    bboxes = [BoundingBox(eos) for eos in s.eoses]
    xmin = minimum([bb.xmin for bb in bboxes])
    xmax = maximum([bb.xmax for bb in bboxes])
    ymin = minimum([bb.ymin for bb in bboxes])
    ymax = maximum([bb.ymax for bb in bboxes])
    BoundingBox(xmin, xmax, ymin, ymax)
end

"Find the appropriate individual EOS in a 2D piecewise `StitchedEOS` at point (`P`,`T`)"
function extracteos(s::StitchedEOS, P, T)
    i = findfirst(e -> (P, T) in e, s.eoses)
    if i == 0
        # no EOS was found for this value of P and T
        OutOfDomainEOS()
    else
        s.eoses[i]
    end
end

# calling a PressurePiecewiseEOS just evaluates the appropriate EOS
(s::StitchedEOS)(P, T) = extracteos(s, P, T)(P, T)

# Full EOS generation

""" Save the full EOS and thermal expansivity to `eos-full.jld`:

    * `raw`: the full `StitchedEOS`, consisting of all H2O tabular and
    functional EOSes
    * `grid`: `GridEOS` version of the above (much faster to evaluate)
    * `gridPlusIdeal`: switches to ideal gas below a certain threshold
    * `thermexp`: the precomputed thermal expansivity for the above """
function save_full_eos!()
    funcs = load_functional_eoses()
    tables = load_tabular_eoses()
    piecewise = load_piecewise_eoses()

    water_idealgas = funcs["misc"]["ideal_gas"]

    # In case of overlapping domains, EOSes listed earlier here have priority
    eos = WaterData.StitchedEOS(tables["iapws"],
                                tables["sugimura"],
                                tables["french"],
                                funcs["choukroungrasset"]["I"],
                                funcs["choukroungrasset"]["III"],
                                funcs["choukroungrasset"]["V"],
                                funcs["choukroungrasset"]["VI"],
                                funcs["misc"]["mgd_iceVII"],
                                funcs["misc"]["iapws_pastfrench"],
                                piecewise["iceX"],
                                piecewise["iceX_beyond"],
                                funcs["misc"]["iapws_highpressure"],
                                funcs["misc"]["iapws_highprestemp"],
                                funcs["misc"]["iapws_hightemp"],
                                funcs["misc"]["fallback"])

    # set the resolution
    Nx = config.grid_resolution
    Ny = config.grid_resolution

    # set the region to be evaluated
    bb = BoundingBox(config.Pmin, config.Pmax, config.Tmin, config.Tmax)
    Ps = logspace(log10(bb.xmin), log10(bb.xmax), Nx)
    Ts = logspace(log10(bb.ymin), log10(bb.ymax), Ny)

    # compute ρ
    ρs = zeros(Nx, Ny)
    meter = Progress(Nx*Ny, "Gridding ρ...")
    for (j, T) in enumerate(Ts), (i, P) in enumerate(Ps)
        ρs[i, j] = eos(P, T)
        next!(meter)
    end

    # compute α
    αs = zeros(Nx, Ny)
    meter = Progress(Nx*Ny, "Calculating α...")
    for (j, T) in enumerate(Ts), (i, P) in enumerate(Ps)
        αs[i, j] = thermalexpansivity(eos, P, T)
        next!(meter)
    end

    # We throw away values less than zero, as we don't physically expect this
    # behaviour (except perhaps for low-pressure ices which have negative
    # expansivity, but these aren't in the region of phase space we're
    # interested in)
    clamp!(αs, 0, Inf)

    # make EOSes and write to file
    # this EOS uses everything except the ideal gas EOS
    grideos = GridEOS(Ps, Ts, ρs)
    # this EOS is the same, but adds the ideal gas below 100bar
    gridPlusIdeal = PressurePiecewiseEOS([water_idealgas, grideos], [0, 100e5, Inf])
    thermexp = GridEOS(Ps, Ts, αs)
    jldopen("$(WaterData.config.datadir)/eos-full.jld", "w") do file
        write(file, "raw", eos)
        write(file, "grid", grideos)
        write(file, "gridPlusIdeal", gridPlusIdeal)
        write(file, "thermexp", thermexp)
    end
end


text:
module CoilFunctions

# import CSV
using Unitful, DataStructures, MappedArrays
using ThreadTools: tmap
using Transducers
import Unitful: 𝐈,𝐋,𝐌,𝐓
using Roots: find_zero

export ideal_fill, AWG_Chart, optimalcoil, estimatetruefill, enclosewinding, resistance

@derived_dimension ResistanceLength dimension(u"Ω/m") #𝐈^-2*𝐋*𝐌*𝐓^-3

include("types.jl")
include("defaults.jl")
include("winding.jl")
include("coil.jl")
include("optimization.jl")

include("drawing.jl")
# using .Drawing

const AWG_Chart=let
    src=readlines(joinpath((@__DIR__),"AWG Chart.csv"))
    data=map(src[2:end]) do x
        x=parse.(Float64,split(x,','))
        AWG(x[1],x[2]u"inch",x[5]u"inch",x[3]*u"Ω/1000ft",x[7]u"ft/lb")
    end
    AWG_Chart=Dict(x.gauge=>x for x in data)
    AWG_Chart=OrderedDict(sort(AWG_Chart, by=first))
    #Filter out half gauages
    filter!(x->isinteger(x[1]),AWG_Chart)
end

const _default_awg=AWG_Chart[30.0]

≲(a,b) = (a ≈ b) || (a < b)
≳(a,b) = (a ≈ b) || (a > b)

function istrianglesorted(a::AbstractVector)
    i = argmax(a)
    issorted(@view a[1:i]) && issorted(@view a[i:end];lt= > )
end

end # module


text:
using PerformanceTestTools

PerformanceTestTools.@include_foreach(
    "test_threaded.jl",
    [nothing,
    ["JULIA_NUM_THREADS" => Threads.nthreads() > 1 ? "1" : "2"],
    ],
)


text:
module PiADS1x15

using PiGPIO

export ADS1015, ADS1115
export set_config, read_conv_result, set_threshold_alert, enable_conv_rdy_alert, is_busy

# I2C Addresses
const addresses = Dict(
  :GND => 0x48, # I2C address if ADDR pin is connected to GND
  :VDD => 0x49, # I2C address if ADDR pin is connected to VCC
  :SDA => 0x4a, # I2C address if ADDR pin is connected to SDA
  :SCL => 0x4b # I2C address if ADDR pin is connected to SCL
)

# Pointer register
const pointers = Dict(
  :CONVERT => 0x00, # Convert register
  :CONFIG => 0x01, # Config register
  :LOWTHRESH => 0x02, # Low threshold register
  :HITHRESH => 0x03 # High threshold register
)

# Config register
const config_os = Dict(
  :MASK_OS => 0x8000,
  :STARTSINGLE => 0x8000, # Write: Start a single-conversion
  :BUSY => 0x0000, # Read: Device is currently performing a conversion
  :NOTBUSY => 0x8000 # Read: Device is in idle mode
)

const config_mux = Dict(
  :MASK_MUX => 0x7000,
  :DIFF_0_1 => 0x0000, # AIN_P = AIN0 and AIN_N = AIN1 (differential, default)
  :DIFF_0_3 => 0x1000, # AIN_P = AIN0 and AIN_N = AIN3 (differential)
  :DIFF_1_3 => 0x2000, # AIN_P = AIN1 and AIN_N = AIN3 (differential)
  :DIFF_2_3 => 0x3000, # AIN_P = AIN2 and AIN_N = AIN3 (differential)
  :SINGLE_0 => 0x4000, # AIN_P = AIN0 and AIN_N = GND (single-ended)
  :SINGLE_1 => 0x5000, # AIN_P = AIN1 and AIN_N = GND (single-ended)
  :SINGLE_2 => 0x6000, # AIN_P = AIN2 and AIN_N = GND (single-ended)
  :SINGLE_3 => 0x7000 # AIN_P = AIN3 and AIN_N = GND (single-ended)
)

const config_pga = Dict(
  :MASK_PGA => 0x0E00,
  :V6_144 => 0x0000, # +-6.144V range
  :V4_096 => 0x0200, # +-4.096V range
  :V2_048 => 0x0400, # +-2.048V range (default)
  :V1_024 => 0x0600, # +-1.024V range
  :V0_512 => 0x0800, # +-0.512V range
  :V0_256 => 0x0A00 # +-0.256V range
)

const config_mode = Dict(
  :MASK_MODE => 0x0100,
  :CONTINOUS => 0x0000, # Continuous-conversion mode
  :SINGLE => 0x0100 # Single-shot mode or power-down state (default)
)

const config_dr = Dict(
  :MASK_DR => 0x00E0,
  :SPS128 => 0x0000, # 128 samples per second
  :SPS250 => 0x0020, # 250 samples per second
  :SPS490 => 0x0040, # 490 samples per second
  :SPS920 => 0x0060, # 920 samples per second
  :SPS1600 => 0x0080, # 1600 samples per second (default)
  :SPS2400 => 0x00A0, # 2400 samples per second
  :SPS3300 => 0x00C0 # 3300 samples per second
)

const config_compmode = Dict(
  :MASK_COMPMODE => 0x0010,
  :TRADITIONAL => 0x0000, # Traditional comparator (default)
  :WINDOW => 0x0010 # Window comparator
)

const config_comppol = Dict(
  :MASK_COMPPOL => 0x0008,
  :ACTIVELOW => 0x0000, # Active low ALERT/RDY pin (default)
  :ACTIVEHI => 0x0008 # Active high ALERT/RDY pin
)

const config_complat = Dict(
  :MASK_COMPLAT => 0x0004,
  :NONLATCH => 0x0000, # Non-latching comparator (default)
  :LATCH => 0x0004 # Latching comparator
)

const config_compque = Dict(
  :MASK_COMPQUE => 0x0003,
  :CONV1 => 0x0000, # Assert ALERT/RDY after one conversions
  :CONV2 => 0x0001, # Assert ALERT/RDY after two conversions
  :CONV4 => 0x0002, # Assert ALERT/RDY after four conversions
  :NONE => 0x0003 # Disable the comparator and put ALERT/RDY to high impedance (default)
)

current_config = nothing # Current config register, needs to be set after startup

abstract type ADS1x15 end

mutable struct ADS1015 <: ADS1x15
  handle          # PiGPIO handle
  i2c_address      # I2C address
  conversion_delay # Conversion delay (s)
  bit_shift        # Bit shift
end

mutable struct ADS1115 <: ADS1x15
  handle          # PiGPIO handle
  i2c_address      # I2C address
  conversion_delay # Conversion delay (s)
  bit_shift        # Bit shift
end

function ADS1015(pi::Pi, i2c_bus::Integer, i2c_address::Symbol=:GND)
  global addresses
  handle = PiGPIO.i2c_open(pi, i2c_bus, addresses[i2c_address])
  return ADS1015(handle, i2c_address, 1//1000, 4)
end

function ADS1115(pi::Pi, i2c_bus::Integer, i2c_address::Symbol=:GND)
  global addresses
  handle = PiGPIO.i2c_open(pi, i2c_bus, addresses[i2c_address])
  return ADS1115(handle, i2c_address, 9//1000, 0)
end

"""
    write_register(pi, ads, register, value)

Write a 16-bits `value` to the specified destination `register`.
"""
function write_register(pi::Pi, ads::ADS1x15, register::Number, value::Number)
  return PiGPIO.i2c_write_word_data(pi, ads.handle, UInt8(register), UInt16(value))
end

"""
  read_register(pi, ads, register)

Read 16-bits from the specified destination `register`.
"""
function read_register(pi::Pi, ads::ADS1x15, register::Number)
  return PiGPIO.i2c_read_word_data(pi, ads.handle, UInt8(register))
end

"""
  set_config(pi, ads [, os; mux, pga, mode, dr, comp_mode, comp_pol, comp_lat, comp_que])

Set the config register of a ADS1x15 device `ads` that is connected to a Raspberry Pi via a PiGPIO.jl object `pi`.
The conversion results can be collected afterwards using [`read_conv_result`](@ref).
"""
function set_config(pi::Pi, ads::ADS1x15;
  os::Symbol=:STARTSINGLE,
  mux::Symbol=:DIFF_0_1,
  pga::Symbol=:V2_048,
  mode::Symbol=:SINGLE,
  dr::Symbol=:SPS1600,
  compmode::Symbol=:TRADITIONAL,
  comppol::Symbol=:ACTIVELOW,
  complat::Symbol=:NONLATCH,
  compque::Symbol=:NONE)

  global pointers
  global config_os
  global config_mux
  global config_pga
  global config_mode
  global config_dr
  global config_compmode
  global config_comppol
  global config_complat
  global config_compque
  global current_config

  config = config_os[os] | config_mux[mux] | config_pga[pga] | config_mode[mode] | config_dr[dr] |
  config_compmode[compmode] | config_comppol[comppol] | config_complat[complat] | config_compque[compque]

  if write_register(pi, ads, pointers[:CONFIG], config) == 0
    current_config = config # Save current config register
    return true
  else
    return false
  end
end

"""
  read_conv_result(pi, ads)

Get a conversion result from an ADS1x15 device connected to a Raspberry Pi via a PiGPIO.jl object `pi`.
The config register needs to be set sepeartely beforehand using [`set_config`](@ref).
"""
function read_conv_result(pi::Pi, ads::ADS1x15)
  global pointers
  global config_mux
  global current_config

  if current_config == nothing
    error("Config register needs to be set before reading values.")
  end

  sleep(ads.conversion_delay)
  result = read_register(pi, ads, pointers[:CONVERT]) >> ads.bit_shift
  if ads.bit_shift != 0 && (current_config & config_mux[:MASK_MUX]) <= config_mux[:DIFF_2_3] && result > 0x07FF
    result |= 0xF000 # Move sign to 16th bit in case of negative differential values
  end
  return result
end

"""
  read_config(pi, ads)

Read the current config register from an ADS1x15 device connected to a Raspberry Pi via a PiGPIO.jl object `pi`.
"""
function read_config(pi::Pi, ads::ADS1x15)
  global pointers
  global config_os
  global config_mux
  global config_pga
  global config_mode
  global config_dr
  global config_compmode
  global config_comppol
  global config_complat
  global config_compque
  global current_config

  res = read_register(pi, ads, pointers[:CONFIG])
  if res != 0
    current_config = res
    if field == :NONE
      return current_config
    elseif field == :OS
      return current_config & config_os[:MASK_OS]
    elseif field == :MUX
      return current_config & config_mux[:MASK_MUX]
    elseif field == :PGA
      return current_config & config_pga[:MASK_PGA]
    elseif field == :MODE
      return current_config & config_mode[:MASK_MODE]
    elseif field == :DR
      return current_config & config_dr[:MASK_DR]
    elseif field == :COMPMODE
      return current_config & config_compmode[:MASK_COMPMODE]
    elseif field == :COMPPOL
      return current_config & config_comppol[:MASK_COMPPOL]
    elseif field == :COMPLAT
      return current_config & config_complat[:MASK_COMPLAT]
    elseif field == :COMPQUE
      return current_config & config_compque[:MASK_COMPQUE]
    else
      error("Unknown field $field.")
    end
  end
end

"""
  set_threshold_ADS1x15(pi, ads, thld_value[, write_low_thld])

Set the low or high (see boolean `write_low_thld`) comperator alert threshold to a specific digital value `thld_value`.
Use function [`enable_conv_rdy_alert`](@ref) in order to enable a alert pin changes on newly available conversion results.
"""
function set_threshold_alert(pi::Pi, ads::ADS1x15, thld_value::Number, write_low_thld::Bool=true)
  global pointers
  register = write_low_thld ? pointers[:LOWTHRESH] : pointers[:HITHRESH]
  return write_register(pi, ads, register, thld_value << ads.bit_shift)
end

"""
  enable_conv_rdy_alert(pi, ads)

Enable alert pin changes on newly available conversion results.
"""
function enable_conv_rdy_alert(pi::Pi, ads::ADS1x15)
  global pointers
  write_register(pi, ads, pointers[:LOWTHRESH], 0x800)
  write_register(pi, ads, pointers[:HITHRESH], 0x7ff)
end

"""
  is_busy(pi, ads)

Check if ADC is currenty performing a measurement.
"""
function is_busy(pi::Pi, ads::ADS1x15)
  global pointers
  global config_os
  return read_register(pi, ads, pointers[:CONFIG]) & config_os[:MASK_OS] == config_os[:BUSY]
end

end #module


text:
"""
    smoptimize(f::Function, search_range::Array{Tuple{Float64,Float64},1}; options=Options())

Optimize the function `f` in the range `search_range` using a Radial Basis Function based surrogate model.
"""
function smoptimize(f::Function, search_range::Array{Tuple{Float64,Float64},1}; options::Options=Options())

    #Load some option values
    @unpack num_start_samples, sampling_plan_opt_gens,
            iterations, trace = options

    #Create sampling plan
    lhc_plan = scaled_LHC_sampling_plan(search_range,num_start_samples,sampling_plan_opt_gens;trace=trace)

    #Evaluate sampling plan
    lhc_samples = f_opt_eval(f,lhc_plan,trace)

    #Initialize variables to be returned
    sm_interpolant = nothing
    infill_type = Array{Symbol,1}(undef,0)
    infill_prediction = Array{Float64,1}(undef,0)
    optres = nothing
    infill_plan = Array{Float64,2}(undef,size(lhc_plan,1),0)
    infill_sample = Array{Float64,2}(undef,1,0)

    #Run the entire optimization iterations number of times
    for i = 1:iterations
        trace && print_iteration(i,iterations)

        #Create the optimized Radial Basis Function interpolant
        samples_all = [lhc_samples infill_sample]
        plan_all = [lhc_plan infill_plan]
        sm_interpolant, optres = surrogate_model(plan_all, samples_all; options=options)

        #Points to add to the sampling plan to improve the interpolant
        infill_plan_new, infill_type_new, infill_prediction_new, options  = model_infill(search_range,plan_all,
                                                                                samples_all,sm_interpolant;options=options)

        #Evaluate the new infill points
        infill_sample_new = f_opt_eval(f,infill_plan_new,samples_all,trace)

        #Add infill points
        infill_plan = [infill_plan infill_plan_new]
        infill_sample = [infill_sample infill_sample_new]
        infill_type = [infill_type; infill_type_new]
        infill_prediction = [infill_prediction; infill_prediction_new]

    end

    return SurrogateResult( lhc_samples, lhc_plan, sm_interpolant,
                            optres, infill_sample, infill_type,
                            infill_plan, infill_prediction,options)
end

"""
    f_opt_eval(f,plan,samples,trace)

Calculate the objective function value(s) and provide intermediate results printing
showing improvements over previous best iteration.
"""
function f_opt_eval(f,plan,samples,trace)

    trace && println("Evaluating function ",size(plan,2)," times ...")

    new_samples = mapslices(f,plan,dims=1)

    if trace
        _, min_loc = findmin(new_samples)
        for i = 1:size(plan,2)
            if i == min_loc[2]
                printstyled(@sprintf("%-15.7g",new_samples[i]); color=:light_green, bold=true)
            else
                printstyled(@sprintf("%-15.7g",new_samples[i]))
            end
        end
        print("\t actual value\n")
        println("---------------------------------------------------------------")

        new_min = minimum(new_samples)
        old_min = minimum(samples)
        new_max = maximum(new_samples)
        old_max = maximum(samples)

        print("Max and min sample value: ")
        printstyled(@sprintf("%.7g",maximum((new_max,old_max))); color=:light_red)
        print("\t")
        printstyled(@sprintf("%.7g",minimum((new_min,old_min))); color=:green, bold=true)

        if isless(new_min,old_min)
            print("\t (Improvement from last iteration ", @sprintf("%.7g",old_min-new_min),")")
        else
            print("\t (Improvement from last iteration N/A)")
        end
        print("\n")
    end

    return new_samples
end

"""
    f_opt_eval(f,plan,trace)

Calculate the objective function value(s) and plot value if trace.
"""
function f_opt_eval(f,plan,trace)

    trace && println("Evaluating function ",size(plan,2)," times ...")

    new_samples = mapslices(f,plan,dims=1)

    if trace
        new_min = minimum(new_samples)
        new_max = maximum(new_samples)

        print("Max and min sample value: ")
        printstyled(@sprintf("%.7g",new_max); color=:light_red)
        print("\t")
        printstyled(@sprintf("%.7g",new_min); color=:light_green, bold=true)
        print("\n")
    end

    return new_samples
end

"""
    print_iteration(i,iterations)

REPL printing of the current iteration if trace.
"""
function print_iteration(i,iterations)
    print("\n \n \n \t Iteration ")
    printstyled(i,bold=true)
    print(" out of ", iterations, "\n")
end


text:
export RepeatingAction, stop!, schedule_repeating!
export SequenceOfActions, ShuffledActions, TentativeAction, NoAction

type RepeatingAction
    interval::Float64
    order::Int64
    actionable

    function RepeatingAction(interval, order, actionable)
        if interval <= 0
            throw(ArgumentError("Interval must be a positive number."))
        end

        new(interval, order, actionable)
    end
end

RepeatingAction(action::Function, interval, order) = RepeatingAction(
    interval, order, action
)

stop!(repeating::RepeatingAction) = repeating.interval = -1

function Base.call(repeating::RepeatingAction, env, schedule::Schedule)
    if repeating.interval > 0
        repeating.actionable(env, schedule)

        schedule_once_in!(
            schedule, repeating, repeating.interval, repeating.order
        )
    end
end

function schedule_repeating!(schedule::Schedule, action,
                             starting_at=1.0, interval=1.0, order=1)
    repeating = RepeatingAction(interval, order, action)
    schedule_once!(
        schedule,
        repeating,
        starting_at,
        order
    )
    repeating
end

function schedule_repeating!(action::Function, schedule::Schedule,
                             starting_at=1.0, interval=1.0, order=1)
    schedule_repeating!(schedule, action, starting_at, interval, order)
end

type SequenceOfActions
    actions

    function SequenceOfActions(actions::AbstractVector)
        new(Vector(actions))
    end
end

function Base.call(seq::SequenceOfActions, env, schedule::Schedule)
    for action in seq.actions
        action(env, schedule)
    end
end

type ShuffledActions
    actions

    function ShuffledActions(actions::AbstractVector)
        new(Vector(actions))
    end
end

function Base.call(seq::ShuffledActions, env, schedule::Schedule)
    shuffle!(seq.actions)
    for action in seq.actions
        action(env, schedule)
    end
end

immutable NoAction
end

function Base.call(::NoAction, env, schedule::Schedule)
end

type TentativeAction
    actionable
end

function Base.call(tentative::TentativeAction, env, schedule::Schedule)
    tentative.actionable(env, schedule)
end

stop!(tentative::TentativeAction) = tentative.actionable = NoAction()


text:
using WiltonInts84

using Base.Test
using FixedSizeArrays

include("num_quad.jl")

if !isdefined(:record)
    record = false
end

p1 = Vec(0.4,0.7000000000000001,0.0)
p2 = Vec(0.5,0.7000000000000001,0.0)
p3 = Vec(0.5,0.6000000000000001,0.0)
z  = Vec(0.0,0.0,1.0)

N = 2
T = eltype(Float64)
P = Vec{3,T}
J = Vector{Vec{N+3,T}}()
L = Vector{Vec{N+3,P}}()
I = similar(J)
K = similar(L)

h = 1.2

# circle through triangle vertex
c = Vec(0.020000000000000004,0.06,0.0)
r, R = 0.0, 0.8
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 0
@test length(q.arcs) == 0
@test length(q.circles) == 0
A,B = wiltonints(q,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end


text:
include("init-namedarrays.jl")

print("enameration, ")

@testset "Enameration" begin
	names = [("one", "a"), ("two", "a"), ("one", "b"), ("two", "b"), ("one", "c"), ("two", "c"), ("one", "d"), ("two", "d")]
	i = 1
	for (nameTuple, val) in enamerate(n)
		@test nameTuple == names[i]
		i += 1
		@test val === n[nameTuple...]
	end
end


text:
function runsim!(sim::Simulation, n)
    sim.step !== 0 && error("Simulation is not at initial state.")
    poiseuille!(sim)
    boundaries!(sim)

    un_1 = copy(sim.u)
    vn_1 = copy(sim.v)
    ti = @elapsed begin
        make_step!(sim)
    end
    with_logger(sim.logger) do
        @info "Step number $step done in $(round(ti * 1.0e3, digits = 2))ms\nSimulation time : $(sim.step * sim.dt)"
    end
    for step in 2:n
        ti = @elapsed begin
            make_step!(sim, un_1, vn_1)
        end
        with_logger(sim.logger) do
            @info """Step number $step done in $(round(ti * 1.0e3, digits = 2))ms\nSimulation time : $(sim.step * sim.dt)
            Diagnosis :
            Re_h = $(reynh(sim))
            Re_h_omega = $(reynhvort(sim))
            """
        end
        un_1 = copy(sim.u)
        vn_1 = copy(sim.v)
    end
end

function make_step!(sim::Simulation)
    # if sim.step == 0
    step_euler!(sim)
    # else

    step_poisson!(sim, sim.params.tol)
    poisson_project!(sim)
    boundaries!(sim)
    sim.step = sim.step + 1
end

function make_step!(sim::Simulation, un_1::UGrid, vn_1::VGrid)
    # if sim.step == 0
    step_AB2!(sim, un_1, vn_1)
    # else

    step_poisson!(sim, sim.params.tol)
    poisson_project!(sim)
    boundaries!(sim)
    sim.step = sim.step + 1
end


text:
module GPUAcceleratedTracking

using
    BenchmarkTools,
    CUDA,
    GNSSSignals,
    StructArrays,
    StaticArrays,
    Parameters,
    Tracking,
    Statistics,
    Dates,
    DataFrames,
    Query,
    PrettyTables,
    CairoMakie,
    ColorSchemes
    # Plots,
    # PGFPlotsX

import Unitful: MHz, kHz, Hz, s, ms, dBHz, ustrip, NoUnits
import Tracking: TrackingState, NumAnts, NumAccumulators

struct KernelAlgorithm{x}
end

KernelAlgorithm(x) = KernelAlgorithm{x}()

struct ReductionAlgorithm{x}
end

ReductionAlgorithm(x) = ReductionAlgorithm{x}()

struct ReplicaAlgorithm{x}
end

ReplicaAlgorithm(x) = ReplicaAlgorithm{x}()

const GNSSDICT = Dict(
    "GPSL1" => GPSL1,
    "GPSL5" => GPSL5,
)

const ALGODICT = Dict(
    # "1_3_pure"                  => 1300,
    # "1_3_pure_textmem"          => 1301,
    # "1_3_cplx"                  => 1320,
    # "1_3_cplx_textmem"          => 1321,
    "1_3_cplx_multi"            => 1330,
    "1_3_cplx_multi_textmem"    => 1331,
    # "1_4_cplx_multi"            => 1430,
    "1_4_cplx_multi_textmem"    => 1431,
    "2_3_cplx_multi"            => 2330,
    "2_3_cplx_multi_textmem"    => 2331,
    "2_4_cplx_multi"            => 2430,
    "2_4_cplx_multi_textmem"    => 2431,
    "3_4_cplx_multi"            => 3430,
    "3_4_cplx_multi_textmem"    => 3431,
    "4_4_cplx_multi_textmem"    => 4431,
    "5_4_cplx_multi_textmem"    => 5431
)

const REDDICT = Dict(
    "pure" => ReductionAlgorithm(1),
    "cplx" => ReductionAlgorithm(2),
    "cplx_multi" => ReductionAlgorithm(3)
)

const MEMDICT = Dict(
    "gmem" => ReplicaAlgorithm(1),
    "textmem" => ReplicaAlgorithm(2)
)

const ALGODICTINV = Dict(
    1300 => "1_3_pure"                   ,
    1301 => "1_3_pure_textmem"           ,
    1320 => "1_3_cplx"                   ,
    1331 => "1_3_cplx_textmem"           ,
    1330 => "1_3_cplx_multi"             ,
    1331 => "1_3_cplx_multi_textmem"     ,
    1430 => "1_4_cplx_multi"             ,
    1431 => "1_4_cplx_multi_textmem"     ,
    2330 => "2_3_cplx_multi"             ,
    2331 => "2_3_cplx_multi_textmem"     ,
    2430 => "2_4_cplx_multi"             ,
    2431 => "2_4_cplx_multi_textmem"     ,
    3430 => "3_4_cplx_multi"             ,
    3431 => "3_4_cplx_multi_textmem"     ,
    4431 => "4_4_cplx_multi_textmem"     ,
    5431 => "5_4_cplx_multi_textmem"
)
include("algorithms.jl")
include("reduction.jl")
include("gen_signal.jl")
include("replica_benchmarks.jl")
include("benchmarks.jl")
include("plots.jl")
include("results.jl")

export
    gen_signal,
    run_track_benchmark,
    run_kernel_benchmark,
    run_reduction_benchmark,
    run_replica_benchmark,
    add_results!,
    add_metadata!,
    eval_results,
    plot_reduction_benchmark,
    plot_replica_benchmark,
    reduce_3,
    reduce_4,
    reduce_5,
    reduce_cplx_3,
    reduce_cplx_4,
    reduce_cplx_5,
    reduce_cplx_multi_3,
    reduce_cplx_multi_31,
    reduce_cplx_multi_4,
    reduce_cplx_multi_5,
    reduce_cplx_multi_nant_5,
    gen_code_replica_kernel!,
    gen_code_replica_strided_kernel!,
    gen_code_replica_texture_mem_kernel!,
    gen_code_replica_texture_mem_strided_kernel!,
    gen_code_replica_texture_mem_strided_nsat_kernel!,
    downconvert_strided_kernel!,
    downconvert_and_accumulate_strided_kernel!,
    downconvert_and_correlate_kernel_1330!,
    downconvert_and_correlate_kernel_1331!,
    downconvert_and_correlate_kernel_1431!,
    downconvert_and_correlate_kernel_3431!,
    downconvert_and_correlate_kernel_4431!,
    downconvert_and_correlate_kernel_5431!,
    downconvert_and_correlate_kernel_3d_4431!,
    cpu_reduce_partial_sum,
    cuda_reduce_partial_sum,
    kernel_algorithm,
    KernelAlgorithm,
    ReductionAlgorithm,
    ReplicaAlgorithm

end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
# make.jl
push!(LOAD_PATH,"../src/")
using Documenter
using ThreeBodyProblem

makedocs(
    sitename = "ThreeBodyProblem.jl",
    # modules = [ThreeBodyProblem],
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Introduction" => "index.md"
        "Functions" => [
             "constants.md",
             "frames.md",
             "orbitalelements.md",
             "parameters.md",
             "util.md",
             "dynamics.md",
             "example.md"
        ]
        #     "costfunctions.md",
        #     "constraints.md",
        #     "creating_problems.md",
        #     "solving.md"
        # ],
        # "Interfaces" => [
        #     "costfunction_interface.md",
        #     "constraint_interface.md",
        #     "solver_interface.md"
        # ],
        # "Documentation" => [
        #     "model_types.md",
        #     "discretization.md",
        #     "cost_api.md",
        #     "constraint_api.md",
        #     "problem.md",
        #     "solvers.md",
        #     "rotations.md"
        # ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/jared711/ThreeBodyProblem.jl.git",
    target = "build"
    deploy_config=Documenter.Travis(),
    push_preview = true,
)

# Let's see if this will work


text:
using MPI

using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Distributed: index2rank, east_halo, west_halo, north_halo, south_halo, top_halo, bottom_halo

# Right now just testing with 4 ranks!
comm = MPI.COMM_WORLD
mpi_ranks = MPI.Comm_size(comm)
@assert mpi_ranks == 4

#####
##### Multi architectures and rank connectivity
#####

function test_triply_periodic_rank_connectivity_with_411_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(4, 1, 1))

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    @test local_rank == index2rank(arch.local_index..., arch.ranks...)

    connectivity = arch.connectivity

    # No communication in y and z.
    @test isnothing(connectivity.south)
    @test isnothing(connectivity.north)
    @test isnothing(connectivity.top)
    @test isnothing(connectivity.bottom)

    # +---+---+---+---+
    # | 0 | 1 | 2 | 3 |
    # +---+---+---+---+

    if local_rank == 0
        @test connectivity.east == 1
        @test connectivity.west == 3
    elseif local_rank == 1
        @test connectivity.east == 2
        @test connectivity.west == 0
    elseif local_rank == 2
        @test connectivity.east == 3
        @test connectivity.west == 1
    elseif local_rank == 3
        @test connectivity.east == 0
        @test connectivity.west == 2
    end

    return nothing
end

function test_triply_periodic_rank_connectivity_with_141_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    @test local_rank == index2rank(arch.local_index..., arch.ranks...)

    connectivity = arch.connectivity

    # No communication in x and z.
    @test isnothing(connectivity.east)
    @test isnothing(connectivity.west)
    @test isnothing(connectivity.top)
    @test isnothing(connectivity.bottom)

    # +---+
    # | 3 |
    # +---+
    # | 2 |
    # +---+
    # | 1 |
    # +---+
    # | 0 |
    # +---+

    if local_rank == 0
        @test connectivity.north == 1
        @test connectivity.south == 3
    elseif local_rank == 1
        @test connectivity.north == 2
        @test connectivity.south == 0
    elseif local_rank == 2
        @test connectivity.north == 3
        @test connectivity.south == 1
    elseif local_rank == 3
        @test connectivity.north == 0
        @test connectivity.south == 2
    end

    return nothing
end

function test_triply_periodic_rank_connectivity_with_114_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 1, 4))

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    @test local_rank == index2rank(arch.local_index..., arch.ranks...)

    connectivity = arch.connectivity

    # No communication in x and y.
    @test isnothing(connectivity.east)
    @test isnothing(connectivity.west)
    @test isnothing(connectivity.north)
    @test isnothing(connectivity.south)

    #   /---/
    #  / 3 /
    # /---/
    #   /---/
    #  / 2 /
    # /---/
    #   /---/
    #  / 1 /
    # /---/
    #   /---/
    #  / 0 /
    # /---/

    if local_rank == 0
        @test connectivity.top == 1
        @test connectivity.bottom == 3
    elseif local_rank == 1
        @test connectivity.top == 2
        @test connectivity.bottom == 0
    elseif local_rank == 2
        @test connectivity.top == 3
        @test connectivity.bottom == 1
    elseif local_rank == 3
        @test connectivity.top == 0
        @test connectivity.bottom == 2
    end

    return nothing
end

function test_triply_periodic_rank_connectivity_with_221_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(2, 2, 1))

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    @test local_rank == index2rank(arch.local_index..., arch.ranks...)

    connectivity = arch.connectivity

    # No communication in z.
    @test isnothing(connectivity.top)
    @test isnothing(connectivity.bottom)

    # +---+---+
    # | 0 | 2 |
    # +---+---+
    # | 1 | 3 |
    # +---+---+

    if local_rank == 0
        @test connectivity.east == 2
        @test connectivity.west == 2
        @test connectivity.north == 1
        @test connectivity.south == 1
    elseif local_rank == 1
        @test connectivity.east == 3
        @test connectivity.west == 3
        @test connectivity.north == 0
        @test connectivity.south == 0
    elseif local_rank == 2
        @test connectivity.east == 0
        @test connectivity.west == 0
        @test connectivity.north == 3
        @test connectivity.south == 3
    elseif local_rank == 3
        @test connectivity.east == 1
        @test connectivity.west == 1
        @test connectivity.north == 2
        @test connectivity.south == 2
    end

    return nothing
end

#####
##### Local grids for distributed models
#####

function test_triply_periodic_local_grid_with_411_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(4, 1, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    local_grid = model.grid
    nx, ny, nz = size(local_grid)

    @test local_grid.xF[1] == 0.25*local_rank
    @test local_grid.xF[nx+1] == 0.25*(local_rank+1)
    @test local_grid.yF[1] == 0
    @test local_grid.yF[ny+1] == 2
    @test local_grid.zF[1] == -3
    @test local_grid.zF[nz+1] == 0

    return nothing
end

function test_triply_periodic_local_grid_with_141_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    local_grid = model.grid
    nx, ny, nz = size(local_grid)

    @test local_grid.xF[1] == 0
    @test local_grid.xF[nx+1] == 1
    @test local_grid.yF[1] == 0.5*local_rank
    @test local_grid.yF[ny+1] == 0.5*(local_rank+1)
    @test local_grid.zF[1] == -3
    @test local_grid.zF[nz+1] == 0

    return nothing
end

function test_triply_periodic_local_grid_with_114_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 1, 4))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    local_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    local_grid = model.grid
    nx, ny, nz = size(local_grid)

    @test local_grid.xF[1] == 0
    @test local_grid.xF[nx+1] == 1
    @test local_grid.yF[1] == 0
    @test local_grid.yF[ny+1] == 2
    @test local_grid.zF[1] == -3 + 0.75*local_rank
    @test local_grid.zF[nz+1] == -3 + 0.75*(local_rank+1)

    return nothing
end

function test_triply_periodic_local_grid_with_221_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(2, 2, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    i, j, k = arch.local_index
    local_grid = model.grid
    nx, ny, nz = size(local_grid)

    @test local_grid.xF[1] == 0.5*(i-1)
    @test local_grid.xF[nx+1] == 0.5*i
    @test local_grid.yF[1] == j-1
    @test local_grid.yF[ny+1] == j
    @test local_grid.zF[1] == -3
    @test local_grid.zF[nz+1] == 0

    return nothing
end

#####
##### Injection of halo communication BCs
#####

function test_triply_periodic_bc_injection_with_411_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(4, 1, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        fbcs = field.boundary_conditions
        @test fbcs.east isa HaloCommunicationBC
        @test fbcs.west isa HaloCommunicationBC
        @test !isa(fbcs.north, HaloCommunicationBC)
        @test !isa(fbcs.south, HaloCommunicationBC)
        @test !isa(fbcs.top, HaloCommunicationBC)
        @test !isa(fbcs.bottom, HaloCommunicationBC)
    end
end

function test_triply_periodic_bc_injection_with_141_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        fbcs = field.boundary_conditions
        @test !isa(fbcs.east, HaloCommunicationBC)
        @test !isa(fbcs.west, HaloCommunicationBC)
        @test fbcs.north isa HaloCommunicationBC
        @test fbcs.south isa HaloCommunicationBC
        @test !isa(fbcs.top, HaloCommunicationBC)
        @test !isa(fbcs.bottom, HaloCommunicationBC)
    end
end

function test_triply_periodic_bc_injection_with_114_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(1, 1, 4))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        fbcs = field.boundary_conditions
        @test !isa(fbcs.east, HaloCommunicationBC)
        @test !isa(fbcs.west, HaloCommunicationBC)
        @test !isa(fbcs.north, HaloCommunicationBC)
        @test !isa(fbcs.south, HaloCommunicationBC)
        @test fbcs.top isa HaloCommunicationBC
        @test fbcs.bottom isa HaloCommunicationBC
    end
end

function test_triply_periodic_bc_injection_with_221_ranks()
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
    arch = MultiCPU(grid=full_grid, ranks=(2, 2, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        fbcs = field.boundary_conditions
        @test fbcs.east isa HaloCommunicationBC
        @test fbcs.west isa HaloCommunicationBC
        @test fbcs.north isa HaloCommunicationBC
        @test fbcs.south isa HaloCommunicationBC
        @test !isa(fbcs.top, HaloCommunicationBC)
        @test !isa(fbcs.bottom, HaloCommunicationBC)
    end
end

#####
##### Halo communication
#####

function test_triply_periodic_halo_communication_with_411_ranks(halo)
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(16, 6, 4), extent=(1, 2, 3), halo=halo)
    arch = MultiCPU(grid=full_grid, ranks=(4, 1, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        interior(field) .= arch.local_rank
        fill_halo_regions!(field, arch)

        @test all(east_halo(field, include_corners=false) .== arch.connectivity.east)
        @test all(west_halo(field, include_corners=false) .== arch.connectivity.west)

        @test all(interior(field) .== arch.local_rank)
        @test all(north_halo(field, include_corners=false) .== arch.local_rank)
        @test all(south_halo(field, include_corners=false) .== arch.local_rank)
        @test all(top_halo(field, include_corners=false) .== arch.local_rank)
        @test all(bottom_halo(field, include_corners=false) .== arch.local_rank)
    end

    return nothing
end

function test_triply_periodic_halo_communication_with_141_ranks(halo)
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(4, 16, 4), extent=(1, 2, 3), halo=halo)
    arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        interior(field) .= arch.local_rank
        fill_halo_regions!(field, arch)

        @test all(north_halo(field, include_corners=false) .== arch.connectivity.north)
        @test all(south_halo(field, include_corners=false) .== arch.connectivity.south)

        @test all(interior(field) .== arch.local_rank)
        @test all(east_halo(field, include_corners=false) .== arch.local_rank)
        @test all(west_halo(field, include_corners=false) .== arch.local_rank)
        @test all(top_halo(field, include_corners=false) .== arch.local_rank)
        @test all(bottom_halo(field, include_corners=false) .== arch.local_rank)
    end

    return nothing
end

function test_triply_periodic_halo_communication_with_114_ranks(halo)
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(4, 4, 16), extent=(1, 2, 3), halo=halo)
    arch = MultiCPU(grid=full_grid, ranks=(1, 1, 4))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        interior(field) .= arch.local_rank
        fill_halo_regions!(field, arch)

        @test all(top_halo(field, include_corners=false) .== arch.connectivity.top)
        @test all(bottom_halo(field, include_corners=false) .== arch.connectivity.bottom)

        @test all(interior(field) .== arch.local_rank)
        @test all(east_halo(field, include_corners=false) .== arch.local_rank)
        @test all(west_halo(field, include_corners=false) .== arch.local_rank)
        @test all(north_halo(field, include_corners=false) .== arch.local_rank)
        @test all(south_halo(field, include_corners=false) .== arch.local_rank)
    end

    return nothing
end

function test_triply_periodic_halo_communication_with_221_ranks(halo)
    topo = (Periodic, Periodic, Periodic)
    full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 3), extent=(1, 2, 3), halo=halo)
    arch = MultiCPU(grid=full_grid, ranks=(2, 2, 1))
    model = DistributedIncompressibleModel(architecture=arch, grid=full_grid, pressure_solver=nothing)

    for field in merge(fields(model), model.pressures)
        interior(field) .= arch.local_rank
        fill_halo_regions!(field, arch)

        @test all(east_halo(field, include_corners=false) .== arch.connectivity.east)
        @test all(west_halo(field, include_corners=false) .== arch.connectivity.west)
        @test all(north_halo(field, include_corners=false) .== arch.connectivity.north)
        @test all(south_halo(field, include_corners=false) .== arch.connectivity.south)

        @test all(interior(field) .== arch.local_rank)
        @test all(top_halo(field, include_corners=false) .== arch.local_rank)
        @test all(bottom_halo(field, include_corners=false) .== arch.local_rank)
    end

    return nothing
end

#####
##### Run tests!
#####

@testset "Distributed MPI Oceananigans" begin

    @info "Testing distributed MPI Oceananigans..."

    @testset "Multi architectures rank connectivity" begin
        @info "  Testing multi architecture rank connectivity..."
        test_triply_periodic_rank_connectivity_with_411_ranks()
        test_triply_periodic_rank_connectivity_with_141_ranks()
        test_triply_periodic_rank_connectivity_with_114_ranks()
        test_triply_periodic_rank_connectivity_with_221_ranks()
    end

    @testset "Local grids for distributed models" begin
        @info "  Testing local grids for distributed models..."
        test_triply_periodic_local_grid_with_411_ranks()
        test_triply_periodic_local_grid_with_141_ranks()
        test_triply_periodic_local_grid_with_114_ranks()
        test_triply_periodic_local_grid_with_221_ranks()
    end

    @testset "Injection of halo communication BCs" begin
        @info "  Testing injection of halo communication BCs..."
        test_triply_periodic_bc_injection_with_411_ranks()
        test_triply_periodic_bc_injection_with_141_ranks()
        test_triply_periodic_bc_injection_with_114_ranks()
        test_triply_periodic_bc_injection_with_221_ranks()
    end

    @testset "Halo communication" begin
        @info "  Testing halo communication..."
        for H in 1:3
            test_triply_periodic_halo_communication_with_411_ranks((H, H, H))
            test_triply_periodic_halo_communication_with_141_ranks((H, H, H))
            test_triply_periodic_halo_communication_with_114_ranks((H, H, H))
            test_triply_periodic_halo_communication_with_221_ranks((H, H, H))
        end
    end

    @testset "Time stepping IncompressibleModel" begin
        topo = (Periodic, Periodic, Periodic)
        full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 8), extent=(1, 2, 3))
        arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))
        model = DistributedIncompressibleModel(architecture=arch, grid=full_grid)

        time_step!(model, 1)
        @test model isa IncompressibleModel
        @test model.clock.time ≈ 1

        simulation = Simulation(model, Δt=1, stop_iteration=2)
        run!(simulation)
        @test model isa IncompressibleModel
        @test model.clock.time ≈ 2
    end

    @testset "Time stepping ShallowWaterModel" begin
        topo = (Periodic, Periodic, Bounded)
        full_grid = RegularRectilinearGrid(topology=topo, size=(8, 8, 1), extent=(1, 2, 3))
        arch = MultiCPU(grid=full_grid, ranks=(1, 4, 1))
        model = DistributedShallowWaterModel(architecture=arch, grid=full_grid, gravitational_acceleration=1)

        set!(model, h=model.grid.Lz)
        time_step!(model, 1)
        @test model isa ShallowWaterModel
        @test model.clock.time ≈ 1

        simulation = Simulation(model, Δt=1, stop_iteration=2)
        run!(simulation)
        @test model isa ShallowWaterModel
        @test model.clock.time ≈ 2
    end
end


text:
@testset "Points" begin
    @test Point() == Point{1}()
    @test Point() == Point{1, Float64}()
    @test Point() == Point(0.0)
    @test Point{2}() == Point(0.0, 0.0)
    @test Point{2, Int}() == Point(0, 0)
    @test Point(1, 2) == Point((1, 2))

    p = Point(1.0)
    @test zero(p) == Point(0.0)
    @test get_coords(p) == (1.0,)
#    @test typeof(Point(1)) == Point{1, Int64, Tuple{Int64}}
    @test typeof(Point(1)) == Point{1, Int64}
    p = Point(2.0, 3.0)
    @test p[1] == 2.0
    @test p[2] == 3.0
end


text:
function write_output(	file_path::String,
						D_empirical::Float64,
						deltat_coarse::Float64,
						number_of_time_points_coarse::Int64,
						msd::Array{Float64, 1},
						msd_x::Array{Float64, 1},
						msd_y::Array{Float64, 1},
						msd_z::Array{Float64, 1},
						t_exec::Float64)
	file_stream::IOStream = open(file_path, "w")

	t::Array{Float64, 1} = 0.0:deltat_coarse:(convert(Float64, number_of_time_points_coarse-1) * deltat_coarse)

	@printf(file_stream, "%s", "<output_diffusion>\n")

	write_key(file_stream, "D_empirical", D_empirical)
	write_key(file_stream, "t", t)
	write_key(file_stream, "msd", msd)
	write_key(file_stream, "msd_x", msd_x)
	write_key(file_stream, "msd_y", msd_y)
	write_key(file_stream, "msd_z", msd_z)
	write_key(file_stream, "t_exec", t_exec)

	@printf(file_stream, "%s", "</output_diffusion>")

	close(file_stream)

	nothing
end