text:
# # Wind and convection-driven mixing in an ocean surface boundary layer
#
# This example simulates mixing by three-dimensional turbulence in an ocean surface
# boundary layer driven by atmospheric winds and convection. It demonstrates:
#
#   * how to use the `SeawaterBuoyancy` model for buoyancy with a linear equation of state;
#   * how to use a turbulence closure for large eddy simulation;
#   * how to use a function to impose a boundary condition;
#   * how to use user-defined `model.parameters` in a boundary condition function.
#
# In addition to `Oceananigans.jl` we need `Plots` for plotting, `Random` for
# generating random initial conditions, and `Printf` for printing progress messages.
# We also need `Oceananigans.OutputWriters` and `Oceananigans.Diagnostics` to access
# some nice features for writing output data to disk.

using Random, Printf, Plots
using Oceananigans, Oceananigans.OutputWriters, Oceananigans.Diagnostics, Oceananigans.Utils

# ## Model parameters
#
# Here we use an isotropic, cubic grid with `Nz` grid points and grid spacing
# `Δz = 1` meter. We specify fluxes of heat, momentum, and salinity via
#
#   1. A temperature flux `Qᵀ` at the top of the domain, which is related to heat flux
#       by `Qᵀ = Qʰ / (ρ₀ * cᴾ)`, where `Qʰ` is the heat flux, `ρ₀` is a reference density,
#       and `cᴾ` is the heat capacity of seawater. With a reference density
#       `ρ₀ = 1026 kg m⁻³`and heat capacity `cᴾ = 3991`, our chosen temperature flux of
#       `Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹` corresponds to a heat flux of `Qʰ = 204.7 W m⁻²`, a
#       relatively powerful cooling rate.
#
#   2. A velocity flux `Qᵘ` at the top of the domain, which is related
#       to the `x` momentum flux `τˣ` via `τˣ = ρ₀ * Qᵘ`, where `ρ₀` is a reference density.
#       Our chosen value of `Qᵘ = -2 × 10⁻⁵ m² s⁻²` roughly corresponds to atmospheric winds
#       of `uᵃ = 2.9 m s⁻¹` in the positive `x`-direction, using the parameterization
#       `τ = 0.0025 * |uᵃ| * uᵃ`.
#
#   3. An evaporation rate `evaporation = 10⁻⁷ m s⁻¹`, or approximately 0.1 millimeter per
#       hour.
#
# Finally, we use an initial temperature gradient of `∂T/∂z = 0.005 K m⁻¹`,
# which implies an iniital buoyancy frequency `N² = α * g * ∂T/∂z = 9.8 × 10⁻⁶ s⁻²`
# with a thermal expansion coefficient `α = 2 × 10⁻⁴ K⁻¹` and gravitational acceleration
# `g = 9.81 s⁻²`. Note that, by default, the `SeawaterBuoyancy` model uses a gravitational
# acceleration `gᴱᵃʳᵗʰ = 9.80665 s⁻²`.

         Nz = 32       # Number of grid points in x, y, z
         Δz = 1.0      # Grid spacing in x, y, z (meters)
         Qᵀ = 5e-5     # Temperature flux at surface
         Qᵘ = -2e-5    # Velocity flux at surface
       ∂T∂z = 0.005    # Initial vertical temperature gradient
evaporation = 1e-7     # Mass-specific evaporation rate [m s⁻¹]
          f = 1e-4     # Coriolis parameter
          α = 2e-4     # Thermal expansion coefficient
          β = 8e-4     # Haline contraction coefficient
nothing # hide

# ## Boundary conditions
#
# Here we define `Flux` boundary conditions at the surface for `u`, `T`, and `S`,
# and a `Gradient` boundary condition on `T` that maintains a constant stratification
# at the bottom. Our flux boundary condition for salinity uses a function that calculates
# the salinity flux in terms of the evaporation rate.

grid = RegularCartesianGrid(size=(Nz, Nz, Nz), length=(Δz*Nz, Δz*Nz, Δz*Nz))

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))

T_bcs = TracerBoundaryConditions(grid,    top = BoundaryCondition(Flux, Qᵀ),
                                       bottom = BoundaryCondition(Gradient, ∂T∂z))

## Salinity flux: Qˢ = - E * S
@inline Qˢ(i, j, grid, time, iter, U, C, p) = @inbounds -p.evaporation * C.S[i, j, 1]

S_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qˢ))
nothing # hide

# ## Model instantiation
#
# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`,
# using a `FPlane` model for rotation (constant rotation rate), a linear equation
# of state for temperature and salinity, the Anisotropic Minimum Dissipation closure
# to model the effects of unresolved turbulence, and the previously defined boundary
# conditions for `u`, `T`, and `S`. We also pass the evaporation rate to the container
# model.parameters for use in the boundary condition function that calculates the salinity
# flux.

model = Model(
         architecture = CPU(),
                 grid = grid,
             coriolis = FPlane(f=f),
             buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
              closure = AnisotropicMinimumDissipation(),
  boundary_conditions = SolutionBoundaryConditions(grid, u=u_bcs, T=T_bcs, S=S_bcs),
           parameters = (evaporation = evaporation,)
)
nothing # hide

# Notes:
#
# * To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than
#   `AnisotropicMinimumDissipation`, use `closure = ConstantSmagorinsky()` in the model constructor.
#
# * To change the `architecture` to `GPU`, replace the `architecture` keyword argument with
#   `architecture = GPU()``

# ## Initial conditions
#
# Our initial condition for temperature consists of a linear stratification superposed with
# random noise damped at the walls, while our initial condition for velocity consists
# only of random noise.

## Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

## Temperature initial condition: a stable density tradient with random noise superposed.
T₀(x, y, z) = 20 + ∂T∂z * z + ∂T∂z * model.grid.Lz * 1e-6 * Ξ(z)

## Velocity initial condition: random noise scaled by the friction velocity.
u₀(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

set!(model, u=u₀, w=u₀, T=T₀, S=35)

# ## Set up output
#
# We set up an output writer that saves all velocity fields, tracer fields, and the subgrid
# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument
# to `JLD2OutputWriter` indicates that output will be saved in
# `ocean_wind_mixing_and_convection.jld2`.

## Create a NamedTuple containing all the fields to be outputted.
fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))

## Instantiate a JLD2OutputWriter to write fields. We will add it to the simulation before
## running it.
field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
                                prefix="ocean_wind_mixing_and_convection", force=true)

# ## Running the simulation
#
# To run the simulation, we instantiate a `TimeStepWizard` to ensure stable time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 0.2.

wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=5.0)
nothing # hide

# A diagnostic that returns the maximum absolute value of `w` by calling
# `wmax(model)`:

wmax = FieldMaximum(abs, model.velocities.w)
nothing # hide

# Finally, we set up and run the the simulation.

simulation = Simulation(model, Δt=wizard, stop_iteration=0, progress_frequency=10)
simulation.output_writers[:fields] = field_writer

anim = @animate for i in 1:100
    ## Run the simulation forward
    simulation.stop_iteration += 10
    walltime = @elapsed run!(simulation)

    ## Print a progress message
    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
            model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
            wmax(model), prettytime(walltime))

    ## Coordinate arrays for plotting
    xC, zF, zC = model.grid.zC, model.grid.zF[1:Nz], model.grid.zC

    ## Slices to plots.
    jhalf = floor(Int, model.grid.Ny/2)
    w = Array(interior(model.velocities.w))[:, jhalf, :]
    T = Array(interior(model.tracers.T))[:, jhalf, :]
    S = Array(interior(model.tracers.S))[:, jhalf, :]

    ## Plot the slices.
    w_plot = heatmap(xC, zF, w', xlabel="x (m)", ylabel="z (m)", color=:balance, clims=(-3e-2, 3e-2))
    T_plot = heatmap(xC, zC, T', xlabel="x (m)", ylabel="z (m)", color=:thermal, clims=(19.75, 20))
    S_plot = heatmap(xC, zC, S', xlabel="x (m)", ylabel="z (m)", color=:haline, clims=(34.99, 35.01))

    ## Arrange the plots side-by-side.
    plot(w_plot, T_plot, S_plot, layout=(1, 3), size=(1600, 400),
         title=["vertical velocity (m/s)" "temperature (C)" "salinity (g/kg)"])
end

mp4(anim, "ocean_wind_mixing_and_convection.mp4", fps = 15) # hide


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

using SparseArrays
Base.getproperty(S::SparseMatrixCSC, ::Symbol) = error("use accessor function")
Base.getproperty(S::SparseVector, ::Symbol) = error("use accessor function")

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Base: Float64
using HTTP: HTTP, request, URI
using JSON3: JSON3, propertynames
using XLSX: XLSX, readxlsx, readtable, readdata, sheetnames
using CSV: CSV, Rows
using DataFrames: DataFrames, DataFrame, select!, combine, groupby, ByRow, dropmissing!,
                  transform!, rename!
const BEA_API_KEY = get(ENV, "api_bea_token", "")

#=
Find the mapping to of industry to NAICS2 (sector) through the APIs
=#
url = URI(scheme = "https", host = "apps.bea.gov", path = "/api/data",
          query = ["UserID" => BEA_API_KEY,
                   "method" => "GetData",
                   "DatasetName" => "GDPbyIndustry",
                   "TableID" => 6,
                   "Year" => join(2008:2019, ','),
                   "Frequency" => "A",
                   "Industry" => "ALL",
                   "ResultFormat" => "JSON"
                  ])
url = URI(scheme = "https", host = "apps.bea.gov", path = "/api/data",
          query = ["UserID" => BEA_API_KEY,
                   "method" => "GetParameterValues",
                   "DatasetName" => "GDPbyIndustry",
                   "ParameterName" => "TableID",
                   "ResultFormat" => "JSON"
                  ])
url = URI(scheme = "https", host = "apps.bea.gov", path = "/api/data",
          query = ["UserID" => BEA_API_KEY,
                   "method" => "GetData",
                   "DatasetName" => "GDPbyIndustry",
                   "TableID" => 25,
                   "Year" => join(2008:2019, ','),
                   "Frequency" => "A",
                   "Industry" => "ALL",
                   "ResultFormat" => "JSON"
                  ])
response = request("GET", url)
@assert response.status == 200
json = JSON3.read(response.body)
json.BEAAPI.Results.Parameter
json.BEAAPI.Results.ParamValue
json.BEAAPI.Results.Dataset
json.BEAAPI.Results.Error
x = DataFrame(json.BEAAPI.Results[1].Data)
select!(x, [:Year, :Industry, :IndustrYDescription, :DataValue])

sort!(combine(groupby(x, :IndustrYDescription), nrow), order(:nrow, rev = true))[1:8,:]

vars = ["Year", "Industry",
        "Compensation of employees", "Energy inputs",
        "Gross operating surplus", "Intermediate inputs", "Materials inputs",
        "Purchased-services inputs", "Taxes on production and imports less subsidies",
        "Value added"]
subset!(x, :IndustrYDescription => ByRow(∈(vars)))
transform!(x, :DataValue => ByRow(x -> parse(Float64, x)), renamecols = false)
y = combine(groupby(x, [:Year, :Industry])) do subdf
    subdf = unstack(subdf, :IndustrYDescription, :DataValue)
    for col in setdiff(vars, names(subdf))
        subdf[!,col] .= 0
    end
    select!(subdf, vars)
end
dropmissing!(y)
select!(subset(y, "Taxes on production and imports less subsidies" => ByRow(<(0))),
        ["Year", "Industry", "Taxes on production and imports less subsidies"])
println(string(names(y)))
transform!(y,
           ["Taxes on production and imports less subsidies", "Compensation of employees"] =>
                ByRow(/) => :top_coe,
           ["Intermediate inputs", "Compensation of employees"] => ByRow(/) => :ii_coe,
           ["Gross operating surplus", "Compensation of employees"] => ByRow(/) => :gos_coe,
           ["Energy inputs", "Purchased-services inputs", "Compensation of employees"] =>
                ByRow((ene, psi, coe) -> (ene + psi) / coe) => :enepsi_coe)
rename!(y, :Year => :year, :Industry => :naics)
select!(y, [:year, :naics, :top_coe, :ii_coe, :gos_coe, :enepsi_coe])
CSV.write(joinpath("data", "klems.csv"), y)

max_stars_repo_name: Nosferican/Estimating-Resource-Costs

max_issues_repo_name: Nosferican/Estimating-Resource-Costs

max_forks_repo_name: Nosferican/Estimating-Resource-Costs


text:
include("CenterOfMass.jl")
include("setup_parameters.jl")

actuator_limit = 4.0
max_robots = 4

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

data = get_data(prior)[:]
@time for ii = 1:max_robots
  for c = combinations(1:length(attachment_ps), ii)
    total_limit = actuator_limit * ii
    action_w = mean(attachment_ws[c])

    field = get_critical_values(circle_ws, action_w, prior, interior_q, total_limit)
    normals = normal_matrix(field[:], sigma^2)

    csqmi = compute_mutual_information(data, normals)
  end
end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
_interpolate(
    ȳ::AbstractVector,
    v::Real,
    stations::AbstractVector;
    interp = LinearInterpolation
) = v .* ones(Float64, length(stations))
_interpolate(
    ȳ::AbstractVector,
    v::AbstractVector,
    stations::AbstractVector;
    interp = LinearInterpolation
) = let A = interp(ȳ, v)
    A.(stations)
end
_interpolate(
    ȳ::AbstractVector,
    v::AbstractMatrix,
    stations::AbstractVector;
    interp = LinearInterpolation
) = permutedims(
    hcat(
        [
            _interpolate(ȳ, v[i, :], stations; interp = interp) for i = 1:size(v, 1)
        ]...
    )
)

_trapz(
    ys::AbstractVector, xs::AbstractVector
) = let dx = xs[2:end] .- xs[1:(end - 1)]
    sum((ys[1:(end - 1)] .+ ys[2:end]) .* dx) / 2
end
_trapz(
    ys::AbstractMatrix, xs::AbstractVector
) = [
    _trapz(ys[i, :], xs) for i = 1:size(ys, 1)
]

_rowcross(
    u::AbstractMatrix,
    v::AbstractMatrix
) = permutedims(
    hcat(
        (u[2, :] .* v[3, :] .- u[3, :] .* v[2, :]),
        (u[3, :] .* v[1, :] .- u[1, :] .* v[3, :]),
        (u[1, :] .* v[2, :] .- u[2, :] .* v[1, :])
    )
)

const default_stations = collect(0.0:0.01:1.0)

"""
    ```
    function simulation(
        ys::AbstractVector,
        cs,
        Ω::Real;
        n_blades::Int64 = 2,
        xs = 0.0,
        CL0 = 0.0,
        CLα = 2 * π,
        Cm0 = 0.0,
        Cmα = 0.0,
        polar = zeros(3),
        e = 0.0,
        γ = 0.0,
        V∞::Real = 0.0,
        a::Real = 340.0,
        ρ::Real = 1.225,
        n_iterations::Int64 = 20,
        h::Real = 1e-3,
        tip_loss::Bool = true,
        stations::AbstractVector = default_stations,
        interp = LinearInterpolation
    )
    ```

Simulate propeller flow

* `ys`: positions of stations in the radial coordinate
    in which properties are defined
* `cs`: chords for each control station in `ys` (or a scalar,
    for a chord applicable to all stations)
* `Ω`: rotational speed (rad/s)
* `n_blades`: number of blades. Defaults to 2
* `xs`: displacements of control section quarter chords from
    radial line. Can be used to apply sweep to the blade, and leads
    to more accurate load predictions
* `CL0, CLα, Cm0, Cmα`: airfoil properties at each control section
* `polar`: polars for airfoils at each control section
    (as a `(3, N)` matrix), or a single polar for
    all sections (as a vector).

    Should follow format: `CD = polar[1] * CL ^ 2 + polar[2] * CL + polar[3]`
* `e`: displacement of elastic axis from quarter-chord, as a fraction
    of the chord. Scalar for the blade, or a vector for each control section
* `γ`: local incidences, in radians. Scalar for the whole blade, or
    vector for each control section
* `V∞`: freestream velocity
* `a`: sound speed
* `ρ`: air density
* `n_iterations`: number of Newton iterations for BEMM
* `h`: step for linearization in Newton method for BEMM
* `tip_loss`: boolean indicating whether or not to apply
    Glauert's tip loss factor
* `stations`: vector ranging from 0 to 1 with positions for
    integration points
* `interp`: interpolation mechanism for obtaining sectional
    properties
"""
function simulation(
    ys::AbstractVector,
    cs,
    Ω::Real;
    n_blades::Int64 = 2,
    xs = 0.0,
    CL0 = 0.0,
    CLα = 2 * π,
    Cm0 = 0.0,
    Cmα = 0.0,
    polar = zeros(3),
    e = 0.0,
    γ = 0.0,
    V∞::Real = 0.0,
    a::Real = 340.0,
    ρ::Real = 1.225,
    n_iterations::Int64 = 20,
    h::Real = 1e-3,
    tip_loss::Bool = true,
    stations::AbstractVector = default_stations,
    interp = LinearInterpolation
)

    polar = (
        isa(polar, AbstractVector) ?
        hcat(
            [polar for _ = 1:length(ys)]...
        ) :
        polar
    )

    R = maximum(ys) - minimum(ys)
    D = 2 * R

    A = R ^ 2 * π

    ȳ = ys ./ R
    ys = _interpolate(ȳ, ys, stations)

    cs = _interpolate(ȳ, cs, stations; interp = interp)
    xs = _interpolate(ȳ, xs, stations; interp = interp)
    CL0 = _interpolate(ȳ, CL0, stations; interp = interp)
    CLα = _interpolate(ȳ, CLα, stations; interp = interp)
    Cm0 = _interpolate(ȳ, Cm0, stations; interp = interp)
    Cmα = _interpolate(ȳ, Cmα, stations; interp = interp)
    γ = _interpolate(ȳ, γ, stations; interp = interp)
    e = _interpolate(ȳ, e, stations; interp = interp)

    polar = _interpolate(ȳ, polar, stations; interp = interp)

    dy = ys[2:end] .- ys[1:(end - 1)]

    Vr = @. ys * Ω

    f = Vi -> begin
        ṁ = A * (Vi + V∞) * ρ

        ϕ = atan.(V∞ + Vi, Vr)

        α = γ .- ϕ

        # Glauert's tip loss factor
        F = @. 2 * acos(exp(- n_blades * ((R - ys + 1e-10) / (ys + 1e-10)) / (2 * sin(ϕ)))) / π

        Vmod = @. sqrt((V∞ + Vi) ^ 2 + Vr .^ 2)
        M = Vmod ./ a

        β = sqrt.(1.0 .- M .^ 2)

        CL = (α .* CLα .+ CL0) .* (tip_loss ? F : 1.0)
        CD = polar[1, :] .* CL .^ 2 .+ polar[2, :] .* CL .+ polar[3, :]

        CL = CL ./ β
        CD = CD ./ β

        qc = @. ((V∞ + Vi) ^ 2 + Vr ^ 2) * ρ * cs / 2

        CX = @. CL * cos(ϕ) - CD * sin(ϕ)
        CY = @. CL * sin(ϕ) + CD * cos(ϕ)

        T = let qcCX = CX .* qc
            sum((qcCX[1:(end - 1)] .+ qcCX[2:end]) .* dy) * n_blades / 2
        end

        T - ṁ * Vi * 2
    end

    h = h * Ω * R

    fp = Vi -> (f(Vi + h) - f(Vi)) / h

    Vi = 0.0

    for _ = 1:n_iterations
        Vi = Vi - f(Vi) / fp(Vi)
    end

    ṁ = A * (Vi + V∞) * ρ

    ϕ = atan.(V∞ + Vi, Vr)

    α = γ .- ϕ

    # Glauert's tip loss factor
    F = @. 2 * acos(exp(- n_blades * ((R - ys + 1e-10) / (ys + 1e-10)) / (2 * sin(ϕ)))) / π

    Vmod = @. sqrt((V∞ + Vi) ^ 2 + Vr .^ 2)
    M = Vmod ./ a

    β = sqrt.(1.0 .- M .^ 2)

    CL = (α .* CLα .+ CL0) .* (tip_loss ? F : 1.0)
    CD = polar[1, :] .* CL .^ 2 .+ polar[2, :] .* CL .+ polar[3, :]

    CL = CL ./ β
    CD = CD ./ β

    qc = @. ((V∞ + Vi) ^ 2 + Vr ^ 2) * ρ * cs / 2

    CX = @. CL * cos(ϕ) - CD * sin(ϕ)
    CY = @. CL * sin(ϕ) + CD * cos(ϕ)

    Cm = (α .* Cmα .+ Cm0) .* (tip_loss ? F : 1.0) ./ β

    T = let qcCX = qc .* CX
        sum(dy .* (qcCX[1:(end - 1)] .+ qcCX[2:end])) * n_blades / 2
    end
    Q = let qcCYy = qc .* CY .* ys
        sum(dy .* (qcCYy[1:(end - 1)] .+ qcCYy[2:end])) * n_blades / 2
    end

    M = Cm .* cs .* qc

    quarter_chords = permutedims(
        [
            xs ys zero(xs)
        ]
    )
    elastic_centers = quarter_chords .+ permutedims(
        [
            (cs .* cos.(γ) .* e) zero(ys) (- cs .* sin.(γ) .* e)
        ]
    )

    dF = permutedims([(CY .* qc) zero(ys) (CX .* qc)])

    forces = hcat(
        [
            _trapz(
                dF[:, i:end], ys[i:end]
            ) for i = 1:length(ys)
        ]...
    )
    moments = let dM = permutedims([zero(ys) M zero(ys)])
        hcat(
            [
                _trapz(
                    dM[:, i:end] .+ _rowcross(
                        quarter_chords[:, i:end] .- elastic_centers[:, i],
                        dF[:, i:end]
                    ),
                    ys[i:end]
                ) for i = 1:length(ys)
            ]...
        )
    end

    My = moments[2, :]
    Mx = moments[1, :] .* cos.(γ) .- moments[3, :] .* sin.(γ)
    Mz = moments[3, :] .* cos.(γ) .+ moments[1, :] .* sin.(γ)

    Vx = - forces[3, :] .* sin.(γ) .+ forces[1, :] .* cos.(γ)
    Vz = forces[3, :] .* cos.(γ) .+ forces[1, :] .* sin.(γ)

    (
        T = T,
        Q = Q,
        CL = CL,
        Cm = Cm,
        CD = CD,
        forces = forces,
        moments = moments,
        Mx = Mx,
        My = My,
        Mz = Mz,
        Vx = Vx,
        Vz = Vz,
        ϕ = ϕ,
        r = ys,
        c = cs,
        Vi = Vi
    )

end

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
using Documenter, DirectDetections

makedocs(
    sitename="DirectDetections.jl",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting-started.md",
        "Tutorials" => [
            "Fit Astrometry" => "modelling.md",
            "Fit Astrometric Acceleration" => "pma.md",
            "Fit Images" => "images.md",
            "Derived Variables" => "derived.md",
            "Connecting Mass with Photometry" => "mass-photometry.md",
            # "Multiple Planets" => "multi-planets.md",
            # "" => "multi-planets.md",
        ],
        "Documentation" => [
            "Priors" => "priors.md",
            "Samplers" => "samplers.md",
            "Chains" => "chains.md",
            "Kepler Solver" => "kepler.md",
            "API" => "api.md"
        ]
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/sefffal/DirectDetections.jl.git",
    devbranch = "main"
)


text:
import PerronFrobenius.TransferOperatorRectangularBinning

@recipe function f(r::TransferOperatorRectangularBinning)
    seriestype  :=  :heatmap
    xlabel --> "state #i"
    ylabel --> "state #j"
    r.transfermatrix
end


text:
export makeGroup,connect,@group
function makeGroup(groupName,groupSize)
  group=Dict()
  group["N"]=sum(groupSize)
  startIndex=cumsum(groupSize)-groupSize+1
  for (idx,name) in enumerate(groupName)
    group[name]=startIndex[idx]:startIndex[idx]+groupSize[idx]-1
  end
  return group
end

macro connect(groups,preGroupName,postGroupName,condition,cmd)
  var=cmd.args[1].args[1]
  lhs=cmd.args[1]
  rhs=cmd.args[2]
  if !isdefined(var)
    eval(:($var=zeros($groups["N"],$groups["N"])))
  end
  ex=quote
    for post=1:length($groups[$postGroupName])
      for pre=1:length($groups[$preGroupName])
        if $condition
          rhsValue=$rhs
          pre=pre+$(groups)[$preGroupName][1]-1
          post=post+$(groups)[$postGroupName][1]-1
          $lhs=rhsValue
        end
      end
    end
  end
  eval(ex)
end

macro group(groups,groupName,cmd)
  var=cmd.args[1].args[1]
  rhs=cmd.args[2]
  lhs=cmd.args[1]
  if !isdefined(var)
    eval(:($var=zeros($groups["N"])))
  end
  ex=quote
    for i=1:length($groups[$groupName])
      rhsValue=$rhs
      i=i+$groups[$groupName][1]-1
      $lhs=rhsValue
    end
  end
  eval(ex)
end


text:
sculpt{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = S
qualia{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = Q
clay{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C})   = C
value{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C})  = (x.lo,x.hi)

isexact{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = qualia(x) == EXACT
isinexact{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = qualia(x) == INEXACT
isclosedclosed{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = sculpt(x) == CLCL
isclosedopen{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = sculpt(x) == CLOP
isopenclosed{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = sculpt(x) == OPCL
isopenopen{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = sculpt(x) == OPOP

# each boundry is closed(false) or Open (true)
# get boundries from Sculpt
boundries(::Type{CLCL}) = (false,false)
boundries(::Type{CLOP}) = (false,true)
boundries(::Type{OPCL}) = (true,false)
boundries(::Type{OPOP}) = (true,true)
# get Sculpt from boundries
boundries{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C}) = boundries(S)
# get Sculpt from boundries
const Boundries = [CLCL, OPCL, CLOP, OPOP]
boundries{B<:Bool}(lo::B,hi::B) = Boundries[ one(Int8)+reinterpret(Int8,lo)+(reinterpret(Int8,hi)<<1) ]

# select type corresponding to negation of type given
@inline negate(::Type{CLCL}) = CLCL
@inline negate(::Type{CLOP}) = OPCL
@inline negate(::Type{OPCL}) = CLOP
@inline negate(::Type{OPOP}) = OPOP

#=
    Closed bounds are reached but Not Exceeded.
    Open bounds are approached but Not Reached.

    closed(Flex) gives the ClCl inclusion of Flex
    opened(Flex) gives the OpOp inclusion of Flex
=#

closed{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = ClCl(x.lo, x.hi)

opened{Q<:Qualia,C<:Clay}(x::Flex{CLCL,Q,C}) = Flex{OPOP,Q,C}(prevFloat(x.lo), nextFloat(x.hi))
opened{Q<:Qualia,C<:Clay}(x::Flex{CLOP,Q,C}) = Flex{OPOP,Q,C}(prevFloat(x.lo), x.hi)
opened{Q<:Qualia,C<:Clay}(x::Flex{OPCL,Q,C}) = Flex{OPOP,Q,C}(x.lo, nextFloat(x.hi))
opened{Q<:Qualia,C<:Clay}(x::Flex{OPOP,Q,C}) = x

clopened{Q<:Qualia,C<:Clay}(x::Flex{CLCL,Q,C}) = Flex{CLOP,Q,C}(x.lo, nextFloat(x.hi))
clopened{Q<:Qualia,C<:Clay}(x::Flex{CLOP,Q,C}) = x
clopened{Q<:Qualia,C<:Clay}(x::Flex{OPCL,Q,C}) = Flex{CLOP,Q,C}(x.lo, nextFloat(x.hi))
clopened{Q<:Qualia,C<:Clay}(x::Flex{OPOP,Q,C}) = Flex{CLOP,Q,C}(x.lo, x.hi)

opclosed{Q<:Qualia,C<:Clay}(x::Flex{CLCL,Q,C}) = Flex{CLOP,Q,C}(prevFloat(x.lo), x.hi)
opclosed{Q<:Qualia,C<:Clay}(x::Flex{CLOP,Q,C}) = Flex{CLOP,Q,C}(prevFloat(x.lo), x.hi)
opclosed{Q<:Qualia,C<:Clay}(x::Flex{OPCL,Q,C}) = x
opclosed{Q<:Qualia,C<:Clay}(x::Flex{OPOP,Q,C}) = Flex{CLOP,Q,C}(x.lo, x.hi)


text:
module MeshAlgorithms

using GeometryTypes
using StaticArrays
using LinearAlgebra

export
    Mesh,
    convert_from_homogenous,
    distance,
    fundamental_error_quadric

include("simplification.jl")

end # module


text:
abstract type AbstractEquationOfState end
abstract type AbstractIdealGas <: AbstractEquationOfState end

struct BarotropicFluid <: AbstractEquationOfState end
struct DryIdealGas <: AbstractIdealGas end
struct MoistIdealGas <: AbstractIdealGas end

@inline function calc_pressure(::BarotropicFluid, state, aux, params)
    ρ = state.ρ
    cₛ = params.cₛ
    ρₒ = params.ρₒ

    return (cₛ * ρ)^2 / (2 * ρₒ)
end

@inline function calc_pressure(eos::DryIdealGas, state, aux, params)
    ρ = state.ρ

    R_d = calc_gas_constant(eos, state, params)
    T = calc_air_temperature(eos, state, aux, params)

    return ρ * R_d * T
end

@inline function calc_pressure(eos::MoistIdealGas, state, aux, params)
    ρ = state.ρ

    R_m = calc_gas_constant(eos, state, params)
    T = calc_air_temperature(eos, state, aux, params)

    return ρ * R_m * T
end

@inline function calc_linear_pressure(eos::DryIdealGas, state, aux, params)
    ρ = state.ρ
    ρe = state.ρe
    Φ = aux.Φ
    T_0 = params.T_0

    γ = calc_heat_capacity_ratio(eos, state, params)
    cv_d = calc_heat_capacity_at_constant_volume(eos, state, params)

    return (γ - 1) * (ρe - ρ * Φ + ρ * cv_d * T_0)
end

@inline function calc_linear_pressure(::MoistIdealGas, state, aux, params)
    ρ = state.ρ
    ρe = state.ρe
    ρ_q_tot = state.ρq
    ρ_q_liq = 0 # zero for now
    ρ_q_ice = 0 # zero for now
    Φ = aux.Φ
    T_0 = params.T_0

    γ = calc_heat_capacity_ratio(DryIdealGas(), state, params)
    cv_d = calc_heat_capacity_at_constant_volume(DryIdealGas(), state, params)

    ρ_e_latent = (ρ_q_tot - ρ_q_liq) * params.e_int_v0 - ρ_q_ice * (params.e_int_v0 + params.e_int_i0)

    return (γ - 1) * (ρe - ρ * Φ - ρ_e_latent + ρ * cv_d * T_0)
end

@inline function calc_very_linear_pressure(eos::DryIdealGas, state, aux, params)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ

    # Reference states
    ρᵣ = aux.ref_state.ρ
    ρuᵣ = aux.ref_state.ρu

    γ = calc_heat_capacity_ratio(eos, state, params)

    return (γ - 1) * (ρe - dot(ρuᵣ, ρu) / ρᵣ + ρ * dot(ρuᵣ, ρuᵣ) / (2 * ρᵣ^2) - ρ * Φ)
end

@inline function calc_sound_speed(::BarotropicFluid, state, aux, params)
    ρ = state.ρ
    cₛ = params.cₛ
    ρₒ = params.ρₒ

    return cₛ * sqrt(ρ / ρₒ)
end

@inline function calc_sound_speed(eos::DryIdealGas, state, aux, params)
    ρ = state.ρ

    γ = calc_heat_capacity_ratio(eos, state, params)
    p = calc_pressure(eos, state, aux, params)

    return sqrt(γ * p / ρ)
end

@inline function calc_sound_speed(eos::MoistIdealGas, state, aux, params)
    ρ = state.ρ

    γ = calc_heat_capacity_ratio(eos, state, params)
    p = calc_pressure(eos, state, aux, params)

    return sqrt(γ * p / ρ)
end

@inline function calc_ref_sound_speed(eos::DryIdealGas, state, aux, params)
    p = aux.ref_state.p
    ρ = aux.ref_state.ρ

    γ = calc_heat_capacity_ratio(eos, state, params)

    return sqrt(γ * p / ρ)
end

@inline function calc_ref_sound_speed(::MoistIdealGas, state, aux, params)
    p = aux.ref_state.p
    ρ = aux.ref_state.ρ

    γ = calc_heat_capacity_ratio(DryIdealGas(), state, params)

    return sqrt(γ * p / ρ)
end

@inline function calc_air_temperature(::DryIdealGas, state, aux, params)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    T_0 = params.T_0

    cv_d = calc_heat_capacity_at_constant_volume(DryIdealGas(), state, params)

    e_int = (ρe - ρu' * ρu / 2ρ - ρ * Φ) / ρ
    T = T_0 + e_int / cv_d

    return T
end

@inline function calc_air_temperature(eos::MoistIdealGas, state, aux, params)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    q_tot = state.ρq / ρ
    q_liq = 0.0 # zero for now
    q_ice = 0.0 # zero for now
    Φ = aux.Φ
    T_0 = params.T_0
    e_int_v0 = params.e_int_v0
    e_int_i0 = params.e_int_i0

    cv_m = calc_heat_capacity_at_constant_volume(eos, state, params)

    e_int = (ρe - ρu' * ρu / 2ρ - ρ * Φ) / ρ
    T = T_0 + (e_int - (q_tot - q_liq) * e_int_v0 + q_ice * (e_int_v0 + e_int_i0)) / cv_m

    return T
end

@inline function calc_total_specific_enthalpy(eos::DryIdealGas, state, aux, params)
    ρ = state.ρ
    ρe = state.ρe

    p = calc_pressure(eos, state, aux, params)

    return (ρe + p) / ρ
end

@inline function calc_total_specific_enthalpy(eos::MoistIdealGas, state, aux, params)
    ρ = state.ρ
    ρe = state.ρe

    p = calc_pressure(eos, state, aux, params)

    return (ρe + p) / ρ
end

@inline function calc_heat_capacity_at_constant_pressure(::DryIdealGas, state, params)
    return params.cp_d
end

@inline function calc_heat_capacity_at_constant_pressure(::MoistIdealGas, state, params)
    q_tot = state.ρq / state.ρ
    q_liq = 0 # zero for now
    q_ice = 0 # zero for now
    cp_d = params.cp_d
    cp_v = params.cp_v
    cp_l = params.cp_l
    cp_i = params.cp_i

    cp_m = cp_d + (cp_v - cp_d) * q_tot + (cp_l - cp_v) * q_liq + (cp_i - cp_v) * q_ice

    return cp_m
end

@inline function calc_heat_capacity_at_constant_volume(::DryIdealGas, state, params)
    return params.cv_d
end

@inline function calc_heat_capacity_at_constant_volume(::MoistIdealGas, state, params)
    q_tot = state.ρq / state.ρ
    q_liq = 0 # zero for now
    q_ice = 0 # zero for now
    cv_d = params.cv_d
    cv_v = params.cv_v
    cv_l = params.cv_l
    cv_i = params.cv_i

    cv_m = cv_d + (cv_v - cv_d) * q_tot + (cv_l - cv_v) * q_liq + (cv_i - cv_v) * q_ice

    return cv_m
end

@inline function calc_gas_constant(::DryIdealGas, state, params)
    return params.R_d
end

@inline function calc_gas_constant(::MoistIdealGas, state, params)
    q_tot = state.ρq / state.ρ
    q_liq = 0 # zero for now
    q_ice = 0 # zero for nov
    R_d = params.R_d
    molmass_ratio = params.molmass_ratio

    R_m = R_d * (1 + (molmass_ratio - 1) * q_tot - molmass_ratio * (q_liq + q_ice))

    return R_m
end

@inline function calc_heat_capacity_ratio(eos::AbstractEquationOfState, state, params)
    cp = calc_heat_capacity_at_constant_pressure(eos, state, params)
    cv = calc_heat_capacity_at_constant_volume(eos, state, params)
    γ = cp / cv

    return γ
end

@inline function calc_saturation_vapor_pressure(T, params)
    pₜᵣ = params.pₜᵣ
    R_v = params.R_v
    Tₜᵣ = params.Tₜᵣ
    T_0 = params.T_0
    cp_v = params.cp_v
    cp_l = params.cp_l
    LH_v0 = params.LH_v0

    Δcp = cp_v - cp_l
    pᵥₛ = pₜᵣ * (T / Tₜᵣ)^(Δcp / R_v) * exp((LH_v0 - Δcp * T_0) / R_v * (1 / Tₜᵣ - 1 / T))

    return pᵥₛ
end

@inline function calc_saturation_specific_humidity(ρ, T, params)
    R_v = params.R_v

    pᵥₛ = calc_saturation_vapor_pressure(T, params)
    qt = pᵥₛ / (ρ * R_v * T)

    return qt
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

"""
    AttributeType{V} <: AbstractAttributeType

Type parameter `V` contains the value type of an attribute.

```
(OBJECT = 0, BOOLEAN = 1, LONG = 2, DOUBLE = 3, STRING = 4, DATETIME = 5)
```
"""
struct AttributeType{V} <: AbstractAttributeType
    label::Label
    is_root::Bool
end

# More convenient for reference
const VALUE_TYPE = Proto.AttributeType_ValueType

# Encode the value type (as defined in VALUE_TYPE) in the type parameter
function AttributeType(label::Label, is_root::Bool, value_type::EnumType)
    AttributeType{value_type}(label, is_root)
end

# Define conversion functions
let
    function _as_something(value_type::EnumType)
        return function(t::AttributeType{V}) where {V}
            t.is_root || throw(TypeDBClientException(INVALID_CONCEPT_CASTING, V))
            return AttributeType(Label("attribute"), true, value_type)
        end
    end

    global as_boolean   = _as_something(VALUE_TYPE.BOOLEAN)
    global as_long      = _as_something(VALUE_TYPE.LONG)
    global as_double    = _as_something(VALUE_TYPE.DOUBLE)
    global as_string    = _as_something(VALUE_TYPE.STRING)
    global as_date_time = _as_something(VALUE_TYPE.DATETIME)
end

proto(::AttributeType{V}) where {V} = V

# TODO What is Object value type in Java? It needs to return `false` for that.
is_writable(::AttributeType) = true
is_writable(::AttributeType{VALUE_TYPE.OBJECT}) = false

is_keyable(::AttributeType) = false
is_keyable(::AttributeType{VALUE_TYPE.LONG}) = true

# ------------------------------------------------------------------------
# Remote functions
# ------------------------------------------------------------------------

# The following functions do not need to be ported since they're already
# defined in ThingType.jl:
#   set_supertype
#   get_instances

# Technically, we can reuse the logic from ThingType. Maybe refactor soon.
function get_subtypes(r::RemoteConcept{C,T}) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    concept = r.concept
    req = TypeRequestBuilder.get_subtypes_req(label(concept))
    res = execute(r.transaction, req)
    typs = res.type_res_part.type_get_subtypes_res_part.types

    if is_root(concept) && proto(concept) !== VALUE_TYPE.OBJECT
        return Iterators.filter(
            t -> proto(t) == proto(concept) || label(t) == label(concept),
            (instantiate(t) for t in typs)
        )
    else
        return instantiate.(typs)
    end
end

function get_owners(r::RemoteConcept{C,T}, only_key = false) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    req = AttributeTypeRequestBuilder.get_owners_req(r.concept.label, only_key)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.attribute_type_get_owners_res_part.owners for r in res)))
end

function get_regex(r::RemoteConcept{C,T}) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    req = AttributeTypeRequestBuilder.get_regex_req(r.concept.label)
    res = execute(r.transaction, req)
    regex = res.type_res.attribute_type_get_regex_res.regex
    return isempty(regex) ? nothing : regex
end

function set_regex(r::RemoteConcept{C,T}, regex::Optional{AbstractString}) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    regex_str = regex === nothing ? "" : regex
    req = AttributeTypeRequestBuilder.set_regex_req(r.concept.label, regex_str)
    res = execute(r.transaction, req)
end

proto_attribute_value(value::Bool) = Proto.Attribute_Value(; boolean = value)
proto_attribute_value(value::AbstractFloat) = Proto.Attribute_Value(; double = value)
proto_attribute_value(value::Integer) = Proto.Attribute_Value(; long = value)
proto_attribute_value(value::AbstractString) = Proto.Attribute_Value(; string = value)

# Java client uses epoch of 1970-01-01T00:00:00Z and value is stored as milliseconds
# https://github.com/vaticle/typedb-client-java/blob/fb535f8c9494ec6ff93421e1a962510c6cb46139/concept/thing/AttributeImpl.java#L491
function proto_attribute_value(value::DateTime)
    milliseconds_since_1970 = round(Int64, datetime2unix(value) * 1000)
    return Proto.Attribute_Value(; datetime = milliseconds_since_1970)
end

function Base.get(r::RemoteConcept{C,T}, value) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    req = AttributeTypeRequestBuilder.get_req(r.concept.label, proto_attribute_value(value))
    res = execute(r.transaction, req)
    if hasproperty(res.type_res.attribute_type_get_res, :attribute)
        return Attribute(res.type_res.attribute_type_get_res.attribute)
    else
        return nothing
    end
end

function put(r::RemoteConcept{C,T}, value) where {
    C <: AbstractAttributeType, T <: AbstractCoreTransaction
}
    req = AttributeTypeRequestBuilder.put_req(r.concept.label, proto_attribute_value(value))
    res = execute(r.transaction, req)
    return Attribute(res.type_res.attribute_type_put_res.attribute)
end


text:
using Test
using Example2018

#=
@show Polynomial([1,2,3])
@show evaluate(Polynomial([1,2,3]), 2)
@show Polynomial([1,2,3]) + Polynomial([4,5])
@show 3 * Polynomial([1,2,4,42])
@show deriv(Polynomial([1,1,1,1,1]))
@show deriv(Polynomial([0,0,1,0,0]))
=#

@testset "Polynomial form a ring" begin

z = Polynomial([0])

  for n in 1:2
    deg = rand(0:100)
    #p = Polynomial([1,2,3])
    p = Polynomial(rand(-1000:1000), deg + 1)

    @test p + z == p
  end

@show 9*9

end


text:
n = 7  # rand(1:10)
m = 17 # rand(1:20)
A = matrixdepot("pei", n, m)
B = ones(n,n)
@test issymmetric(A)
@test matrixdepot("pei", n) - triu(B) - tril(B) ≈ zeros(n,n)
println("'pei' passed test...")


text:
using Static

as_bool(b::Bool) = b
as_bool(b::StaticBool) = Bool(b)

function stable_tuple_prod(@nospecialize(x::Tuple))
    @assert length(x) > 0
    x1 = first(x)
    Base.@simd for xi in Base.tail(x)
        x1 *= xi
    end
    return x1
end

as_char(c::Char) = c
as_char(c::StaticInt) = Char(c)


text:
# ------------------------------------------------------------------------------------------
# # Optimizing DiffEq Code
#
# In this notebook we will walk through some of the main tools for optimizing your code in
# order to efficiently solve DifferentialEquations.jl. User-side optimizations are important
# because, for sufficiently difficult problems, most of the time will be spent inside of
# your `f` function, the function you are trying to solve. "Efficient" integrators are those
# that reduce the required number of `f` calls to hit the error tolerance. The main ideas
# for optimizing your DiffEq code, or any Julia function, are the following:
#
# - Make it non-allocating
# - Use StaticArrays for small arrays
# - Use broadcast fusion
# - Make it type-stable
# - Reduce redundant calculations
# - Make use of BLAS calls
# - Optimize algorithm choice
#
# We'll discuss these strategies in the context of small and large systems. Let's start with
# small systems.
#
# ## Optimizing Small Systems (<100 DEs)
#
# Let's take the classic Lorenz system from before. Let's start by naively writing the
# system in its out-of-place form:
# ------------------------------------------------------------------------------------------

function lorenz(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 [dx,dy,dz]
end

# ------------------------------------------------------------------------------------------
# Here, `lorenz` returns an object, `[dx,dy,dz]`, which is created within the body of
# `lorenz`.
#
# This is a common code pattern from high-level languages like MATLAB, SciPy, or R's
# deSolve. However, the issue with this form is that it allocates a vector, `[dx,dy,dz]`, at
# each step. Let's benchmark the solution process with this choice of function:
# ------------------------------------------------------------------------------------------

using DifferentialEquations, BenchmarkTools
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# The BenchmarkTools package's `@benchmark` runs the code multiple times to get an accurate
# measurement. The minimum time is the time it takes when your OS and other background
# processes aren't getting in the way. Notice that in this case it takes about 5ms to solve
# and allocates around 11.11 MiB. However, if we were to use this inside of a real user code
# we'd see a lot of time spent doing garbage collection (GC) to clean up all of the arrays
# we made. Even if we turn off saving we have these allocations.
# ------------------------------------------------------------------------------------------

@benchmark solve(prob,Tsit5(),save_everystep=false)

# ------------------------------------------------------------------------------------------
# The problem of course is that arrays are created every time our derivative function is
# called. This function is called multiple times per step and is thus the main source of
# memory usage. To fix this, we can use the in-place form to ***make our code non-
# allocating***:
# ------------------------------------------------------------------------------------------

function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

# ------------------------------------------------------------------------------------------
# Here, instead of creating an array each time, we utilized the cache array `du`. When the
# inplace form is used, DifferentialEquations.jl takes a different internal route that
# minimizes the internal allocations as well. When we benchmark this function, we will see
# quite a difference.
# ------------------------------------------------------------------------------------------

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,Tsit5(),save_everystep=false)

# ------------------------------------------------------------------------------------------
# There is a 4x time difference just from that change! Notice there are still some
# allocations and this is due to the construction of the integration cache. But this doesn't
# scale with the problem size:
# ------------------------------------------------------------------------------------------

tspan = (0.0,500.0) # 5x longer than before
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5(),save_everystep=false)

# ------------------------------------------------------------------------------------------
# since that's all just setup allocations.
#
# #### But if the system is small we can optimize even more.
#
# Allocations are only expensive if they are "heap allocations". For a more in-depth
# definition of heap allocations, [there are a lot of sources online](http://net-
# informations.com/faq/net/stack-heap.htm). But a good working definition is that heap
# allocations are variable-sized slabs of memory which have to be pointed to, and this
# pointer indirection costs time. Additionally, the heap has to be managed and the garbage
# controllers has to actively keep track of what's on the heap.
#
# However, there's an alternative to heap allocations, known as stack allocations. The stack
# is statically-sized (known at compile time) and thus its accesses are quick. Additionally,
# the exact block of memory is known in advance by the compiler, and thus re-using the
# memory is cheap. This means that allocating on the stack has essentially no cost!
#
# Arrays have to be heap allocated because their size (and thus the amount of memory they
# take up) is determined at runtime. But there are structures in Julia which are stack-
# allocated. `struct`s for example are stack-allocated "value-type"s. `Tuple`s are a stack-
# allocated collection. The most useful data structure for DiffEq though is the
# `StaticArray` from the package
# [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). These arrays have their
# length determined at compile-time. They are created using macros attached to normal array
# expressions, for example:
# ------------------------------------------------------------------------------------------

using StaticArrays
A = @SVector [2.0,3.0,5.0]

# ------------------------------------------------------------------------------------------
# Notice that the `3` after `SVector` gives the size of the `SVector`. It cannot be changed.
# Additionally, `SVector`s are immutable, so we have to create a new `SVector` to change
# values. But remember, we don't have to worry about allocations because this data structure
# is stack-allocated. `SArray`s have a lot of extra optimizations as well: they have fast
# matrix multiplication, fast QR factorizations, etc. which directly make use of the
# information about the size of the array. Thus, when possible they should be used.
#
# Unfortunately static arrays can only be used for sufficiently small arrays. After a
# certain size, they are forced to heap allocate after some instructions and their compile
# time balloons. Thus static arrays shouldn't be used if your system has more than 100
# variables. Additionally, only the native Julia algorithms can fully utilize static arrays.
#
# Let's ***optimize `lorenz` using static arrays***. Note that in this case, we want to use
# the out-of-place allocating form, but this time we want to output a static array:
# ------------------------------------------------------------------------------------------

function lorenz_static(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 @SVector [dx,dy,dz]
end

# ------------------------------------------------------------------------------------------
# To make the solver internally use static arrays, we simply give it a static array as the
# initial condition:
# ------------------------------------------------------------------------------------------

u0 = @SVector [1.0,0.0,0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,Tsit5(),save_everystep=false)

# ------------------------------------------------------------------------------------------
# And that's pretty much all there is to it. With static arrays you don't have to worry
# about allocating, so use operations like `*` and don't worry about fusing operations
# (discussed in the next section). Do "the vectorized code" of R/MATLAB/Python and your code
# in this case will be fast, or directly use the numbers/values.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Implement the out-of-place array, in-place array, and out-of-place static array forms for
# the [Henon-Heiles System](https://en.wikipedia.org/wiki/H%C3%A9non%E2%80%93Heiles_system)
# and time the results.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Optimizing Large Systems
#
# ### Interlude: Managing Allocations with Broadcast Fusion
#
# When your system is sufficiently large, or you have to make use of a non-native Julia
# algorithm, you have to make use of `Array`s. In order to use arrays in the most efficient
# manner, you need to be careful about temporary allocations. Vectorized calculations
# naturally have plenty of temporary array allocations. This is because a vectorized
# calculation outputs a vector. Thus:
# ------------------------------------------------------------------------------------------

A = rand(1000,1000); B = rand(1000,1000); C = rand(1000,1000)
test(A,B,C) = A + B + C
@benchmark test(A,B,C)

# ------------------------------------------------------------------------------------------
# That expression `A + B + C` creates 2 arrays. It first creates one for the output of `A +
# B`, then uses that result array to `+ C` to get the final result. 2 arrays! We don't want
# that! The first thing to do to fix this is to use broadcast fusion. [Broadcast
# fusion](https://julialang.org/blog/2017/01/moredots) puts expressions together. For
# example, instead of doing the `+` operations separately, if we were to add them all at the
# same time, then we would only have a single array that's created. For example:
# ------------------------------------------------------------------------------------------

test2(A,B,C) = map((a,b,c)->a+b+c,A,B,C)
@benchmark test2(A,B,C)

# ------------------------------------------------------------------------------------------
# Puts the whole expression into a single function call, and thus only one array is required
# to store output. This is the same as writing the loop:
# ------------------------------------------------------------------------------------------

function test3(A,B,C)
    D = similar(A)
    @inbounds for i in eachindex(A)
        D[i] = A[i] + B[i] + C[i]
    end
    D
end
@benchmark test3(A,B,C)

# ------------------------------------------------------------------------------------------
# However, Julia's broadcast is syntactic sugar for this. If multiple expressions have a
# `.`, then it will put those vectorized operations together. Thus:
# ------------------------------------------------------------------------------------------

test4(A,B,C) = A .+ B .+ C
@benchmark test4(A,B,C)

# ------------------------------------------------------------------------------------------
# is a version with only 1 array created (the output). Note that `.`s can be used with
# function calls as well:
# ------------------------------------------------------------------------------------------

sin.(A) .+ sin.(B)

# ------------------------------------------------------------------------------------------
# Also, the `@.` macro applys a dot to every operator:
# ------------------------------------------------------------------------------------------

test5(A,B,C) = @. A + B + C #only one array allocated
@benchmark test5(A,B,C)

# ------------------------------------------------------------------------------------------
# Using these tools we can get rid of our intermediate array allocations for many vectorized
# function calls. But we are still allocating the output array. To get rid of that
# allocation, we can instead use mutation. Mutating broadcast is done via `.=`. For example,
# if we pre-allocate the output:
# ------------------------------------------------------------------------------------------

D = zeros(1000,1000);

# ------------------------------------------------------------------------------------------
# Then we can keep re-using this cache for subsequent calculations. The mutating
# broadcasting form is:
# ------------------------------------------------------------------------------------------

test6!(D,A,B,C) = D .= A .+ B .+ C #only one array allocated
@benchmark test6!(D,A,B,C)

# ------------------------------------------------------------------------------------------
# If we use `@.` before the `=`, then it will turn it into `.=`:
# ------------------------------------------------------------------------------------------

test7!(D,A,B,C) = @. D = A + B + C #only one array allocated
@benchmark test7!(D,A,B,C)

# ------------------------------------------------------------------------------------------
# Notice that in this case, there is no "output", and instead the values inside of `D` are
# what are changed (like with the DiffEq inplace function). Many Julia functions have a
# mutating form which is denoted with a `!`. For example, the mutating form of the `map` is
# `map!`:
# ------------------------------------------------------------------------------------------

test8!(D,A,B,C) = map!((a,b,c)->a+b+c,D,A,B,C)
@benchmark test8!(D,A,B,C)

# ------------------------------------------------------------------------------------------
# Some operations require using an alternate mutating form in order to be fast. For example,
# matrix multiplication via `*` allocates a temporary:
# ------------------------------------------------------------------------------------------

@benchmark A*B

# ------------------------------------------------------------------------------------------
# Instead, we can use the mutating form `A_mul_B!` into a cache array to avoid allocating
# the output:
# ------------------------------------------------------------------------------------------

@benchmark A_mul_B!(D,A,B) # same as D = A * B

# ------------------------------------------------------------------------------------------
# For repeated calculations this reduced allocation can stop GC cycles and thus lead to more
# efficient code. Additionally, ***we can fuse together higher level linear algebra
# operations using BLAS***. The package
# [SugarBLAS.jl](https://github.com/lopezm94/SugarBLAS.jl) makes it easy to write higher
# level operations like `alpha*B*A + beta*C` as mutating BLAS calls.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Example Optimization: Gierer-Meinhardt Reaction-Diffusion PDE Discretization
#
# Let's optimize the solution of a Reaction-Diffusion PDE's discretization. In its
# discretized form, this is the ODE:
#
# $$ du = D_1 (A_y u + u A_x) + \frac{au^2}{v} + \bar{u} - \alpha u $$
# $$ dv = D_2 (A_y v + v A_x) + a u^2 + \beta v $$
#
# where $u$, $v$, and $A$ are matrices. Here, we will use the simplified version where $A$
# is the tridiagonal stencil $[1,-2,1]$, i.e. it's the 2D discretization of the LaPlacian.
# The native code would be something along the lines of:
# ------------------------------------------------------------------------------------------

# Generate the constants
p = (1.0,1.0,1.0,10.0,0.001,100.0) # a,α,ubar,β,D1,D2
N = 100
Ax = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
Ax[2,1] = 2.0
Ax[end-1,end] = 2.0
Ay[1,2] = 2.0
Ay[end,end-1] = 2.0

function basic_version!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = r[:,:,1]
  v = r[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  dr[:,:,1] = Du + a*u*u./v + ubar - α*u
  dr[:,:,2] = Dv + a*u*u - β*v
end

a,α,ubar,β,D1,D2 = p
uss = (ubar+β)/α
vss = (a/β)*uss^2
r0 = zeros(100,100,2)
r0[:,:,1] = uss+0.1*rand()
r0[:,:,2] = vss

prob = ODEProblem(basic_version!,r0,(0.0,0.1),p)

# ------------------------------------------------------------------------------------------
# In this version we have encoded our initial condition to be a 3-dimensional array, with
# `u[:,:,1]` being the `A` part and `u[:,:,2]` being the `B` part.
# ------------------------------------------------------------------------------------------

@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# While this version isn't very efficient,
#
# #### We recommend writing the "high-level" code first, and iteratively optimizing it!
#
# The first thing that we can do is get rid of the slicing allocations. The operation
# `r[:,:,1]` creates a temporary array instead of a "view", i.e. a pointer to the already
# existing memory. To make it a view, add `@view`. Note that we have to be careful with
# views because they point to the same memory, and thus changing a view changes the original
# values:
# ------------------------------------------------------------------------------------------

A = rand(4)
@show A
B = @view A[1:3]
B[2] = 2
@show A

# ------------------------------------------------------------------------------------------
# Notice that changing `B` changed `A`. This is something to be careful of, but at the same
# time we want to use this since we want to modify the output `dr`. Additionally, the last
# statement is a purely element-wise operation, and thus we can make use of broadcast fusion
# there. Let's rewrite `basic_version!` to ***avoid slicing allocations*** and to ***use
# broadcast fusion***:
# ------------------------------------------------------------------------------------------

function gm2!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm2!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# Now, most of the allocations are taking place in `Du = D1*(Ay*u + u*Ax)` since those
# operations are vectorized and not mutating. We should instead replace the matrix
# multiplications with `A_mul_B!`. When doing so, we will need to have cache variables to
# write into. This looks like:
# ------------------------------------------------------------------------------------------

Ayu = zeros(N,N)
uAx = zeros(N,N)
Du = zeros(N,N)
Ayv = zeros(N,N)
vAx = zeros(N,N)
Dv = zeros(N,N)
function gm3!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  A_mul_B!(Ayu,Ay,u)
  A_mul_B!(uAx,u,Ax)
  A_mul_B!(Ayv,Ay,v)
  A_mul_B!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm3!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# But our temporary variables are global variables. We need to either declare the caches as
# `const` or localize them. We can localize them by adding them to the parameters, `p`. It's
# easier for the compiler to reason about local variables than global variables.
# ***Localizing variables helps to ensure type stability***.
# ------------------------------------------------------------------------------------------

p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2
function gm4!(dr,r,p,t)
  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  A_mul_B!(Ayu,Ay,u)
  A_mul_B!(uAx,u,Ax)
  A_mul_B!(Ayv,Ay,v)
  A_mul_B!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm4!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# We could then use the BLAS `gemmv` to optimize the matrix multiplications some more, but
# instead let's devectorize the stencil.
# ------------------------------------------------------------------------------------------

p = (1.0,1.0,1.0,10.0,0.001,100.0,N)
function fast_gm!(du,u,p,t)
  a,α,ubar,β,D1,D2,N = p

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for j in 2:N-1
    i = 1
    du[1,j,1] = D1*(2u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = 1
    du[1,j,2] = D2*(2u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,1] = D1*(2u[i-1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
           a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,2] = D2*(2u[i-1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
           a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for i in 2:N-1
    j = 1
    du[i,1,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = 1
    du[i,1,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1,1] = D1*(2u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,1,2] = D2*(2u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]

    i = 1; j = N
    du[1,N,1] = D1*(2u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,N,2] = D2*(2u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = 1
    du[N,1,1] = D1*(2u[i-1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[N,1,2] = D2*(2u[i-1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = N
    du[end,end,1] = D1*(2u[i-1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[end,end,2] = D2*(2u[i-1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
   end
end
prob = ODEProblem(fast_gm!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

# ------------------------------------------------------------------------------------------
# Lastly, we can do other things like multithread the main loops, but these optimizations
# get the last 2x-3x out. The main optimizations which apply everywhere are the ones we just
# performed (though the last one only works if your matrix is a stencil. This is known as a
# matrix-free implementation of the PDE discretization).
#
# This gets us to about 20,000x faster than our original MATLAB/SciPy/R vectorized style
# code!
#
# The last thing to do is then ***optimize our algorithm choice***. We have been using
# `Tsit5()` as our test algorithm, but in reality this problem is a stiff PDE discretization
# and thus one recommendation is to use `CVODE_BDF()`. However, instead of using the default
# dense Jacobian, we should make use of the sparse Jacobian afforded by the problem. The
# Jacobian is the matrix $\frac{df_i}{dr_j}$, where $r$ is read by the linear index (i.e.
# down columns). But since the $u$ variables depend on the $v$, the band size here is large,
# and thus this will not do well with a Banded Jacobian solver. Instead, we utilize sparse
# Jacobian algorithms. `CVODE_BDF` allows us to use a sparse Newton-Krylov solver by setting
# `linear_solver = :GMRES` (see [the solver
# documentation](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Sundials.jl-1),
# and thus we can solve this problem efficiently. Let's see how this scales as we increase
# the integration time.
# ------------------------------------------------------------------------------------------

prob = ODEProblem(fast_gm!,r0,(0.0,10.0),p)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))

prob = ODEProblem(fast_gm!,r0,(0.0,100.0),p)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))

# ------------------------------------------------------------------------------------------
# What's happening is that, because the problem is stiff, the number of steps required by
# the explicit Runge-Kutta method grows rapidly, whereas `CVODE_BDF` is taking large steps.
# Additionally, the `GMRES` linear solver form is quite an efficient way to solve the
# implicit system in this case. This is problem-dependent, and in many cases using a Krylov
# method effectively requires a preconditioner, so you need to play around with testing
# other algorithms and linear solvers to find out what works best with your problem.
#
# ## Conclusion
#
# Julia gives you the tools to optimize the solver "all the way", but you need to make use
# of it. The main thing to avoid is temporary allocations. For small systems, this is
# effectively done via static arrays. For large systems, this is done via in-place
# operations and cache arrays. Either way, the resulting solution can be immensely sped up
# over vectorized formulations by using these principles.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module Interpolation

end # module


text:
import Adapt

"""
    struct BoundaryCondition{C<:BCType, T}

Container for boundary conditions.
"""
struct BoundaryCondition{C<:BCType, T}
    condition :: T
end

"""
    BoundaryCondition(BC, condition)

Construct a boundary condition of type `BC` with a number or array as a `condition`.

Boundary condition types include `Periodic`, `Flux`, `Value`, `Gradient`, and `NormalFlow`.
"""
BoundaryCondition(BC, condition) = BoundaryCondition{BC, typeof(condition)}(condition)

"""
    BoundaryCondition(BC, condition::Function; parameters=nothing, discrete_form=false)

Construct a boundary condition of type `BC` with a function boundary `condition`.

By default, the function boudnary `condition` is assumed to have the 'continuous form'
`condition(ξ, η, t)`, where `t` is time and `ξ` and `η` vary along the boundary.
In particular:

- On `x`-boundaries, `condition(y, z, t)`.
- On `y`-boundaries, `condition(x, z, t)`.
- On `z`-boundaries, `condition(x, y, t)`.

If `parameters` is not `nothing`, then function boundary conditions have the form
`func(ξ, η, t, parameters)`, where `ξ` and `η` are spatial coordinates varying along
the boundary as explained above.

If `discrete_form=true`, the function `condition` is assumed to have the "discrete form",

    `condition(i, j, grid, clock, model_fields)`,

where `i`, and `j` are indices that vary along the boundary. If `discrete_form=true` and
`parameters` is not `nothing`, the function `condition` is called with

    `condition(i, j, grid, clock, model_fields, parameters)`.
"""
function BoundaryCondition(TBC, condition::Function;
                           parameters = nothing,
                           discrete_form = false,
                           field_dependencies=())

    if discrete_form
        field_dependencies != () && error("Cannot set `field_dependencies` when `discrete_form=true`!")
        condition = DiscreteBoundaryFunction(condition, parameters)
    else
        # Note that the boundary :x and location Center, Center are in general incorrect.
        # These are corrected in the FieldBoundaryConditions constructor.
        condition = ContinuousBoundaryFunction(condition, parameters, field_dependencies)
    end

    return BoundaryCondition{TBC, typeof(condition)}(condition)
end

bctype(bc::BoundaryCondition{C}) where C = C

# Adapt boundary condition struct to be GPU friendly and passable to GPU kernels.
Adapt.adapt_structure(to, b::BoundaryCondition{C, A}) where {C<:BCType, A<:AbstractArray} =
    BoundaryCondition(C, Adapt.adapt(to, parent(b.condition)))

#####
##### Some abbreviations to make life easier.
#####

# These type aliases make dispatching on BCs easier (not exported).
const BC   = BoundaryCondition
const FBC  = BoundaryCondition{<:Flux}
const PBC  = BoundaryCondition{<:Periodic}
const NFBC = BoundaryCondition{<:NormalFlow}
const VBC  = BoundaryCondition{<:Value}
const GBC  = BoundaryCondition{<:Gradient}
const ZFBC = BoundaryCondition{Flux, Nothing} # "zero" flux

# More readable BC constructors for the public API.
    PeriodicBoundaryCondition() = BoundaryCondition(Periodic,   nothing)
      NoFluxBoundaryCondition() = BoundaryCondition(Flux,       nothing)
ImpenetrableBoundaryCondition() = BoundaryCondition(NormalFlow, nothing)

      FluxBoundaryCondition(val; kwargs...) = BoundaryCondition(Flux, val; kwargs...)
     ValueBoundaryCondition(val; kwargs...) = BoundaryCondition(Value, val; kwargs...)
  GradientBoundaryCondition(val; kwargs...) = BoundaryCondition(Gradient, val; kwargs...)
NormalFlowBoundaryCondition(val; kwargs...) = BoundaryCondition(NormalFlow, val; kwargs...)

# Support for various types of boundary conditions
@inline getbc(bc::BC{<:NormalFlow, Nothing}, args...) = 0
@inline getbc(bc::BC{C, <:Number},        args...)                  where C = bc.condition
@inline getbc(bc::BC{C, <:AbstractArray}, i, j, grid, clock, model_fields) where C = @inbounds bc.condition[i, j]
@inline getbc(bc::BC{C, <:Function},      i, j, grid, clock, model_fields) where C = bc.condition(i, j, grid, clock, model_fields)

@inline Base.getindex(bc::BC{C, <:AbstractArray}, i, j) where C = getindex(bc.condition, i, j)


text:
# ENV["PYCALL_JL_RUNTIME_PYTHON"] = "venv/bin/python"
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../../"))

include(joinpath(@__DIR__, "../../../PlotUtils/imports.jl"))
using Pandas; const pd = Pandas

using PyCall
TSNE = PyCall.pyimport("sklearn.manifold").TSNE

function get_best_lam(samples)
  n = length(samples)
  Ws = [s[:theta__W] for s in samples]
  Zs = [s[:theta__Z] for s in samples]
  lams = [s[:theta__lam] for s in samples]
  I = length(lams[1])

  idx_best = [estimate_ZWi_index(Zs, Ws, i) for i in 1:I]
  return [Int64.(lams[idx_best[i]][i]) for i in 1:I]
end

function compute_tsne(path_to_simdat; use_complete_data=true, seed=0, verbose=2)
  simdat = BSON.load(path_to_simdat)[:simdat]

  if use_complete_data
    y = Matrix{Float64}.(simdat[:y_complete])
  else
    # TODO
  end

  # number of samples in y
  I = length(y)

  # Create TSNE objects
  tsne = [TSNE(verbose=verbose, random_state=seed) for _ in 1:I]

  # Fit TSNE for each sample
  foreach(i -> tsne[i].fit(y[i]), 1:I)

  return tsne
end

function plot_tsne(tsne, path_to_output, imgdir; method, markersize=10)
  output = BSON.load(path_to_output)
  mkpath(imgdir)

  # number of samples in y
  I = length(tsne)

  # Best lambdas
  best_lambda = get_best_lam(output[:samples][1])

  # Noisy class couunts (shouold be 0)
  num_noisy = zeros(Int, I)

  for i in 1:I
    # This shouold be zero
    num_noisy[i] = sum(best_lambda[i] .== 0)

    tsne_df = pd.DataFrame(tsne[i].embedding_, columns=["comp1", "comp2"])
    tsne_df[Symbol("rfam")] = best_lambda[i]

    sns.pairplot(x_vars=:comp1, y_vars=:comp2, data=tsne_df, hue=:rfam,
                 plot_kws=edgecolor=Dict(:linewidth => 0, :s=> markersize),
                 aspect=1, height=5);
    plt.savefig("$(imgdir)/rfam$(i).pdf", bbox_inches="tight");
    plt.close();
  end

  open("$(imgdir)/num_noisy.txt", "w") do io
    println(io, num_noisy)
  end
end

function gen_imgdir(path_to_output, suffix)
  d = splitdir(path_to_output)[1]
  d1, d2 = splitdir(d)
  return joinpath(d1, "tsne", d2, suffix)
end

## MAIN ###
get_simdat_path(path_to_output) = joinpath(splitdir(path_to_output)[1], "simdat.bson")
results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-8-1"

tsne = [let
          path = joinpath(results_dir, "pmiss0.0-phi0.0-zind$(zind)/simdat.bson")
          compute_tsne(path, use_complete_data=true, seed=0, verbose=2)
        end for zind in (1, 2, 3)]

path_to_outputs = let
  path = [joinpath(results_dir, "pmiss$(pmiss)-phi$(phi)-zind$(zind)/output.bson")
          for phi in (0.0, 1.0, 10.0, 25.0),
              zind in (1, 2, 3),
              pmiss in (0.6, 0.0)]  # pmiss in (0.0, 0.6)
  vec(path)
end

for path_to_output in path_to_outputs
  for method in ["rfam"]
    imgdir = gen_imgdir(path_to_output, method)
    mkpath(imgdir)

    println("Now on: ", imgdir)

    zind = let
      m = match(r"(?<=zind)\d+", path_to_output).match
      parse(Int, m)
    end

    plot_tsne(tsne[zind], path_to_output, imgdir;
              markersize=10, method=method)
  end
end


text:
# interface for skipping null entries

function skipnull_init(f, op, X::NullableArray,
                       ifirst::Int, ilast::Int)
    # Get first non-null element
    ifirst = Base.findnext(x -> x == false, X.isnull, ifirst)
    @inbounds v1 = X.values[ifirst]

    # Get next non-null element
    ifirst = Base.findnext(x -> x == false, X.isnull, ifirst + 1)
    @inbounds v2 = X.values[ifirst]

    # Reduce first two elements
    return op(f(v1), f(v2)), ifirst
end

# sequential map-reduce
function mapreduce_seq_impl_skipnull(f, op, X::NullableArray,
                                     ifirst::Int, ilast::Int)
    # initialize first reduction
    v, i = skipnull_init(f, op, X, ifirst, ilast)

    while i < ilast
        i += 1
        @inbounds isnull = X.isnull[i]
        isnull && continue
        @inbounds entry = X.values[i]
        v = op(v, f(entry))
    end
    return Nullable(v)
end

# pairwise map-reduce
function mapreduce_pairwise_impl_skipnull{T}(f, op, X::NullableArray{T},
                                             ifirst::Int, ilast::Int,
                                            #  n_notnull::Int, blksize::Int)
                                            blksize::Int)
    if ifirst + blksize > ilast
        # fall back to Base implementation if no nulls in block
        # if any(isnull, slice(X, ifirst:ilast))
            return mapreduce_seq_impl_skipnull(f, op, X, ifirst, ilast)
        # else
            # Nullable(Base.mapreduce_seq_impl(f, op, X.values, ifirst, ilast))
        # end
    else
        imid = (ifirst + ilast) >>> 1
        # n_notnull1 = imid - ifirst + 1 - countnz(X.isnull[ifirst:imid])
        # n_notnull2 = ilast - imid - countnz(X.isnull[imid+1:ilast])
        v1 = mapreduce_pairwise_impl_skipnull(f, op, X, ifirst, imid,
                                              blksize)
        v2 = mapreduce_pairwise_impl_skipnull(f, op, X, imid+1, ilast,
                                              blksize)
        return op(v1, v2)
    end
end

# from comment: https://github.com/JuliaLang/julia/pull/16217#issuecomment-223768129
if VERSION < v"0.5.0-dev+4441"
    sum_pairwise_blocksize = Base.sum_pairwise_blocksize
else
    sum_pairwise_blocksize(T) = Base.pairwise_blocksize(T, +)
end

mapreduce_impl_skipnull{T}(f, op, X::NullableArray{T}) =
    mapreduce_seq_impl_skipnull(f, op, X, 1, length(X.values))
mapreduce_impl_skipnull(f, op::typeof(+), X::NullableArray) =
    mapreduce_pairwise_impl_skipnull(f, op, X, 1, length(X.values),
                                   max(128, sum_pairwise_blocksize(f)))

# general mapreduce interface

function _mapreduce_skipnull{T}(f, op, X::NullableArray{T}, missingdata::Bool)
    n = length(X)
    !missingdata && return Nullable(Base.mapreduce_impl(f, op, X.values, 1, n))

    nnull = countnz(X.isnull)
    nnull == n && return Base.mr_empty(f, op, T)
    nnull == n - 1 && return Base.r_promote(op, f(X.values[findnext(x -> x == false), X, 1]))
    # nnull == 0 && return Base.mapreduce_impl(f, op, X, 1, n)

    return mapreduce_impl_skipnull(f, op, X)
end

function Base._mapreduce(f, op, X::NullableArray, missingdata)
    missingdata && return Base._mapreduce(f, op, X)
    Nullable(Base._mapreduce(f, op, X.values))
end

# to fix ambiguity warnings
function Base.mapreduce(f, op::Union{typeof(&), typeof(|)},
                        X::NullableArray, skipnull::Bool = false)
    missingdata = any(isnull, X)
    if skipnull
        return _mapreduce_skipnull(f, op, X, missingdata)
    else
        return Base._mapreduce(f, op, X, missingdata)
    end
end

if VERSION >= v"0.5.0-dev+3701"
    const specialized_binary = identity
else
    const specialized_binary = Base.specialized_binary
end

"""
    mapreduce(f, op::Function, X::NullableArray; [skipnull::Bool=false])

Map a function `f` over the elements of `X` and reduce the result under the
operation `op`. One can set the behavior of this method to skip the null entries
of `X` by setting the keyword argument `skipnull` equal to true. If `skipnull`
behavior is enabled, `f` will be automatically lifted over the elements of `X`.
Note that, in general, mapreducing over a `NullableArray` will return a
`Nullable` object regardless of whether `skipnull` is set to `true` or `false`.
"""
function Base.mapreduce(f, op::Function, X::NullableArray;
                        skipnull::Bool = false)
    missingdata = any(isnull, X)
    if skipnull
        return _mapreduce_skipnull(f, specialized_binary(op),
                                   X, missingdata)
    else
        return Base._mapreduce(f, specialized_binary(op), X, missingdata)
    end
end

function Base.mapreduce(f, op, X::NullableArray; skipnull::Bool = false)
    missingdata = any(isnull, X)
    if skipnull
        return _mapreduce_skipnull(f, op, X, missingdata)
    else
        return Base._mapreduce(f, op, X, missingdata)
    end
end

"""
    mapreduce(f, op::Function, X::NullableArray; [skipnull::Bool=false])

Reduce `X`under the operation `op`. One can set the behavior of this method to
skip the null entries of `X` by setting the keyword argument `skipnull` equal
to true. If `skipnull` behavior is enabled, `f` will be automatically lifted
over the elements of `X`. Note that, in general, mapreducing over a
`NullableArray` will return a `Nullable` object regardless of whether `skipnull`
is set to `true` or `false`.
"""
Base.reduce(op, X::NullableArray; skipnull::Bool = false) =
    mapreduce(identity, op, X; skipnull = skipnull)

# standard reductions

for (fn, op) in ((:(Base.sum), +),
                 (:(Base.prod), *),
                 (:(Base.minimum), scalarmin),
                 (:(Base.maximum), scalarmax))
    @eval begin
        $fn(f::Function, X::NullableArray; skipnull::Bool = false) =
                mapreduce(f, $op, X; skipnull = skipnull)
        $fn(X::NullableArray; skipnull::Bool = false) =
            mapreduce(identity, $op, X; skipnull = skipnull)
    end
end

for (fn, f, op) in ((:(Base.sumabs), abs, +),
                    (:(Base.sumabs2), abs2, +))
    @eval $fn(X::NullableArray; skipnull::Bool = false) =
        mapreduce($f, $op, X; skipnull=skipnull)
end

# internal methods for Base.minimum and Base.maximum
for op in (scalarmin, scalarmax)
    @eval begin
        function Base._mapreduce{T}(::typeof(identity), ::$(typeof(op)),
                                    X::NullableArray{T}, missingdata)
            missingdata && return Nullable{T}()
            Nullable(Base._mapreduce(identity, $op, X.values))
        end
    end
end

function Base.mapreduce_impl{T}(f, op::typeof(scalarmin), X::NullableArray{T},
                                first::Int, last::Int)
    i = first
    v = f(X[i])
    i += 1
    while i <= last
        @inbounds x = f(X[i])
        if isnull(x) | isnull(v)
            return Nullable{eltype(x)}()
        elseif x.value < v.value
            v = x
        end
        i += 1
    end
    return v
end

function Base.mapreduce_impl{T}(f, op::typeof(scalarmax), X::NullableArray{T},
                                first::Int, last::Int)
    i = first
    v = f(X[i])
    i += 1
    while i <= last
        @inbounds x = f(X[i])
        if isnull(x) | isnull(v)
            return Nullable{eltype(x)}()
        elseif x.value > v.value
            v = x
        end
        i += 1
    end
    return v
end

function Base.extrema{T}(X::NullableArray{T}; skipnull::Bool = false)
    length(X) > 0 || throw(ArgumentError("collection must be non-empty"))
    vmin = Nullable{T}()
    vmax = Nullable{T}()
    @inbounds for i in 1:length(X)
        x = X.values[i]
        null = X.isnull[i]
        if skipnull && null
            continue
        elseif null
            return (Nullable{T}(), Nullable{T}())
        elseif isnull(vmax) # Equivalent to isnull(vmin)
            vmax = vmin = Nullable(x)
        elseif x > vmax.value
            vmax = Nullable(x)
        elseif x < vmin.value
            vmin = Nullable(x)
        end
    end
    return (vmin, vmax)
end


text:
module SIMDTypes

# using Static: StaticInt

struct Bit; data::Bool; end # Dummy for Ptr
# @inline Base.convert(::Type{Bool}, b::Bit) = getfield(b, :data)

const FloatingTypes = Union{Float16,Float32,Float64}
const SignedHW = Union{Int8,Int16,Int32,Int64}
const UnsignedHW = Union{UInt8,UInt16,UInt32,UInt64}
const IntegerTypesHW = Union{SignedHW,UnsignedHW}
# const IntegerTypes = Union{StaticInt,IntegerTypesHW}

const NativeTypesExceptBitandFloat16 = Union{Bool,Base.HWReal}
const NativeTypesExceptBit = Union{Bool,Base.HWReal,Float16}
const NativeTypesExceptFloat16 = Union{Bool,Base.HWReal,Bit}
const NativeTypes = Union{NativeTypesExceptBit, Bit}

const _Vec{W,T<:Number} = NTuple{W,Core.VecElement{T}}

end


text:
#using Distributions

function jfpGMM{T1<:Number, T2<:Number}(x::Array{Float64, 2},
  prop::Vector{Float64},
  mu::Array{Float64, 2},
  sigma::Array{Array, 1},
  k::Int64,
  df::T1,
  lambda::T2,
  itermax::Int64 = 300,
  tol::Float64 = 1e-06
  )
  # mu a k by d matrix
  # sigma a [d,d,k] array
  const n::Int64 = size(x, 1)
  const d::Int64 = size(x, 2)

  delta::Float64 = 1.0
  prop_old = copy(prop)
  mu_old = deepcopy(mu)
  sigma_old = deepcopy(sigma)

  stp::Int64 = 0
  pdf_est = Array{Float64,2}(n,k)
  prob0 = Array{Float64,2}(n,k)
  while stp < itermax
    # E step
    for i in 1:k
      @inbounds tmp_mu = mu_old[i,:]
      @inbounds tmp_sigma = sigma_old[i]
      @inbounds pdf_est[:,i] = pdf(MvNormal(tmp_mu, PDMat(cholfact(Hermitian(tmp_sigma)))),x')
      @inbounds prob0[:,i] = pdf_est[:,i].*prop_old[i]
    end
    h_est = prob0./sum(prob0,2)

    # M step
    mu_new = *(h_est',x)./sum(h_est,1)'
    sigma_new = deepcopy(sigma_old)
    for i in 1:k
      @inbounds sigma_new[i] = *((x-reshape(kron(mu[i,:], rep(1,times = n)), (n,d)))',
                      diagm(h_est[:,i]),
                      (x-reshape(kron(mu[i,:], rep(1,times = n)), (n,d))))./sum(h_est[:,i])
    end
    prop_new = vec((sum(h_est,1).-lambda*df)./(n-k*lambda*df))

    if any(prop_new .< 1e-04)
      prop_new[prop_new .< 1e-04] = 0
    end
    prop_new = prop_new./sum(prop_new)

    delta = sum(abs(prop_new .- prop_old))
    if delta < tol
      break
    end

    if any(prop_new .== 0)
        keepidx = find(x -> x > 0, prop_new)
        k = length(keepidx)

        prop_old = prop_new[keepidx]
        mu_old = mu_new[keepidx,:]
        sigma_old = sigma_new[keepidx]
        pdf_est = pdf_est[:,keepidx]
        prob0 = h_est[:,keepidx]
        h_est = h_est[:,keepidx]
        delta = 1
      else
        prop_old = prop_new
        mu_old = mu_new
        sigma_old = sigma_new
    end

    if k <= 1
      break
    end
    stp += 1
  end
  tag = mapslices(indmax, h_est, 2)
  return Dict("k" => k,
              "prop" => prop_old,
              "mu" => mu_old,
              "sigma" => sigma_old,
              "pdf_est" => pdf_est,
              "ll" => sum(log.(sum(prob0,2))),
              "cluster" => tag)
end

max_stars_repo_name: hillarykoch/pGMCM

max_issues_repo_name: hillarykoch/pGMCM

max_forks_repo_name: hillarykoch/pGMCM


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -13.110901270484293 <= q <= 11.380819522292354)

                     add_NL_constraint(m, :(swish(-0.7773677229963925 + -0.31823111196481735*swish(0.13898320039278467 + -0.4573031975765818*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.05869662098445039*swish(0.7177153834386609 + 0.046197114005355644*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.31920840200501477*swish(-0.05621054225465549 + -0.4685096496776473*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) + swish(-0.2800668423712529 + 0.5241807567635033*swish(0.13898320039278467 + -0.4573031975765818*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.9218858277065758*swish(0.7177153834386609 + 0.046197114005355644*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.0068491484591501894*swish(-0.05621054225465549 + -0.4685096496776473*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) + swish(-0.0923150493931284 + -0.18688964829603716*swish(0.13898320039278467 + -0.4573031975765818*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.027363793963717686*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.7118578003946046*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + 0.4699561510424055*swish(0.7177153834386609 + 0.046197114005355644*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + 0.8341427099559482*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + -0.5961814570241306*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3]))) + -0.26284635568255643*swish(-0.05621054225465549 + -0.4685096496776473*swish(-0.1658975405753922 + 0.8948668185465118*$(x[1]) + 0.5677868700589963*$(x[2]) + -0.004745367697401193*$(x[3])) + -0.723697578017108*swish(0.8084147176169236 + 0.6887886991760568*$(x[1]) + 0.9976411528938542*$(x[2]) + 0.9554800615635051*$(x[3])) + 0.643520408356284*swish(-0.23665048660081967 + 0.22941639672835912*$(x[1]) + 0.6694807543091978*$(x[2]) + 0.08612932953468189*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Atomic solution for an SDE.

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `p`: current solution of p
* `p̅`: previous solution of p
* `p̃`: compensated summation error of p
* `ΔW`: Wiener process driving the stochastic process q
* `ΔZ`: Wiener process driving the stochastic process q
* `K`:  integer parameter defining the truncation of the increments of the Wiener process (for strong solutions)
"""
mutable struct AtomicSolutionPSDE{DT,TT} <: AtomicSolution{DT,TT}
    t::TT
    t̅::TT

    q::Vector{DT}
    q̅::Vector{DT}
    q̃::Vector{DT}

    p::Vector{DT}
    p̅::Vector{DT}
    p̃::Vector{DT}

    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    K::Int

    function AtomicSolutionPSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real}
        new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nm), zeros(DT, nm))
    end
end

AtomicSolutionPSDE(DT, TT, nd, nm) = AtomicSolutionPSDE{DT, TT}(nd, nm)

function set_solution!(asol::AtomicSolutionPSDE, sol)
    t, q, p = sol
    asol.t  = t
    asol.q .= q
    asol.p .= p
end

function set_increments!(asol::AtomicSolutionPSDE, incs)
    ΔW, ΔZ = incs
    asol.ΔW .= ΔW
    asol.ΔZ .= ΔZ
end

function get_solution(asol::AtomicSolutionPSDE)
    (asol.t, asol.q, asol.p)
end

function get_increments(asol::AtomicSolutionPSDE)
    (asol.ΔW, asol.ΔZ)
end

function get_increments!(asol::AtomicSolutionPSDE, ΔW, ΔZ)
    ΔW .= asol.ΔW
    ΔZ .= asol.ΔZ
end

function CommonFunctions.reset!(asol::AtomicSolutionPSDE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.p̅ .= asol.p
    asol.t += Δt
end

function update!(asol::AtomicSolutionPSDE{DT}, y::Vector{DT}, z::Vector{DT}) where {DT}
    for k in eachindex(y,z)
        update!(asol, y[k], z[k], k)
    end
end

function update!(asol::AtomicSolutionPSDE{DT}, y::DT, z::DT, k::Int) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(y, asol.q[k], asol.q̃[k])
    asol.p[k], asol.p̃[k] = compensated_summation(z, asol.p[k], asol.p̃[k])
end