text:
"""Transit function to call the specified bottom drag function."""
function bottom_drag!(  u::AbstractMatrix,
                        v::AbstractMatrix,
                        η::AbstractMatrix,
                        Diag::DiagnosticVars,
                        S::ModelSetup)

    @unpack bottom_drag = S.parameters

    if bottom_drag == "quadratic"
        bottom_drag_quadratic!(u,v,η,Diag,S)
    elseif bottom_drag == "linear"
        bottom_drag_linear!(u,v,η,Diag,S)
    elseif bottom_drag == "none"
        no_bottom_drag!(u,v,η,Diag,S)
    end
end

"""Quadratic bottom drag Bu,Bv = cD/h * | u⃗ | * u⃗"""
function bottom_drag_quadratic!(u::AbstractMatrix,
                                v::AbstractMatrix,
                                η::AbstractMatrix,
                                Diag::DiagnosticVars,
                                S::ModelSetup)

    @unpack h,h_u,h_v = Diag.VolumeFluxes
    @unpack u²,v²,KEu,KEv = Diag.Bernoulli
    @unpack Bu,Bv,sqrtKE,sqrtKE_u,sqrtKE_v = Diag.Bottomdrag
    @unpack ep = S.grid
    @unpack cD = S.constants
    @unpack H = S.forcing

    thickness!(h,η,H)
    Ix!(h_u,h)
    Iy!(h_v,h)

    speed!(u²,v²,u,v)
    Ix!(KEu,u²)
    Iy!(KEv,v²)

    # sqrt of KE, which is actually the kinetic energy without the 0.5 factor
    m,n = size(sqrtKE)
    @boundscheck (m+ep,n+2) == size(KEu) || throw(BoundsError())
    @boundscheck (m+2,n) == size(KEv) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            sqrtKE[i,j] = sqrt(KEu[i+ep,j+1] + KEv[i+1,j])
        end
    end

    Ix!(sqrtKE_u,sqrtKE)
    Iy!(sqrtKE_v,sqrtKE)

    m,n = size(Bu)
    @boundscheck (m,n) == size(sqrtKE_u) || throw(BoundsError())
    @boundscheck (m,n) == size(h_u) || throw(BoundsError())
    @boundscheck (m+2+ep,n+2) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            Bu[i,j] = cD*sqrtKE_u[i,j] * u[i+1+ep,j+1] / h_u[i,j]
        end
    end

    m,n = size(Bv)
    @boundscheck (m,n) == size(sqrtKE_v) || throw(BoundsError())
    @boundscheck (m,n) == size(h_v) || throw(BoundsError())
    @boundscheck (m+2,n+2) == size(v) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            Bv[i,j] = cD*sqrtKE_v[i,j] * v[i+1,j+1] / h_v[i,j]
        end
    end
end

"""Linear bottom drag computed as r_B*(u,v). r_B is negative and contains the
grid spacing Δ as gradient operators are dimensionless."""
function bottom_drag_linear!(   u::AbstractMatrix,
                                v::AbstractMatrix,
                                η::AbstractMatrix,
                                Diag::DiagnosticVars,
                                S::ModelSetup)

    @unpack Bu,Bv = Diag.Bottomdrag
    @unpack ep = S.grid
    @unpack rD = S.constants

    m,n = size(Bu)
    @boundscheck (m+2+ep,n+2) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            Bu[i,j] = rD * u[i+1+ep,j+1]
        end
    end

    m,n = size(Bv)
    @boundscheck (m+2,n+2) == size(v) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            Bv[i,j] = rD * v[i+1,j+1]
        end
    end
end

"""No bottom drag."""
function no_bottom_drag!(   u::AbstractMatrix,
                            v::AbstractMatrix,
                            η::AbstractMatrix,
                            Diag::DiagnosticVars,
                            S::ModelSetup)

    @unpack Bu,Bv = Diag.Bottomdrag
    Bu .= zero(eltype(Bu))
    Bv .= zero(eltype(Bu))
end


text: const OUTAGE = "outage"

ext: jl


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using Surrogates
using Random
using Optim
using Plots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,γ) = p
    @inbounds begin
        du[1] = -β*S*I
        du[2] = β*S*I - γ*I
        du[3] = γ*I
    end
    nothing
end;

tmax = 10000.0
tspan = (0.0,tmax);

cb_ss = TerminateSteadyState();

n_train = 25 # Number of samples
n_test = 1000; # Number of samples

# Parameters are β
lb = 0.00005
ub = 0.001;

N = 1000.0
u0 = [990.0,10.0,0.0]
p = [0.0005,0.25]
prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

Random.seed!(123)
sampler = LatinHypercubeSample();

β = Surrogates.sample(n_train,lb,ub,sampler)
β₀ = copy(β);

final_size = function(β)
  p = prob_ode.p
  p[1] = β
  prob = remake(prob_ode;p=p)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  fsp = sol[end][3]/N
  fsp
end;

fs = final_size.(β);

β_grid = lb:0.00001:ub
plot(β_grid,final_size.(β_grid),label="True function",xlab="Infectivity parameter, β",ylab="Final size")
scatter!(β,fs,legend=:right,label="Sampled points")

fs_surrogate = RadialBasis(β, fs, lb, ub, rad=cubicRadial);

scatter(fs,
        fs_surrogate.(β),
        legend=false,
        xlabel="True final size",
        ylabel="Surrogate final size")
Plots.abline!(1.0,0.0)

surrogate_optimize(final_size,
                   SRBF(),
                   lb,
                   ub,
                   fs_surrogate,
                   sampler;
                   maxiters=1000,
                   num_new_samples=1000);

β₁ = copy(β)
length(β₀),length(β)

fs_opt = final_size.(β)
scatter(fs_opt,fs_surrogate.(β), xlabel="Final size", ylabel="Surrogate final size",legend=false)
Plots.abline!(1.0, 0.0)

β_test = sample(n_test,lb,ub,UniformSample())
fs_test = final_size.(β_test)
fs_test_pred = fs_surrogate.(β_test);

l = @layout [a b]
pl1 = plot(β_grid,final_size.(β_grid),color=:red,label="Model",xlabel="β",ylabel="Final size",legend=:right)
scatter!(pl1, β_test,fs_test_pred,color=:blue,label="Surrogate")
pl2 = scatter(fs_test,fs_test_pred,color=:blue,legend=false,xlabel="Final size",ylabel="Surrogate final size")
Plots.abline!(pl2,1.0,0.0)
plot(pl1,pl2,layout=l)

@benchmark final_size(0.0005)

@benchmark fs_surrogate(0.0005)

peak_infected = function(β)
  p = prob_ode.p
  p[1] = β
  prob = remake(prob_ode;p=p)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  tss = sol.t[end]
  f = (t) -> -sol(t,idxs=2)
  opt = Optim.optimize(f,0.0,tss)
  pk = -opt.minimum
  pk
end

plot(β_grid, peak_infected,xlabel="β",ylabel="Peak infected",legend=false)

pk = peak_infected.(β)
pk_surrogate = RadialBasis(β,pk,lb,ub,rad = cubicRadial)
surrogate_optimize(peak_infected,
                   SRBF(),
                   lb,
                   ub,
                   pk_surrogate,
                   sampler;
                   maxiters=1000,
                   num_new_samples=1000)
pk_test = peak_infected.(β_test)
pk_test_pred = pk_surrogate.(β_test);

l = @layout [a b]
pl1 = plot(β_grid,peak_infected.(β_grid),color=:red,label="Model",xlabel="β",ylabel="Peak infected",legend=:right)
scatter!(pl1, β_test,pk_test_pred,color=:blue,label="Surrogate")
pl2 = scatter(pk_test,pk_test_pred,color=:blue,legend=false,xlabel="Peak infected",ylabel="Surrogate peak infected")
Plots.abline!(pl2,1.0,0.0)
plot(pl1,pl2,layout=l)

@benchmark peak_infected(0.0005)

@benchmark pk_surrogate(0.0005)

peak_time = function(β)
  p = prob_ode.p
  p[1] = β
  prob = remake(prob_ode;p=p)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  tss = sol.t[end]
  f = (t) -> -sol(t,idxs=2)
  opt = Optim.optimize(f,0.0,tss)
  pt = opt.minimizer
  pt
end

plot(β_grid,peak_time,xlabel="β",ylabel="Peak time",legend=false)

pt = peak_time.(β)
pt_surrogate = RadialBasis(β,pt,lb,ub,rad = cubicRadial)
surrogate_optimize(peak_time,
                   SRBF(),
                   lb,
                   ub,
                   pt_surrogate,
                   sampler;
                   maxiters=1000,
                   num_new_samples=1000)
pt_test = peak_time.(β_test)
pt_test_pred = pt_surrogate.(β_test);

l = @layout [a b]
pl1 = plot(β_grid,peak_time.(β_grid),color=:red,label="Model",xlabel="β",ylabel="Peak time",legend=:right)
scatter!(pl1, β_test,pt_test_pred,color=:blue,label="Surrogate")
pl2 = scatter(pt_test,pt_test_pred,color=:blue,legend=false,xlabel="Peak time",ylabel="Surrogate peak time")
Plots.abline!(pl2,1.0,0.0)
plot(pl1,pl2,layout=l)

@benchmark peak_time(0.0005)

@benchmark pt_surrogate(0.0005)


text:
# # Dry Atmosphere GCM diagnostics
#
# This file computes selected diagnostics for the GCM and outputs them
# on the spherical interpolated diagnostic grid.
#
# Use it by calling `Diagnostics.setup_atmos_default_diagnostics()`.
#
# TODO:
# - enable zonal means and calculation of covariances using those means
#   - ds.T_zm = mean(.*1., ds.T; dims = 3)
#   - ds.u_zm = mean((ds.u); dims = 3 )
#   - v_zm = mean(ds.v; dims = 3)
#   - w_zm = mean(ds.w; dims = 3)
#   - ds.uvcovariance = (ds.u .- ds.u_zm) * (ds.v .- v_zm)
#   - ds.vTcovariance = (ds.v .- v_zm) * (ds.T .- ds.T_zm)
# - add more variables, including horiz streamfunction from laplacial of vorticity (LN)
# - density weighting
# - maybe change thermo/dyn separation to local/nonlocal vars?

using LinearAlgebra
using Printf
using Statistics

using ..Atmos
using ..Atmos: thermo_state, turbulence_tensors

include("diagnostic_fields.jl")

function atmos_gcm_default_init(dgngrp::DiagnosticsGroup, currtime)
    if !(dgngrp.interpol isa InterpolationCubedSphere)
        @warn """
            Diagnostics ($dgngrp.name): currently requires `InterpolationCubedSphere`!
            """
    end
end

# 3D variables
function vars_atmos_gcm_default_simple_3d(atmos::AtmosModel, FT)
    @vars begin
        u::FT
        v::FT
        w::FT
        rho::FT
        temp::FT
        pres::FT
        thd::FT                 # θ_dry
        et::FT                  # e_tot
        ei::FT                  # e_int
        ht::FT
        hi::FT
        vort::FT                # Ω₃

        moisture::vars_atmos_gcm_default_simple_3d(atmos.moisture, FT)
    end
end
vars_atmos_gcm_default_simple_3d(::MoistureModel, FT) = @vars()
function vars_atmos_gcm_default_simple_3d(m::EquilMoist, FT)
    @vars begin
        qt::FT                  # q_tot
        ql::FT                  # q_liq
        qv::FT                  # q_vap
        qi::FT                  # q_ice
        thv::FT                 # θ_vir
        thl::FT                 # θ_liq

    end
end
num_atmos_gcm_default_simple_3d_vars(m, FT) =
    varsize(vars_atmos_gcm_default_simple_3d(m, FT))
atmos_gcm_default_simple_3d_vars(m, array) =
    Vars{vars_atmos_gcm_default_simple_3d(m, eltype(array))}(array)

function atmos_gcm_default_simple_3d_vars!(
    atmos::AtmosModel,
    state_conservative,
    thermo,
    dyni,
    vars,
)
    vars.u = state_conservative.ρu[1] / state_conservative.ρ
    vars.v = state_conservative.ρu[2] / state_conservative.ρ
    vars.w = state_conservative.ρu[3] / state_conservative.ρ
    vars.rho = state_conservative.ρ
    vars.temp = thermo.temp
    vars.pres = thermo.pres
    vars.thd = thermo.θ_dry
    vars.et = state_conservative.ρe / state_conservative.ρ
    vars.ei = thermo.e_int
    vars.ht = thermo.h_tot
    vars.hi = thermo.h_int

    vars.vort = dyni.Ω₃

    atmos_gcm_default_simple_3d_vars!(
        atmos.moisture,
        state_conservative,
        thermo,
        vars,
    )

    return nothing
end
function atmos_gcm_default_simple_3d_vars!(
    ::MoistureModel,
    state_conservative,
    thermo,
    vars,
)
    return nothing
end
function atmos_gcm_default_simple_3d_vars!(
    moist::EquilMoist,
    state_conservative,
    thermo,
    vars,
)
    vars.moisture.qt = state_conservative.moisture.ρq_tot / state_conservative.ρ
    vars.moisture.ql = thermo.moisture.q_liq
    vars.moisture.qv = thermo.moisture.q_vap
    vars.moisture.qi = thermo.moisture.q_ice
    vars.moisture.thv = thermo.moisutre.θ_vir
    vars.moisture.thl = thermo.moisture.θ_liq_ice

    return nothing
end

# Dynamic variables
function vars_dyn(FT)
    @vars begin
        Ω₁::FT
        Ω₂::FT
        Ω₃::FT
    end
end
dyn_vars(array) = Vars{vars_dyn(eltype(array))}(array)

"""
    atmos_gcm_default_collect(bl, currtime)

    Master function that performs a global grid traversal to compute various
    diagnostics using the above functions.
"""
function atmos_gcm_default_collect(dgngrp::DiagnosticsGroup, currtime)
    interpol = dgngrp.interpol
    if !(interpol isa InterpolationCubedSphere)
        @warn """
            Diagnostics ($dgngrp.name): currently requires `InterpolationCubedSphere`!
            """
        return nothing
    end

    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    bl = dg.balance_law
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get needed arrays onto the CPU
    device = array_device(Q)
    if device isa CPU
        ArrayType = Array
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
    else
        ArrayType = CuArray
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
    end
    FT = eltype(state_data)

    # TODO: can this be done in one pass?
    #
    # Non-local vars, e.g. relative vorticity
    vgrad = VectorGradients(dg, Q)
    vort = Vorticity(dg, vgrad)

    # Compute thermo variables
    thermo_array = Array{FT}(undef, npoints, num_thermo(bl, FT), nrealelem)
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state_conservative = extract_state_conservative(dg, state_data, ijk, e)
        state_auxiliary = extract_state_auxiliary(dg, aux_data, ijk, e)

        thermo = thermo_vars(bl, view(thermo_array, ijk, :, e))
        compute_thermo!(bl, state_conservative, state_auxiliary, thermo)
    end

    # Interpolate the state, thermo and dyn vars to sphere (u and vorticity
    # need projection to zonal, merid). All this may happen on the GPU.
    istate =
        ArrayType{FT}(undef, interpol.Npl, number_state_conservative(bl, FT))
    interpolate_local!(interpol, Q.realdata, istate)

    ithermo = ArrayType{FT}(undef, interpol.Npl, num_thermo(bl, FT))
    interpolate_local!(interpol, ArrayType(thermo_array), ithermo)

    idyn = ArrayType{FT}(undef, interpol.Npl, size(vort.data, 2))
    interpolate_local!(interpol, vort.data, idyn)

    # TODO: get indices here without hard-coding them
    _ρu, _ρv, _ρw = 2, 3, 4
    project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))
    _Ω₁, _Ω₂, _Ω₃ = 1, 2, 3
    project_cubed_sphere!(interpol, idyn, (_Ω₁, _Ω₂, _Ω₃))

    # FIXME: accumulating to rank 0 is not scalable
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    all_thermo_data = accumulate_interpolated_data(mpicomm, interpol, ithermo)
    all_dyn_data = accumulate_interpolated_data(mpicomm, interpol, idyn)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # adjust the level dimension for `planet_radius`
        level_val = dims["level"]
        dims["level"] =
            (level_val[1] .- FT(planet_radius(param_set)), level_val[2])

        # set up the array for the diagnostic variables based on the interpolated grid
        nlong = length(dims["long"][1])
        nlat = length(dims["lat"][1])
        nlevel = length(dims["level"][1])

        simple_3d_vars_array = Array{FT}(
            undef,
            1,
            nlong,
            nlat,
            nlevel,
            num_atmos_gcm_default_simple_3d_vars(bl, FT),
        )

        @visitI nlong nlat nlevel begin
            statei = Vars{vars_state_conservative(bl, FT)}(view(
                all_state_data,
                lo,
                la,
                le,
                :,
            ))
            thermoi = thermo_vars(bl, view(all_thermo_data, lo, la, le, :))
            dyni = dyn_vars(view(all_dyn_data, lo, la, le, :))
            simple_3d_vars = atmos_gcm_default_simple_3d_vars(
                bl,
                view(simple_3d_vars_array, 1, lo, la, le, :),
            )

            atmos_gcm_default_simple_3d_vars!(
                bl,
                statei,
                thermoi,
                dyni,
                simple_3d_vars,
            )
        end

        # assemble the diagnostics for writing
        dim_names = tuple("time", collect(keys(dims))...)
        varvals = OrderedDict()
        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_gcm_default_simple_3d(bl, FT)),
        )
        for (vari, varname) in enumerate(varnames)
            var = Variables[varname]
            varvals[varname] = (
                dim_names,
                simple_3d_vars_array[:, :, :, :, vari],
                OrderedDict(
                    "units" => var.units,
                    "long_name" => var.long,
                    "standard_name" => var.standard,
                ),
            )
        end

        # write output
        dprefix = @sprintf(
            "%s_%s_%s_num%04d",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
            dgngrp.num
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        write_data(dgngrp.writer, dfilename, dims, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_gcm_default_fini(dgngrp::DiagnosticsGroup, currtime) end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
md"""
# Evolution Strategies

Evolution strategies are the functions used to evolve the population in the genetic algorithm from one generation to the next. The following is performed by evolution strategies:

* Select which candidates to use for the next generation
* Produce new candidates, e.g by mutating the selected candidates

Important to note about evolution strategies is that they generally expect candidates which can provide a precomputed fitness value,
e.g. [`FittedCandidate`](@ref)s. This is because the fitness value is used by things like sorting where it is not only impractical
to recompute it, but is also might lead to undefined behaviour if it is not always the same. Use [`Population`](@ref) to get some
help with computing fitness for all candidates before passing them on to evolution.

Note that there is no general requirement on an evolution strategy to return the same population size as it was given. It is also
free to create completely new candidates without basing anything on any given candidate.
"""

@testset "Evolution strategies" begin #src
# For controlled randomness in the examples.
struct FakeRng end
Base.rand(::FakeRng) = 0.7

# Dummy candidate for brevity.
struct Cand <: AbstractCandidate
    fitness
end
NaiveGAflux.fitness(d::Cand) = d.fitness

# [`EliteSelection`](@ref) selects the n best candidates.
elitesel = EliteSelection(2)
@test evolve(elitesel, Cand.(1:10)) == Cand.([10, 9])

# [`EvolveCandidates`](@ref) maps candidates to new candidates (e.g. through mutation).
evocands = EvolveCandidates(c -> Cand(fitness(c) + 0.1))
@test evolve(evocands, Cand.(1:10)) == Cand.(1.1:10.1)

# [`SusSelection`](@ref) selects candidates randomly using stochastic uniform sampling.
# Selected candidates will be forwarded to the wrapped evolution strategy before returned.
sussel = SusSelection(5, evocands, FakeRng())
@test evolve(sussel, Cand.(1:10)) == Cand.([4.1, 6.1, 8.1, 9.1, 10.1])

# [`CombinedEvolution`](@ref) combines the populations from several evolution strategies.
comb = CombinedEvolution(elitesel, sussel)
@test evolve(comb, Cand.(1:10)) == Cand.(Any[10, 9, 4.1, 6.1, 8.1, 9.1, 10.1])
end #src


text:
# __precompile__()

module LightDarkPOMDPs

importall POMDPs

using StaticArrays
using Plots
using POMDPToolbox
using Parameters # for @with_kw
using ParticleFilters # for AbstractParticleBelief

export
    AbstractLD2,
    LightDark2D,
    LightDark2DTarget,
    LightDark2DKalman,
    SymmetricNormal2,
    Vec2

include("lightdark2d.jl")
include("lightdark2dtarget.jl")
include("lightdark2dfilter.jl")
include("lightdark2dvis.jl")

end # module


text:
incx = abs(incx)
@inbounds for i = ix:incx:ix-1+n*incx
    x[i] = a
end


text:
module Devito

using LinearAlgebra, MPI, PyCall, Strided

const numpy = PyNULL()
const sympy = PyNULL()
const devito = PyNULL()
const seismic = PyNULL()

function __init__()
    copy!(numpy, pyimport("numpy"))
    copy!(sympy, pyimport("sympy"))
    copy!(devito, pyimport("devito"))
    copy!(seismic, pyimport("examples.seismic"))
end

numpy_eltype(dtype) = dtype == numpy.float32 ? Float32 : Float64

PyCall.PyObject(::Type{Float32}) = numpy.float32
PyCall.PyObject(::Type{Float64}) = numpy.float64

"""
    configuration!(key, value)

Configure Devito.  Examples include
```julia
configuration!("log-level", "DEBUG")
configuration!("language", "openmp")
configuration!("mpi", false)
```
"""
function configuration!(key, value)
    c = PyDict(devito."configuration")
    c[key] = value
    c[key]
end
configuration(key) = PyDict(devito."configuration")[key]
configuration() = PyDict(devito."configuration")

_reverse(argument::Tuple) = reverse(argument)
_reverse(argument) = argument

function reversedims(arguments)
    _arguments = collect(arguments)
    keys = first.(_arguments)
    values = @. _reverse(last(_arguments))
    (; zip(keys, values)...)
 end

struct DevitoArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    o::PyObject # Python object for the numpy array
    p::A # copy-free
end

function DevitoArray{T,N}(o) where {T,N}
    p = unsafe_wrap(Array{T,N}, Ptr{T}(o.__array_interface__["data"][1]), reverse(o.shape); own=false)
    DevitoArray{T,N,Array{T,N}}(o, p)
end

function DevitoArray(o)
    T = numpy_eltype(o.dtype)
    N = length(o.shape)
    DevitoArray{T,N}(o)
end

Base.size(x::DevitoArray{T,N}) where {T,N} = size(x.p)
Base.parent(x::DevitoArray) = x.p

Base.getindex(x::DevitoArray{T,N,A}, i) where {T,N,A<:Array} = getindex(parent(x), i)
Base.getindex(x::DevitoArray{T,N,A}, I::Vararg{Int,N}) where {T,N,A<:StridedView} = getindex(parent(x), I...)
Base.setindex!(x::DevitoArray{T,N,A}, v, i) where {T,N,A<:Array} = setindex!(parent(x), v, i)
Base.setindex!(x::DevitoArray{T,N,A}, v, I::Vararg{Int,N}) where {T,N,A<:StridedView} = setindex!(parent(x), v, I...)
Base.IndexStyle(::Type{<:DevitoArray{<:Any,<:Any,<:Array}}) = IndexLinear()
Base.IndexStyle(::Type{<:DevitoArray{<:Any,<:Any,<:StridedView}}) = IndexCartesian()

Base.view(x::DevitoArray{T,N,Array{T,N}}, I::Vararg{Any}) where {T,N} = DevitoArray(x.o, sview(x.p, I...))

struct DevitoMPIArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    o::PyObject
    p::A
    local_indices::NTuple{N,UnitRange{Int}}
end

function DevitoMPIArray{T,N}(o, idxs) where {T,N}
    p = unsafe_wrap(Array{T,N}, Ptr{T}(o.__array_interface__["data"][1]), length.(idxs); own=false)
    DevitoMPIArray{T,N,Array{T,N}}(o, p, idxs)
end

function Base.size(x::DevitoMPIArray{T,N}) where {T,N}
    MPI.Initialized() || MPI.Init()
    n = ntuple(i->( localsize(x)[i] > 0 ? x.local_indices[i][end] : 0), N)
    MPI.Allreduce(n, max, MPI.COMM_WORLD)
end

Base.parent(x::DevitoMPIArray) = x.p

localsize(x::DevitoMPIArray{T,N}) where {T,N} = ntuple(i->size(x.local_indices[i])[1], N)

localindices(x::DevitoMPIArray{T,N}) where {T,N} = x.local_indices

function Base.convert(::Type{Array}, x::DevitoMPIArray{T}) where {T}
    MPI.Initialized() || MPI.Init()
    y = zeros(T, size(x))
    y[localindices(x)...] .= parent(x)
    MPI.Reduce!(y, +, 0, MPI.COMM_WORLD)
    y
end

function Base.copyto!(dst::DevitoMPIArray, src::AbstractArray)
    MPI.Initialized() || MPI.Init()
    MPI.Bcast!(src, 0, MPI.COMM_WORLD)
    parent(dst) .= src[localindices(dst)...]
    MPI.Barrier(MPI.COMM_WORLD)
    dst
end

function Base.fill!(x::DevitoMPIArray, v)
    MPI.Initialized() || MPI.Init()
    parent(x) .= v
    MPI.Barrier(MPI.COMM_WORLD)
    x
end

function Base.getindex(x::DevitoMPIArray{T,N}, I::Vararg{Int,N}) where {T,N}
    if all(ntuple(idim->I[idim] ∈ x.local_indices[idim], N))
        J = ntuple(idim->I[idim]-x.local_indices[idim]+1, N)
        v = getindex(x.p, J...)
    end
    v
end
Base.setindex!(x::DevitoMPIArray{T,N}, v, i) where {T,N} = @warn "not implemented"
Base.IndexStyle(::Type{<:DevitoMPIArray}) = IndexCartesian()

# TODO -- need to implement broadcasting interface for DevitoMPIArray

struct DevitoMPISparseArray{T,N,NM1} <: AbstractArray{T,N}
    o::PyObject
    p::Array
    local_indices::Array{Int,NM1}
end

function DevitoMPISparseArray{T,N}(o, idxs) where {T,N}
    local p
    if length(idxs) == 0
        p = Array{T,N}(undef, ntuple(_->0, N))
    else
        p = unsafe_wrap(Array{T,N}, Ptr{T}(o.__array_interface__["data"][1]), reverse(o.shape); own=false)
    end
    DevitoMPISparseArray{T,N,N-1}(o, p, idxs)
end
Base.IndexStyle(::Type{<:DevitoMPISparseArray}) = IndexCartesian()
Base.getindex(x::DevitoMPISparseArray{T,N}, I::Vararg{Int,N}) where {T,N} = error("not implemented")
Base.setindex!(x::DevitoMPISparseArray{T,N}, v, I::Vararg{Int,N}) where {T,N} = error("not implemented")

Base.parent(x::DevitoMPISparseArray) = x.p

function Base.size(x::DevitoMPISparseArray{T,N}) where {T,N}
    MPI.Initialized() || MPI.Init()
    n = MPI.Allreduce(x.o.shape[2], +, MPI.COMM_WORLD)
    (n,x.o.shape[1])
end

function Base.convert(::Type{Array}, x::DevitoMPISparseArray{T,N}) where {T,N}
    MPI.Initialized() || MPI.Init()
    y = zeros(T, size(x))
    _x = parent(x)
    n = x.o.shape[1]
    for j = 1:n
        for (i,idx) in enumerate(x.local_indices)
            y[idx,j] = _x[i,j]
        end
    end
    MPI.Reduce!(y, +, 0, MPI.COMM_WORLD)
    y
end
Base.convert(::Type{Array}, x::Transpose{T, <:DevitoMPISparseArray{T,2}}) where {T} = (convert(Array, parent(x)))'

function Base.copyto!(dst::DevitoMPISparseArray, src::Array, transposeflag = false)
    MPI.Initialized() || MPI.Init()
    MPI.Bcast!(src, 0, MPI.COMM_WORLD)
    n = transposeflag ? size(src, 1) : size(src, 2)
    dst_parent = parent(dst)
    for (i,idx) in enumerate(dst.local_indices)
        for j = 1:n
            dst_parent[i,j] = transposeflag ? src[j,idx] : src[idx,j]
        end
    end
    MPI.Barrier(MPI.COMM_WORLD)
    dst
end
Base.copyto!(dst::Transpose{T, <:DevitoMPISparseArray{T,2}}, src::Matrix) where {T} = copyto!(parent(dst), src, true)

#
# Dimension
#
abstract type AbstractDimension end

# here is Devito's dimension type hierarchy:
# https://github.com/devitocodes/devito/blob/02bbefb7e380d299a2508fef2923c1a4fbd5c59d/devito/types/dimension.py

# Python <-> Julia quick-and-dirty type/struct for dimensions
for (M,F) in ((:devito,:SpaceDimension), (:devito,:SteppingDimension), (:devito,:TimeDimension), (:devito,:ConditionalDimension),(:devito,:Dimension),(:devito,:DefaultDimension))
    @eval begin
        struct $F <: AbstractDimension
            o::PyObject
        end
        PyCall.PyObject(x::$F) = x.o
        Base.convert(::Type{$F}, x::PyObject) = $F(x)
        $F(args...; kwargs...) = pycall($M.$F, $F, args...; kwargs...)
        export $F
    end
end

# Python <-> Julia quick-and-dirty type/struct mappings
for (M,F) in ((:devito,:Constant), (:devito,:Eq), (:devito,:Injection), (:devito,:Operator))
    @eval begin
        struct $F
            o::PyObject
        end
        PyCall.PyObject(x::$F) = x.o
        Base.convert(::Type{$F}, x::PyObject) = $F(x)
        $F(args...; kwargs...) = pycall($M.$F, $F, args...; kwargs...)
        export $F
    end
end

function SpaceDimension end
"""
    SpaceDimension(;kwargs...)

Construct a space dimension that defines the extend of a physical grid.

See https://www.devitoproject.org/devito/dimension.html?highlight=spacedimension#devito.types.dimension.SpaceDimension.

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
````
"""

Base.:(==)(x::AbstractDimension,y::AbstractDimension) = x.o == y.o

function Operator end
"""
    Opertor(expressions; kwargs...)

Generate, JIT-compile and run C code starting from an ordered sequence of symbolic expressions.

See: https://www.devitoproject.org/devito/operator.html?highlight=operator#devito.operator.operator.Operator""

# Example
Assuming that one has constructed the following Devito expressions: `stencil_p`, `src_term` and `rec_term`,
```julia
op = Operator([stencil_p, src_term, rec_term]; name="opIso")
```
"""

function ConditionalDimension end
"""
    ConditionalDimension(;kwargs)

Symbol defining a non-convex iteration sub-space derived from a parent Dimension, implemented by the compiler generating conditional “if-then” code within the parent Dimension’s iteration space.

See: https://www.devitoproject.org/devito/dimension.html?highlight=conditional#devito.types.dimension.ConditionalDimension

# Example
```julia
size, factor = 16, 4
i  = SpaceDimension(name="i")
grid = Grid(shape=(size,),dimensions=(i,))
ci = ConditionalDimension(name="ci", parent=i, factor=factor)
g  = Devito.Function(name="g", grid=grid, shape=(size,), dimensions=(i,))
f  = Devito.Function(name="f", grid=grid, shape=(div(size,factor),), dimensions=(ci,))
op = Operator([Eq(g, 1), Eq(f, g)],name="Cond")
```
"""

factor(x::ConditionalDimension) = x.o.factor
export factor
Base.parent(x::Union{ConditionalDimension,SteppingDimension}) = x.o.parent

#
# Grid
#
struct Grid{T,N}
    o::PyObject
end

"""
    Grid(; shape[, optional key-word arguments...])

Construct a grid that can be used in the construction of a Devito Functions.

See: https://www.devitoproject.org/devito/grid.html?highlight=grid#devito.types.Grid

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (x,z), # z is fast (row-major)
    shape = (251,501),
    origin = (0.0,0.0),
    extent = (1250.0,2500.0),
    dtype = Float32)
```
"""
function Grid(args...; kwargs...)
    o = pycall(devito.Grid, PyObject, args...; reversedims(kwargs)...)
    T = numpy_eltype(o.dtype)
    N = length(o.shape)
    Grid{T,N}(o)
end

PyCall.PyObject(x::Grid) = x.o

Base.:(==)(x::Grid{T,N},y::Grid{T,N}) where{T,N} = x.o == y.o
Base.size(grid::Grid{T,N}) where {T,N} = reverse((grid.o.shape)::NTuple{N,Int})
extent(grid::Grid{T,N}) where {T,N} = reverse((grid.o.extent)::NTuple{N,Float64})
"""
    origin(grid)
returns the tuple corresponding to the grid's origin
"""
origin(grid::Grid{T,N}) where {T,N} = reverse((grid.o.origin)::NTuple{N,Float64})
size_with_halo(grid::Grid{T,N}, h) where {T,N} = ntuple(i->grid.o.shape[N-i+1] + h[i][1] + h[i][2], N)
Base.size(grid::Grid, i) = size(grid)[i]
Base.ndims(grid::Grid{T,N}) where {T,N} = N
Base.eltype(grid::Grid{T}) where {T} = T

spacing(x::Grid{T,N}) where {T,N} = reverse(x.o.spacing)
spacing_map(x::Grid) = PyDict(x.o."spacing_map")

#
# SubDomain
#

struct SubDomain{N}
    o::PyObject
end

"""
    name(subdomain)

returns the name of the subdomain
"""
function name(x::SubDomain)
    x.o.name
end

PyCall.PyObject(x::SubDomain) = x.o

"""
    subdomains(grid)

returns subdomains associated with a Devito grid
"""
function subdomains(x::Grid{T,N}) where {T,N}
    dictpre =  PyDict(x.o."subdomains")
    dict = Dict()
    for key in keys(dictpre)
        dict[key] = SubDomain{N}(dictpre[key])
    end
    return dict
end

#
# Functions
#
abstract type DevitoMPI end
struct DevitoMPITrue <: DevitoMPI end
struct DevitoMPIFalse <: DevitoMPI end

abstract type DiscreteFunction{T,N,M} end

struct Function{T,N,M} <: DiscreteFunction{T,N,M}
    o::PyObject
end

ismpi_distributed(o::PyObject) = o._distributor.nprocs == 1 ? DevitoMPIFalse : DevitoMPITrue  # TODO - when should should o._distributed == None ??

"""
    Devito.Function(; kwargs...)

Tensor symbol representing a discrete function in symbolic equations.

See: https://www.devitoproject.org/devito/function.html?highlight=function#devito.types.Function

# Example
```
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (x,z),
    shape = (251,501), # assume x is first, z is second (i.e. z is fast in python)
    origin = (0.0,0.0),
    extent = (1250.0,2500.0),
    dtype = Float32)

b = Devito.Function(name="b", grid=grid, space_order=8)
```
"""
function Function(args...; kwargs...)
    o = pycall(devito.Function, PyObject, args...; kwargs...)
    T = numpy_eltype(o.dtype)
    N = length(o.shape)
    M = ismpi_distributed(o)
    Function{T,N,M}(o)
end

struct SubFunction{T,N,M} <: DiscreteFunction{T,N,M}
    o::PyObject
end

struct TimeFunction{T,N,M} <: DiscreteFunction{T,N,M}
    o::PyObject
end

"""
    TimeFunction(; kwargs...)

Tensor symbol representing a discrete function in symbolic equations.

See https://www.devitoproject.org/devito/timefunction.html?highlight=timefunction#devito.types.TimeFunction.

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (x,z),
    shape = (251,501), # assume x is first, z is second (i.e. z is fast in python)
    origin = (0.0,0.0),
    extent = (1250.0,2500.0),
    dtype = Float32)

p = TimeFunction(name="p", grid=grid, time_order=2, space_order=8)
```
"""
function TimeFunction(args...; kwargs...)
    o = pycall(devito.TimeFunction, PyObject, args...; kwargs...)
    T = numpy_eltype(o.dtype)
    N = length(o.shape)
    M = ismpi_distributed(o)
    TimeFunction{T,N,M}(o)
end

struct SparseTimeFunction{T,N,M} <: DiscreteFunction{T,N,M}
    o::PyObject
end

"""
    SparseTimeFunction(; kwargs...)

Tensor symbol representing a space- and time-varying sparse array in symbolic equations.

See: https://www.devitoproject.org/devito/sparsetimefunction.html?highlight=sparsetimefunction

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (x,z),
    shape = (251,501), # assume x is first, z is second (i.e. z is fast in python)
    origin = (0.0,0.0),
    extent = (1250.0,2500.0),
    dtype = Float32)

time_range = 0.0f0:0.5f0:1000.0f0
src = SparseTimeFunction(name="src", grid=grid, npoint=1, nt=length(time_range))
```
"""
function SparseTimeFunction(args...; kwargs...)
    o = pycall(devito.SparseTimeFunction, PyObject, args...; kwargs...)
    T = numpy_eltype(o.dtype)
    N = length(o.shape)
    M = ismpi_distributed(o)
    SparseTimeFunction{T,N,M}(o)
end

PyCall.PyObject(x::DiscreteFunction) = x.o

"""
    grid(f::DiscreteFunction)

Return the grid corresponding to the discrete function `f`.
"""
grid(x::Function{T,N}) where {T,N} = Grid{T,N}(x.o.grid)
grid(x::TimeFunction{T,N}) where {T,N} = Grid{T,N-1}(x.o.grid)

function grid(x::SparseTimeFunction{T}) where {T}
    N = length(x.o.grid.shape)
    Grid{T,N}(x.o.grid)
end

"""
    halo(x::DiscreteFunction)

Return the Devito "outer" halo size corresponding to the discrete function `f`.
"""
halo(x::DiscreteFunction{T,N}) where {T,N} = reverse(x.o.halo)::NTuple{N,Tuple{Int,Int}}

"""
    inhalo(x::DiscreteFunction)

Return the Devito "inner" halo size used for domain decomposition, and corresponding to
the discrete function `f`.
"""
inhalo(x::DiscreteFunction{T,N}) where {T,N} = reverse(x.o._size_inhalo)::NTuple{N,Tuple{Int,Int}}

"""
    size(x::DiscreteFunction)

Return the shape of the grid for the discrete function `x`.
"""
Base.size(x::DiscreteFunction{T,N}) where {T,N} = reverse(x.o.shape)::NTuple{N,Int}

"""
    ndims(x::DiscreteFunction)

Return the number of dimensions corresponding to the discrete function `x`.
"""
Base.ndims(x::DiscreteFunction{T,N}) where {T,N} = N

"""
    size_with_halo(x::DiscreteFunction)

Return the size of the grid associated with `x`, inclusive of the Devito "outer" halo.
"""
size_with_halo(x::DiscreteFunction{T,N}) where{T,N} = reverse(convert.(Int, x.o.shape_with_halo))::NTuple{N,Int}

"""
    size_with_inhalo(x::DiscreteFunction)

Return the size of the grid associated with `z`, inclusive the the Devito "inner" and "outer" halos.
"""
size_with_inhalo(x::DiscreteFunction{T,N}) where {T,N} = reverse(x.o._shape_with_inhalo)::NTuple{N,Int}

function Base.size(x::SparseTimeFunction{T,N,DevitoMPITrue}) where {T,N}
    MPI.Initialized() || MPI.Init()

    _shape = zeros(Int, N)
    if x.o._decomposition[1] == nothing
        if reduce(*, x.o.shape) != 0
            map(i->_shape[i] = x.o.shape[N-i+1], 1:N)
            _shape[1] = x.o.shape[N]
        end
        MPI.Reduce!(_shape, +, 0, MPI.COMM_WORLD)
    else
        error("not implemented")
    end
    ntuple(i->_shape[i], N)
end

size_with_halo(x::SparseTimeFunction) = size(x)

localmask(x::DiscreteFunction{T,N}) where {T,N} = ntuple(i->convert(Int,x.o._mask_domain[N-i+1].start)+1:convert(Int,x.o._mask_domain[N-i+1].stop), N)::NTuple{N,UnitRange{Int}}
localmask_with_halo(x::DiscreteFunction{T,N}) where {T,N} = ntuple(i->convert(Int,x.o._mask_outhalo[N-i+1].start)+1:convert(Int,x.o._mask_outhalo[N-i+1].stop), N)::NTuple{N,UnitRange{Int}}
localmask_with_inhalo(x::DiscreteFunction{T,N}) where {T,N} = ntuple(i->convert(Int,x.o._mask_inhalo[N-i+1].start)+1:convert(Int,x.o._mask_inhalo[N-i+1].stop), N)::NTuple{N,UnitRange{Int}}

localindices(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = localmask(x)
localindices_with_halo(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = localmask_with_halo(x)
localindices_with_inhalo(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = localmask_with_inhalo(x)

"""
    forward(x::TimeFunction)

Returns the symbol for the time-forward state of the `TimeFunction`.

See: https://www.devitoproject.org/devito/timefunction.html?highlight=forward#devito.types.TimeFunction.forward
"""
forward(x::TimeFunction) = x.o.forward

"""
    backward(x::TimeFunction)

Returns the symbol for the time-backward state of the `TimeFunction`.

See: https://www.devitoproject.org/devito/timefunction.html?highlight=forward#devito.types.TimeFunction.backward
"""
backward(x::TimeFunction) = x.o.backward

"""
    time_dim(x::Union{Grid,TimeFunction})

Returns the time dimension for the associated object.
"""
time_dim(x::Union{Grid,TimeFunction}) = dimension(x.o.time_dim)

"""
    stepping_dim(x::Grid)

Returns the stepping dimension for the associated grid.
"""
stepping_dim(x::Grid) = dimension(x.o.stepping_dim)
export time_dim, stepping_dim

"""
    data(x::DiscreteFunction)

Return the data associated with the grid that corresponds to the discrete function `x`.  This is the
portion of the grid that excludes the halo.  In the case of non-MPI Devito, this returns an array
of type `DevitoArray`.  In the case of the MPI Devito, this returns an array of type `DevitoMPIArray`.

The `data` can be converted to an `Array` via `convert(Array, data(x))`.  In the case where `data(x)::DevitoMPIArray`,
this also *collects* the data onto MPI rank 0.
"""
data(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = view(DevitoArray{T,N}(x.o."_data_allocated"), localindices(x)...)

"""
    data_with_halo(x::DiscreteFunction)

Return the data associated with the grid that corresponds to the discrete function `x`.  This is the
portion of the grid that excludes the inner halo and includes the outer halo.  In the case of non-MPI
Devito, this returns an array of type `DevitoArray`.  In the case of the MPI Devito, this returns an
array of type `DevitoMPIArray`.

The `data` can be converted to an `Array` via `convert(Array, data(x))`.  In the case where `data(x)::DevitoMPIArray`,
this also *collects* the data onto MPI rank 0.
"""
data_with_halo(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = view(DevitoArray{T,N}(x.o."_data_allocated"), localindices_with_halo(x)...)

"""
    data_with_inhalo(x::DiscreteFunction)

Return the data associated with the grid that corresponds to the discrete function `x`.  This is the
portion of the grid that includes the inner halo and includes the outer halo.  In the case of non-MPI
Devito, this returns an array of type `DevitoArray`.  In the case of the MPI Devito, this returns an
array of type `DevitoMPIArray`.

The `data` can be converted to an `Array` via `convert(Array, data(x))`.  In the case where `data(x)::DevitoMPIArray`,
this also *collects* the data onto MPI rank 0.
"""
data_with_inhalo(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = view(DevitoArray{T,N}(x.o."_data_allocated"), localindices_with_inhalo(x)...)

"""
    data_allocated(x::DiscreteFunction)

Return the data associated with the grid that corresponds to the discrete function `x`.  This is the
portion of the grid that includes the inner halo and includes the outer halo.  We expect this to be
equivalent to `data_with_inhalo`.

The `data` can be converted to an `Array` via `convert(Array, data(x))`.  In the case where `data(x)::DevitoMPIArray`,
this also *collects* the data onto MPI rank 0.
"""
data_allocated(x::DiscreteFunction{T,N,DevitoMPIFalse}) where {T,N} = DevitoArray{T,N}(x.o."_data_allocated")

function localindices(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    localinds = PyCall.trygetproperty(x.o,"local_indices",nothing)
    if localinds === nothing
        return ntuple(i -> 0:-1, N)
    else
        return ntuple(i->convert(Int,localinds[N-i+1].start)+1:convert(Int,localinds[N-i+1].stop), N)
    end
end

topology(x::DiscreteFunction) = reverse(x.o._distributor.topology)
mycoords(x::DiscreteFunction) = reverse(x.o._distributor.mycoords) .+ 1
decomposition(x::DiscreteFunction) = reverse(x.o._decomposition)

function localindices_with_inhalo(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    h = inhalo(x)
    localidxs = localindices(x)
    n = size_with_inhalo(x)

    _mycoords = mycoords(x)
    _decomposition = decomposition(x)

    ntuple(idim->begin
            local strt,stop
            if _decomposition[idim] == nothing
                strt = 1
                stop = n[idim]
            else
                strt = localidxs[idim][1] + (_mycoords[idim]-1)*(h[idim][1] + h[idim][2])
                stop = strt + length(localidxs[idim]) - 1 + h[idim][1] + h[idim][2]
            end
            strt:stop
        end, N)
end

function localindices_with_halo(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    h = halo(x)
    localidxs = localindices(x)
    n = size_with_halo(x)

    _mycoords = mycoords(x)
    _topology = topology(x)
    _decomposition = decomposition(x)

    ntuple(idim->begin
            local strt,stop
            if _decomposition[idim] == nothing
                strt = 1
                stop = n[idim]
            else
                strt = _mycoords[idim] == 1 ? localidxs[idim][1] : localidxs[idim][1] + h[idim][1]
                stop = _mycoords[idim] == _topology[idim] ? localidxs[idim][end] + h[idim][1] + h[idim][2] : localidxs[idim][end] + h[idim][1]
            end
            strt:stop
        end, N)
end

function data(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    p = sview(parent(data_allocated(x)), localmask(x)...)
    DevitoMPIArray{T,N,typeof(p)}(x.o."_data_allocated", p, localindices(x))
end

function data_with_halo(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    p = sview(parent(data_allocated(x)), localmask_with_halo(x)...)
    DevitoMPIArray{T,N,typeof(p)}(x.o."_data_allocated", p, localindices_with_halo(x))
end

function data_with_inhalo(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    p = sview(parent(data_allocated(x)), localmask_with_inhalo(x)...)
    DevitoMPIArray{T,N,typeof(p)}(x.o."_data_allocated", p, localindices_with_inhalo(x))
end

function data_allocated(x::DiscreteFunction{T,N,DevitoMPITrue}) where {T,N}
    DevitoMPIArray{T,N}(x.o."_data_allocated", localindices_with_inhalo(x))
end

# TODO - needed? <--
function data_with_inhalo(x::SparseTimeFunction{T,N,DevitoMPITrue}) where {T,N}
    MPI.Initialized() || MPI.Init()
    rnk = MPI.Comm_rank(MPI.COMM_WORLD)
    x.o._decomposition[1] == nothing || error("Sam does not know what he is doing!")
    idxs = x.o._decomposition[2][rnk+1] .+ 1
    d = DevitoMPISparseArray{T,N}(x.o."_data_allocated", idxs)
    MPI.Barrier(MPI.COMM_WORLD)
    transpose(d)
end

function data_with_inhalo(x::SparseTimeFunction{T,N,DevitoMPIFalse}) where {T,N}
    d = DevitoArray{T,N}(x.o."_data_allocated")
    transpose(d)
end

data_with_halo(x::SparseTimeFunction{T,N,M}) where {T,N,M} = data_with_inhalo(x)
data(x::SparseTimeFunction{T,N,M}) where {T,N,M} = data_with_inhalo(x)
# -->

function coordinates(x::SparseTimeFunction{T,N,DevitoMPIFalse}) where {T,N}
    n = ndims(grid(x))
    view(DevitoArray{T,N}(x.o.coordinates."_data_allocated"), n:-1:1, :)
end
#coordinates(x::SparseTimeFunction{T,N,DevitoMPITrue}) where {T,N} = DevitoMPIArray{T,N}(x.o.coordinates."_data_allocated", localindices(SubFunction{T,N,DevitoMPITrue}(x.o.coordinates)))
coordinates(x::SparseTimeFunction{T,N,DevitoMPITrue}) where {T,N} = DevitoMPIArray{T,N}(x.o.coordinates."_data_allocated", localindices(SubFunction{T,N,DevitoMPITrue}(x.o.coordinates)))

export DevitoArray, localindices, SubFunction
function dimension(o::PyObject)
    if o.is_Conditional
        return ConditionalDimension(o)
    elseif o.is_Stepping
        return SteppingDimension(o)
    elseif o.is_Space
        return SpaceDimension(o)
    elseif o.is_Time
        return TimeDimension(o)
    elseif o.is_Default
        return DefaultDimension(o)
    elseif o.is_Dimension
        return Dimension(o)
    else
        error("not implemented")
    end
end

"""
    dimensions(x::Union{Grid,DiscreteFunction})

Returns a tuple with the dimensions associated with the Devito grid.
"""
function dimensions(x::Union{Grid{T,N},DiscreteFunction{T,N},SubDomain{N}}) where {T,N}
    ntuple(i->dimension(x.o.dimensions[N-i+1]), N)
end

"""
    inject(x::SparseTimeFunction; kwargs...)

Generate equations injecting an arbitrary expression into a field.

See: Generate equations injecting an arbitrary expression into a field.

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (x,z),
    shape = (251,501), # assume x is first, z is second (i.e. z is fast in python)
    origin = (0.0,0.0),
    extent = (1250.0,2500.0),
    dtype = Float32)

time_range = 0.0f0:0.5f0:1000.0f0
src = SparseTimeFunction(name="src", grid=grid, npoint=1, nt=length(time_range))
src_term = inject(src; field=forward(p), expr=2*src)
```
"""
inject(x::SparseTimeFunction, args...; kwargs...) = pycall(PyObject(x).inject, Injection, args...; kwargs...)

"""
    interpolate(x::SparseTimeFunction; kwargs...)

Generate equations interpolating an arbitrary expression into self.

See: https://www.devitoproject.org/devito/sparsetimefunction.html#devito.types.SparseTimeFunction.interpolate

# Example
```julia
x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))
grid = Grid(
    dimensions = (z,x),
    shape = (501,251), # assume z is first, x is second
    origin = (0.0,0.0),
    extent = (2500.0,1250.0),
    dtype = Float32)

p = TimeFunction(name="p", grid=grid, time_order=2, space_order=8)

time_range = 0.0f0:0.5f0:1000.0f0
nz,nx,δz,δx = size(grid)...,spacing(grid)...
rec = SparseTimeFunction(name="rec", grid=grid, npoint=nx, nt=length(time_range))
rec_coords = coordinates(rec)
rec_coords[1,:] .= 10.0
rec_coords[2,:] .= δx*(0:nx-1)

rec_term = interpolate(rec, expr=p)
```
"""
interpolate(x::SparseTimeFunction; kwargs...) = pycall(PyObject(x).interpolate, PyObject; kwargs...)

"""
apply(    operator::Operator; kwargs...)

Execute the Devito operator, `Operator`.

See: https://www.devitoproject.org/devito/operator.html?highlight=apply#devito.operator.operator.Operator.apply

Note that this returns a `summary::Dict` of the action of applying the operator.  This contains information
such as the number of floating point operations executed per second.
"""
function apply(x::Operator, args...; kwargs...)
    _summary = pycall(PyObject(x).apply, PyObject, args...; kwargs...)

    summary = Dict()
    for (k,v) in _summary.items()
        summary[k] = Dict(
            "time"=>v[1],
            "gflopss"=>v[2],
            "gpointss"=>v[3],
            "oi"=>v[4],
            "ops"=>v[5],
            "itershape"=>v[6])
    end

    summary["globals"] = Dict()
    if haskey(_summary.globals, "fdlike")
        summary["globals"]["fdlike"] = Dict(
            "time"=>_summary.globals["fdlike"][1],
            "gflopss"=>_summary.globals["fdlike"][2],
            "gpointss"=>_summary.globals["fdlike"][3],
            "oi"=>_summary.globals["fdlike"][4],
            "ops"=>_summary.globals["fdlike"][5],
            "itershape"=>_summary.globals["fdlike"][6])
    end

    if haskey(_summary.globals, "vanilla")
        summary["globals"]["vanilla"] = Dict(
            "time"=>_summary.globals["vanilla"][1],
            "gflopss"=>_summary.globals["vanilla"][2],
            "gpointss"=>_summary.globals["vanilla"][3],
            "oi"=>_summary.globals["vanilla"][4],
            "ops"=>_summary.globals["vanilla"][5],
            "itershape"=>_summary.globals["vanilla"][6])
    end
    summary
end

# metaprograming for various derivatives
for F in (:dx,:dy,:dz,:dxr,:dyr,:dzr,:dxl,:dyl,:dzl)
    @eval begin
        $F(x::Union{DiscreteFunction,PyObject}, args...; kwargs...) = pycall(PyObject(x).$F, PyObject, args...; kwargs...)
        export $F
    end
end
"""
    dx(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first derivative with respect to x.
"""
function dx end

"""
    dy(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first derivative with respect to y.
"""
function dy end

"""
    dz(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first derivative with respect to z.
"""

function dz end

"""
    dxl(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first backward one-sided derivative with respect to x.
"""
function dxl end

"""
    dyl(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first backward one-sided derivative with respect to y.
"""
function dyl end

"""
    dzl(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first backward one-sided derivative with respect to z.
"""
function dzl end

"""
    dxr(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first forward one-sided derivative with respect to x.
"""
function dxr end

"""
    dyr(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first forward one-sided derivative with respect to y.
"""
function dyr end

"""
    dzr(f::DiscreteFunction, args...; kwargs...)

Returns the symbol for the first forward one-sided derivative with respect to z.
"""
function dz end

# metaprograming for various derivatives
for F in (:dt,:dt2)
    @eval begin
        $F(x::Union{TimeFunction,PyObject}, args...; kwargs...) = pycall(PyObject(x).$F, PyObject, args...; kwargs...)
        export $F
    end
end

"""
    dt(f::TimeFunction, args...; kwargs...)
Returns the symbol for the first time derivative of a time function
"""
function dt end

"""
    dt2(f::TimeFunction, args...; kwargs...)
Returns the symbol for the second time derivative of a time function
"""
function dt2 end

# metaprogramming for basic operations
for F in ( :+, :-, :*, :/, :^)
    @eval begin
        Base.$F(x::Real,y::DiscreteFunction) = $F(PyObject(x),PyObject(y))
        Base.$F(x::DiscreteFunction, y::DiscreteFunction) = $F(PyObject(x),PyObject(y))
        Base.$F(x::DiscreteFunction, y::PyObject) = $F(x.o,y)
        Base.$F(x::PyObject, y::DiscreteFunction) = $F(x,y.o)
        Base.$F(x::DiscreteFunction, y::Real) = $F(PyObject(x),PyObject(y))
        Base.$F(x::AbstractDimension, y::Real) = $F(PyObject(x),PyObject(y))
        Base.$F(x::Real, y::AbstractDimension) = $F(PyObject(x),PyObject(y))
        Base.$F(x::AbstractDimension, y::DiscreteFunction) = $F(PyObject(x),PyObject(y))
        Base.$F(x::DiscreteFunction, y::AbstractDimension) = $F(PyObject(x),PyObject(y))
        Base.$F(x::AbstractDimension, y::PyObject) = $F(PyObject(x),y)
        Base.$F(x::PyObject, y::AbstractDimension) = $F(x,PyObject(y))
    end
end

Base.:(-)(x::Union{AbstractDimension,DiscreteFunction,PyObject}) = -1*x

# metaprogramming to access Devito dimension boolean attributes
for F in (:is_Dimension, :is_Space, :is_Time, :is_Default, :is_Custom, :is_Derived, :is_NonlinearDerived, :is_Sub, :is_Conditional, :is_Stepping, :is_Modulo, :is_Incr)
    @eval begin
        $F(x::AbstractDimension) = x.o.$F::Bool
        export $F
    end
end
# metaprogramming for devito conditionals
for (M,F) in ((:devito,:Ne),(:devito,:Gt),(:devito,:Ge),(:devito,:Lt),(:devito,:Le))
    @eval begin
        $F(x::Union{Real,DiscreteFunction,PyObject,AbstractDimension},y::Union{Real,DiscreteFunction,PyObject,AbstractDimension}) = $M.$F(PyObject(x),PyObject(y))
        export $F
    end
end

# metaprogramming for symbolic operations on Devito dimensions
for F in (:symbolic_min, :symbolic_max, :spacing, :symbolic_size)
    @eval begin
        $F(x::AbstractDimension) = PyObject(x).$F
        export $F
    end
end

"""
    symbolic_min(x::Dimension)
Symbol defining the minimum point of the Dimension
"""
function symbolic_min end

"""
    symbolic_max(x::Dimension)
Symbol defining the maximum point of the Dimension
"""
function symbolic_max end

"""
    spacing(x::Dimension)
Symbol representing the physical spacing along the Dimension.
"""
function spacing end

"""
    is_Derived(x::Dimension)
Returns true when dimension is derived, false when it is not
"""
function is_Derived end

"""
    symbolic_size(x::Dimension)
Symbol defining the size of the Dimension
"""
function symbolic_size end

# metaprograming for Devito functions taking variable number of arguments
for (M,F) in ((:devito,:Min), (:devito,:Max),(:sympy,:And))
    @eval begin
        $F(args...) = $M.$F((PyObject.(args))...)
        export $F
    end
end

"""
    Min(args...)

Can be used in a Devito.Eq to return the minimum of a collection of arguments
Example:
```julia
    eqmin = Eq(f,Min(g,1))
```
Is equivalent to f = Min(g,1) for Devito functions f,g
"""
function Min end

"""
    Max(args...)

Can be used in a Devito.Eq to return the minimum of a collection of arguments
Example:
```julia
    eqmax = Eq(f,Max(g,1))
```
Is equivalent to f = Max(g,1) for Devito functions f,g
"""
function Max end

# metaprograming for Devito mathematical operations ( more exist and may be added as required, find them at https://github.com/devitocodes/devito/blob/a8a33dc55ac3be008644c58a76b671028625679a/devito/finite_differences/elementary.py )

# these are broken out into four groups to help keep track of how they behave for unit testing

# functions defined on real numbers with equivalent in base
for F in (:cos, :sin, :tan, :sinh, :cosh, :tanh, :exp, :floor)
    @eval begin
        Base.$F(x::Union{AbstractDimension,DiscreteFunction,PyObject}) = devito.$F(PyObject(x))
    end
end
# functions defined on real numbers who are written differently in base
for F in (:Abs,:ceiling)
    @eval begin
        $F(x::Union{AbstractDimension,DiscreteFunction,PyObject}) = devito.$F(PyObject(x))
        export $F
    end
end
# functions defined on positive numbers with equivalent in base
for F in (:sqrt,)
    @eval begin
        Base.$F(x::Union{AbstractDimension,DiscreteFunction,PyObject}) = devito.$F(PyObject(x))
    end
end
# functions defined on positive numbers who are written differently in base
for F in (:ln,)
    @eval begin
        $F(x::Union{AbstractDimension,DiscreteFunction,PyObject}) = devito.$F(PyObject(x))
        export $F
    end
end

"""
    Mod(x::AbstractDimension,y::Int)
Perform Modular division on a dimension
"""
Mod(x::AbstractDimension,y::Int) = sympy.Mod(PyObject(x),PyObject(y))
export Mod

"""Get symbolic representation for function index object"""
function Base.getindex(x::Union{TimeFunction,Function},args...)
    py"""
    def indexobj(x,*args):
        return x[args]
    """
   return py"indexobj"(x,reverse(args)...)
end

"""
    ```
    ccode(x::Operator; filename="")
    ```
Print the ccode associated with a devito operator.
If filename is provided, writes ccode to disk using that filename
"""
function ccode(x::Operator; filename="")
    py"""
    def ccode(x, filename):
        if filename is "":
            return print(x)
        else:
            with open(filename, 'w') as f:
                print(x,file=f)
    """
   py"ccode"(x.o,filename)
   return nothing
end

"""
    SubDomain(name, instructions)

Create a subdomain by passing a list of instructions for each dimension.
Using an instruction with (nothing,) implies that the whole dimension should be used for that subdomain, as will ("middle",0,0)

Examples:
```julia
instructions = ("left",2),("middle",3,3)
SubDomain("subdomain_name",instructions)
```
or
```julia
instructions = [("right",4),("middle",1,2)]
SubDomain("subdomain_name",instructions)
```
or
```julia
SubDomain("subdomain_name",("right",2),("left",1))
```
"""
SubDomain(name::String, instructions::Vector) = SubDomain(name, instructions...)
SubDomain(name::String, instructions::Tuple{Vararg{Tuple}}) = SubDomain(name, instructions...)
function SubDomain(name::String, instructions...)
    # copy and reverse instructions
    instructions = reverse(instructions)
    N = length(instructions)
    @pydef mutable struct subdom <: devito.SubDomain
        function __init__(self, name, instructions)
            self.name = name
            self.instructions = instructions
        end
        function define(self, dimensions)
            dims = PyDict()
            for idim = 1:length(dimensions)
                dims[dimensions[idim]] = self.instructions[idim] in ((nothing,),("middle",0,0)) ? dimensions[idim] : self.instructions[idim]
            end
            return dims
        end
    end
    return SubDomain{N}(subdom(name,instructions))
end

export DiscreteFunction, Grid, Function, SparseTimeFunction, SubDomain,  TimeFunction, apply, backward, ccode, configuration, configuration!, coordinates, data, data_allocated, data_with_halo, data_with_inhalo, dimension, dimensions, dx, dy, dz, extent, forward, grid, halo, inject, interpolate, localindices, localindices_with_halo, localindices_with_inhalo, localsize, name, origin, size_with_halo, spacing, spacing_map, step, subdomains

end


text:
using DynamicEnergyBudgets, Test

@testset "ParallelComplementarySU" begin
    @test synthesizing_unit(ParallelComplementarySU(), 2.0, 2.0) ≈ 1.333333333333333333
    @test synthesizing_unit(ParallelComplementarySU(), 1000000000000000000000000.0, 2.0) == 2.0
    @test synthesizing_unit(ParallelComplementarySU(), 0.0, 2.0) == 0.0
end

@testset "MinimumRuleSU" begin
    @test synthesizing_unit(MinimumRuleSU(), 2.0, 2.0) == 2.0
    @test synthesizing_unit(MinimumRuleSU(), 1000000000000000000000000.0, 2.0) == 2.0
    @test synthesizing_unit(MinimumRuleSU(), 0.0, 2.0) == 0.0
end

@testset "KfamilySU" begin
    @test synthesizing_unit(KfamilySU(1.0), 2.0, 2.0) == 1.0
    @test synthesizing_unit(KfamilySU(2.0), 2.0, 2.0) > 1.0
    @test synthesizing_unit(KfamilySU(0.5), 2.0, 2.0) < 1.0

    @test synthesizing_unit(KfamilySU(2.0), 1000000000000000000000000.0, 2.0) ≈ 2.0
    @test synthesizing_unit(KfamilySU(1.0), 1000000000000000000000000.0, 2.0) ≈ 2.0
    @test synthesizing_unit(KfamilySU(0.5), 1000000000000000000000000.0, 2.0) ≈ 2.0

    @test synthesizing_unit(KfamilySU(2.0), 0.0, 2.0) == 0.0
    @test synthesizing_unit(KfamilySU(1.0), 0.0, 2.0) == 0.0
    @test synthesizing_unit(KfamilySU(0.5), 0.0, 2.0) == 0.0
end

@testset "SU comparisions" begin

    # Kfamiy ≈ ParallelComplementary when k ≈ 1.71. See Ledder et al, k family section. TODO: use exact page when published.
    @test synthesizing_unit(KfamilySU(1.71), 2.0, 2.0) ≈ synthesizing_unit(ParallelComplementarySU(), 2.0, 2.0) atol=0.001
    @test synthesizing_unit(KfamilySU(1.71), 3.0, 2.0) ≈ synthesizing_unit(ParallelComplementarySU(), 3.0, 2.0) atol=0.01

    # Kfamily cant handle large K falues because of floating point power limitations.
    # So we can't test the k → ∞ case where Kfamily == MinimumRule
    # @test synthesizing_unit(KfamilySU(10000.0), 2.0, 2.0) == synthesizing_unit(MinimumRuleSU(), 2.0, 2.0) == 2.0
end

@testset "Stoichiometric reserve merging" begin
    (a, b, c) = stoich_merge(ParallelComplementarySU(), 1.0, 1.0, 2.0, 2.0)

    @test a ≈ 0.3333333333333333333
    @test b ≈ 0.3333333333333333333
    @test c ≈ 1.3333333333333333333

    @test sum((a, b, c)) ≈ 1.0 + 1.0
end


text:
immutable UTF16String <: String
    data::Array{Uint16,1}
end

utf16_is_lead(c::Uint16) = (c & 0xfc00) == 0xd800
utf16_is_trail(c::Uint16) = (c & 0xfc00) == 0xdc00
utf16_is_surrogate(c::Uint16) = (c & 0xf800) == 0xd800
utf16_get_supplementary(lead::Uint16, trail::Uint16) = char((lead-0xd7f7)<<10 + trail)

function endof(s::UTF16String)
    d = s.data
    i = length(d)
    i == 0 && return i
    utf16_is_surrogate(d[i]) ? i-1 : i
end

function next(s::UTF16String, i::Int)
    if !utf16_is_surrogate(s.data[i])
        return char(s.data[i]), i+1
    elseif length(s.data) > i && utf16_is_lead(s.data[i]) && utf16_is_trail(s.data[i+1])
        return utf16_get_supplementary(s.data[i], s.data[i+1]), i+2
    end
    error("invalid UTF-16 character index")
end

# TODO: optmize this
function encode16(s::String)
    buf = Uint16[]
    for c in s
        if c < 0x10000
            push!(buf, uint16(c))
        else
            push!(buf, uint16(0xd7c0 + (c>>10) & 0x3ff))
            push!(buf, uint16(0xdc00 + c & 0x3ff))
        end
    end
    UTF16String(buf)
end

utf16(x) = convert(UTF16String, x)
convert(::Type{UTF16String}, s::UTF16String) = s
convert(::Type{UTF16String}, s::String) = encode16(s)
convert(::Type{UTF8String}, s::UTF16String) =
    sprint(length(s.data), io->for c in s; write(io,c::Char); end)
convert(::Type{Array{Uint16,1}}, s::UTF16String) = s.data
convert(::Type{Array{Uint16}}, s::UTF16String) = s.data

sizeof(s::UTF16String) = sizeof(s.data)
convert{T<:Union(Int16,Uint16)}(::Type{Ptr{T}}, s::UTF16String) =
    convert(Ptr{T}, s.data)

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using Printf, Test, Plots, LinearAlgebra, Distributions

include("../../hybmc/termstructures/YieldCurve.jl")
include("../../hybmc/models/QuasiGaussianModel.jl")

function setUp()
    yc = YieldCurve(0.03)
    d     = 3

    times = [ 1.0,    2.0,    5.0,    10.0     ]

    sigma = [  0.0060 0.0070 0.0080  0.0100 ;
               0.0030 0.0035 0.0040  0.0050 ;
               0.0025 0.0025 0.0025  0.0025 ]

    slope = [  0.10   0.15   0.20    0.25   ;
               0.20   0.35   0.40    0.55   ;
               0.10   0.10   0.10    0.10   ]

    curve = [  0.05   0.05   0.05    0.05   ;
               0.10   0.10   0.10    0.10   ;
               0.20   0.15   0.10    0.00   ]

    delta = [ 1.0,  5.0,  20.0 ]
    chi   = [ 0.01, 0.05, 0.15 ]

    Gamma = [ 1.0  0.8  0.6 ;
              0.8  1.0  0.8 ;
              0.6  0.8  1.0 ]

    # test matricies
    return QuasiGaussianModel(yc,d,times,sigma,slope,curve,delta,chi,Gamma)
end

@testset "test_ModelSetup" begin
    model = setUp()
    @test norm(model.DfT * transpose(model.DfT) - model.Gamma) == 0.0
    # Hf_H_inv from manual calculation
    Hf_H_inv = [ 0.99004983 0.95122942 0.86070798 ;
                 0.95122942 0.77880078 0.47236655 ;
                 0.81873075 0.36787944 0.04978707 ]
    @test norm(model.HHfInv * Hf_H_inv - Matrix{Float64}(I, 3, 3)) < 1.0e-7
    @test stateSize(model) == 13
    @test factors(model) == 3
    # test parameter functions
    times = Array(range(0.0, stop=11.0, length=23))

    sigma_ = [ sigma(model,i,t) for t in times, i = 1:model.d ]
    sigma_ = hcat(times, sigma_)
    sigmaRef = [
        0.00e+00  6.00e-03  3.00e-03  2.50e-03 ;
        5.00e-01  6.00e-03  3.00e-03  2.50e-03 ;
        1.00e+00  6.00e-03  3.00e-03  2.50e-03 ;
        1.50e+00  7.00e-03  3.50e-03  2.50e-03 ;
        2.00e+00  7.00e-03  3.50e-03  2.50e-03 ;
        2.50e+00  8.00e-03  4.00e-03  2.50e-03 ;
        3.00e+00  8.00e-03  4.00e-03  2.50e-03 ;
        3.50e+00  8.00e-03  4.00e-03  2.50e-03 ;
        4.00e+00  8.00e-03  4.00e-03  2.50e-03 ;
        4.50e+00  8.00e-03  4.00e-03  2.50e-03 ;
        5.00e+00  8.00e-03  4.00e-03  2.50e-03 ;
        5.50e+00  1.00e-02  5.00e-03  2.50e-03 ;
        6.00e+00  1.00e-02  5.00e-03  2.50e-03 ;
        6.50e+00  1.00e-02  5.00e-03  2.50e-03 ;
        7.00e+00  1.00e-02  5.00e-03  2.50e-03 ;
        7.50e+00  1.00e-02  5.00e-03  2.50e-03 ;
        8.00e+00  1.00e-02  5.00e-03  2.50e-03 ;
        8.50e+00  1.00e-02  5.00e-03  2.50e-03 ;
        9.00e+00  1.00e-02  5.00e-03  2.50e-03 ;
        9.50e+00  1.00e-02  5.00e-03  2.50e-03 ;
        1.00e+01  1.00e-02  5.00e-03  2.50e-03 ;
        1.05e+01  1.00e-02  5.00e-03  2.50e-03 ;
        1.10e+01  1.00e-02  5.00e-03  2.50e-03 ]
    @test norm(sigma_-sigmaRef) == 0.0

    slope_ = [ slope(model,i,t) for t in times, i = 1:model.d ]
    slope_ = hcat(times, slope_)
    slopeRef = [
         0.     0.1    0.2    0.1 ;
         0.5    0.1    0.2    0.1 ;
         1.     0.1    0.2    0.1 ;
         1.5    0.15   0.35   0.1 ;
         2.     0.15   0.35   0.1 ;
         2.5    0.2    0.4    0.1 ;
         3.     0.2    0.4    0.1 ;
         3.5    0.2    0.4    0.1 ;
         4.     0.2    0.4    0.1 ;
         4.5    0.2    0.4    0.1 ;
         5.     0.2    0.4    0.1 ;
         5.5    0.25   0.55   0.1 ;
         6.     0.25   0.55   0.1 ;
         6.5    0.25   0.55   0.1 ;
         7.     0.25   0.55   0.1 ;
         7.5    0.25   0.55   0.1 ;
         8.     0.25   0.55   0.1 ;
         8.5    0.25   0.55   0.1 ;
         9.     0.25   0.55   0.1 ;
         9.5    0.25   0.55   0.1 ;
        10.     0.25   0.55   0.1 ;
        10.5    0.25   0.55   0.1 ;
        11.     0.25   0.55   0.1 ]
    @test norm(slope_-slopeRef) == 0.0

    curve_ = [ curve(model,i,t) for t in times, i = 1:model.d ]
    curve_ = hcat(times, curve_)
    curveRef = [
         0.     0.05   0.1    0.2  ;
         0.5    0.05   0.1    0.2  ;
         1.     0.05   0.1    0.2  ;
         1.5    0.05   0.1    0.15 ;
         2.     0.05   0.1    0.15 ;
         2.5    0.05   0.1    0.1  ;
         3.     0.05   0.1    0.1  ;
         3.5    0.05   0.1    0.1  ;
         4.     0.05   0.1    0.1  ;
         4.5    0.05   0.1    0.1  ;
         5.     0.05   0.1    0.1  ;
         5.5    0.05   0.1    0.   ;
         6.     0.05   0.1    0.   ;
         6.5    0.05   0.1    0.   ;
         7.     0.05   0.1    0.   ;
         7.5    0.05   0.1    0.   ;
         8.     0.05   0.1    0.   ;
         8.5    0.05   0.1    0.   ;
         9.     0.05   0.1    0.   ;
         9.5    0.05   0.1    0.   ;
        10.     0.05   0.1    0.   ;
        10.5    0.05   0.1    0.   ;
        11.     0.05   0.1    0.   ]
        @test norm(curve_-curveRef) == 0.0
end

@testset "test_ZeroBond" begin
    yts = YieldCurve(0.03)
    d     = 3
    times = [  10.0   ]
    sigma = [  0.0060 ;
               0.0050 ;
               0.0040 ]
    slope = [  0.10   ;
               0.20   ;
               0.30   ]
    curve = [  0.05   ;
               0.10   ;
               0.20   ]
    delta = [ 1.0, 5.0, 20.0 ]
    chi   = [ 0.01, 0.05, 0.15 ]
    Gamma = [ 1.0  0.8  0.6 ;
              0.8  1.0  0.8 ;
              0.6  0.8  1.0 ]
    model = QuasiGaussianModel(yts,d,times,sigma,slope,curve,delta,chi,Gamma)
    uni = Uniform(0.0,0.05)
    x = rand(uni, d)
    y = rand(uni, (d,d))
    y = y * transpose(y)
    t = 5.0
    T = 10.0
    X = cat(x,vec(y),dims=1)
    X = cat(X,[0],dims=1)
    zcb = zeroBond(model,t,T,X,nothing)
    #
    G_ = (1.0 .- exp.(-chi*(T-t))) ./ chi
    zcbRef = discount(yts,T) / discount(yts,t) * exp(-transpose(G_) * x - 0.5 * transpose(G_) * y * G_)
    @test zcb == zcbRef
end

@testset "test_Sigma_xT" begin
    yts = YieldCurve(0.03)
    d     = 3
    times = [  10.0   ]
    sigma = [  0.0060 ;
               0.0050 ;
               0.0040 ]
    slope = [  0.10   ;
               0.20   ;
               0.30   ]
    curve = [  0.05   ;
               0.10   ;
               0.20   ]
    delta = [ 1.0, 5.0, 20.0 ]
    chi   = [ 0.01, 0.05, 0.15 ]
    Gamma = [ 1.0  0.8  0.6 ;
              0.8  1.0  0.8 ;
              0.6  0.8  1.0 ]
    model = QuasiGaussianModel(yts,d,times,sigma,slope,curve,delta,chi,Gamma)
    x = [ 0.00206853, 0.01093424, 0.00696041 ]  # np.random.uniform(0.0,0.05,d)
    y = [ 0.02980810  0.00656372 0.0041344  ;
          0.04282753  0.02406725 0.0001211  ;
          0.02906687  0.01078535 0.02328022 ] # np.random.uniform(0.0,0.05,[d,d])
    y = y * transpose(y)
    t = 5.0
    X = cat(x,vec(y),dims=1)
    X = cat(X,[0],dims=1)
    sigma_xT_ = sigma_xT(model,t,X)
    # reference result from QuantLib
    sigma_xT_Ref = [
         0.031658475910   0.017375011335   0.084573118461 ;
        -0.007058576846   0.021654204395  -0.134352049767 ;
        -0.016417240779  -0.043917617245   0.051199735933 ]
    @test norm(sigma_xT_ - sigma_xT_Ref) < 1.0e-12
end

@testset "test_Evolve" begin
    yts = YieldCurve(0.03)
    d     = 3
    times = [  10.0   ]
    sigma = [  0.0060 ;
               0.0050 ;
               0.0040 ]
    slope = [  0.10   ;
               0.20   ;
               0.30   ]
    curve = [  0.05   ;
               0.10   ;
               0.20   ]
    delta = [ 1.0, 5.0, 20.0 ]
    chi   = [ 0.01, 0.05, 0.15 ]
    Gamma = [ 1.0  0.8  0.6 ;
              0.8  1.0  0.8 ;
              0.6  0.8  1.0 ]
    model = QuasiGaussianModel(yts,d,times,sigma,slope,curve,delta,chi,Gamma)
    times = [0.0, 5.0, 10.0]
    dW = [ -0.8723107  -0.00585635  0.31102388 ;
           -0.15673275  0.28482762  0.79041947 ]
    Xref = [ 0.          0.           0.           0.           0.           0.           0.           0.           0.           0.           0.           0.           0.         ;
            -0.02182255  0.0397456   -0.01717905   0.0006412   -0.00116167   0.00059747  -0.00116167   0.00259077  -0.00141925   0.00059747  -0.00141925   0.00088815   0.00185998 ;
            -0.02476015  0.04705765  -0.0083194    0.00127306  -0.00232704   0.00102049  -0.00232704   0.00518216  -0.00243878   0.00102049  -0.00243878   0.00133699   0.03866521 ]
   #
   X = zeros((3,13))
   evolve(model, 0.0, X[1,:], 5.0, dW[1,:], @view X[2,:])
   evolve(model, 5.0, X[2,:], 5.0, dW[2,:], @view X[3,:])
   @test norm(X-Xref) < 1.0e-7
   # test zero coupon bond (again)
   zcbRef = 0.851676232405887
   zcb = zeroBond(model,5.0,10.0,X[2,:],nothing)
   @test abs(zcb - zcbRef) < 2.0e-10
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
ranges, my_ticket, all_tickets = split(read("day16/input.txt", String), "\n\n")
# ranges, my_ticket, all_tickets = split(read("day16/small.txt", String), "\n\n")

ranges = split(ranges, "\n")
my_ticket = parse.(Int, split(split(my_ticket, "\n")[2], ","))
all_tickets = split(all_tickets, "\n")[2:end]

all_ranges = []

for r in ranges
    s1, e1, s2, e2 = parse.(Int, match(r"(\d+)-(\d+) or (\d+)-(\d+)", r).captures)
    push!(all_ranges, Set([s1:e1..., s2:e2...]))
end

all_valid_nums = union(all_ranges...)

filter!(r -> all(n -> (n in all_valid_nums), parse.(Int, split(r, ","))), all_tickets)

all_tickets = parse.(Int, vcat(permutedims.(split.(all_tickets, ","))...))

matches_col(c, valid_set) = all(n -> n in valid_set, all_tickets[:, c])

matches_per_rule = [Set(filter(c -> matches_col(c, all_ranges[i]), 1:size(all_tickets, 2))) for i in 1:length(all_ranges)]

println(matches_per_rule)

function range_solver(start_col, taken_cols)
    if start_col == length(all_ranges)+1
        return []
    end
    for c in 1:size(all_tickets, 2)
        if c in taken_cols
            continue
        end
        if c in matches_per_rule[start_col]
            result = range_solver(start_col+1, Set([taken_cols..., c]))
            if result !== nothing
                return vcat([c], result)
            end
        end
    end
end

cols = range_solver(1, Set())

println(prod(i -> my_ticket[i], cols[1:6]))

# not proud of this guy, but it does the job


text:
using QDXML
using Base.Test

@test XML("<n><k a=b c=d g=\" h i\">g</k><y / ></n>s")["k"].attrs["g"] == " h i"


text:
DYNARE_ROOT = "/data/projects/dynare/git/preprocessor/src"
cd("models/example1")
run(`$DYNARE_ROOT/dynare_m example1.mod language=julia output=first json=compute`)
#cd("../example2")
#run(`$DYNARE_ROOT/dynare_m example2.mod language=julia output=first json=compute`)
cd("../..")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Experimental

!!! warning
    Types, methods, or macros defined in this module are experimental and subject
    to change and will not have deprecations. Caveat emptor.
"""
module Experimental

"""
    Const(A::Array)

Mark an Array as constant/read-only. The invariant guaranteed is that you will not
modify an Array (through another reference) within an `@aliasscope` scope.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N} <: DenseArray{T,N}
    a::Array{T,N}
end

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
@eval Base.getindex(A::Const, i1::Int) =
    (Base.@_inline_meta; Core.const_arrayref($(Expr(:boundscheck)), A.a, i1))
@eval Base.getindex(A::Const, i1::Int, i2::Int, I::Int...) =
  (Base.@_inline_meta; Core.const_arrayref($(Expr(:boundscheck)), A.a, i1, i2, I...))

"""
    @aliasscope expr

Allows the compiler to assume that all `Const`s are not being modified through stores
within this scope, even if the compiler can't prove this to be the case.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
macro aliasscope(body)
    sym = gensym()
    quote
        $(Expr(:aliasscope))
        $sym = $(esc(body))
        $(Expr(:popaliasscope))
        $sym
    end
end

end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
import LinearAlgebra: ldiv!, dot

export Preconditioner, InvPreconditioner, precond!

"""
$(TYPEDEF)

An immutable structure to hold a preconditioner and its Cholesky factorization.

## Fields
$(TYPEDFIELDS)

"""
struct Preconditioner
    P::Symmetric{Float64}

    F::Cholesky{Float64, Matrix{Float64}}
end

function Preconditioner(P::Matrix{Float64})
    return Preconditioner(Symmetric(P), cholesky(Symmetric(P)))
end

ldiv!(x, P::Preconditioner, b) = copyto!(x, P.F \ b)
dot(A::Array, P::Preconditioner, B::Vector) = dot(A, P.P, B)

"""
    InvPreconditioner

An immutable structure to hold the inverse of the preconditioner.
For instance, this structure can be used to hold the estimate of the inverse of the Hessian in the BFGS algorithm.
"""
struct InvPreconditioner
    InvP::Symmetric{Float64}

    F::Cholesky{Float64, Matrix{Float64}}
end

function InvPreconditioner(InvP::Matrix{Float64})
    return InvPreconditioner(Symmetric(InvP), cholesky(Symmetric(InvP)))
end

ldiv!(x, P::InvPreconditioner, b) = copyto!(x, P.InvP * b)
dot(A::Array, P::InvPreconditioner, B::Vector) = ldiv!(A, P.F, B)

"""
$(TYPEDSIGNATURES)

Computes in-place a Hessian preconditioner based on the Gauss-Newton approximation (outer-product of gradient of the cost function).
"""
function precond!(P, coeff, S::Storage, C::HermiteMapComponent, X)
    Nψ = C.Nψ
    NxX, Ne = size(X)
    @assert NxX == C.Nx "Wrong dimension of the sample X"
    @assert size(S.ψoff, 1) == Ne
    @assert size(S.ψoff, 2) == Nψ

    # Output objective, gradient
    xlast = view(X,NxX,:)#)

    fill!(S.cache_integral, 0)

    # Integrate at the same time for the objective, gradient
    function integrand!(v::Vector{Float64}, t::Float64)
        repeated_grad_xk_basis!(S.cache_dcψxdt, S.cache_gradxd, C.I.f, t*xlast)

        # @avx @. S.cache_dψxd = (S.cache_dcψxdt .* S.ψoff) *ˡ coeff
        @avx @. S.cache_dcψxdt *= S.ψoff
        mul!(S.cache_dψxd, S.cache_dcψxdt, coeff)

        # Integration for J
        vJ = view(v,1:Ne)
        evaluate!(vJ, C.I.g, S.cache_dψxd)

        # Integration for dcJ

        grad_x!(S.cache_dψxd, C.I.g, S.cache_dψxd)

        # @avx for j=2:Nψ+1
        #     for i=1:Ne
        #         v[(j-1)*Ne+i] = S.cache_dψxd[i]*S.cache_dcψxdt[i,j-1]
        #     end
        # end

        v[Ne+1:Ne+Ne*Nψ] .= reshape(S.cache_dψxd .* S.cache_dcψxdt , (Ne*Nψ))
        # v[Ne+1:Ne+Ne*Nψ] .= reshape(grad_x(C.I.g, S.cache_dψxd) .* S.cache_dcψxdt , (Ne*Nψ))
    end

    quadgk!(integrand!, S.cache_integral, 0.0, 1.0)#; rtol = 1e-3)#; order = 9, rtol = 1e-10)

    # Multiply integral by xk (change of variable in the integration)
    @inbounds for j=1:Nψ+1
        @. S.cache_integral[(j-1)*Ne+1:j*Ne] *= xlast
    end
    # @avx for j=1:Nψ+1
    #     for i=1:Ne
    #         S.cache_integral[(j-1)*Ne+i] *= xlast[i]
    #     end
    # end

    # Add f(x_{1:d-1},0) i.e. (S.ψoff .* S.ψd0)*coeff to S.cache_integral
    @avx for i=1:Ne
        f0i = zero(Float64)
        for j=1:Nψ
            f0i += S.ψoffψd0[i,j]*coeff[j]
        end
        S.cache_integral[i] += f0i
    end

    # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
    @avx for i=1:Ne
        prelogJi = zero(Float64)
        for j=1:Nψ
            prelogJi += S.ψoffdψxd[i,j]*coeff[j]
        end
        S.cache_g[i] = prelogJi
    end

    reshape_cacheintegral = reshape(S.cache_integral[Ne+1:Ne+Ne*Nψ], (Ne, Nψ))

    fill!(P, 0.0)
    @inbounds for l=1:Ne
        # Exploit symmetry of the Hessian
        for i=1:Nψ
            for j=i:Nψ
            # P[i,j] +=  reshape2_cacheintegral[l,i,j]*S.cache_integral[l]
            P[i,j] +=  (reshape_cacheintegral[l,i] + S.ψoffψd0[l,i]) * (reshape_cacheintegral[l,j] + S.ψoffψd0[l,j])
            P[i,j] -=  ( (S.ψoffdψxd[l,i]) * (S.ψoffdψxd[l,j])*(
                            hess_x(C.I.g, S.cache_g[l]) * C.I.g(S.cache_g[l]) -
                            grad_x(C.I.g, S.cache_g[l])^2))/C.I.g(S.cache_g[l])^2

            P[j,i] = P[i,j]
            end
        end
    end
    rmul!(P, 1/Ne)
    # Add derivative of the L2 penalty term ∂^2_c α ||c||^2 = 2 *α *I
    @inbounds for i=1:Nψ
        P[i,i] += 2*C.α*I
    end
    return P
end

"""
$(TYPEDSIGNATURES)

Computes in-place a Hessian preconditioner based on the Gauss-Newton approximation (outer-product of gradient of the cost function).
"""
precond!(S::Storage, C::HermiteMapComponent, X) = (P, coeff) -> precond!(P, coeff, S, C, X)

"""
$(TYPEDSIGNATURES)

Computes in-place a Hessian preconditioner based on the diagonal entries of the Gauss-Newton approximation (outer-product of gradient of the cost function).
"""
function diagprecond!(P, coeff, S::Storage, C::HermiteMapComponent, X::Array{Float64,2})
    Nψ = C.Nψ
    Nx = C.Nx
    NxX, Ne = size(X)
    @assert NxX == Nx "Wrong dimension of the sample X"
    @assert size(S.ψoff, 1) == Ne
    @assert size(S.ψoff, 2) == Nψ

    # Output objective, gradient
    xlast = view(X, Nx,:)#)

    fill!(S.cache_integral, 0)

    # Integrate at the same time for the objective, gradient
    function integrand!(v::Vector{Float64}, t::Float64)
        S.cache_dcψxdt .= repeated_grad_xk_basis(C.I.f, t*xlast)

        # @avx @. S.cache_dψxd = (S.cache_dcψxdt .* S.ψoff) *ˡ coeff
        S.cache_dcψxdt .*= S.ψoff
        mul!(S.cache_dψxd, S.cache_dcψxdt, coeff)

        # Integration for J
        vJ = view(v,1:Ne)
        evaluate!(vJ, C.I.g, S.cache_dψxd)
        # v[1:Ne] .= C.I.g(S.cache_dψxd)

        # Integration for dcJ
        v[Ne+1:Ne+Ne*Nψ] .= reshape(grad_x(C.I.g, S.cache_dψxd) .* S.cache_dcψxdt , (Ne*Nψ))
    end

    quadgk!(integrand!, S.cache_integral, 0.0, 1.0)#; rtol = 1e-3)#; order = 9, rtol = 1e-10)

    # Multiply integral by xk (change of variable in the integration)
    @inbounds for j=1:Nψ+1
        @. S.cache_integral[(j-1)*Ne+1:j*Ne] *= xlast
    end

    # Add f(x_{1:d-1},0) i.e. (S.ψoff .* S.ψd0)*coeff to S.cache_integral
    @avx for i=1:Ne
        f0i = zero(Float64)
        for j=1:Nψ
            f0i += S.ψoffψd0[i,j]*coeff[j]
        end
        S.cache_integral[i] += f0i
    end

    # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
    @avx for i=1:Ne
        prelogJi = zero(Float64)
        for j=1:Nψ
            prelogJi += S.ψoffdψxd[i,j]*coeff[j]
        end
        S.cache_g[i] = prelogJi
    end

    reshape_cacheintegral = reshape(S.cache_integral[Ne+1:Ne+Ne*Nψ], (Ne, Nψ))
    # reshape2_cacheintegral = reshape(S.cache_integral[Ne + Ne*Nψ + 1: Ne + Ne*Nψ + Ne*Nψ*Nψ], (Ne, Nψ, Nψ))
    # @show reshape2_cacheintegral
    fill!(P, 0.0)
    @inbounds for l=1:Ne
        # Exploit symmetry of the Hessian
        for i=1:Nψ
            # P[i,j] +=  reshape2_cacheintegral[l,i,j]*S.cache_integral[l]
            P[i] +=  (reshape_cacheintegral[l,i] + S.ψoffψd0[l,i])^2# * (reshape_cacheintegral[l,j] + S.ψoff[l,j]*S.ψd0[l,j])
            P[i] -=  ( (S.ψoffdψxd[l,i])^2*(
                            hess_x(C.I.g, S.cache_g[l]) * C.I.g(S.cache_g[l]) -
                            grad_x(C.I.g, S.cache_g[l])^2))/C.I.g(S.cache_g[l])^2
        end
    end
    rmul!(P, 1/Ne)
    # Add derivative of the L2 penalty term ∂^2_c α ||c||^2 = 2 *α *I
    @inbounds for i=1:Nψ
        P[i] += 2*C.α
    end
    return P
end

"""
$(TYPEDSIGNATURES)

Computes in-place a Hessian preconditioner based on the diagonal entries of the Gauss-Newton approximation (outer-product of gradient of the cost function).
"""
diagprecond!(S::Storage, C::HermiteMapComponent, X::Array{Float64,2}) = (P, coeff) -> diagprecond!(P, coeff, S, C, X)


text:
# Assumes the space dimension and field dimension are the same.  In other words, when the
# space coordinate indices are (i,j,k), then the field has three vector components.
# Therefore, for the input field array F[i,j,k,w], we assume w = 1:3.

# The functions add the calculated values to the existing values of the output array.
# Therefore, if the derivative values themselves are desired, pass the output array
# initialized with zeros.

export apply_divg!

# Wrapper to apply the discrete divergence by default
#
# Note that when the positional arguments are SVector's, this function may call the concrete
# implementation directly instead of the next wrapper function taking ∆l⁻¹ even if the
# keyword arguments are not SVector's, because methods are dispatched based only on the
# types of positional arguments.  Therefore, we convert the keyword arguments into SVector's.
apply_divg!(g::AbsArrNumber,  # output array of scalar; in 3D, g[i,j,k] is g at (i,j,k)
            F::AbsArrNumber,  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: create ∂w by forward|backward difference
            ∆l⁻¹::Tuple{Vararg{Number}}=ntuple(x->1.0,length(isfwd)),  # ∆l⁻¹[w]: inverse of uniform distance between grid planes in w-direction
            isbloch::AbsVecBool=fill(true,length(isfwd)),  # boundary conditions in K dimensions
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(length(isfwd));  # Bloch phase factors in K dimensions
            α::Number=1.0  # scale factor to multiply to result before adding it to g: g += α ∇⋅F
            ) where {OP} =
    (N = size(g); K = length(isfwd); apply_divg!(g, F, Val(OP), isfwd, fill.(∆l⁻¹,N), isbloch, e⁻ⁱᵏᴸ; α))

# Wrapper for converting AbstractVector's to SVec's
apply_divg!(g::AbsArrNumber,  # output array of scalar; in 3D, g[i,j,k] is g at (i,j,k)
            F::AbsArrNumber,  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: create ∂w by forward|backward difference
            ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l⁻¹[w]: inverse of distances between grid planes in x-direction
            isbloch::AbsVecBool=fill(true,K),  # boundary conditions in K dimensions
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # Bloch phase factors in K dimensions
            α::Number=1.0  # scale factor to multiply to result before adding it to g: g += α ∇⋅F
            ) where {K,OP} =
    apply_divg!(g, F, Val(OP), SBool{K}(isfwd), ∆l⁻¹, SBool{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); α)

# Concrete implementation
function apply_divg!(g::AbsArrNumber{K},  # output array of scalar; in 3D, g[i,j,k] is g at (i,j,k)
                     F::AbsArrNumber{K₊₁},  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
                     ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                     isfwd::SBool{K},  # isfwd[w] = true|false: create ∂w by forward|backward difference
                     ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l[w]: inverse of distances between grid planes in x-direction
                     isbloch::SBool{K},  # boundary conditions in K dimensions
                     e⁻ⁱᵏᴸ::SNumber{K};  # Bloch phase factors in K dimensions
                     α::Number=1.0  # scale factor to multiply to result before adding it to g: g += α ∇⋅F
                     ) where {K,K₊₁,OP}
    @assert K₊₁==K+1

    # First partial derivative
    nw = 1
    nu = nw  # component of input field to feed to ∂w
    Fu = selectdim(F, K₊₁, nu)  # nu-th component of input field
    apply_∂!(g, Fu, Val(OP), nw, isfwd[nw], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]; α)  # g += α ∂Fu/∂w

    # Remaining partial derivatives
    for nw = 2:K  # direction of differentiation
        nu = nw  # component of input field to feed to ∂w
        Fu = selectdim(F, K₊₁, nu)  # nu-th component of input field
        apply_∂!(g, Fu, Val(:(+=)), nw, isfwd[nw], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]; α)  # g += α ∂Fu/∂w
    end

    return nothing
end


text:
# using RoMEPlotting
# using IncrementalInference, RoME
using Test

## See tests in DFG, IIF, and RoME for more complete coverage of features

include("pcl/testPointCloud2.jl")

include("testPose2AprilTag4Corner.jl")

using ZMQ
using Caesar, Caesar.ZmqCaesar

@testset "ZMQ Interface" begin
    # Unit tests
    include("multilangzmq/callbackCompatibilityTest.jl")
    # Integration tests
    include("multilangzmq/zmqInternal.jl")
    # End to end tests
    # Do this manually...
    # include("multilangzmq/serverTest.jl")
end

# include("testScatterAlignPose2.jl")


text:
module VisionModel

using Flux, Images

# Models
include("alexnet.jl")
include("vgg.jl")

end


text:
using DDF

using DifferentialForms: bitsign
using LinearAlgebra
using StaticArrays
using Test

const manifolds = Dict{Any,Any}()
@testset "Create manifolds D=$D" for D in 0:Dmax
    S = Float64
    manifolds[D] = [empty_manifold(Val(D), S), simplex_manifold(Val(D), S),
                    hypercube_manifold(Val(D), S),
                    delaunay_hypercube_manifold(Val(D), S),
                    large_hypercube_manifold(Val(D), S; optimize_mesh=false),
                    large_delaunay_hypercube_manifold(Val(D), S),
                    refined_simplex_manifold(Val(D), S),
                    boundary_simplex_manifold(Val(D), S)]
end

@testset "Manifold ops D=$D P=$P R=$R mfd=$(mfd.name)" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D,
mfd in manifolds[D]

    C = length(eltype(get_coords(mfd)))
    S = eltype(eltype(get_coords(mfd)))
    mfd::Manifold{D,C,S}

    funs = Fun{D,P,R,C,S}[]
    push!(funs, zero(Fun{D,P,R,C,S,S}, mfd))
    e = id(Fun{D,P,0,C,S,SVector{C,S}}, mfd)
    if R == 0
        push!(funs, e)
    end

    dR = P == Pr ? +1 : -1

    # Boundary
    R2 = R - 2dR
    if 0 ≤ R2 ≤ D
        b1 = boundary(Val(P), Val(R), mfd)
        b2 = boundary(Val(P), Val(R - dR), mfd)
        b21 = b2 * b1
        b21::Op{D,P,R2,P,R,Int8}
        @test iszero(b21)
    end

    R1 = R - dR
    if 0 ≤ R1 ≤ D
        b = boundary(Val(P), Val(R), mfd)
        for f in funs
            bf = b * f
            T = eltype(f)
            bf::Fun{D,P,R1,C,S,T}
            bf′ = boundary(f)
            @test bf′ == bf
        end
    end

    # Derivative
    R2 = R + 2dR
    if 0 ≤ R2 ≤ D
        d1 = deriv(Val(P), Val(R), mfd)
        d2 = deriv(Val(P), Val(R + dR), mfd)
        d21 = d2 * d1
        d21::Op{D,P,R2,P,R,Int8}
        @test iszero(d21)
    end

    R1 = R + dR
    if 0 ≤ R1 ≤ D
        d = deriv(Val(P), Val(R), mfd)
        for f in funs
            df = d * f
            T = eltype(f)
            df::Fun{D,P,R1,C,S,T}
            df′ = deriv(f)
            @test df′ == df
        end
    end

    # TODO: derivatives are linear, product rule; same for hodge,
    # coderivatives etc.

    # Only if completely well-centred [arXiv:0802.2108 [cs.CG]]
    if mfd.dualkind == BarycentricDuals ||
       mfd.name ∉ ["hypercube manifold", "delaunay hypercube manifold",
        "large delaunay hypercube manifold", "refined simplex manifold"]

        # Hodge
        h = hodge(Val(P), Val(R), mfd)
        h′ = hodge(Val(!P), Val(R), mfd)
        @test isempty(h) || all(x -> x != 0 && isfinite(x), h.values.diag)
        @test isempty(h′) || all(x -> x != 0 && isfinite(x), h′.values.diag)
        @test invhodge(Val(P), Val(R), mfd) ==
              bitsign(isodd(D) ? 0 : P == Pr ? R : D - R) * h
        @test invhodge(Val(!P), Val(R), mfd) ==
              bitsign(isodd(D) ? 0 : !P == Pr ? R : D - R) * h′
        h21 = h′ * h
        h21::Op{D,P,R,P,R,S}
        @test isempty(h21) || all(x -> x != 0 && isfinite(x), h21.values.diag)
        if !isempty(h21)
            s = bitsign(R * (D - R))
            @test norm((h21 - s * one(h21)).values, Inf) ≤ 10eps(S)
        end

        for f in funs
            hf = h * f
            T = eltype(f)
            hf::Fun{D,!P,R,C,S,T}
            hf′ = ⋆f
            @test hf′ == hf
        end

        # Coderivative
        R2 = R - 2dR
        if 0 ≤ R2 ≤ D
            δ1 = coderiv(Val(P), Val(R), mfd)
            δ2 = coderiv(Val(P), Val(R - dR), mfd)
            δ21 = δ2 * δ1
            δ21::Op{D,P,R2,P,R,S}
            nδ1 = norm(δ1, Inf)
            nδ2 = norm(δ2, Inf)
            maxerr = nδ1 * nδ2 * 10eps(S)
            @test norm(δ21, Inf) ≤ maxerr
        end

        R1 = R - dR
        if 0 ≤ R1 ≤ D
            δ = coderiv(Val(P), Val(R), mfd)
            @test isempty(δ) ||
                  all(x -> x != 0 && isfinite(x), nonzeros(δ.values))
            for f in funs
                δf = δ * f
                T = eltype(f)
                δf::Fun{D,P,R1,C,S,T}
                δf′ = coderiv(f)
                @test δf′ == δf
            end
        end

        # Laplace
        Δ = laplace(Val(P), Val(R), mfd)
        @test isempty(Δ) || all(x -> x != 0 && isfinite(x), nonzeros(Δ.values))
        for f in funs
            Δf = Δ * f
            T = eltype(f)
            Δf::Fun{D,P,R,C,S,T}
            Δf′ = laplace(f)
            @test Δf′ == Δf
        end
    end
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_transmission_flows(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	# Transmission related values
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	L = inputs["L"]     # Number of transmission lines
	# Power flows on transmission lines at each time step
	dfFlow = DataFrame(Line = 1:L, Sum = Array{Union{Missing,Float32}}(undef, L))
	if setup["ParameterScale"] == 1
		for i in 1:L
			dfFlow[!,:Sum][i] = sum(inputs["omega"].* (value.(EP[:vFLOW])[i,:])) * ModelScalingFactor
		end
		dfFlow = hcat(dfFlow, DataFrame((value.(EP[:vFLOW])) * ModelScalingFactor, :auto))
	else
		for i in 1:L
			dfFlow[!,:Sum][i] = sum(inputs["omega"].* (value.(EP[:vFLOW])[i,:]))
		end
		dfFlow = hcat(dfFlow, DataFrame(value.(EP[:vFLOW]), :auto))
	end
	auxNew_Names=[Symbol("Line");Symbol("Sum");[Symbol("t$t") for t in 1:T]]
	rename!(dfFlow,auxNew_Names)
	total = DataFrame(["Total" sum(dfFlow[!,:Sum]) fill(0.0, (1,T))], :auto)
	for t in 1:T
		total[:,t+2] .= sum(dfFlow[:,Symbol("t$t")][1:L])
	end
	rename!(total,auxNew_Names)
	dfFlow = vcat(dfFlow, total)

	CSV.write(joinpath(path, "flow.csv"), dftranspose(dfFlow, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
#######################################################################
#
# A demo of an various Excel chart data tools that are available via
# an XlsxWriter chart.
#
# These include, Trendlines, Data Labels, Error Bars, Drop Lines,
# High-Low Lines and Up-Down Bars.
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter
#

using XlsxWriter

function test()

	wb = Workbook("chart_data_table.xlsx")
	ws = add_worksheet!(wb)

	# Add a format for the headings.
	bold = add_format!(wb, Dict("bold"  => 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Number" "Data 1" "Data 2"]
	data = [
		[2 3 4 5 6 7];
		[10 40 50 20 10 50];
		[30 60 70 50 40 30];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])
	write_column!(ws, "C2", data[3,:])

	#######################################################################
	#
	# Trendline example.
	#
	# Create a Line chart.
	chart1 = add_chart!(wb, Dict("type"=> "line"))

	# Configure the first series with a polynomial trendline.
	add_series!(chart1, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
		"trendline"=> Dict(
			"type"=> "polynomial",
			"order"=> 3,
		),
	))

	# Configure the second series with a moving average trendline.
	add_series!(chart1, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
		"trendline"=> Dict("type"=> "linear"),
	))

	# Add a chart title. and some axis labels.
	set_title!(chart1, Dict("name"=> "Chart with Trendlines"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D2", chart1, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Data Labels and Markers example.
	#
	# Create a Line chart.
	chart2 = add_chart!(wb, Dict("type"=> "line"))

	# Configure the first series.
	add_series!(chart2, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
		"data_labels"=> Dict("value"=> 1),
		"marker"=> Dict("type"=> "automatic"),
	))

	# Configure the second series.
	add_series!(chart2, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title. and some axis labels.
	set_title!(chart2, Dict("name"=> "Chart with Data Labels and Markers"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D18", chart2, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Error Bars example.
	#
	# Create a Line chart.
	chart3 = add_chart!(wb, Dict("type"=> "line"))

	# Configure the first series.
	add_series!(chart3, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
		"y_error_bars"=> Dict("type"=> "standard_error"),
	))

	# Configure the second series.
	add_series!(chart3, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=> "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title. and some axis labels.
	set_title!(chart3, Dict("name"=> "Chart with Error Bars"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D34", chart3, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Up-Down Bars example.
	#
	# Create a Line chart.
	chart4 = add_chart!(wb, Dict("type"=> "line"))

	# Add the Up-Down Bars.
	set_up_down_bars!(chart4)

	# Configure the first series.
	add_series!(chart4, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure the second series.
	add_series!(chart4, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title. and some axis labels.
	set_title!(chart4, Dict("name"=> "Chart with Up-Down Bars"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D50", chart4, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# High-Low Lines example.
	#
	# Create a Line chart.
	chart5 = add_chart!(wb, Dict("type"=> "line"))

	# Add the High-Low lines.
	set_high_low_lines!(chart5)

	# Configure the first series.
	add_series!(chart5, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure the second series.
	add_series!(chart5, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title. and some axis labels.
	set_title!(chart5, Dict("name"=> "Chart with High-Low Lines"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D66", chart5, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Drop Lines example.
	#
	# Create a Line chart.
	chart6 = add_chart!(wb, Dict("type"=> "line"))

	# Add Drop Lines.
	set_drop_lines!(chart6)

	# Configure the first series.
	add_series!(chart6, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure the second series.
	add_series!(chart6, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title. and some axis labels.
	set_title!(chart6, Dict("name"=> "Chart with Drop Lines"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D82", chart6, Dict("x_offset"=> 25, "y_offset"=> 10))

	close(wb)
	isfile("chart_data_table.xlsx")
end

test()