text:
# Lorenz.jl

"""
   lorenz(v)

Compute the relative Lorenz Curve of a vector `v` .

Returns two vectors. The first one contains the cumulative proportion of people.
The second contains the cumulative share of income earned.

# Examples
```julia
julia> using Inequality
julia> lorenz_curve([8, 5, 1, 3, 5, 6, 7, 6, 3])
([0.0, 0.1111111111111111, 0.2222222222222222, 0.3333333333333333, 0.4444444444444444, 0.5555555555555556, 0.6666666666666666, 0.7777777777777778, 0.8888888888888888, 1.0],
│ [0.0, 0.022727272727272728, 0.09090909090909091, 0.1590909090909091, 0.2727272727272727, 0.38636363636363635, 0.5227272727272727, 0.6590909090909091, 0.8181818181818182, 1.0])
```
"""
function lorenz_curve(v::AbstractVector{<:Real})
   v = sort(v)
   cumsum!(v,v)
   cumsum_x = Vector{Float64}(undef, length(v))
   cumsum_y = Vector{Float64}(undef, length(v))

   @inbounds for (i, val) in enumerate(v)
    cumsum_x[i] = i/length(v)
    cumsum_y[i] = val/v[end]
   end

   pushfirst!(cumsum_x, 0.0)
   pushfirst!(cumsum_y, 0.0)

   return cumsum_x, cumsum_y

end

"""
   lorenz(v, w)
Compute the weighted Lorenz Curve of a vector `v` using weights given by a weight vector `w`.

Weights must not be negative, missing or NaN. The weights and data vectors must have the same length.

Returns two vectors. The first one contains the cumulative proportion of weighted people.
The second contains the cumulative share of income earned.

# Examples
```julia
julia> using Inequality
julia> lorenz_curve([8, 5, 1, 3, 5, 6, 7, 6, 3], collect(0.1:0.1:0.9))
([0.0, 0.06666666666666667, 0.08888888888888889, 0.13333333333333333, 0.2222222222222222, 0.3333333333333333, 0.5333333333333333, 0.6666666666666666, 0.8444444444444444, 1.0],
[0.0, 0.013761467889908256, 0.05045871559633028, 0.0963302752293578, 0.1513761467889908, 0.2660550458715596, 0.38990825688073394, 0.555045871559633, 0.7752293577981653, 1.0])
```
"""
function lorenz_curve(v::AbstractVector{<:Real}, w::AbstractVector{<:Real})

    checks_weights(v, w)

    v = v .* w

    w = w[sortperm(v)]/sum(w)
    sort!(v)

    cumsum!(v,v)
    cumsum!(w,w)

    cumsum_x = Vector{Float64}(undef, length(v))
    cumsum_y = Vector{Float64}(undef, length(v))

    @inbounds for (i, (val, wgt)) in enumerate(zip(v, w))

     cumsum_x[i] = wgt
     cumsum_y[i] = val/v[end]
    end

    pushfirst!(cumsum_x, 0.0)
    pushfirst!(cumsum_y, 0.0)

    return cumsum_x, cumsum_y

 end

 function lorenz_curve(v::AbstractVector{<:Real}, w::AbstractWeights)

   checks_weights(v, w)

   v = v .* w

   w = w[sortperm(v)]/w.sum
   sort!(v)

   cumsum!(v,v)
   cumsum!(w,w)

   cumsum_x = Vector{Float64}(undef, length(v))
   cumsum_y = Vector{Float64}(undef, length(v))

   @inbounds for (i, (val, wgt)) in enumerate(zip(v, w))

    cumsum_x[i] = wgt
    cumsum_y[i] = val/v[end]
   end

   pushfirst!(cumsum_x, 0.0)
   pushfirst!(cumsum_y, 0.0)

   return cumsum_x, cumsum_y

end

 """
   wlorenz_curve(v, w)

 Compute the weighted relative Lorenz Curve of `v` with weights `w`. See also [`lorenz_curve`](lorenz_curve)
 """
 wlorenz_curve(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}) = lorenz_curve(v, w)

 wlorenz_curve(v::AbstractVector{<:Real}, w::AbstractWeights) = lorenz_curve(v, w)


text:
mutable struct Kriging{N,M,XS<:Scaler,YS<:Scaler}
    xscaler::XS
    yscaler::YS
    gps::Vector{GP.GPE}
end

include("tune.jl")
include("predict.jl")
include("append.jl")
function Kriging(
    x::AbstractMatrix,
    y::AbstractMatrix,
    mean::GP.Mean,
    kernel::GP.Kernel;
    xscaler::Scaler = StandardScaler(x),
    yscaler::Scaler = StandardScaler(y),
    capacity::Integer = 1024,
    stepsize::Integer = 1024,
)
    nx, mx = size(x)
    ny, my = size(y)
    if mx != my
        throw(ArgumentError("`x` and `y` must have same number of observations"))
    end
    xs = xscaler * x
    ys = yscaler * y
    logNoise = log((1000 + mx) * eps())
    if isa(mean, GP.Mean)
        mean = [deepcopy(mean) for i in 1:ny]
    elseif length(mean) != ny || !all(isa(mean, GP.Mean))
        throw(ArgumentError("`mean` must be a GP.Mean object OR a list/tuple of $ny GP.Mean objects"))
    end

    if isa(kernel, GP.Kernel)
        kernel = [deepcopy(kernel) for i in 1:ny]
    elseif length(kernel) != ny || !all(isa.(kernel, GP.Kernel))
        throw(ArgumentError("`kernel` must be a GP.Kernel object OR a list/tuple of $ny GP.Kernel objects"))
    end
    gps = Vector{GP.GPE}(undef, ny)
    for i in 1:ny
        gps[i] = GP.ElasticGPE(
            nx;
            mean = mean[i],
            kernel = kernel[i],
            logNoise = logNoise,
            capacity = capacity,
            stepsize = stepsize,
        )
        append!(gps[i], xs, ys[i, :])
    end
    Kriging{nx,ny,typeof(xscaler),typeof(yscaler)}(xscaler, yscaler, gps)
end

getx(krg::Kriging, i::Int) = krg.xscaler \ view(krg.gps[1].x, :, i)
getx(krg::Kriging) = krg.xscaler \ krg.gps[1].x

function gety(krg::Kriging{N,M}, i::Int) where {N,M}
    y = similar(krg.gps[1].y, M)
    @inbounds for j in 1:M
        y[j] = krg.gps[j].y[i]
    end
    inverse!(krg.yscaler, y)
    y
end

function gety(krg::Kriging{N,M}) where {N,M}
    nsamples = length(krg.gps[1].y)
    y = similar(krg.gps[1].y, M, nsamples)
    @inbounds for j in 1:nsamples
        for i in 1:M
            y[i, j] = krg.gps[i].y[j]
        end
    end
    inverse!(krg.yscaler, y)
    y
end

Base.getindex(krg::Kriging, i::Integer) = gety(krg, i), getx(krg, i)

function get_samples(krg::Kriging{N,M}) where {N,M}
    getx(krg), gety(krg)
end

ext: jl


text:
### RAM holds the fields of the robust adaptive Metropolis (RAM) sampler

immutable RAM <: MHSampler
  jumpscale::Float64 # Scale factor for adapting the jump size
  targetrate::Float64 # Target acceptance rate

  function RAM(s::Float64, r::Float64)
    @assert s > 0 "Scale factor for adapting the jump size is not positive."
    @assert 0 < r < 1 "Target acceptance rate is not between 0 and 1"
    new(s, r)
  end
end

RAM() = RAM(1., 0.234)
RAM(s::Float64) = RAM(s, 0.234)

RAM(; jumpscale::Float64=1.0, targetrate::Float64=0.234) = RAM(jumpscale, targetrate)

### RAMHeap type holds the internal state ("local variables") of the RAM sampler

type RAMHeap <: MCHeap{MCBaseSample}
  instate::MCState{MCBaseSample} # Monte Carlo state used internally by the sampler
  outstate::MCState{MCBaseSample} # Monte Carlo state outputted by the sampler
  tune::MCTune
  count::Int # Current number of iterations
  ratio::Float64
  eta::Float64
  SS::Matrix{Float64}
  S::Matrix{Float64}
  rnormal::Vector{Float64}
end

RAMHeap() =
  RAMHeap(MCState(MCBaseSample(), MCBaseSample()), MCState(MCBaseSample(), MCBaseSample()), VanillaMCTune(), 0, NaN,
  NaN, Array(Float64, 0, 0), Array(Float64, 0, 0), Float64[])

RAMHeap(l::Int) =
  RAMHeap(MCState(MCBaseSample(l), MCBaseSample(l)), MCState(MCBaseSample(l), MCBaseSample(l)), VanillaMCTune(), 0,
  NaN, NaN, fill(NaN, l, l), fill(NaN, l, l), fill(NaN, l))

### Initialize RAM sampler

function initialize_heap(m::MCModel, s::RAM, r::MCRunner, t::MCTuner)
  heap::RAMHeap = RAMHeap(m.size)

  heap.instate.current = MCBaseSample(copy(m.init))
  logtarget!(heap.instate.current, m.eval)
  @assert isfinite(heap.instate.current.logtarget) "Initial values out of model support."

  heap.S = Float64[i == j ? (m.scale.*s.jumpscale)[i] : 0. for i in 1:m.size, j in 1:m.size]

  heap.tune = VanillaMCTune()

  heap.count = 1

  heap
end

function reset!(heap::RAMHeap, x::Vector{Float64}, m::MCModel)
  heap.instate.current = MCBaseSample(copy(x))
  logtarget!(heap.instate.current, m.eval)
end

function initialize_task!(heap::RAMHeap, m::MCModel, s::RAM, r::MCRunner, t::MCTuner)
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, (x::Vector{Float64})->reset!(heap, x, m))

  while true
    iterate!(heap, m, s, r, t, produce)
  end
end

### Perform iteration for RAM sampler

function iterate!(heap::RAMHeap, m::MCModel, s::RAM, r::MCRunner, t::MCTuner, send::Function)
  if t.verbose
    heap.tune.proposed += 1
  end

  heap.rnormal = randn(m.size)
  heap.instate.successive = MCBaseSample(heap.instate.current.sample+heap.S*heap.rnormal)
  logtarget!(heap.instate.successive, m.eval)

  heap.ratio = heap.instate.successive.logtarget-heap.instate.current.logtarget
  if heap.ratio > 0 || (heap.ratio > log(rand()))
    heap.outstate = MCState(heap.instate.successive, heap.instate.current,
    Dict{Any, Any}("accept" => true, "scale" => trace(heap.S)))
    heap.instate.current = deepcopy(heap.instate.successive)

    if t.verbose
      heap.tune.accepted += 1
    end
  else
    heap.outstate = MCState(heap.instate.current, heap.instate.current,
    Dict{Any, Any}("accept" => false, "scale" => trace(heap.S)))
  end

  if t.verbose && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    rate!(heap.tune)
    println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
  end

  heap.eta = min(1, m.size*heap.count^(-2/3))

  heap.SS =
    (heap.rnormal*heap.rnormal')/Base.dot(heap.rnormal, heap.rnormal)*heap.eta*(min(1, exp(heap.ratio))-s.targetrate)
  heap.SS = heap.S*(eye(m.size)+heap.SS)*heap.S'
  heap.S = chol(heap.SS)'

  heap.count += 1

  send(heap.outstate)
end


text:
@test mysignif(SpecificG(20.0, 0.697_992_849E3), TestDigits) ≈ mysignif(SpecificG_Ps(20.0, 5.75), TestDigits)
@test mysignif(SpecificG(20.0u"MPa", 0.697_992_849E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(20.0u"MPa", 5.75u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificG(80.0, 0.854_011_484E3), TestDigits) ≈ mysignif(SpecificG_Ps(80.0, 5.25), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.854_011_484E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(80.0u"MPa", 5.25u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificG(80.0, 0.949_017_998E3), TestDigits) ≈ mysignif(SpecificG_Ps(80.0, 5.75), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.949_017_998E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(80.0u"MPa", 5.75u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(20.0, 0.697_992_849E3), TestDigits) ≈ mysignif(SpecificF_Ps(20.0, 5.75), TestDigits)
@test mysignif(SpecificF(20.0u"MPa", 0.697_992_849E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(20.0u"MPa", 5.75u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(80.0, 0.854_011_484E3), TestDigits) ≈ mysignif(SpecificF_Ps(80.0, 5.25), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.854_011_484E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(80.0u"MPa", 5.25u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(80.0, 0.949_017_998E3), TestDigits) ≈ mysignif(SpecificF_Ps(80.0, 5.75), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.949_017_998E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(80.0u"MPa", 5.75u"kJ/kg/K"), TestDigits)

@test SpecificV(20.0, 0.697_992_849E3) ≈ SpecificV_Ps(20.0, 5.75)
@test SpecificV(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpecificV_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificV(80.0, 0.854_011_484E3) ≈ SpecificV_Ps(80.0, 5.25)
@test SpecificV(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpecificV_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpecificV(80.0, 0.949_017_998E3) ≈ SpecificV_Ps(80.0, 5.75)
@test SpecificV(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpecificV_Ps(80.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificU(20.0, 0.697_992_849E3) ≈ SpecificU_Ps(20.0, 5.75)
@test SpecificU(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpecificU_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificU(80.0, 0.854_011_484E3) ≈ SpecificU_Ps(80.0, 5.25)
@test SpecificU(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpecificU_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpecificU(80.0, 0.949_017_998E3) ≈ SpecificU_Ps(80.0, 5.75)
@test SpecificU(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpecificU_Ps(80.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificH(20.0, 0.697_992_849E3) ≈ SpecificH_Ps(20.0, 5.75)
@test SpecificH(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpecificH_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificH(80.0, 0.854_011_484E3) ≈ SpecificH_Ps(80.0, 5.25)
@test SpecificH(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpecificH_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpecificH(80.0, 0.949_017_998E3) ≈ SpecificH_Ps(80.0, 5.75)
@test SpecificH(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpecificH_Ps(80.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificCP(20.0, 0.697_992_849E3) ≈ SpecificCP_Ps(20.0, 5.75)
@test SpecificCP(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpecificCP_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificCP(80.0, 0.854_011_484E3) ≈ SpecificCP_Ps(80.0, 5.25)
@test SpecificCP(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpecificCP_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpecificCP(80.0, 0.949_017_998E3) ≈ SpecificCP_Ps(80.0, 5.75)
@test SpecificCP(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpecificCP_Ps(80.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificCV(20.0, 0.697_992_849E3) ≈ SpecificCV_Ps(20.0, 5.75)
@test SpecificCV(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpecificCV_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpecificCV(80.0, 0.854_011_484E3) ≈ SpecificCV_Ps(80.0, 5.25)
@test SpecificCV(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpecificCV_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpecificCV(80.0, 0.949_017_998E3) ≈ SpecificCV_Ps(80.0, 5.75)
@test SpecificCV(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpecificCV_Ps(80.0u"MPa", 5.75u"kJ/kg/K")

@test SpeedOfSound(20.0, 0.697_992_849E3) ≈ SpeedOfSound_Ps(20.0, 5.75)
@test SpeedOfSound(20.0u"MPa", 0.697_992_849E3u"K") ≈ SpeedOfSound_Ps(20.0u"MPa", 5.75u"kJ/kg/K")

@test SpeedOfSound(80.0, 0.854_011_484E3) ≈ SpeedOfSound_Ps(80.0, 5.25)
@test SpeedOfSound(80.0u"MPa", 0.854_011_484E3u"K") ≈ SpeedOfSound_Ps(80.0u"MPa", 5.25u"kJ/kg/K")

@test SpeedOfSound(80.0, 0.949_017_998E3) ≈ SpeedOfSound_Ps(80.0, 5.75)
@test SpeedOfSound(80.0u"MPa", 0.949_017_998E3u"K") ≈ SpeedOfSound_Ps(80.0u"MPa", 5.75u"kJ/kg/K")


text:
type Graph{T}
    V::Vector{T}
    E::Vector{Tuple{T, T}}
end

vertices(G::Graph) = G.V
edges(G::Graph) = G.E

function acyclic_path_table{T}(G::Graph{T}, max_path_length::Int = 20)
    paths = [Dict{Vector{T}, Vector{T}}()]

    for E in edges(G)
        push!(get!(paths[1], [E[1]], []), E[2])
        push!(get!(paths[1], [E[2]], []), E[1])
    end

    for i in 1:max_path_length
        push!(paths, Dict{Vector{T}, Vector{T}}())

        for (stem, leaves) in paths[i]
            for leaf in leaves
                new_stem = [stem ; leaf]
                paths[i + 1][new_stem] = T[]
                prefix, suffix = new_stem[1], new_stem[2:end]
                for p in paths[i][suffix]
                    p != prefix && push!(paths[i + 1][new_stem], p)
                end
            end
        end

        isempty(paths[i + 1]) && break
    end

    return paths
end

reversal_invariant(v) = hash(v) ⊻ hash(reverse(v))

unique_paths{T}(paths::Vector{Dict{Vector{T}, Vector{T}}}) = (p -> unique(reversal_invariant, p)).(collect.(keys.(paths)))

acyclic_paths{T}(G::Graph{T}, max_path_length::Int = 20) = unique_paths(acyclic_path_table(G, max_path_length))

function generate_bonds{T}(atoms::Atoms, bond_types::Dict{Tuple{Symbol, Symbol}, Tuple{T, T}})
    vertices = Set{Atom}()
    edges = Vector{Tuple{Atom, Atom}}()

    coords = atoms[:coords]
    types = atoms[:type]

    n = size(coords, 2)::Int

    for i in 2:n, j in 1:i-1
        type_i, type_j = types[i]::Symbol, types[j]::Symbol

        dist_min, dist_max = get(bond_types, (type_i, type_j), get(bond_types, (type_j, type_i), (Inf, -Inf)))
        if dist_min != Inf
            if dist_min^2 <= distancesq(Inf, coords, i, j) <= dist_max^2
                push!(vertices, atoms[i], atoms[j])
                push!(edges, (atoms[i], atoms[j]))
            end
        end
    end

    return Graph(collect(vertices), edges)
end

function unique_upto_reversal{T}(v::Vector{Vector{T}})
    unique = Set{Vector{T}}()
    for a in v
        if a ∉ unique && reverse(a) ∉ unique
            push!(unique, a)
        end
    end
    collect(unique)
end

function unique_neighborhoods(path_table::Vector{Dict{Vector{Atom}, Vector{Atom}}})
    neighborhoods = Dict{Symbol, Dict{Vector{Symbol}, Vector{Atom}}}()
    for (atom, neighbors) in path_table[1]
        atom_type = atom[1][:type]
        neighborhood = sort(Atoms(neighbors)[:type])
        push!(get!(get!(neighborhoods, atom_type, Dict{Vector{Symbol}, Vector{Atom}}()), neighborhood, Symbol[]), atom[1])
    end
    neighborhoods
end

function neighborhood_type_ids(neighborhoods)
    Dict(
        atom_type => Dict(
            first(neighborhood_type) => i
        for (i, neighborhood_type) in enumerate(neighborhood_types))
    for (atom_type, neighborhood_types) in neighborhoods)
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
#
println("Ha) Tests of the HighHarmonic module to calculate a high-harmonic spectrum in single-electron approximation.")
#
setDefaults("print summary: open", "zzz-HHG.sum")
grid       = JAC.Radial.Grid(true)
setDefaults("standard grid", grid::Radial.Grid; printout=true)

observable = HighHarmonic.HhgSpectrum()
approach   = HighHarmonic.HhgHydrogenicSaddlePoint()
pulse      = HighHarmonic.HhgLinearPlaneWaveLaser(2.0, intensity=3.51e16)   ## omega = 2.0 eV; intensity = 10^14 W/cm^2
timeMesh   = HighHarmonic.computeTimeMesh(2, 50)                                    ## 2 cycles with 50 time points / cycle
initialOrb = Radial.Orbital(Subshell("1s_1/2"), 0.5)
target     = HighHarmonic.TargetCloud( HighHarmonic.LocalizedTargetModel( CartesianVector(0., 0., 0.)), WeightedCartesian[])
detector   = HighHarmonic.DetectorScreen( HighHarmonic.LocalizedDetectorModel( CartesianVector(0., 0., 0.)), WeightedCartesian[])

wa = HighHarmonic.Computation( observable, approach, pulse, timeMesh, Nuclear.Model(2.0), initialOrb, target, detector)

wb = HighHarmonic.perform(wa; output=true)
setDefaults("print summary: close", "")


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks
import SciMLBase: RightRootFind

using FMIImport: fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2NewDiscreteStates, fmi2GetContinuousStates, fmi2GetNominalsOfContinuousStates, fmi2SetContinuousStates, fmi2GetDerivatives!
using FMIImport.FMICore: fmi2StatusOK, fmi2TypeCoSimulation, fmi2TypeModelExchange
using FMIImport.FMICore: fmi2ComponentState, fmi2ComponentStateInstantiated, fmi2ComponentStateInitializationMode, fmi2ComponentStateEventMode, fmi2ComponentStateContinuousTimeMode, fmi2ComponentStateTerminated, fmi2ComponentStateError, fmi2ComponentStateFatal
using FMIImport: FMU2Solution, FMU2Event

using ChainRulesCore
import ForwardDiff

import ProgressMeter

############ Model-Exchange ############

# Read next time event from fmu and provide it to the integrator
function time_choice(c::FMU2Component, integrator, tStart, tStop)

    if c.eventInfo.nextEventTimeDefined == fmi2True

        if c.eventInfo.nextEventTime >= tStart && c.eventInfo.nextEventTime <= tStop
            return c.eventInfo.nextEventTime
        else
            # the time event is outside the simulation range!
            @debug "Next time event @$(c.eventInfo.nextEventTime)s is outside simulation time range ($(tStart), $(tStop)), skipping."
            return nothing
        end
    else
        return nothing
    end

end

# Handles events and returns the values and nominals of the changed continuous states.
function handleEvents(c::FMU2Component)

    @assert c.state == fmi2ComponentStateEventMode "handleEvents(...): Must be in event mode!"

    # trigger the loop
    c.eventInfo.newDiscreteStatesNeeded = fmi2True

    valuesOfContinuousStatesChanged = fmi2False
    nominalsOfContinuousStatesChanged = fmi2False
    nextEventTimeDefined = fmi2False
    nextEventTime = 0.0

    while c.eventInfo.newDiscreteStatesNeeded == fmi2True
        fmi2NewDiscreteStates!(c, c.eventInfo)

        if c.eventInfo.valuesOfContinuousStatesChanged == fmi2True
            valuesOfContinuousStatesChanged = fmi2True
        end

        if c.eventInfo.nominalsOfContinuousStatesChanged == fmi2True
            nominalsOfContinuousStatesChanged = fmi2True
        end

        if c.eventInfo.nextEventTimeDefined == fmi2True
            nextEventTimeDefined = fmi2True
            nextEventTime = c.eventInfo.nextEventTime
        end

        if c.eventInfo.terminateSimulation == fmi2True
            @error "handleEvents(...): FMU throws `terminateSimulation`!"
        end
    end

    c.eventInfo.valuesOfContinuousStatesChanged = valuesOfContinuousStatesChanged
    c.eventInfo.nominalsOfContinuousStatesChanged = nominalsOfContinuousStatesChanged
    c.eventInfo.nextEventTimeDefined = nextEventTimeDefined
    c.eventInfo.nextEventTime = nextEventTime

    fmi2EnterContinuousTimeMode(c)

    return nothing
end

# Returns the event indicators for an FMU.
function condition(c::FMU2Component, out::SubArray{<:Real}, x, t, integrator, inputFunction, inputValues::Array{fmi2ValueReference})

    @assert c.state == fmi2ComponentStateContinuousTimeMode "condition(...): Must be called in mode continuous time."

    fmi2SetContinuousStates(c, x)
    fmi2SetTime(c, integrator.t)
    if inputFunction !== nothing
        fmi2SetReal(c, inputValues, inputFunction(integrator.t))
    end
    fmi2GetEventIndicators!(c, out)

    return nothing
end

# Handles the upcoming events.
# Sets a new state for the solver from the FMU (if needed).
function affectFMU!(c::FMU2Component, integrator, idx, inputFunction, inputValues::Array{fmi2ValueReference}, solution::FMU2Solution)

    @assert c.state == fmi2ComponentStateContinuousTimeMode "affectFMU!(...): Must be in continuous time mode!"

    # there are fx-evaluations before the event is handled, reset the FMU state to the current integrator step
    fmi2SetContinuousStates(c, integrator.u)
    fmi2SetTime(c, integrator.t)
    if inputFunction !== nothing
        fmi2SetReal(c, inputValues, inputFunction(integrator.t))
    end

    fmi2EnterEventMode(c)

    # Event found - handle it
    handleEvents(c)

    left_x = nothing
    right_x = nothing

    if c.eventInfo.valuesOfContinuousStatesChanged == fmi2True
        left_x = integrator.u
        right_x = fmi2GetContinuousStates(c)
        @debug "affectFMU!(...): Handled event at t=$(integrator.t), new state is $(new_u)"
        integrator.u = right_x
    else
        @debug "affectFMU!(...): Handled event at t=$(integrator.t), no new state."
    end

    if c.eventInfo.nominalsOfContinuousStatesChanged == fmi2True
        x_nom = fmi2GetNominalsOfContinuousStates(c)
    end

    ignore_derivatives() do
        if idx != -1 # -1 no event, 0, time event, >=1 state event with indicator
            e = FMU2Event(integrator.t, UInt64(idx), left_x, right_x)
            push!(solution.events, e)
        end
    end

    #fmi2EnterContinuousTimeMode(c)
end

# This callback is called every time the integrator finishes an (accpeted) integration step.
function stepCompleted(c::FMU2Component, x, t, integrator, inputFunction, inputValues::Array{fmi2ValueReference}, progressMeter, tStart, tStop, solution::FMU2Solution)

    @assert c.state == fmi2ComponentStateContinuousTimeMode "stepCompleted(...): Must be in continuous time mode."

    if progressMeter !== nothing
        ProgressMeter.update!(progressMeter, floor(Integer, 1000.0*(t-tStart)/(tStop-tStart)) )
    end

    (status, enterEventMode, terminateSimulation) = fmi2CompletedIntegratorStep(c, fmi2True)

    if terminateSimulation == fmi2True
        @error "stepCompleted(...): FMU requested termination!"
    end

    if enterEventMode == fmi2True
        affectFMU!(c, integrator, -1, inputFunction, inputValues, solution)
    else
        if inputFunction != nothing
            fmi2SetReal(c, inputValues, inputFunction(integrator.t))
        end
    end
end

# save FMU values
function saveValues(c::FMU2Component, recordValues, x, t, integrator)

    @assert c.state == fmi2ComponentStateContinuousTimeMode "saveValues(...): Must be in continuous time mode."

    fmi2SetContinuousStates(c, x)
    fmi2SetTime(c, t)

    return (fmiGetReal(c, recordValues)...,)
end

function fx(c::FMU2Component,
    dx::Array{<:Real},
    x::Array{<:Real},
    p::Array,
    t::Real)

    if isa(t, ForwardDiff.Dual)
        t = ForwardDiff.value(t)
    end

    fmi2SetContinuousStates(c, x)

    if all(isa.(dx, ForwardDiff.Dual))
        dx_tmp = collect(ForwardDiff.value(e) for e in dx)
        fmi2GetDerivatives!(c, dx_tmp)
        T, V, N = fd_eltypes(dx)
        dx[:] = collect(ForwardDiff.Dual{T, V, N}(dx_tmp[i], ForwardDiff.partials(dx[i])    ) for i in 1:length(dx))
    else
        fmi2GetDerivatives!(c, dx)
    end

    return dx
end

# ForwardDiff-Dispatch for fx
function fx(comp::FMU2Component,
            dx::Array{<:Real},
            x::Array{<:ForwardDiff.Dual{Tx, Vx, Nx}},
            p::Array,
            t::Real) where {Tx, Vx, Nx}

    return _fx_fd(comp, dx, x, p, t)
end

# function _fx_fd(TVNx, comp, dx, x, p, t)

#     Tx, Vx, Nx = TVNx

#     ȧrgs = [NoTangent(), NoTangent(), collect(ForwardDiff.partials(e) for e in dx), collect(ForwardDiff.partials(e) for e in x), collect(ForwardDiff.partials(e) for e in p), ForwardDiff.partials(t)]
#     args = [fx,          comp,        collect(ForwardDiff.value(e) for e in dx), collect(ForwardDiff.value(e) for e in x),    collect(ForwardDiff.value(e) for e in p),    ForwardDiff.value(t),  ]

#     ȧrgs = (ȧrgs...,)
#     args = (args...,)

#     y, _, sdx, sx, _, _ = ChainRulesCore.frule(ȧrgs, args...)

#     if Vx != Float64
#         Vx = Float64
#     end

#     [collect( ForwardDiff.Dual{Tx, Vx, Nx}(y[i], sx[i]) for i in 1:length(sx) )...]
# end

function fd_eltypes(e::ForwardDiff.Dual{T, V, N}) where {T, V, N}
    return (T, V, N)
end
function fd_eltypes(e::Array{<:ForwardDiff.Dual{T, V, N}}) where {T, V, N}
    return (T, V, N)
end
function _fx_fd(comp, dx, x, p, t)

    ȧrgs = []
    args = []

    push!(ȧrgs, NoTangent())
    push!(args, fx)

    push!(ȧrgs, NoTangent())
    push!(args, comp)

    T = nothing
    V = nothing

    dx_set = length(dx) > 0 && all(isa.(dx, ForwardDiff.Dual))
    x_set = length(x) > 0 && all(isa.(x, ForwardDiff.Dual))
    p_set = length(p) > 0 && all(isa.(p, ForwardDiff.Dual))
    t_set = isa(t, ForwardDiff.Dual)

    if dx_set
        T, V, N = fd_eltypes(dx)
        push!(ȧrgs, collect(ForwardDiff.partials(e) for e in dx))
        push!(args, collect(ForwardDiff.value(e) for e in dx))
        #@info "dx_set num=$(length(dx)) partials=$(length(ForwardDiff.partials(dx[1])))"
    else
        push!(ȧrgs, NoTangent())
        push!(args, dx)
    end

    if x_set
        T, V, N = fd_eltypes(x)
        push!(ȧrgs, collect(ForwardDiff.partials(e) for e in x))
        push!(args, collect(ForwardDiff.value(e) for e in x))
        #@info "x_set num=$(length(x)) partials=$(length(ForwardDiff.partials(x[1])))"
    else
        push!(ȧrgs, NoTangent())
        push!(args, x)
    end

    if p_set
        T, V, N = fd_eltypes(p)
        push!(ȧrgs, collect(ForwardDiff.partials(e) for e in p))
        push!(args, collect(ForwardDiff.value(e) for e in p))
    else
        push!(ȧrgs, NoTangent())
        push!(args, p)
    end

    if t_set
        T, V, N = fd_eltypes(t)
        push!(ȧrgs, ForwardDiff.partials(t))
        push!(args, ForwardDiff.value(t))
    else
        push!(ȧrgs, NoTangent())
        push!(args, t)
    end

    ȧrgs = (ȧrgs...,)
    args = (args...,)

    y, _, sdx, sx, sp, st = ChainRulesCore.frule(ȧrgs, args...)

    ys = []

    #[collect( ForwardDiff.Dual{Tx, Vx, Nx}(y[i], ForwardDiff.partials(x_partials[i], t_partials[i])) for i in 1:length(y) )...]
    for i in 1:length(y)
        is = NoTangent()

        if dx_set
            is = sdx[i]#.values
        end
        if x_set
            is = sx[i]#.values
        end

        if p_set
            is = sp[i]#.values
        end
        if t_set
            is = st[i]#.values
        end

        #display("dx: $dx")
        #display("sdx: $sdx")

        #partials = (isdx, isx, isp, ist)

        #display(partials)

        #V = Float64
        #N = length(partials)
        #display("$T $V $N")

        #display(is)

        @assert is != ZeroTangent() && is != NoTangent() "is: $(is)"

        push!(ys, ForwardDiff.Dual{T, V, N}(y[i], is    )   ) #  ForwardDiff.Partials{N, V}(partials)
    end

    ys
end

# frule for fx
function ChainRulesCore.frule((Δself, Δcomp, Δdx, Δx, Δp, Δt),
                              ::typeof(fx),
                              comp, #::FMU2Component,
                              dx,
                              x,#::Array{<:Real},
                              p,
                              t)

    y = fx(comp, dx, x, p, t)
    function fx_pullforward(Δdx, Δx, Δt)

        # if t >= 0.0
        #     fmi2SetTime(comp, t)
        # end

        # if all(isa.(x, ForwardDiff.Dual))
        #     xf = collect(ForwardDiff.value(e) for e in x)
        #     fmi2SetContinuousStates(comp, xf)
        # else
        #     fmi2SetContinuousStates(comp, x)
        # end

        c̄omp = ZeroTangent()
        d̄x = ZeroTangent()
        x̄ = ZeroTangent()
        p̄ = ZeroTangent()
        t̄ = ZeroTangent()

        if Δdx != NoTangent()
            d̄x = Δdx
        end

        if Δx != NoTangent()
            if comp.A == nothing || size(comp.A) != (length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
                comp.A = zeros(length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
            end
            comp.jacobianUpdate!(comp.A, comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)
            x̄ = comp.A * Δx
        end

        if Δt != NoTangent()
            dt = 1e-6
            dx1 = fmi2GetDerivatives(comp)
            fmi2SetTime(comp, t + dt)
            dx2 = fmi2GetDerivatives(comp)
            ∂t = (dx2-dx1)/dt
            t̄ = ∂t * Δt
        end

        return (c̄omp, d̄x, x̄, p̄, t̄)
    end
    return (y, fx_pullforward(Δdx, Δx, Δt)...)
end

# rrule for fx
function ChainRulesCore.rrule(::typeof(fx),
                              comp::FMU2Component,
                              dx,
                              x,
                              p,
                              t)

    y = fx(comp, dx, x, p, t)
    function fx_pullback(ȳ)

        if t >= 0.0
            fmi2SetTime(comp, t)
        end

        fmi2SetContinuousStates(comp, x)

        if comp.A == nothing || size(comp.A) != (length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
            comp.A = zeros(length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
        end
        comp.jacobianUpdate!(comp.A, comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)

        n_dx_x = @thunk(comp.A' * ȳ)

        f̄ = NoTangent()
        c̄omp = ZeroTangent()
        d̄x = ZeroTangent()
        x̄ = n_dx_x
        p̄ = ZeroTangent()
        t̄ = ZeroTangent()

        return f̄, c̄omp, d̄x, x̄, p̄, t̄
    end
    return (y, fx_pullback)
end

# frule for fx
# function ChainRulesCore.frule((Δself, Δcomp, Δdx, Δx, Δp, Δt),
#                               ::typeof(fx),
#                               comp, #::FMU2Component,
#                               dx,
#                               x,#::Array{<:Real},
#                               p,
#                               t)

#     y = fx(comp, dx, x, p, t)
#     function fx_pullforward(Δx)

#         if t >= 0.0
#             fmi2SetTime(comp, t)
#         end

#         if all(isa.(x, ForwardDiff.Dual))
#             xf = collect(ForwardDiff.value(e) for e in x)
#             fmi2SetContinuousStates(comp, xf)
#         else
#             fmi2SetContinuousStates(comp, x)
#         end

#         if comp.A == nothing || size(comp.A) != (length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
#             comp.A = zeros(length(comp.fmu.modelDescription.derivativeValueReferences), length(comp.fmu.modelDescription.stateValueReferences))
#         end
#         comp.jacobianUpdate!(comp.A, comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)

#         n_dx_x = comp.A * Δx

#         c̄omp = ZeroTangent()
#         d̄x = ZeroTangent()
#         x̄ = n_dx_x
#         p̄ = ZeroTangent()
#         t̄ = ZeroTangent()

#         return (c̄omp, d̄x, x̄, p̄, t̄)
#     end
#     return (y, fx_pullforward(Δx)...)
# end

function prepareFMU(fmu::FMU2, c::Union{Nothing, FMU2Component}, instantiate::Union{Nothing, Bool}, terminate::Union{Nothing, Bool}, reset::Union{Nothing, Bool}, setup::Union{Nothing, Bool}, parameters::Union{Dict{<:Any, <:Any}, Nothing}, t_start, t_stop, tolerance;
    x0::Union{Array{<:Real}, Nothing}=nothing)

    if instantiate === nothing
        instantiate = fmu.executionConfig.instantiate
    end

    if terminate === nothing
        terminate = fmu.executionConfig.terminate
    end

    if reset === nothing
        reset = fmu.executionConfig.reset
    end

    if setup === nothing
        setup = fmu.executionConfig.setup
    end

    c = nothing

    # instantiate (hard)
    if instantiate
        c = fmi2Instantiate!(fmu)
    else
        if c === nothing
            c = fmu.components[end]
        end
    end

    # soft terminate (if necessary)
    if terminate
        retcode = fmi2Terminate(c; soft=true)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Termination failed with return code $(retcode)."
    end

    # soft reset (if necessary)
    if reset
        retcode = fmi2Reset(c; soft=true)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Reset failed with return code $(retcode)."
    end

    # enter setup (hard)
    if setup
        retcode = fmi2SetupExperiment(c, t_start, t_stop; tolerance=tolerance)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Setting up experiment failed with return code $(retcode)."

        retcode = fmi2EnterInitializationMode(c)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Entering initialization mode failed with return code $(retcode)."
    end

    if x0 !== nothing
        retcode = fmi2SetContinuousStates(c, x0)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Setting initial state failed with return code $(retcode)."
    end

    if parameters !== nothing
        retcodes = fmi2Set(c, collect(keys(parameters)), collect(values(parameters)) )
        @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial parameters failed with return code $(retcode)."
    end

    # exit setup (hard)
    if setup
        retcode = fmi2ExitInitializationMode(c)
        @assert retcode == fmi2StatusOK "fmi2Simulate(...): Exiting initialization mode failed with return code $(retcode)."
    end

    return c
end

function finishFMU(fmu::FMU2, c::FMU2Component, freeInstance::Union{Nothing, Bool})
    if freeInstance === nothing
        freeInstance = fmu.executionConfig.freeInstance
    end

    # freeInstance (hard)
    if freeInstance
        fmi2FreeInstance!(c)
    end
end

# wrapper
function fmi2SimulateME(c::FMU2Component, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing; kwargs...)
    fmi2SimulateME(c.fmu, c, t_start, t_stop; kwargs...)
end

"""
Simulates a FMU instance for the given simulation time interval.
State- and Time-Events are handled correctly.

Via the optional keyword arguemnts `inputValues` and `inputFunction`, a custom input function of the time `t` can be defined, that should return a array of values for `fmi2SetReal(..., inputValues, inputFunction(t))`.

Keywords:
    - solver: Any Julia-supported ODE-solver (default is Tsit5)
    - customFx: [deperecated] Ability to give a custom state derivative function ẋ=f(x,t)
    - recordValues: Array of variables (strings or variableIdentifiers) to record. Results are returned as `DiffEqCallbacks.SavedValues`
    - saveat: Time points to save values at (interpolated)
    - setup: Boolean, if FMU should be setup (default=true)
    - reset: Union{Bool, :auto}, if FMU should be reset before simulation (default reset=:auto)
    - inputValueReferences: Array of input variables (strings or variableIdentifiers) to set at every simulation step
    - inputFunction: Function to retrieve the values to set the inputs to
    - parameters: Dictionary of parameter variables (strings or variableIdentifiers) and values (Real, Integer, Boolean, String) to set parameters during initialization
    - `callbacks`: custom callbacks to add

Returns:
    - If keyword `recordValues` is not set, a struct of type `ODESolution`.
    - If keyword `recordValues` is set, a tuple of type (ODESolution, DiffEqCallbacks.SavedValues).
"""
function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing;
    tolerance::Union{Real, Nothing} = nothing,
    dt::Union{Real, Nothing} = nothing,
    solver = nothing,
    customFx = nothing,
    recordValues::fmi2ValueReferenceFormat = nothing,
    saveat = [],
    x0::Union{Array{<:Real}, Nothing} = nothing,
    setup::Union{Bool, Nothing} = nothing,
    reset::Union{Bool, Nothing} = nothing,
    instantiate::Union{Bool, Nothing} = nothing,
    freeInstance::Union{Bool, Nothing} = nothing,
    terminate::Union{Bool, Nothing} = nothing,
    inputValueReferences::fmi2ValueReferenceFormat = nothing,
    inputFunction = nothing,
    parameters::Union{Dict{<:Any, <:Any}, Nothing} = nothing,
    dtmax::Union{Real, Nothing} = nothing,
    callbacks = [],
    showProgress::Bool = true,
    kwargs...)

    @assert fmi2IsModelExchange(fmu) "fmi2SimulateME(...): This function supports Model Excahnge FMUs only."

    recordValues = prepareValueReference(fmu, recordValues)
    inputValueReferences = prepareValueReference(fmu, inputValueReferences)

    fmusol = FMU2Solution(fmu)

    savingValues = (length(recordValues) > 0)
    hasInputs = (length(inputValueReferences) > 0)
    hasParameters = (parameters !== nothing)
    hasStartState = (x0 !== nothing)

    cbs = []

    for cb in callbacks
        push!(cbs, cb)
    end

    t_start = t_start === nothing ? fmi2GetDefaultStartTime(fmu.modelDescription) : t_start
    t_start = t_start === nothing ? 0.0 : t_start
    t_stop = t_stop === nothing ? fmi2GetDefaultStopTime(fmu.modelDescription) : t_stop
    t_stop = t_stop === nothing ? 1.0 : t_stop
    tolerance = tolerance === nothing ? fmi2GetDefaultTolerance(fmu.modelDescription) : tolerance
    tolerance = tolerance === nothing ? 1e-4 : tolerance
    dt = dt === nothing ? fmi2GetDefaultStepSize(fmu.modelDescription) : dt
    dt = dt === nothing ? 1e-5 : dt

    if dtmax === nothing
        dtmax = (t_stop-t_start)/100.0
    end

    if solver === nothing
        solver = Tsit5()
    end

    c = prepareFMU(fmu, c, instantiate, terminate, reset, setup, parameters, t_start, t_stop, tolerance; x0=x0)

    # from here on, we are in event mode, if `setup=false` this is the job of the user
    @assert c.state == fmi2ComponentStateEventMode "FMU needs to be in event mode after setup."

    if x0 === nothing
        x0 = fmi2GetContinuousStates(c)
        x0_nom = fmi2GetNominalsOfContinuousStates(c)
    end

    # initial event handling
    handleEvents(c)
    #fmi2EnterContinuousTimeMode(c)

    c.fmu.hasStateEvents = (c.fmu.modelDescription.numberOfEventIndicators > 0)
    c.fmu.hasTimeEvents = (c.eventInfo.nextEventTimeDefined == fmi2True)

    if customFx === nothing
        customFx = (dx, x, p, t) -> fx(c, dx, x, p, t)
    end

    p = []
    problem = ODEProblem(customFx, x0, (t_start, t_stop), p,)

    progressMeter = nothing
    if showProgress
        progressMeter = ProgressMeter.Progress(1000; desc="Simulating ME-FMU ...", color=:blue, dt=1.0) #, barglyphs=ProgressMeter.BarGlyphs("[=> ]"))
        ProgressMeter.update!(progressMeter, 0) # show it!
    end

    # callback functions

    if savingValues
        fmusol.values = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})
        fmusol.valueReferences = copy(recordValues)

        savingCB = SavingCallback((u,t,integrator) -> saveValues(c, recordValues, u, t, integrator),
                                  fmusol.values,
                                  saveat=saveat)
        push!(cbs, savingCB)
    end

    # use step callback always if we have inputs or need event handling (or just want to see our simulation progress)
    if hasInputs || c.fmu.hasStateEvents || c.fmu.hasTimeEvents || showProgress
        stepCb = FunctionCallingCallback((x, t, integrator) -> stepCompleted(c, x, t, integrator, inputFunction, inputValueReferences, progressMeter, t_start, t_stop, fmusol);
                                            func_everystep = true,
                                            func_start = true)
        push!(cbs, stepCb)
    end

    if c.fmu.hasStateEvents

        eventCb = VectorContinuousCallback((out, x, t, integrator) -> condition(c, out, x, t, integrator, inputFunction, inputValueReferences),
                                           (integrator, idx) -> affectFMU!(c, integrator, idx, inputFunction, inputValueReferences, fmusol),
                                           Int64(c.fmu.modelDescription.numberOfEventIndicators);
                                           rootfind = RightRootFind,
                                           save_positions=(false,false))
        push!(cbs, eventCb)
    end

    if c.fmu.hasTimeEvents
        timeEventCb = IterativeCallback((integrator) -> time_choice(c, integrator, t_start, t_stop),
                                        (integrator) -> affectFMU!(c, integrator, 0, inputFunction, inputValueReferences, fmusol), Float64;
                                        initial_affect = (c.eventInfo.nextEventTime == t_start),
                                        save_positions=(false,false))
        push!(cbs, timeEventCb)
    end

    fmusol.states = solve(problem, solver; callback = CallbackSet(cbs...), saveat = saveat, reltol=tolerance, dt=dt, dtmax=dtmax, kwargs...)
    fmusol.success = (fmusol.states.retcode == :Success)

    # cleanup progress meter
    if showProgress
        ProgressMeter.finish!(progressMeter)
    end

    finishFMU(fmu, c, freeInstance)

    return fmusol
end

# wrapper
function fmi2SimulateCS(c::FMU2Component, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing; kwargs...)
    fmi2SimulateCS(c.fmu, c, t_start, t_stop; kwargs...)
end

############ Co-Simulation ############

"""
Starts a simulation of the Co-Simulation FMU instance.

Via the optional keyword arguments `inputValues` and `inputFunction`, a custom input function of the time `t` can be defined, that should return a array of values for `fmi2SetReal(..., inputValues, inputFunction(t))`.

Keywords:
    - recordValues: Array of variables (strings or variableIdentifiers) to record. Results are returned as `DiffEqCallbacks.SavedValues`
    - saveat: Time points to save values at (interpolated)
    - setup: Boolean, if FMU should be setup (default=true)
    - reset: Boolean, if FMU should be reset before simulation (default reset=setup)
    - inputValueReferences: Array of input variables (strings or variableIdentifiers) to set at every simulation step
    - inputFunction: Function to retrieve the values to set the inputs to
    - parameters: Dictionary of parameter variables (strings or variableIdentifiers) and values (Real, Integer, Boolean, String) to set parameters during initialization
Returns:
    - If keyword `recordValues` is not set, a boolean `success` is returned (simulation success).
    - If keyword `recordValues` is set, a tuple of type (true, DiffEqCallbacks.SavedValues) or (false, nothing).
"""
function fmi2SimulateCS(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing;
                        tolerance::Union{Real, Nothing} = nothing,
                        dt::Union{Real, Nothing} = nothing,
                        recordValues::fmi2ValueReferenceFormat = nothing,
                        saveat = [],
                        setup::Union{Bool, Nothing} = nothing,
                        reset::Union{Bool, Nothing} = nothing,
                        instantiate::Union{Bool, Nothing} = nothing,
                        freeInstance::Union{Bool, Nothing} = nothing,
                        terminate::Union{Bool, Nothing} = nothing,
                        inputValueReferences::fmi2ValueReferenceFormat = nothing,
                        inputFunction = nothing,
                        showProgress::Bool=true,
                        parameters::Union{Dict{<:Any, <:Any}, Nothing} = nothing)

    fmusol = FMU2Solution(fmu)

    recordValues = prepareValueReference(fmu, recordValues)
    inputValueReferences = prepareValueReference(fmu, inputValueReferences)

    variableSteps = fmi2IsCoSimulation(fmu) && fmu.modelDescription.coSimulation.canHandleVariableCommunicationStepSize

    t_start = t_start === nothing ? fmi2GetDefaultStartTime(fmu.modelDescription) : t_start
    t_start = t_start === nothing ? 0.0 : t_start
    t_stop = t_stop === nothing ? fmi2GetDefaultStopTime(fmu.modelDescription) : t_stop
    t_stop = t_stop === nothing ? 1.0 : t_stop
    tolerance = tolerance === nothing ? fmi2GetDefaultTolerance(fmu.modelDescription) : tolerance
    tolerance = tolerance === nothing ? 0.0 : tolerance
    dt = dt === nothing ? fmi2GetDefaultStepSize(fmu.modelDescription) : dt
    dt = dt === nothing ? 1e-3 : dt

    c = prepareFMU(fmu, c, instantiate, terminate, reset, setup, parameters, t_start, t_stop, tolerance)

    # default setup
    if length(saveat) == 0
        saveat = t_start:dt:t_stop
    end

    # setup if no variable steps
    if variableSteps == false
        if length(saveat) >= 2
            dt = saveat[2] - saveat[1]
        end
    end

    t = t_start

    record = length(recordValues) > 0

    progressMeter = nothing
    if showProgress
        progressMeter = ProgressMeter.Progress(1000; desc="Simulating CS-FMU ...", color=:blue, dt=1.0) #, barglyphs=ProgressMeter.BarGlyphs("[=> ]"))
        ProgressMeter.update!(progressMeter, 0) # show it!
    end

    if record
        fmusol.values = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...} )
        fmusol.valueReferences = copy(recordValues)

        i = 1

        svalues = (fmi2GetReal(c, recordValues)...,)
        DiffEqCallbacks.copyat_or_push!(fmusol.values.t, i, t)
        DiffEqCallbacks.copyat_or_push!(fmusol.values.saveval, i, svalues, Val{false})

        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            if inputFunction != nothing
                fmi2SetReal(c, inputValueReferences, inputFunction(t))
            end

            fmi2DoStep(c, dt; currentCommunicationPoint=t)
            t = t + dt
            i += 1

            svalues = (fmi2GetReal(c, recordValues)...,)
            DiffEqCallbacks.copyat_or_push!(fmusol.values.t, i, t)
            DiffEqCallbacks.copyat_or_push!(fmusol.values.saveval, i, svalues, Val{false})

            if progressMeter !== nothing
                ProgressMeter.update!(progressMeter, floor(Integer, 1000.0*(t-t_start)/(t_stop-t_start)) )
            end
        end

        if progressMeter !== nothing
            ProgressMeter.finish!(progressMeter)
        end

        fmusol.success = true

    else
        i = 1
        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            if inputFunction != nothing
                fmi2SetReal(c, inputValueReferences, inputFunction(t))
            end

            fmi2DoStep(c, dt; currentCommunicationPoint=t)
            t = t + dt
            i += 1

            if progressMeter !== nothing
                ProgressMeter.update!(progressMeter, floor(Integer, 1000.0*(t-t_start)/(t_stop-t_start)) )
            end
        end

        if progressMeter !== nothing
            ProgressMeter.finish!(progressMeter)
        end

        fmusol.success = true
    end

    finishFMU(fmu, c, freeInstance)

    return fmusol
end

##### CS & ME #####

# wrapper
function fmi2Simulate(c::FMU2Component, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing; kwargs...)
    fmi2Simulate(c.fmu, c, t_start, t_stop; kwargs...)
end

"""
Starts a simulation of the FMU instance for the matching FMU type, if both types are available, CS is preferred.

Keywords:
    - recordValues: Array of variables (strings or variableIdentifiers) to record. Results are returned as `DiffEqCallbacks.SavedValues`
    - setup: Boolean, if FMU should be setup (default=true)
    - reset: Boolean, if FMU should be reset before simulation (default reset=setup)
    - inputValues: Array of input variables (strings or variableIdentifiers) to set at every simulation step
    - inputFunction: Function to retrieve the values to set the inputs to
    - saveat: [ME only] Time points to save values at (interpolated)
    - solver: [ME only] Any Julia-supported ODE-solver (default is Tsit5)
    - customFx: [ME only, deperecated] Ability to give a custom state derivative function ẋ=f(x,t)

Returns:
    - `success::Bool` for CS-FMUs
    - `ODESolution` for ME-FMUs
    - if keyword `recordValues` is set, a tuple of type (success::Bool, DiffEqCallbacks.SavedValues) for CS-FMUs
    - if keyword `recordValues` is set, a tuple of type (ODESolution, DiffEqCallbacks.SavedValues) for ME-FMUs
"""
function fmi2Simulate(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing; kwargs...)

    if fmu.type == fmi2TypeCoSimulation
        return fmi2SimulateCS(fmu, c, t_start, t_stop; kwargs...)
    elseif fmu.type == fmi2TypeModelExchange
        return fmi2SimulateME(fmu, c, t_start, t_stop; kwargs...)
    else
        error(unknownFMUType)
    end
end


text:
using AdversarialAttacks
using Test

@testset "AdversarialAttacks.jl" begin
    # Write your tests here.
end


text:
using AWS

@service S3

struct S3Store <: AbstractStore
    bucket::String
    listversion::Int
    aws::AWS.AbstractAWSConfig
end

function S3Store(bucket::String;
  listversion = 2,
  aws = nothing,
  )
  if aws === nothing
    aws = AWS.global_aws_config()
  end
  S3Store(bucket, listversion, aws)
end

Base.show(io::IO,::S3Store) = print(io,"S3 Object Storage")

function error_is_ignorable(e)
  #isa(e,AWSException) && (e.code=="NoSuchKey" || e.code=="404")
  true
end

function Base.getindex(s::S3Store, i::String)
  try
    return S3.get_object(s.bucket,i,aws_config=s.aws)
  catch e
    if error_is_ignorable(e)
      return nothing
    else
      throw(e)
    end
  end
end

function Base.setindex!(s::S3Store, v, i::String)
  return S3.put_object(s.bucket,i,Dict("body"=>v),aws_config=s.aws)
end

Base.delete!(s::S3Store, d::String) = S3.delete_object(s.bucket,d, aws_config=s.aws)

function storagesize(s::S3Store,p)
  r = cloud_list_objects(s,p)
  haskey(r,"Contents") || return 0
  contents = r["Contents"]
  datafiles = filter(entry -> !any(filename -> endswith(entry["Key"], filename), [".zattrs",".zarray",".zgroup"]), contents)
  if isempty(datafiles)
    0
  else
    sum(datafiles) do f
      parse(Int, f["Size"])
    end
  end
end

function isinitialized(s::S3Store, i::String)
  try
    S3.head_object(s.bucket,i,aws_config=s.aws)
    return true
  catch e
    if error_is_ignorable(e)
      return false
    else
      println(i)
      throw(e)
    end
  end
end

function cloud_list_objects(s::S3Store,p)
  prefix = (isempty(p) || endswith(p,"/")) ? p : string(p,"/")
  listfun = s.listversion==2 ? S3.list_objects_v2 : S3.list_objects
  listfun(s.bucket, Dict("prefix"=>prefix, "delimiter" => "/"), aws_config = s.aws)
end
function subdirs(s::S3Store, p)
  s3_resp = cloud_list_objects(s, p)
  !haskey(s3_resp,"CommonPrefixes") && return String[]
  allstrings(s3_resp["CommonPrefixes"],"Prefix")
end
function subkeys(s::S3Store, p)
  s3_resp = cloud_list_objects(s, p)
  !haskey(s3_resp,"Contents") && return String[]
  r = allstrings(s3_resp["Contents"],"Key")
  map(i->splitdir(i)[2],r)
end
allstrings(v::AbstractArray,prefixkey) = map(i -> rstrip(String(i[prefixkey]),'/'), v)
allstrings(v,prefixkey) = [rstrip(String(v[prefixkey]),'/')]

# Some special AWS configs
struct AnonymousGCS <:AbstractAWSConfig end
#struct NoCredentials end
#AWS.region(::AnonymousGCS) = "" # No region
AWS.credentials(::AnonymousGCS) = nothing # No credentials
#AWS.check_credentials(c::NoCredentials) = c # Skip credentials check
AWS.sign!(::AnonymousGCS, ::AWS.Request) = nothing # Don't sign request
function AWS.generate_service_url(::AnonymousGCS, service::String, resource::String)
    service == "s3" || throw(ArgumentError("Can only handle s3 requests to GCS"))
    return string("https://storage.googleapis.com.", resource)
end

push!(storageregexlist,r"^gs://"=>AnonymousGCS)
push!(storageregexlist,r"^s3://"=>S3Store)

function storefromstring(::Type{<:S3Store}, s, _)
  decomp = split(s,"/",keepempty=false)
  bucket = decomp[2]
  path = join(decomp[3:end],"/")
  S3Store(String(bucket),aws=AWS.global_aws_config()),path
end

function storefromstring(::Type{<:AnonymousGCS}, s, _)
    decomp = split(s,"/",keepempty=false)
    bucket = decomp[2]
    path = join(decomp[3:end],"/")
    S3Store(String(bucket), aws=AnonymousGCS(), listversion=1), path
end


text:
function listavailablecropondate(date)
    #transform date to weeks from beginning of year
    wdate = Week(Date(date)).value

    #check table of range of harvest dates
    possible_harvest_weeks = CSV.read("data/possible_harvest_weeks.csv", DataFrame)
    subdf = @where(possible_harvest_weeks, :week .== wdate)

    #check price table
    price_table = CSV.read("data/crop plan - jm_yields_and_sales.csv", DataFrame)
    cropinfo = @pipe leftjoin(subdf, price_table, on = :vegetable) |>
        @select(_, :vegetable, :week, :date, :our_yield_kg_or_unit, :desired_selling_price,
        :unit_of_selling_price, :recommended_quantity, :roughly_equivalent_to, :price_of_recommended_quantity, :predicted_revenu_per_our_bed)

    #check standard quantity table

    #calculate price per standard quantity

    #format table
    println(sum(skipmissing(cropinfo[:, :price_of_recommended_quantity])))
    return  cropinfo# vegetable, price, unit, price/m2, standard_quantity, price_per_standard_quantity
end


text:
## matrix class stuff
## Work with Array{Sym}, not python array objects
## requires conversion from SymMatrix -> Array{Sym} in outputs, as appropriate

## covert back to Array{Sym}
subs(ex::Array{Sym}, args...; kwargs...) = Sym[subs(u, args...; kwargs...) for u in ex]

getindex(s::SymMatrix, i::Integer...) = pyeval("x[i]", Sym, x=project(s), i= tuple(([i...].-1)...))
getindex(s::SymMatrix, i::Integer) = pyeval("x[i]", Sym, x=project(s), i=map(x->x-1, ind2sub(size(s), i)))
getindex(s::SymMatrix, i::Symbol) = project(s)[i] # is_nilpotent, ... many such predicates
getindex(s::Array{Sym}, i::Symbol) = project(s)[i] # digaonalize..

## size
function size(x::SymMatrix)
    a = x[:shape]               # a PyObject tuple
    isa(a, PyObject) ? (a[:__getitem__](0), a[:__getitem__](1)) : a
end

function size(x::SymMatrix, dim::Integer)
    if dim <= 2 && pyisinstance(x.x, matrixtype)
        return x[:shape][dim]
    else
        return 1
    end
end

## we want our matrices to be arrays of Sym objects, not symbolic matrices
## so that julia manages them
## it is convenient (for printing, say) to convert to a sympy matrix
convert(::Type{SymMatrix}, a::Array{Sym}) = Sym(sympy[:Matrix](a))
convert(::Type{Sym}, a::Array{Sym}) = Sym(sympy[:Matrix](a))
function convert(::Type{Array{Sym}}, a::SymMatrix)
    sz = size(a)
    ndims = length(sz)
    if ndims == 0
        a
    elseif ndims == 1
        Sym[a[i] for i in 1:length(a)]
    elseif ndims == 2
        Sym[a[i,j] for i in 1:size(a)[1], j in 1:size(a)[2]]
    else
        ## need something else for arrays... XXX -- can't linear index a
        b = Sym[a[i] for i in 1:length(a)]
        reshape(b, sz)
    end
end

## when projecting, we convert to a symbolic matrix then project
project(x::Array{Sym}) = convert(SymMatrix, x) |> project

## linear algebra functions that are methods of sympy.Matrix
## return a "scalar"
for meth in (:condition_number,
           )

    cmd = "x." * string(meth) * "()"
    meth_name=string(meth)
    @eval begin
           @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(a::SymMatrix) = pycall(a[($cmd)], Sym, x) #pyeval(($cmd), x=project(a)))
        ($meth)(a::Matrix{Sym}) = ($meth)(convert(SymMatrix, a))
    end
    eval(Expr(:export, meth))
end

## can be any dimension
for meth in (:norm,
             )
    meth_name = string(meth)
    @eval begin
        @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(a::SymMatrix, args...; kwargs...) = object_meth(a, $meth_name, args...;kwargs...)
        ($meth)(a::Array{Sym}, args...; kwargs...) = ($meth)(convert(SymMatrix, a), args...;kwargs...)
    end
    eval(Expr(:export, meth))
end

for meth in (
           :has,
           )
    meth_name = string(meth)
    @eval begin
      @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(a::SymMatrix, args...; kwargs...) = object_meth(a, $meth_name, args...;kwargs...)
        ($meth)(a::Matrix{Sym}, args...; kwargs...) = ($meth)(convert(SymMatrix, a), args...;kwargs...)
    end
    eval(Expr(:export, meth))
end

## dont' define inv -- it has amiguity with base inv
inverse(ex::SymMatrix) = call_matrix_meth(ex, :inv)
inverse(ex::Matrix{Sym}) = call_matrix_meth(convert(SymMatrix, ex),:inv)
export inverse

## But
## is_symmetric <-> issym
## istriu, istril,
for meth in  (:is_anti_symmetric, :is_diagonal, :is_diagonalizable,:is_nilpotent,
                :is_symbolic, :is_symmetric)

    meth_name = string(meth)
    @eval begin
      @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(ex::SymMatrix, args...; kwargs...) = ex[@compat(Symbol($meth_name))]()
        ($meth)(ex::Matrix{Sym}, args...; kwargs...) = ex[@compat(Symbol($meth_name))]()
    end
    eval(Expr(:export, meth))
end

## methods called as properties

## is_upper <-> itriu
## is_lower <-> istril
matrix_operators = (:H, :C,
                    :is_lower, :is_lower_hessenberg, :is_square, :is_upper,  :is_upper_hessenberg, :is_zero
)

for meth in matrix_operators
     meth_name = string(meth)
    @eval begin
      @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(ex::SymMatrix, args...; kwargs...) = ex[@compat(Symbol($meth_name))]
        ($meth)(ex::Matrix{Sym}, args...; kwargs...) = ex[@compat(Symbol($meth_name))]
    end
    eval(Expr(:export, meth))
end

## These take a matrix, return a container of symmatrices. Here we convert these to arrays of sym
## Could use
## - qr generically instead of QRdecomposition
## - lu instead of
## - chol will call cholesky
map_matrix_methods = (:LDLsolve,
                      :LDLdecomposition, :LDLdecompositionFF,
                      :LUdecomposition_Simple,
                      :LUsolve,
                      :QRdecomposition, :QRsolve,
                      :adjoint, :adjugate,
                      :cholesky, :cholesky_solve, :cofactor, :conjugate,
                      :diagaonal_solve, :diagonalize, :dual,
                      :expand,
                      :integrate,
                      :inverse_ADJ, :inverse_GE, :inverse_LU,
                      :jordan_cells, :jordan_form,
                      :limit,
                      :lower_triangular_solve,
                      :minorEntry, :minorMatrix,
                      :normalized, :nullspace,
                      :permuteBkwd, :permuteFwd,
                      :print_nonzero,
                      :singular_values,
                      :upper_triangular_solve,
                      :vech
                      )

for meth in map_matrix_methods
    meth_name = string(meth)
    @eval begin
      @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(ex::SymMatrix, args...; kwargs...) = call_matrix_meth(ex, @compat(Symbol($meth_name)), args...; kwargs...)
        ($meth)(ex::Matrix{Sym}, args...; kwargs...) = call_matrix_meth(convert(SymMatrix, ex), @compat(Symbol($meth_name)), args...; kwargs...)
    end
    eval(Expr(:export, meth))
end

### Some special functions
Base.chol(a::Matrix{Sym}) = cholesky(a)

expm(ex::Matrix{Sym}) = convert(Array{Sym}, object_meth(convert(SymMatrix, ex), :exp))
expm(a::SymMatrix) = a[:exp]()

Base.conj(a::SymMatrix) = conjugate(a)
Base.conj(a::Sym) = conjugate(a)

## :eigenvals, returns {val => mult, val=> mult}
## we return an array of eigen values, as eigvals does
function eigvals(a::Matrix{Sym})
    ## this is a hack, as  d = a[:eigenvals]() may not convert to a Julia dict (Ubuntu...)
    ds = a[:eigenvects]()
    [d[1] for d in ds]
end
eigvals(a::SymMatrix) = eigvals(convert(Array{Sym}, a))

## :eigenvects ## returns list of triples (eigenval, multiplicity, basis).
"""

The `eigvecs` function returns a list of triples (eigenval, multiplicity, basis) for an `Matrix{Sym}`.

"""
function eigvecs(a::Matrix{Sym})
    ds = a[:eigenvects]()

    hcat([hcat([convert(Array{Sym}, v) for v in d[3]]...) for d in ds]...)

end
eigvecs(a::SymMatrix) = eigvecs(convert(Array{Sym}, a))

## Take any matrix and return reduced row-echelon form and indices of pivot vars
## To simplify elements before finding nonzero pivots set simplified=True
## We return only the rref form -- not the pivot vars.
"""

Return reduced row echelon form. This function does not return the
pivot variables generated by SymPy. For those, something like `m,
pivots = convert(SymMatrix, a)[:rref]()` would work.

As `rref` is no longer in base, as of v"0.4.0", we extend the
interface to `Matrix{Int}` and `Matrix{Rational}`. The return value
will be `Matrix{Int}` if possible, otherwise a matrix of type
rational.

"""
function rref(a::Matrix{Sym}; kwargs...)
    rref(convert(SymMatrix, a); kwargs...)
end

## rref. The sympy method returns
function rref(a::SymMatrix)
  d = a[:rref]()
  convert(Array{Sym}, d[1]) ## return Array{Sym}, not SymMatrix
end

rref{T <: Integer}(a::Matrix{T}) = N(rref(convert(Matrix{Sym}, a)))
rref{T <: Integer}(a::Matrix{Rational{T}}) = N(rref(convert(Matrix{Sym}, a)))

## call with a (A,b), return array
for fn in (:cross,
           :LUSolve)
    cmd = "x." * string(fn) * "()"
    meth_name = string(fn)
    @eval begin
      @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($fn)(A::SymMatrix, b::Sym) = object_meth(A, fn, b) #convert(Array{Sym}, pyeval(($cmd), A=project(A), b=project(b)))
        ($fn)(A::Array{Sym, 2}, b::Vector{Sym}) = $(fn)(convert(SymMatrix,A), convert(SymMatrix, b))
    end
end

## GramSchmidt -- how to call?
## call with a (A,b), return scalar
# for fn in ()
#     @eval ($fn)(A::Sym, b::Sym) = convert(Array{Sym}, pyeval("A.($fn)(b)", A=project(A), b=project(b)))
#     @eval ($fn)(A::Array{Sym, 2}, b::Vector{Sym}) = ($fn)(convert(SymMatrix, A), convert(SymMatrix, b))
# end

### Higher dimensional derivatives

## For gradient we have [diff(ex,var) for var in free_symbols(ex)]... but order is of importance...

"""

Jacobian of a symbolic matrix

Example:

```
rho, phi = symbols("rho, phi")
M = [rho*cos(phi), rho*sin(phi), rho^2]
Y = [rho, phi]
jacobian(M, Y)
```
"""
function jacobian(x::Array{Sym}, y::Array{Sym})
    X = convert(SymMatrix, x)
    Y = convert(SymMatrix, y)
    call_matrix_meth(X, :jacobian, Y)
end

"""

Find Hessian of a symbolic expression

Example:
```
x, y = symbols("x y")
f = x^2 - 2x*y
hessian(f, [x,y])
```
"""
function hessian(f::Sym, x::Vector{Sym})
    out = sympy_meth(:hessian, f, x)
    convert(SymMatrix, out) |> u -> convert(Array{Sym}, u)
end
hessian(ex::Sym) = hessian(ex, free_symbols(ex))


text:
module CoffeeBeetles

export main

using Serialization, DungBase

using DataStructures, CoordinateTransformations, Rotations, DataFrames, Missings, Distributions, AngleBetweenVectors, LinearAlgebra, StatsBase, OnlineStats, Colors, PrettyTables, Measurements, HypothesisTests, GLM, DelimitedFiles, Printf

using CairoMakie, MakieLayout, FileIO, AbstractPlotting
import AbstractPlotting:px
CairoMakie.activate!()

include("generate_artifacts.jl")
include("preparedata.jl")
include("stats.jl")
include("plot.jl")

"""
    main()
Create all the tables and figures included in the manuscript.
"""
function main()
    data = deserialize(datafile)
    df = getdf(data)
    speeds!(df)
    descriptive_stats(df)
    displaced_stats(df)
    save_figures(df)
end

end # module


text:
function checkmt!(
    latt_vars::LatticeVars, atm_vars::AtomicVars, spsymb,
    mtr_vars::MuffinTinRadialVars
)

    rmt = mtr_vars.rmt
    rmtdelta = mtr_vars.rmtdelta
    epslat = latt_vars.epslat
    nspecies = atm_vars.nspecies

    rmt0 = zeros(Float64,nspecies)
    rmt0[1:nspecies] = rmt[1:nspecies]

    println("rmtdelta = ", rmtdelta)

    while true
        # find the minimum distance between muffin-tin surfaces
        dmin, is, js = mtdmin(latt_vars, atm_vars, rmt)
        # adjust muffin-tin radii if required1
        if dmin < (rmtdelta - epslat)
            println("Adjusting rmt")
            t1 = rmt[is] + rmt[js]
            t2 = (t1 + dmin - rmtdelta)/t1
            rmt[is] = rmt[is]*t2
            if is != js
                rmt[js] = rmt[js]*t2
            end
        else
            break
        end
    end

    for is in 1:nspecies
        if rmt[is] < 0.25
            @printf("Error(checkmt): muffin-tin radius too small for species %d %s\n", is, spsymb[is])
            @printf("Radius : %18.10f\n", rmt[is])
        end
        if rmt[is] < rmt0[is]
            @printf("Info(checkmt): reduced muffin-tin radius of species %3d %s", is, spsymb[is])
            @printf(" is reduced from %f to %f\n", rmt0[is], rmt[is])
        end
    end

    return
end


text:
function prime_factors(number)
    factors = []
    divisor = 2
    while divisor <= number
        quotient, reminder = divrem(number, divisor)

        if iszero(reminder)
            push!(factors, divisor)
            number = quotient
            divisor = 1
        end

        divisor += 1
    end
    return factors
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
code_container = (borderLeft = "thin #C8D4E3 solid",)

container_padding_bottom = (
    padding = 0,
    paddingBottom = "15px",
)

notebox = (
    width = "90%",
    margin = "0 auto",
    border = "1px solid lightgrey",
    textAlign = "center",
    borderRadius = "5px",
    backgroundColor = "lightgrey",
)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using JLD
using PyPlot
using Distributed
using SharedArrays
function plot_optimization_path()
		#filename = string("../data/",
		#				  "optimization/",
		#				  "beta_dJds",
		#				  ".jld")

		#filename1 = string("../data/",
		#				  "optimization/",
		#				  "beta_dJds_1",
		#				  ".jld")

		#X = load(filename)
		#dJds = X["dJds_path"]
		#beta_path = X["beta_path"]

		#Y = load(filename1)
		#dJds1 = Y["dJds_path"]
		#beta1_path = Y["beta_path"]

		X = load("../data/optimization/Eac_along_path_05.jld")
		beta = X["beta"]
		#beta1 = X["beta1"]
		Eac = X["Eac"]
		#Eac1 = X["Eac1"]

		X = load("../data/attractor/more_les/Eac.jld")
		Y = load("../data/attractor/more_les/Eac1.jld")
		Z = load("../data/attractor/Eac.jld")
		W = load("../data/attractor/Eac_and_beta.jld")
		J1 = X["Eac"]
		J2 = Y["Eac"]
		s1 = X["beta"]
		s2 = Y["beta"]
		J3 = Z["Eac"]
		J4 = W["Eac"]
		s3 = Z["beta"]
		s4 = W["beta"]

		ll = 4.5
		ul = 7.5

		J1 = J1[ll .<= s1 .<= ul]
		J2 = J2[ll .<= s2 .<= ul]
		J3 = J3[ll .<= s3 .<= ul]
		J4 = J4[ll .<= s4 .<= ul]

		s1 = s1[ll .<= s1 .<= ul]
		s2 = s2[ll .<= s2 .<= ul]
		s3 = s3[ll .<= s3 .<= ul]
		s4 = s4[ll .<= s4 .<= ul]

		fig, ax = subplots(1,1)
		ax.xaxis.set_tick_params(labelsize=36)
		ax.yaxis.set_tick_params(labelsize=36)

		ax.plot(s1, J1, "r.",ms=6.0)
		ax.plot(s2, J2, "r.",ms=6.0)
		ax.plot(s3, J3, "r.",ms=6.0)
		ax.plot(s4, J4, "r.",ms=6.0)

		ax.plot(beta[1:10], Eac[1:10], "bP", ms=8.0)
		#ax.plot(beta1[1:6], Eac1[1:6], "kP", ms=8.0)
		#ax.plot(6.5, 23.1, "bP", ms=8.0)

		ax.grid(true)
		ax.set_xlabel(L"$\beta$", fontsize=36)
		ax.set_ylabel(L"$\langle J_{\rm ac}\rangle$",
					  rotation=0,fontsize=36)

		ll = 6.0
		ul = 7.3

		J1 = J1[ll .<= s1 .<= ul]
		J2 = J2[ll .<= s2 .<= ul]
		J3 = J3[ll .<= s3 .<= ul]
		J4 = J4[ll .<= s4 .<= ul]

		s1 = s1[ll .<= s1 .<= ul]
		s2 = s2[ll .<= s2 .<= ul]
		s3 = s3[ll .<= s3 .<= ul]
		s4 = s4[ll .<= s4 .<= ul]

		ax_in = fig.add_axes([0.3, 0.5, 0.3, 0.35])
		ax_in.plot(s1, J1, "r.",ms=6.0)
		ax_in.plot(s2, J2, "r.",ms=6.0)
		ax_in.plot(s3, J3, "r.",ms=6.0)
		ax_in.plot(s4, J4, "r.",ms=6.0)
		ax_in.plot(beta[1:5], Eac[1:5], "bP", ms=8.0)
		#ax_in.plot(beta[1], Eac[1], "bP", ms=8.0)
		#ax_in.plot(beta1[1], Eac1[1], "kP", ms=8.0)
		#ax_in.plot(6.5, 23.1, "bP", ms=8.0)

		ax_in.grid(true)
		ax_in.xaxis.set_tick_params(labelsize=36)
		ax_in.yaxis.set_tick_params(labelsize=36)
		@show beta
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
module InfiniteOptSymbolics

# Import the ncessary packages
import InfiniteOpt, ModelingToolkit

const _Symbolics = ModelingToolkit.Symbolics
const _JuMP = InfiniteOpt.JuMP

# Include the code
include("utilities.jl")
include("conversion.jl")

# Define additional stuff that should not be exported
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

# Following JuMP, export everything that doesn't start with a _
for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") && Base.isidentifier(sym_string[2:end])))
        continue
    end
    @eval export $sym
end

end # end of module


text:
#=
The power set of a set is the set of all its subsets. Write a function that, given a set, generates its power set.

For example, given the set {1, 2, 3}, it should return {{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}.

You may also use a list or array to represent a set.
=#

function get_powerset(set::Array{Int,1})
    powerset = [[]]
    for number in set
        push!(powerset, [number])
    end
    set_length = 1
    while set_length < length(set)
        newpowerset = copy(powerset)
        for arr_index in eachindex(powerset)
            arr = powerset[arr_index]
            for number in set
                if length(arr) == set_length && !in(number, arr) && number > arr[end]
                    new_array = union(arr, number)
                    push!(newpowerset, new_array)
                end
            end
        end
        powerset = newpowerset
        set_length = set_length + 1
    end
    return powerset
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using HTTP, JSON

"""
    query

Fire a query to the server.
Display the search results.
"""
function query(pattern)
    resp = HTTP.request("POST",
            "http://127.0.0.1:30005/search",
            [],
            "{\"pattern\": \"$pattern\"}"
        )
    jsonresp = JSON.parse(String(resp.body))
    @assert jsonresp["success"]
    for x in jsonresp["data"]
        println(x["file"], ":", x["line"])
    end
end

max_stars_repo_name: tanmaykm/KubernetesIntegratedJulia

max_issues_repo_name: tanmaykm/KubernetesIntegratedJulia

max_forks_repo_name: tanmaykm/KubernetesIntegratedJulia


text:
import brainflow

brainflow.enable_dev_board_logger()

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))
sampling_rate = brainflow.get_sampling_rate(Integer(brainflow.SYNTHETIC_BOARD))

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
downsampled_data = brainflow.perform_downsampling(data_first_channel, 3, Integer(brainflow.EACH))
println("After Downsampling")
println(downsampled_data)

max_stars_repo_name: kusumikakd/Brainflow

max_issues_repo_name: kusumikakd/Brainflow

max_forks_repo_name: kusumikakd/Brainflow