text:
using Oceananigans.Architectures: architecture, device_event
using Oceananigans.BoundaryConditions: DefaultBoundaryCondition
using Oceananigans.BuoyancyModels: ∂z_b, top_buoyancy_flux
using Oceananigans.Operators: ℑzᵃᵃᶜ

function hydrostatic_turbulent_kinetic_energy_tendency end

struct TKEBasedVerticalDiffusivity{TD, CK, CD, CL, CQ} <: AbstractTurbulenceClosure{TD}
    diffusivity_scaling :: CK
    dissipation_parameter :: CD
    mixing_length_parameter :: CL
    surface_model :: CQ

    function TKEBasedVerticalDiffusivity{TD}(
        diffusivity_scaling :: CK,
        dissipation_parameter :: CD,
        mixing_length_parameter :: CL,
        surface_model :: CQ) where {TD, CK, CD, CL, CQ}

        return new{TD, CK, CD, CL, CQ}(diffusivity_scaling, dissipation_parameter, mixing_length_parameter, surface_model)
    end
end

"""
    TKEBasedVerticalDiffusivity(FT=Float64;
                                diffusivity_scaling = RiDependentDiffusivityScaling{FT}(),
                                dissipation_parameter = 2.91,
                                mixing_length_parameter = 1.16,
                                surface_model = TKESurfaceFlux{FT}(),
                                time_discretization::TD = ExplicitTimeDiscretization())

Returns the `TKEBasedVerticalDiffusivity` turbulence closure for vertical mixing by
small-scale ocean turbulence based on the prognostic evolution of subgrid
Turbulent Kinetic Energy (TKE).

`TKEBasedVerticalDiffusivity` is a downgradient, diffusive
closure formulated with three different eddy diffusivities for momentum, tracers, and TKE.
Each eddy diffusivity is the product of a diffusivity "scaling", a mixing length, and a turbulent
velocity scale which is the square root of the local TKE, such that

```math
Kᵠ = Cᵠ ℓ √e
```

where `Kᵠ` is the eddy diffusivity of `ϕ` where `ϕ` is either `u` (for momentum) `c` (for tracers), or
`e` (for TKE). `Cᵠ` is the diffusivity scaling for `ϕ`, `ℓ` is the mixing length
and `√e` is the turbulent velocity scale. The mixing length `ℓ` is modeled as

```math
ℓ = min(ℓᵇ, ℓᶻ)
```

where `ℓᵇ = Cᵇ * √e / N` and `ℓᶻ` is the distance to the nearest boundary.
`TKEBasedVerticalDiffusivity` also invokes a model for the flux of TKE across the numerical
ocean surface due to unstable buoyancy forcing and wind stress.

The `TKEBasedVerticalDiffusivity` is formulated in terms of 12 free parameters. These parameters
are _experimentally_ calibrated against large eddy simulations of ocean surface boundary layer turbulence
in idealized scenarios involving monotonic boundary layer deepening into variable stratification
due to constant surface momentum fluxes and/or destabilizing surface buoyancy flux.
This calibration has not been peer-reviewed, may be inaccurate and imperfect, and may not
be appropriate for three-dimensional ocean simulations.

See https://github.com/CliMA/LESbrary.jl for more information about the large eddy simulations.

The calibration procedure is not documented and is part of ongoing research.
The calibration was performed using a combination of Markov Chain Monte Carlo (MCMC)-based simulated
annealing and noisy Ensemble Kalman Inversion methods.

The one positional argument determines the floating point type of the free parameters
of `TKEBasedVerticalDiffusivity`. The default is `Float64`.

Keyword arguments
=================

* `diffusivity_scaling` : A group of parameters that scale the eddy diffusivity for momentum, tracers, and TKE.
                          The default is `RiDependentDiffusivityScaling{FT}()`, which represents a group of
                          parameters that implement a "smoothed step function" scaling that varies with the
                          local gradient Richardson number `Ri = ∂z(b) / (∂z(u)² + ∂z(v)²)`.

* `dissipation_parameter` : Parameter `Cᴰ` in the closure `ϵ = Cᴰ * e^3/2 / ℓ` that models the dissipation of TKE,
                            `ϵ`, appearing in the TKE evolution equation. The default is 2.91 via calibration
                            against large eddy simulations.

* `mixing_length_parameter` : Parameter `Cᵇ` that multiplies the "buoyancy mixing length" `ℓᵇ = Cᵇ * √e / N`,
                            that appears in `TKEBasedVerticalDiffusivity`'s mixing length model.
                            The default is 1.16 via calibration against large eddy simulations.

* `time_discretization` : Either `ExplicitTimeDiscretization` or `VerticallyImplicitTimeDiscretization`.

"""
function TKEBasedVerticalDiffusivity(FT=Float64;
                                     diffusivity_scaling = RiDependentDiffusivityScaling{FT}(),
                                     dissipation_parameter = 2.91,
                                     mixing_length_parameter = 1.16,
                                     surface_model = TKESurfaceFlux{FT}(),
                                     time_discretization::TD = ExplicitTimeDiscretization()) where TD

    @warn "TKEBasedVerticalDiffusivity is an experimental turbulence closure that \n" *
          "is unvalidated and whose default parameters are not calibrated for \n" *
          "realistic ocean conditions or for use in a three-dimensional \n" *
          "simulation. Use with caution and report bugs and problems with physics \n" *
          "to https://github.com/CliMA/Oceananigans.jl/issues. \n\n" *

          "Note: `TKEBasedVerticalDiffusivity` will generally produce \n" *
          "incorrect results if its parameters are altered _after_ \n" *
          "model instantiation."

    dissipation_parameter = convert(FT, dissipation_parameter)
    mixing_length_parameter = convert(FT, mixing_length_parameter)
    diffusivity_scaling = convert_eltype(FT, diffusivity_scaling)
    surface_model = convert_eltype(FT, surface_model)

    return TKEBasedVerticalDiffusivity{TD}(diffusivity_scaling,
                                           dissipation_parameter,
                                           mixing_length_parameter,
                                           surface_model)
end

const TKEVD = TKEBasedVerticalDiffusivity

"""
    struct RiDependentDiffusivityScaling{FT}

A diffusivity model in which momentum, tracers, and TKE
each have Richardson-number-dependent diffusivities with
free parameter of type `FT`.

The Richardson number is

    ``Ri = ∂z B / ( (∂z U)² + (∂z V)² )`` ,

where ``B`` is buoyancy and ``∂z`` denotes a vertical derviative.
The Richardson-number dependent diffusivities are multiplied by the stability
function

    ``σ(Ri) = σ⁻ + (σ⁺ - σ⁻) * step(Ri, Riᶜ, Riʷ)``

where ``σ⁰``, ``σᵟ``, ``Riᶜ``, and ``Riʷ`` are free parameters,
and ``step`` is a smooth step function defined by

    ``step(x, c, w) = (1 + tanh((x - c) / w)) / 2``.

The 8 free parameters in `RiDependentDiffusivityScaling` have been _experimentally_ calibrated
against large eddy simulations of ocean surface boundary layer turbulence in idealized
scenarios involving monotonic boundary layer deepening into variable stratification
due to constant surface momentum fluxes and/or destabilizing surface buoyancy flux.
See https://github.com/CliMA/LESbrary.jl for more information about the large eddy simulations.
The calibration was performed using a combination of Markov Chain Monte Carlo (MCMC)-based simulated
annealing and noisy Ensemble Kalman Inversion methods.
"""
Base.@kwdef struct RiDependentDiffusivityScaling{FT}
    Cᴷu⁻  :: FT = 0.15
    Cᴷu⁺  :: FT = 0.73
    Cᴷc⁻  :: FT = 0.40
    Cᴷc⁺  :: FT = 1.77
    Cᴷe⁻  :: FT = 0.13
    Cᴷe⁺  :: FT = 1.22
    CᴷRiʷ :: FT = 0.72
    CᴷRiᶜ :: FT = 0.76
end

"""
    struct TKESurfaceFlux{FT}

A model for the flux of TKE across the numerical ocean surface with
free parameters of type `FT`, parameterized in
terms of the kinematic surface stress and buoyancy flux:

```math
Qᵉ = - Cᴰ * (Cᵂu★ * u★³ + CᵂwΔ * w★³)
```

where `Qᵉ` is the surface flux of TKE, `Cᴰ = TKEBasedVerticalDiffusivity.dissipation_parameter`,
`u★ = (Qᵘ^2 + Qᵛ^2)^(1/4)` is the friction velocity and `w★ = (Qᵇ * Δz)^(1/3)` is the
turbulent velocity scale associated with the surface vertical grid spacing `Δz` and the
surface buoyancy flux `Qᵇ`.

The 2 free parameters in `TKESurfaceFlux` have been _experimentally_ calibrated
against large eddy simulations of ocean surface boundary layer turbulence in idealized
scenarios involving monotonic boundary layer deepening into variable stratification
due to constant surface momentum fluxes and/or destabilizing surface buoyancy flux.
See https://github.com/CliMA/LESbrary.jl for more information about the large eddy simulations.
The calibration was performed using a combination of Markov Chain Monte Carlo (MCMC)-based simulated
annealing and noisy Ensemble Kalman Inversion methods.
"""
Base.@kwdef struct TKESurfaceFlux{FT}
    Cᵂu★ :: FT = 3.62
    CᵂwΔ :: FT = 1.31
end

@inline function top_tke_flux(i, j, grid, surface_model::TKESurfaceFlux, closure,
                              buoyancy, fields, top_tracer_bcs, top_velocity_bcs, clock)

    wΔ³ = top_convective_turbulent_velocity³(i, j, grid, clock, fields, buoyancy, top_tracer_bcs)
    u★ = friction_velocity(i, j, grid, clock, fields, top_velocity_bcs)

    Cᴰ = closure.dissipation_parameter
    Cᵂu★ = surface_model.Cᵂu★
    CᵂwΔ = surface_model.CᵂwΔ

    return - Cᴰ * (Cᵂu★ * u★^3 + CᵂwΔ * wΔ³)
end

for S in (:RiDependentDiffusivityScaling, :TKESurfaceFlux)
    @eval @inline convert_eltype(::Type{FT}, s::$S) where FT = $S{FT}(; Dict(p => getproperty(s, p) for p in propertynames(s))...)
    @eval @inline convert_eltype(::Type{FT}, s::$S{FT}) where FT = s
end

#####
##### TKE top boundary condition
#####

""" Computes the friction velocity based on fluxes of u and v. """
@inline function friction_velocity(i, j, grid, clock, fields, velocity_bcs)
    FT = eltype(grid)
    Qᵘ = getbc(velocity_bcs.u, i, j, grid, clock, fields)
    Qᵛ = getbc(velocity_bcs.v, i, j, grid, clock, fields)
    return sqrt(sqrt(Qᵘ^2 + Qᵛ^2))
end

@inline function top_convective_turbulent_velocity³(i, j, grid, clock, fields, buoyancy, tracer_bcs)
    FT = eltype(grid)
    Qᵇ = top_buoyancy_flux(i, j, grid, buoyancy, tracer_bcs, clock, fields)
    Δz = Δzᵃᵃᶜ(i, j, grid.Nz, grid)
    return max(zero(FT), Qᵇ) * Δz
end

@inline function top_tke_flux(i, j, grid, clock, fields, parameters)
    buoyancy = parameters.buoyancy
    top_tracer_bcs = parameters.top_tracer_boundary_conditions
    top_velocity_bcs = parameters.top_velocity_boundary_conditions
    closure = parameters.closure # problematic because model.closure can be changed.

    return top_tke_flux(i, j, grid, closure.surface_model, closure,
                        buoyancy, fields, top_tracer_bcs, top_velocity_bcs, clock)
end

#####
##### Utilities for model constructors
#####

""" Infer tracer boundary conditions from user_bcs and tracer_names. """
function top_tracer_boundary_conditions(grid, tracer_names, user_bcs)
    user_bc_names = keys(user_bcs)
    default_top_bc = DefaultBoundaryCondition(topology(grid, 3), Center)

    tracer_bcs = Tuple(name ∈ user_bc_names ? user_bcs[name].top : default_top_bc
                       for name in tracer_names)

    return NamedTuple{tracer_names}(tracer_bcs)
end

""" Infer velocity boundary conditions from user_bcs and tracer_names. """
function top_velocity_boundary_conditions(grid, user_bcs)

    user_bc_names = keys(user_bcs)

    u_top_bc = :u ∈ user_bc_names ? user_bcs.u.top : DefaultBoundaryCondition(topology(grid, 3), Center)
    v_top_bc = :v ∈ user_bc_names ? user_bcs.v.top : DefaultBoundaryCondition(topology(grid, 3), Center)

    return (u=u_top_bc, v=v_top_bc)
end

""" Add TKE boundary conditions specific to TKEBasedVerticalDiffusivity. """
function add_closure_specific_boundary_conditions(closure::TKEVD,
                                                  user_bcs,
                                                  grid,
                                                  tracer_names,
                                                  buoyancy)

    top_tracer_bcs = top_tracer_boundary_conditions(grid, tracer_names, user_bcs)
    top_velocity_bcs = top_velocity_boundary_conditions(grid, user_bcs)

    parameters = (buoyancy = buoyancy,
                  top_tracer_boundary_conditions = top_tracer_bcs,
                  top_velocity_boundary_conditions = top_velocity_bcs,
                  closure = closure)

    top_tke_bc = FluxBoundaryCondition(top_tke_flux, discrete_form=true, parameters=parameters)

    if :e ∈ keys(user_bcs)
        @warn "Replacing top boundary conditions for tracer `e` with " *
              "boundary condition specific to $(typeof(closure).name.wrapper)"

        e_bcs = user_bcs[:e]

        tke_bcs = TracerBoundaryConditions(grid,
                                           top = top_tke_bc,
                                           bottom = e_bcs.bottom,
                                           north = e_bcs.north,
                                           south = e_bcs.south,
                                           east = e_bcs.east,
                                           west = e_bcs.west)

    else
        tke_bcs = TracerBoundaryConditions(grid, top=top_tke_bc)
    end

    new_boundary_conditions = merge(user_bcs, (e = tke_bcs,))

    return new_boundary_conditions
end

function DiffusivityFields(arch, grid, tracer_names, bcs, closure::TKEVD)

    Kᵘ_bcs = :Kᵘ ∈ keys(bcs) ? bcs[:Kᵘ] : DiffusivityBoundaryConditions(grid)
    Kᶜ_bcs = :Kᶜ ∈ keys(bcs) ? bcs[:Kᶜ] : DiffusivityBoundaryConditions(grid)
    Kᵉ_bcs = :Kᵉ ∈ keys(bcs) ? bcs[:Kᵉ] : DiffusivityBoundaryConditions(grid)

    Kᵘ = CenterField(arch, grid, Kᵘ_bcs)
    Kᶜ = CenterField(arch, grid, Kᶜ_bcs)
    Kᵉ = CenterField(arch, grid, Kᵉ_bcs)

    return (; Kᵘ, Kᶜ, Kᵉ)
end

function with_tracers(tracer_names, closure::TKEVD)
    :e ∈ tracer_names || error("Tracers must contain :e to represent turbulent kinetic energy for `TKEBasedVerticalDiffusivity`.")
    return closure
end

#####
##### Diffusivity field utilities
#####

function calculate_diffusivities!(diffusivities, arch, grid, closure::TKEVD, buoyancy, velocities, tracers)

    e = tracers.e

    event = launch!(arch, grid, :xyz,
                    calculate_tke_diffusivities!, diffusivities, grid, closure, e, velocities, tracers, buoyancy,
                    dependencies=device_event(arch))

    wait(device(arch), event)

    return nothing
end

@kernel function calculate_tke_diffusivities!(diffusivities, grid, closure, e, velocities, tracers, buoyancy)
    i, j, k, = @index(Global, NTuple)
    @inbounds begin
        diffusivities.Kᵘ[i, j, k] = Kuᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
        diffusivities.Kᶜ[i, j, k] = Kcᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
        diffusivities.Kᵉ[i, j, k] = Keᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
    end
end

#####
##### Mixing length
#####

@inline surface(i, j, k, grid)                = znode(Center(), Center(), Face(), i, j, grid.Nz+1, grid)
@inline bottom(i, j, k, grid)                 = znode(Center(), Center(), Face(), i, j, 1, grid)
@inline depthᶜᶜᶜ(i, j, k, grid)               = surface(i, j, k, grid) - znode(Center(), Center(), Center(), i, j, k, grid)
@inline height_above_bottomᶜᶜᶜ(i, j, k, grid) = znode(Center(), Center(), Center(), i, j, k, grid) - bottom(i, j, k, grid)

@inline wall_vertical_distanceᶜᶜᶜ(i, j, k, grid) = min(depthᶜᶜᶜ(i, j, k, grid), height_above_bottomᶜᶜᶜ(i, j, k, grid))

@inline function sqrt_∂z_b(i, j, k, grid, buoyancy, tracers)
    FT = eltype(grid)
    N² = ∂z_b(i, j, k, grid, buoyancy, tracers)
    N²⁺ = max(zero(FT), N²)
    return sqrt(N²⁺)
end

@inline function buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    FT = eltype(grid)
    Cᵇ = closure.mixing_length_parameter
    N⁺ = ℑzᵃᵃᶜ(i, j, k, grid, sqrt_∂z_b, buoyancy, tracers)

    @inbounds e⁺ = max(zero(FT), e[i, j, k])

    return @inbounds ifelse(N⁺ == 0, FT(Inf), Cᵇ * sqrt(e⁺) / N⁺)
end

@inline function dissipation_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    ℓᶻ = wall_vertical_distanceᶜᶜᶜ(i, j, k, grid)
    ℓᵇ = buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    ℓ = min(ℓᶻ, ℓᵇ)
    ℓ_min = Δzᵃᵃᶜ(i, j, k, grid) / 2 # minimum mixing length...
    return max(ℓ_min, ℓ)
end

#####
##### Diffusivities
#####

@inline ϕ²(i, j, k, grid, ϕ, args...) = ϕ(i, j, k, grid, args...)^2

@inline function Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    FT = eltype(grid)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.v)
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return ifelse(N² == 0, zero(FT), N² / (∂z_u² + ∂z_v²))
end

@inline step(x, c, w) = (1 + tanh((x - c) / w)) / 2

@inline scale(Ri, σ⁻, σ⁺, c, w) = σ⁻ + (σ⁺ - σ⁻) * step(Ri, c, w)

@inline function momentum_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.diffusivity_scaling.Cᴷu⁻,
                 closure.diffusivity_scaling.Cᴷu⁺,
                 closure.diffusivity_scaling.CᴷRiᶜ,
                 closure.diffusivity_scaling.CᴷRiʷ)
end

@inline function tracer_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.diffusivity_scaling.Cᴷc⁻,
                 closure.diffusivity_scaling.Cᴷc⁺,
                 closure.diffusivity_scaling.CᴷRiᶜ,
                 closure.diffusivity_scaling.CᴷRiʷ)
end

@inline function TKE_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.diffusivity_scaling.Cᴷe⁻,
                 closure.diffusivity_scaling.Cᴷe⁺,
                 closure.diffusivity_scaling.CᴷRiᶜ,
                 closure.diffusivity_scaling.CᴷRiʷ)
end

@inline turbulent_velocity(i, j, k, grid, e) = @inbounds sqrt(max(zero(eltype(grid)), e[i, j, k]))

@inline function unscaled_eddy_diffusivityᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    ℓ = dissipation_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    u★ = turbulent_velocity(i, j, k, grid, e)
    return ℓ * u★
end

@inline function Kuᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
    K = unscaled_eddy_diffusivityᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    σu = momentum_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    return σu * K
end

@inline function Kcᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
    K = unscaled_eddy_diffusivityᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    σc = tracer_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    return σc * K
end

@inline function Keᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy)
    K = unscaled_eddy_diffusivityᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    σe = TKE_diffusivity_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    return σe * K
end

#####
##### Terms in the turbulent kinetic energy equation, all at cell centers
#####

@inline ϕ²(i, j, k, grid, ϕ) = ϕ(i, j, k, grid)^2

@inline function shear_production(i, j, k, grid, closure, clock, velocities, tracers, buoyancy, diffusivities)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.v)
    Ku = Kuᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy)
    return Ku * (∂z_u² + ∂z_v²)
end

@inline function buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Kc = Kcᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy)
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return - Kc * N²
end

@inline function dissipation(i, j, k, grid, closure, tracers, buoyancy)
    e = tracers.e
    FT = eltype(grid)
    three_halves = FT(3/2)
    @inbounds ẽ³² = abs(e[i, j, k])^three_halves

    ℓ = dissipation_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    Cᴰ = closure.dissipation_parameter

    return Cᴰ * ẽ³² / ℓ
end

#####
##### Viscous flux, diffusive fluxes, plus shenanigans for diffusive fluxes of TKE (eg TKE "transport")
#####

# Special "index type" alternative to Val for dispatch
struct TKETracerIndex{N} end
@inline TKETracerIndex(N) = TKETracerIndex{N}()

@inline function viscous_flux_uz(i, j, k, grid, closure::TKEVD, clock, velocities, diffusivities, tracers, buoyancy)
    Ku = ℑxzᶠᵃᶠ(i, j, k, grid, diffusivities.Kᵘ)
    return - Ku * ∂zᵃᵃᶠ(i, j, k, grid, velocities.u)
end

@inline function viscous_flux_vz(i, j, k, grid, closure::TKEVD, clock, velocities, diffusivities, tracers, buoyancy)
    Kv = ℑyzᵃᶠᶠ(i, j, k, grid, diffusivities.Kᵘ)
    return - Kv * ∂zᵃᵃᶠ(i, j, k, grid, velocities.v)
end

@inline function viscous_flux_wz(i, j, k, grid, closure::TKEVD, clock, velocities, diffusivities, tracers, buoyancy)
    @inbounds Kw = diffusivities.Kᵘ[i, j, k]
    return - Kw * ∂zᵃᵃᶜ(i, j, k, grid, velocities.w)
end

@inline function diffusive_flux_z(i, j, k, grid, closure::TKEVD, c, tracer_index, clock, diffusivities, tracers, buoyancy, velocities)
    Kcᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, diffusivities.Kᶜ)
    return - Kcᶜᶜᶠ * ∂zᵃᵃᶠ(i, j, k, grid, c)
end

# Diffusive flux of TKE!
@inline function diffusive_flux_z(i, j, k, grid, closure::TKEVD, e, ::TKETracerIndex, clock, diffusivities, tracers, buoyancy, velocities)
    Keᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, diffusivities.Kᵉ)
    return - Keᶜᶜᶠ * ∂zᵃᵃᶠ(i, j, k, grid, e)
end

# "Translations" for diffusive transport by non-TKEVD closures
@inline diffusive_flux_x(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_x(i, j, k, grid, closure, e, Val(N), args...)
@inline diffusive_flux_y(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_y(i, j, k, grid, closure, e, Val(N), args...)
@inline diffusive_flux_z(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_z(i, j, k, grid, closure, e, Val(N), args...)

# Shortcuts --- TKEVD incurs no horizontal transport
@inline diffusive_flux_x(i, j, k, grid, ::TKEVD, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, ::TKEVD, args...) = zero(eltype(grid))

# Disambiguate
@inline diffusive_flux_x(i, j, k, grid, ::TKEVD, e, ::TKETracerIndex, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, ::TKEVD, e, ::TKETracerIndex, args...) = zero(eltype(grid))

#####
##### Support for VerticallyImplicitTimeDiscretization
#####

const VITD = VerticallyImplicitTimeDiscretization

@inline z_viscosity(closure::TKEVD, diffusivities, velocities, tracers, buoyancy) = diffusivities.Kᵘ

@inline function z_diffusivity(closure::TKEVD, ::Val{tracer_index}, diffusivities, velocities, tracers, buoyancy) where tracer_index
    tke_index = findfirst(name -> name === :e, keys(tracers))

    if tracer_index === tke_index
        return diffusivities.Kᵉ
    else
        return diffusivities.Kᶜ
    end
end

const VerticallyBoundedGrid{FT} = AbstractPrimaryGrid{FT, <:Any, <:Any, <:Bounded}

@inline diffusive_flux_z(i, j, k, grid::APG{FT}, ::VITD, closure::TKEVD, args...) where FT = zero(FT)
@inline viscous_flux_uz(i, j, k, grid::APG{FT}, ::VITD, closure::TKEVD, args...) where FT = zero(FT)
@inline viscous_flux_vz(i, j, k, grid::APG{FT}, ::VITD, closure::TKEVD, args...) where FT = zero(FT)

@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::TKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), # on boundaries, calculate fluxes explicitly
                  zero(FT))
end

@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::TKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), # on boundaries, calculate fluxes explicitly
                  zero(FT))
end

@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::TKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), # on boundaries, calculate fluxes explicitly
                  zero(FT))
end

@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::TKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), # on boundaries, calculate fluxes explicitly
                  zero(FT))
end


text:
# Playing Cards: Standard Card
import Base: *, ==, <, >

@enum Suit ♣=1 ♢ ♡ ♠
const clubs, diamonds, hearts, spades = ♣, ♢, ♡, ♠
const Clubs, Diamonds, Hearts, Spades = ♣, ♢, ♡, ♠

@enum Rank ace=1 two three four five six seven eight nine ten jack queen king
const Ace, Two, Three, Four = ace, two, three, four
const Five, Six, Seven, Eight = five, six, seven, eight
const Nine, Ten, Jack, Queen, King = nine, ten, jack, queen, king

const suits = ["Clubs", "Diamonds", "Hearts", "Spades"]
const ranks = ["A", "2", "3", "4", "5", "6", "7", "8", "9", "10", "J", "Q", "K"]

struct Card
    rank::Rank
    suit::Suit
    function Card(rank::Int, suit::Int)
        @assert(1 ≤ rank ≤ 13, "Invalid rank!")
        @assert(1 ≤ suit ≤ 4, "Invalid suit!")
        return new(Rank(rank), Suit(suit));
    end
end

@inline Card(rank::Rank, suit::Suit) = Card(Int(rank), Int(suit));
@inline Card(rank::Rank, suit::Int) = Card(Int(rank), suit);
@inline Card(rank::Int, suit::Suit) = Card(rank, Int(suit));
@inline rank::Int * suit::Suit = Card(rank, suit);

@inline Base.show(io::IO, card::Card) = print(io, ranks[Int(card.rank)], card.suit);
@inline Base.in(suit::Suit, card::Card) = suit == card.suit;
@inline Base.in(rank::Rank, card::Card) = rank == card.rank;
@inline Base.in(rank::Int, card::Card) = Rank(rank) == card.rank;

@inline ==(x::Card, y::Card) = x.rank == y.rank;
@inline ===(x::Card, y::Card) = (x.rank == y.rank) && (x.suit == y.suit);
@inline Base.isequal(x::Card, y::Card) = x == y;

@inline <(x::Card, y::Card) = isless(x, y);
@inline >(x::Card, y::Card) = isless(y, x);
function Base.isless(x::Card, y::Card)
    xval = (x.rank == ace) ? Int(x.rank) + 13 : Int(x.rank)
    yval = (y.rank == ace) ? Int(y.rank) + 13 : Int(y.rank)
    return xval < yval;
end

@inline suit(card::Card) = card.suit;
@inline rank(card::Card) = card.rank;
@inline capitalize(s::String, i::Int = 1) = s[1:i-1] * uppercase(s[i:i]) * s[i+1:end];
@inline name(card::Card) = string(capitalize(string(rank(card))), " of ", suits[Int(suit(card))]);

ext: jl

max_stars_repo_name: Joshua-Robison/PlayingCards

max_issues_repo_name: Joshua-Robison/PlayingCards

max_forks_repo_name: Joshua-Robison/PlayingCards


text:
# Loading config for OpenCV 3.2.0 shared libraries (.dylib, .so)
# Supports OSX (10.12.3) and Linux
# see README for details

cvlibdir = "/usr/local/lib/"
cvheaderdir = "/usr/local/include/"

version = "3.2.0"

libprefix = "libopencv_"

libNames = [
              "shape",
              "stitching",
              "objdetect",
              "superres",
              "videostab",
              "calib3d",
              "features2d",
              "highgui",
              "videoio",
              "imgcodecs",
              "video",
              "photo",
              "ml",
              "imgproc",
              #"flann",   // TODO: resolve typeid error due to rtti flag
              "viz"
            ]

function getFullLibNames()
    return map(x -> "$(libprefix)$(x).$(version)", libNames)
end


text:
export @vectorize

# Error thrown from ill-formed uses of @vec
struct VectorizeError <: Exception
    msg::String
end

julia_ivdep() = Symbol("julia.ivdep")
vectorize_enable(flag) = (Symbol("llvm.loop.vectorize.enable"), convert(Bool, flag))
vectorize_width(n) = (Symbol("llvm.loop.vectorize.width"), convert(Int, n))

function loopinfo(name, expr, nodes...)
    if expr.head != :for
        error("Syntax error: pragma $name needs a for loop")
    end
    push!(expr.args[2].args, Expr(:loopinfo, nodes...))
    return expr
end

"""
    @vectorize [ivdep] [length=512] for ...
        ...
    end

Instrument "for" loop IR with metadata controling its vectorization. Without further
options the llvm.loop.vectorize.enable(true) marker is set. The optional "ivdep" adds
metadata to mark memory accesses inside the loop as independent of each other. The
"length" option can set a different vectorization width than the default (256 on VE).
Loops using only 32 bit entities might benefit of packed vectorization by setting
`length=512`.
"""
macro vectorize(expr...)
    loop = expr[end]
    args = expr[1:end-1]
    N = 0
    ivdep = false
    for arg in args
        if arg.head == :(=) && arg.args[1] == :length
            if arg.args[2] isa Integer
                N = arg.args[2]
            else
                throw(VectorizeError("@vectorize length value must be an integer!"))
            end
        elseif arg === :ivdep
            ivdep = true
        else
            throw(VectorizeError("Only 'ivdep' and 'length=...' are "*
                "valid as the first argument to @vectorize"))
        end
    end
    nodes = (vectorize_enable(true),)
    if N > 0
        nodes = (vectorize_width(N), nodes...)
    end
    if ivdep
        nodes = (julia_ivdep(), nodes...)
    end
    return esc(loopinfo("@vectorize", loop, nodes...))
end


text:
using DataFrames
#using StatsBase
using Statistics
using Distributions
using CSV
using StatsBase
using Printf
export spearman_cor, consolidate_dataframe, read_dataframe, write_dataframe_with_comments, write_dataframe
export goal_lookup_from_df, add_edf_col_from_df, frequencies_from_counts, variance1, variance2

# Returns a dataframe by reading a CSV file with comments that start with "#"
function read_dataframe( df_filename::String )
  df = CSV.read(df_filename,DataFrame,header=true,comment="#")
  df
end

function scor_df( filename::String )
  df = read_dataframe( filename )
  rnames = vcat( [:logsteps], names(df)[12:20])  # result column names
  cdf = DataFrame()
  cdf.names = String.(rnames)
  for i = 1:length(rnames)
    cdf[!,rnames[i]] = zeros(Float64,length(rnames))
    for j = (i+1):length(rnames)
      cdf[j,rnames[i]] = StatsBase.corspearman( df[!,rnames[i]], df[!,rnames[j]] )
    end
  end
  cdf
end

function t_value( v1::Vector{Float64}, v2::Vector{Float64} )
  r = corspearman( v1, v2 )
  r*sqrt((length(v1)-2)/(1-r^2))
end

# Returns a pair of the Spearman correlation between two columns of df and a p-value for this correlation
function spearman_cor( df::DataFrame, name1::Symbol, name2::Symbol )
  v1 = Vector{Float64}(df[!,name1])
  v2 = Vector{Float64}(df[!,name2])
  r = StatsBase.corspearman( v1, v2 )
  t_value = r*sqrt((length(v1)-2)/(1-r^2))
  # p_value is always positive
  p_value = r>=0.0 ? Distributions.ccdf( TDist(length(v1)-2), t_value ) : Distributions.cdf( TDist(length(v1)-2), t_value )
  (r, p_value)
end

function spearman_cor( v1::AbstractVector, v2::AbstractVector )
  r = StatsBase.corspearman( v1, v2 )
  t_value = r*sqrt((length(v1)-2)/(1-r^2))
  # p_value is always positive
  p_value = r>=0.0 ? Distributions.ccdf( TDist(length(v1)-2), t_value ) : Distributions.cdf( TDist(length(v1)-2), t_value )
  (r, p_value)
end

# Consolidates a dataframe with the columns specified below by averaging columns with the same parameter values
# Revised 10/27/20
# Similar to R aggregate function.
function consolidate_dataframe( in_filename::String, out_filename::String )
  df = read_dataframe( in_filename )
  cdf = consolidate_dataframe( df )
  open( out_filename, "w" ) do outfile
    open( in_filename, "r" ) do infile
      line = readline(infile)
      while line[1] == '#'
        #println("line: ",line)
        #if line[1:6] != "# cor("
          write(outfile,string(line,"\n"))
        #end
        line = readline(infile)
      end
    end
   #write(outfile,"# consolidate increement: $increment \n")
   CSV.write( outfile, cdf, append=true, writeheader=true )
  end
end

# Consolidates a dataframe with the columns specified below by averaging rows with the same parameter values over all columns except goals
# Assumes that the first column is the "goal" column which may either be Vectors of MyInts or strings.
# The :evo_count column is not averaged since these are cummulative over the rows corresponding to a goal.
# Thus, for the :evo_count column, the last row corresponding to the goal is the consolidated value.
# Revised 10/18/21   Previous values for evo_count were not correct
# Similar to R aggregate function.
function consolidate_dataframe( df::DataFrame )
  new_df = DataFrame()
  #df.n_combined = map(Float64,df.n_combined)  # Only for neutral walk files, remove otherwise
  for n in names(df)
    if n!= "ntries" && n != "evo_count"
      new_df[!,n] = typeof(df[1,n])[]
    else
      new_df[!,n] = Float64[]
    end
  end
  #println("size(new_df): ",size(new_df))
  #println("names(new_df): ",names(new_df))
  new_names = vcat(["goal"],names(df)[2:end])
  #println("new_names: ",new_names)
  select!(df,new_names)
  df = sort(df,[:goal])
  ff = findfirst(x->x!=df.goal[1],df.goal)
  increment = ff != nothing ? ff - 1 : size(df)[1]
  #println("increment: ",increment)
  if Int(floor(size(df)[1]/increment)) != Int(ceil(size(df)[1]/increment))
    error("The number of rows in the dataframe must be a multiple of increment")
  end
  for r = 1:increment:Int(floor(size(df)[1]))
    ndf_row = Any[df.goal[r]]   # first entry of row is the current goal
    isum = zeros(Int64,size(df)[2])
    fsum = zeros(Float64,size(df)[2])
    for row = r:(r+increment-1)
      #println("r: ",r,"  row: ",row)
      for i = 2:size(df)[2]
        if names(df)[i] == "evo_count"
          #println("evo_count skip value: ",df[row,:evo_count])
          continue
        elseif typeof(df[row,i]) <: Int64
          isum[i] += df[row,i]
        elseif typeof(df[row,i]) <: Float64
          fsum[i] += df[row,i]
        else
          error("Illegal type of dataframe element")
        end
      end
    end
    for i = 2:size(df)[2]
      if names(df)[i] == "evo_count"
        #println("evo_count push value: ",df[r+increment-1,:evo_count])
        push!(ndf_row,df[r+increment-1,:evo_count])  # last row for goal has correct evo_count
      elseif typeof(df[r,i]) <: Int64
        #print(" ",names(df)[i]," ",isum[i])
        push!(ndf_row,isum[i]/increment)
      elseif typeof(df[r,i]) <: Float64
        #print(" ",names(df)[i]," ",fsum[i])
        push!(ndf_row,fsum[i]/increment)
      end
      #println()
    end
    #println("size(new_df): ",new_df)
    #println("length(ndf_row): ",length(ndf_row))
    #println("r:",r," ",ndf_row)
    push!(new_df,ndf_row)
  end
  new_df
end

function write_dataframe( df::DataFrame, out_filename::String )
  CSV.write( out_filename, df, append=true, writeheader=true )
end

# Writes the datarame df to the file out_filename with comments taken from file comments_filename.
# No dataframe is read from comments_filename.
function write_dataframe_with_comments( df::DataFrame, comments_filename::String, out_filename::String )
  open( out_filename, "w" ) do outfile
    open( comments_filename, "r" ) do infile
      line = readline(infile)
      while line[1] == '#'
        #println("line: ",line)
        if line[1:6] != "# cor("
          write(outfile,string(line,"\n"))
        end
        line = readline(infile)
      end
    end
   #write(outfile,"# increement: $increment \n")
   CSV.write( outfile, df, append=true, writeheader=true )
  end
end

function write_comments( in_filename::String, out_filename::String )
  open( out_filename, "w" ) do outfile
    open( in_filename, "r" ) do infile
      line = readline(infile)
      while line[1] == '#'
        println("line: ",line)
        write(outfile,string(line,"\n"))
        line = readline(infile)
      end
    end
  end
end

# Given a goal represented as a string, extract the df values of the columns fields for that goal.
# Example:  goal_lookup_from_df(df,"0x4208",[:counts10ints]) returns 1192.
function goal_lookup_from_df( df::DataFrame, goal::String, fields::Vector{Symbol} )
  [df[findfirst(x->x==goal,df.goal),f] for f in fields]
end

# Adds a column to dataframe edf whose values are obtained by goal_lookup_from_df( df, goal_string, [field] )
function add_edf_col_from_df( df::DataFrame, edf::DataFrame, field::Symbol )
  edf[!,field] =  zeros(Float64,size(edf)[1])
  for i = 1:size(edf)[1]
    # Extract MyInt goal from edf.goal[i] by parsing.  See Julia documentation for Metaprogramming.
    goal = Meta.parse(edf[i,:goal]).args[2]
    goal_string = @sprintf("0x%x",goal)   # Convert to string
    edf[i,field] = goal_lookup_from_df(df,goal_string,[field])[1]
  end
  edf
end

# Return the frequencies of the unique elements of V
function frequencies_from_counts( V::Vector{Int64} )
  sort!(V)
  #println("V: ",V)
  freq = Int64[]
  i = 1
  while i <= length(V)
    push!(freq,1)
    while i < length(V) && V[i+1] == V[i]
      freq[end] += 1
      i += 1
    end
    i += 1
  end
  #println("sum(freq): ",sum(freq),"  freq: ",freq)
  # Two assertions to check correctness
  @assert sum(freq) == length(V)
  @assert length(freq) == length(unique(V))
  freq
end

function pad_missing( V::Vector, new_length::Int64 )
  result = Union{Missing,Int64}[]
    for v in V push!(result,v)
  end
  for i = length(V)+1:new_length
    push!(result,missing)
  end
  result
end

# Simple sample variance
function variance1( X::Vector )
  N = length(X)
  avg = sum(X)/N
  sum( (x - avg)^2 for x in X )/(N-1)
end

# Simple sample variance
function variance2( X::Vector )
  N = length(X)
  ( sum( x^2 for x in X) - sum(X)^2/N )/(N-1)
end


text:
mutable struct Figure{B<:Backend}
    # NOTE: couldn't use @kw_args here, clashed with the format of the base constructor.
    g ::B      # description buffer
    id::String # id of the figure
    # ---
    axes     ::Vector{Axes{B}}   # subplots
    size     ::T2F               # (width, heigth)
    textstyle::TextStyle         # parent font etc
    bgcolor  ::Option{Colorant}  # background col, nothing=transparent
    # ---
    texlabels   ::Option{Bool}   # true if has tex
    texscale    ::Option{String} # scale latex * hei (def=1)
    texpreamble ::Option{String} # latex preamble
    transparency::Option{Bool}   # if true, use cairo device
    # ---
    subroutines::Dict{String,String}
end
Figure(g::B, id::String) where B<:Backend =
    Figure(g, id, Vector{Axes{B}}(), (12., 9.), TextStyle(font="texcmss", hei=0.35),
           c"white", ∅, ∅, ∅, ∅, Dict{String,String}())

"""
    reset!(f)

Internal function to completely refresh the figure `f` only keeping its id and size.
"""
function reset!(f::Figure{B}) where B
    take!(f.g) # empty the buffer
    f.axes         = Vector{Axes{B}}() # clean axes
    f.textstyle    = TextStyle(font="texcmss", hei=0.35) # default fontstyle
    f.bgcolor      = c"white" # default bg color
    f.texlabels    = ∅
    f.texscale     = ∅
    f.texpreamble  = ∅
    f.transparency = ∅
    f.subroutines  = Dict{String,String}()
    GP_ENV["CURFIG"]  = f
    GP_ENV["CURAXES"] = nothing
    return f
end

function Base.show(io::IO, ::MIME"text/plain", f::Figure{GLE})
    tbg = (f.bgcolor === nothing)
    wbg = (f.bgcolor == RGB(1,1,1))
    s = "GPlot.Figure{GLE}" *
        "\n\t"*rpad("Name:", 15) * (f.id == "_fig_" ? "default (\"_fig_\")" : "\"$(f.id)\"") *
        "\n\t"*rpad("Size:", 15) * "$(round.(f.size, digits=1))" *
        "\n\t"*rpad("Bg. color:", 15) * (tbg ? "none" : ifelse(wbg, "white", col2str(f.bgcolor))) *
        "\n\t"*rpad("N. axes:", 15) * "$(length(f.axes))" *
        "\n\t"*rpad("LaTeX:", 15) * (isdef(f.texlabels) ? "$(f.texlabels)" : "false") *
        "\n\t"*rpad("Transparent:", 15) * (isdef(f.transparency) ? "$(f.transparency)" : "false")
    write(io, s)
    return nothing
end


text:
export TwoStageCost, two_stage_method

struct TwoStageCost{F,F2,D} <: Function
  cost_function::F
  cost_function2::F2
  estimated_solution::D
  estimated_derivative::D
end

(f::TwoStageCost)(p) = f.cost_function(p)
(f::TwoStageCost)(p,g) = f.cost_function2(p,g)

function decide_kernel(kernel)
    if kernel == :Epanechnikov
        return Epanechnikov_kernel
    elseif kernel == :Uniform
        return Uniform_kernel
    elseif kernel == :Triangular
        return Triangular_kernel
    elseif kernel == :Quartic
      return Quartic_Kernel
    elseif kernel == :Triweight
      return Triweight_Kernel
    elseif kernel == :Tricube
      return Tricube_Kernel
    elseif kernel == :Gaussian
      return Gaussian_Kernel
    elseif kernel == :Cosine
      return Cosine_Kernel
    elseif kernel == :Logistic
      return Logistic_Kernel
    elseif kernel == :Sigmoid
      return Sigmoid_Kernel
    else
      return Silverman_Kernel
    end
end

function construct_t1(t,tpoints)
    T1 = []
    for i in 1:length(tpoints)
        push!(T1,[1 tpoints[i]-t])
    end
    foldl(vcat,T1)
end
function construct_t2(t,tpoints)
    T2 = []
    for i in 1:length(tpoints)
        push!(T2,[1 tpoints[i]-t (tpoints[i]-t)^2])
    end
    foldl(vcat,T2)
end
function construct_w(t,tpoints,h,kernel_function)
    n = length(tpoints)
    W = zeros(n)
    for i in 1:n
        W[i] = kernel_function((tpoints[i]-t)/h)/h
    end
    Matrix(Diagonal(W))
end
function construct_estimated_solution_and_derivative!(estimated_solution,estimated_derivative,e1,e2,data,kernel_function,tpoints,h,n)
  for i in 1:n
      T1 = construct_t1(tpoints[i],tpoints)
      T2 = construct_t2(tpoints[i],tpoints)
      W = construct_w(tpoints[i],tpoints,h,kernel_function)
      estimated_solution[:,i] = e1'*inv(T1'*W*T1)*T1'*W*data'
      estimated_derivative[:,i] = e2'*inv(T2'*W*T2)T2'*W*data'
  end
end

function two_stage_method(prob::DiffEqBase.DEProblem,tpoints,data;kernel= :Epanechnikov,
                          loss_func = L2Loss,mpg_autodiff = false,
                          verbose = false,verbose_steps = 100,
                          autodiff_prototype = mpg_autodiff ? zeros(length(prob.p)) : nothing,
                          autodiff_chunk = mpg_autodiff ? ForwardDiff.Chunk(autodiff_prototype) : nothing)
    f = prob.f
    n = length(tpoints)
    h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16))
    estimated_solution = zeros(size(data)[1],n)
    estimated_derivative = zeros(size(data)[1],n)
    kernel_function = decide_kernel(kernel)
    e1 = [1;0]
    e2 = [0;1;0]
    construct_estimated_solution_and_derivative!(estimated_solution,estimated_derivative,e1,e2,data,kernel_function,tpoints,h,n)
    # Step - 2
    cost_function = function (p)
        du = similar(prob.u0, promote_type(eltype(prob.u0), eltype(p)))
        sol = Vector{typeof(du)}(undef,n)
        f = prob.f
        for i in 1:n
          est_sol = @view estimated_solution[:,i]
          f(du,est_sol,p,tpoints[i])
          sol[i] = copy(du)
        end
        sqrt(sum(abs2,vec(estimated_derivative)[i] - vec(VectorOfArray(sol))[i] for i in 1:length(vec(estimated_derivative))))
    end

    if mpg_autodiff
      gcfg = ForwardDiff.GradientConfig(cost_function, autodiff_prototype, autodiff_chunk)
      g! = (x, out) -> ForwardDiff.gradient!(out, cost_function, x, gcfg)
    else
      g! = (x, out) -> Calculus.finite_difference!(cost_function,x,out,:central)
    end
    if verbose
      count = 0 # keep track of # function evaluations
    end
    cost_function2 = function (p,grad)
      if length(grad)>0
        g!(p,grad)
      end
      loss_val = cost_function(p)
      if verbose
        count::Int += 1
        if mod(count,verbose_steps) == 0
          println("Iteration: $count")
          println("Current Cost: $loss_val")
          println("Parameters: $p")
        end
      end
      loss_val
    end

    return TwoStageCost(cost_function, cost_function2, estimated_solution, estimated_derivative)
end


text:
# to do: fix this markdown.
"""
    setupcubicitp(ϕ::Array{T,D}, x_ranges::Vector{LinRange{T,L}}, amplification_factor::T)

Separable cubic spline interpolation on compact interval, with gain parameter.

Inputs: `x_ranges` is a vector of LinRange that specifies the sampling coordinates.
        `Φ is the set of sampling values associated with x_ranges.
        `amplification_factor` is the gain multiplier applied to the interpolated function.

Returns:    `f` is the interpolated function.
            `df` is the gradient of `f`.
            `d2f` is the Hessian of `f`. It throws an error if evaluated at any x in the compact interval that is specified by `x_ranges`.
"""
function setupcubicitp(ϕ::Array{T,D}, x_ranges::Vector{LinRange{T,L}}, amplification_factor::T) where {T,D,L}

    @assert D == length(x_ranges)
    N_array = collect( length(x_ranges[d]) for d = 1:D )

    itp_ϕ = Interpolations.interpolate(ϕ,
                Interpolations.BSpline(Interpolations.Cubic(
                    Interpolations.Flat(Interpolations.OnGrid()))))

    etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())
    #etp_ϕ = Interpolations.extrapolate(itp_ϕ, 0)

    st = collect( x_ranges[d][1] for d = 1:D )
    fin = collect( x_ranges[d][end] for d = 1:D )
    f = xx->etp_ϕ(interval2itpindex(xx,
                st,
                fin,
                N_array)...)*amplification_factor

    # chain rule for first derivatives.
    df = xx->( Interpolations.gradient(etp_ϕ,interval2itpindex(xx,
                st,
                fin,
                N_array)...) .* derivativeinterval2itpindex(st,fin,N_array) .*amplification_factor )

    # chain rule for second derivatives.
    d2f = xx->( Interpolations.hessian(etp_ϕ,interval2itpindex(xx,
                st,
                fin,
                N_array)...) .* (derivativeinterval2itpindex(st,fin,N_array).^2) .*amplification_factor )

    return f, df, d2f
end

function interval2itpindex(x::T, a::T, b::T, N::Int)::T where T <: Real
    return (x-a)/(b-a)*(N-1) + 1
end

# converts an x ∈ [a[d],b[d]]^D, x ∈ ℝ^D, to an index i ∈ [N]^D.
function interval2itpindex(x::Vector{T}, a::Vector{T}, b::Vector{T}, N::Vector{Int})::Vector{T} where T <: Real
    return collect( interval2itpindex(x[d], a[d], b[d], N[d]) for d = 1:length(N) )
end

# do for each entry in X.
function interval2itpindex(X::Vector{Vector{T}}, a::Vector{T}, b::Vector{T}, N::Vector{Int})::Vector{Vector{T}} where T <: Real
    return collect( collect( interval2itpindex(X[n][d], a[d], b[d], N[d]) for d = 1:length(N) ) for n = 1:length(X) )
end

# for HCubature.
function interval2itpindex(x::AbstractArray{T,D}, a::Vector{T}, b::Vector{T}, N::Vector{Int})::Vector{T} where {T,D}
    return collect( interval2itpindex(x[d], a[d], b[d], N[d]) for d = 1:length(N) )
end
# function interval2itpindex(x::StaticArrays.SArray{Tuple{D},T,1,D}, a::Vector{T}, b::Vector{T}, N::Vector{Int})::Vector{T} where {T,D}
#     return collect( interval2itpindex(x[d], a[d], b[d], N[d]) for d = 1:length(N) )
# end

### for AD.
function interval2itpindex(x::T, a, b, N::Int) where T <: Real
    return (x-a)/(b-a)*(N-1) + 1
end

# converts an x ∈ [a[d],b[d]]^D, x ∈ ℝ^D, to an index i ∈ [N]^D.
function interval2itpindex(x::Vector{T}, a, b, N::Vector{Int}) where T <: Real
    return collect( interval2itpindex(x[d], a[d], b[d], N[d]) for d = 1:length(N) )
end

# do for each entry in X.
function interval2itpindex(X::Vector{Vector{T}}, a, b, N::Vector{Int}) where T <: Real
    return collect( collect( interval2itpindex(X[n][d], a[d], b[d], N[d]) for d = 1:length(N) ) for n = 1:length(X) )
end

function derivativeinterval2itpindex(a::T, b::T, N::Int)::T where T <: Real
    return (N-1)/(b-a)
end

function derivativeinterval2itpindex(a::Vector{T}, b::Vector{T}, N::Vector{Int})::Vector{T} where T <: Real
    return collect( derivativeinterval2itpindex(a[d], b[d], N[d]) for d = 1:length(N) )
end

# # WIP: what is with the hessian?
# function setupcubicitpmutator(ϕ::Array{T,D}, x_ranges::Vector{LinRange{T,L}}, amplification_factor::T) where {T,D,L}
#
#     @assert D == length(x_ranges)
#     N_array = collect( length(x_ranges[d]) for d = 1:D )
#
#     itp_ϕ = Interpolations.interpolate(ϕ,
#                 Interpolations.BSpline(Interpolations.Cubic(
#                     Interpolations.Flat(Interpolations.OnGrid()))))
#
#     etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())
#     #etp_ϕ = Interpolations.extrapolate(itp_ϕ, 0)
#
#     st = collect( x_ranges[d][1] for d = 1:D )
#     fin = collect( x_ranges[d][end] for d = 1:D )
#     f = xx->etp_ϕ(interval2itpindex(xx,
#                 st,
#                 fin,
#                 N_array)...)*amplification_factor
#
#     # chain rule for first derivatives.
#     df = xx->( Interpolations.gradient(etp_ϕ,interval2itpindex(xx,
#                 st,
#                 fin,
#                 N_array)...) .* derivativeinterval2itpindex(st,fin,N_array) .*amplification_factor )
#
#     # chain rule for second derivatives.
#     d2f = xx->( Interpolations.hessian(etp_ϕ,interval2itpindex(xx,
#                 st,
#                 fin,
#                 N_array)...) .* (derivativeinterval2itpindex(st,fin,N_array).^2) .*amplification_factor )
#
#
#     return f, df, d2f
# end
#
#
# """
# Chain rule for first derivatives. Implements:
#
# df = xx->( Interpolations.gradient(etp_ϕ,interval2itpindex(xx,
#             st,
#             fin,
#             N_array)...) .* derivativeinterval2itpindex(st,fin,N_array) .*amplification_factor )
# """
# function firstderivativechainrulewrapper!( g::Vector{T},
#                                             x,
#                                             N_array::Vector{Int},
#                                             etp_ϕ,
#                                             a::T)::Nothing where T
#     #
#     @assert length(x) == length(g)
#
#     Interpolations.gradient!(g, etp_ϕ, interval2itpindex(  x,
#                                                         st,
#                                                         fin,
#                                                         N_array)...)
#
#     for d = 1:D
#         g[d] = g[d] * derivativeinterval2itpindex(a[d], b[d], N[d]) * a
#     end
#
#     return nothing
# end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
using Sundials

## f routine. Compute function f(t,y).

function f(t, y_nv, ydot_nv, user_data)
    y = convert(Vector, y_nv)
    ydot = convert(Vector, ydot_nv)
    ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
    ydot[3] = 3.0e7*y[2]*y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return Sundials.CV_SUCCESS
end

## g routine. Compute functions g_i(t,y) for i = 0,1.

function g(t, y_nv, gout_ptr, user_data)
    y = convert(Vector, y_nv)
    gout = Sundials.asarray(gout_ptr, (2,))
    gout[1] = y[1] - 0.0001
    gout[2] = y[3] - 0.01
    return Sundials.CV_SUCCESS
end

## Jacobian routine. Compute J(t,y) = df/dy.
# broken -- needs a wrapper from Sundials._DlsMat to Matrix and Jac user function wrapper
function Jac(N, t, ny, fy, Jptr, user_data, tmp1, tmp2, tmp3)
    y = convert(Vector, ny)
    dlsmat = unpack(IOString(unsafe_wrap(convert(Ptr{UInt8}, Jptr),
                                                 (sum(map(sizeof, Sundials._DlsMat))+10,), false)),
                    Sundials._DlsMat)
    J = unsafe_wrap(unsafe_ref(dlsmat.cols), (Int(neq), Int(neq)), false)
    J[1,1] = -0.04
    J[1,2] = 1.0e4*y[3]
    J[1,3] = 1.0e4*y[2]
    J[2,1] = 0.04
    J[2,2] = -1.0e4*y[3] - 6.0e7*y[2]
    J[2,3] = -1.0e4*y[2]
    J[3,2] = 6.0e7*y[2]
    return Sundials.CV_SUCCESS
end

neq = 3

t0 = 0.0
t1 = 0.4
tmult = 10.0
nout = 12
y0 = [1.0, 0.0, 0.0]
reltol = 1e-4
abstol = [1e-8, 1e-14, 1e-6]
userdata = nothing
mem_ptr = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
cvode_mem = Sundials.Handle(mem_ptr)
userfun = Sundials.UserFunctionAndData(f, userdata)
Sundials.CVodeSetUserData(cvode_mem, userfun)

Sundials.CVodeInit(cvode_mem, Sundials.old_cfunction(Sundials.cvodefun,
                    Cint, Tuple{Sundials.realtype, Sundials.N_Vector,
                    Sundials.N_Vector, Ref{typeof(userfun)}}), t1,
                    convert(Sundials.N_Vector, y0))
Sundials.@checkflag Sundials.CVodeInit(cvode_mem, f, t0, y0)
Sundials.@checkflag Sundials.CVodeSVtolerances(cvode_mem, reltol, abstol)
Sundials.@checkflag Sundials.CVodeRootInit(cvode_mem, 2, g)
A = Sundials.SUNDenseMatrix(neq,neq)
mat_handle = Sundials.MatrixHandle(A,Sundials.DenseMatrix())
LS = Sundials.SUNDenseLinearSolver(convert(Sundials.N_Vector,y0),A)
LS_handle = Sundials.LinSolHandle(LS,Sundials.Dense())
Sundials.@checkflag Sundials.CVDlsSetLinearSolver(cvode_mem, LS, A)
#Sundials.@checkflag Sundials.CVDlsSetDenseJacFn(cvode_mem, Jac)

iout = 0
tout = t1
t = [t0]

while iout < nout
    y = similar(y0)
    flag = Sundials.CVode(cvode_mem, tout, y, t, Sundials.CV_NORMAL)
    println("T=", tout, ", Y=", y)
    if flag == Sundials.CV_ROOT_RETURN
        rootsfound = zeros(Cint, 2)
        Sundials.@checkflag Sundials.CVodeGetRootInfo(cvode_mem, rootsfound)
        println("roots=", rootsfound)
    elseif flag == Sundials.CV_SUCCESS
      global iout += 1
      global tout *= tmult
    end
end

empty!(cvode_mem)
empty!(mat_handle)
empty!(LS_handle)


text:
"""
    read_sn_data(folder::String)

Reads a network folder and builds the mn_data dictionary.

# Arguments
- `folder::String`: the path to the folder containing all the network data
"""
function read_sn_data(folder::String; kwargs...)
   directory_path = "$folder/"
   subnet_file = "$(directory_path)/subnetworks.csv"

   subnets = CSV.read(subnet_file, DataFrame)

   if nrow(subnets) == 1
      interfaces = DataFrame(index=Int64[],file=String[],variable_f=Bool[],f_base=Float64[],f_min=Float64[],fmax=Float64[])
   else
      interface_file = "$(directory_path)/interfaces.csv"
      interfaces = CSV.read(interface_file, DataFrame)
   end

   networks = Dict{String,Any}()

   for subnet in eachrow(subnets)
      file_name = directory_path*subnet.file
      networks[subnet.file] = PowerModels.parse_file(file_name)

      # Add zeros to turn linear objective functions into quadratic ones
      # so that additional parameter checks are not required
      PowerModels.standardize_cost_terms!(networks[subnet.file], order=2)

   end

   mn_data = make_mn_data(subnets, interfaces, networks; kwargs...)

   return mn_data
end

"""
    function make_mn_data(
        subnetworks,
        interfaces,
        networks::Dict
    )

Builds the mn_data dictionary from the specifications of the subnetworks and interfaces DataFrames and the network data in the networks Dict.

# Arguments
- `subnetworks`: a DataFrame in the format of *subnetworks.csv*, and *interfaces.csv*
- `interfaces`: a DataFrame in the format of *interfaces.csv*
- `networks::Dict`: a Dict of all subnetworks, as PowerModels networks
"""
function make_mn_data(
   subnetworks,
   interfaces,
   networks::Dict;
   no_converter_loss::Bool=false
   )

   mn_data = Dict{String,Any}(
      "sn" => Dict{String,Any}()
   )

   for subnet in eachrow(subnetworks)
      file_name = subnet.file
      # println("reading $(file_name)")

      mn_data["sn"]["$(subnet.index)"] = networks[file_name]
      mn_data["sn"]["$(subnet.index)"]["variable_f"] = typeof(subnet.variable_f)==Bool ? subnet.variable_f : (lowercase(subnet.variable_f) == "true")
      mn_data["sn"]["$(subnet.index)"]["f_base"] = subnet.f_base
      if mn_data["sn"]["$(subnet.index)"]["variable_f"]
         mn_data["sn"]["$(subnet.index)"]["f_min"] = subnet.f_min
         mn_data["sn"]["$(subnet.index)"]["f_max"] = subnet.f_max
      end

      # REMOVING BUSES OF TYPE 4 (ISOLATED)
      # mn_data["sn"]["$(subnet.index)"]["bus"] = filter(p->(last(p)["bus_type"]!=4), mn_data["sn"]["$(subnet.index)"]["bus"])
      # Disable loads at type 4 buses
      for (idx,load) in mn_data["sn"]["$(subnet.index)"]["load"]
         # println("Bus $(load["load_bus"])")
         # println(mn_data["sn"]["$(subnet.index)"]["bus"]["$(load["load_bus"])"])
         if mn_data["sn"]["$(subnet.index)"]["bus"]["$(load["load_bus"])"]["bus_type"]==4 && load["status"] == 1
            println("Bus $(load["load_bus"]) is type 4 (isolated) and has a load. The load will be disabled.")
            load["status"] = 0
         end
      end
      for (idx,gen) in mn_data["sn"]["$(subnet.index)"]["gen"]
         if mn_data["sn"]["$(subnet.index)"]["bus"]["$(gen["gen_bus"])"]["bus_type"]==4 && gen["gen_status"] == 1
            println("Bus $(gen["gen_bus"]) is type 4 (isolated) and has a generator. The generator will be disabled.")
            gen["gen_status"] = 0
         end
      end
      for (idx,shunt) in mn_data["sn"]["$(subnet.index)"]["shunt"]
         if mn_data["sn"]["$(subnet.index)"]["bus"]["$(shunt["shunt_bus"])"]["bus_type"]==4 && shunt["status"] == 1
            println("Bus $(shunt["shunt_bus"]) is type 4 (isolated) and has a shunt. The shunt will be disabled.")
            shunt["status"] = 0
         end
      end
   end

   mn_data["converter"] = Dict{Int64, Any}()
   interface_bus = Array{Tuple{Int64,Int64},1}()
   for interface in eachrow(interfaces)
      # Check that we have transformer/filter branch data in the interfaces file and for this terminal specifically
      if !no_converter_loss && all(indexin(["R","X","G","B"], names(interfaces)) .!= nothing) && !(any(x -> ismissing(x), interface[["R","X","G","B"]])) && !(interface.R == 0 && interface.X == 0)
         # First create the converter terminal buses. New bus index is '90X',
         # where X is the existing bus, e.g. bus 334 gets terminal bus 90334.
         # If 90X is already a bus in the network, try 900X, and so on.
         busmult=0
         new_bus_idx = Int(interface.bus + 900*10^(busmult+floor(log10(interface.bus))))
         while new_bus_idx in parse.(Int,keys(mn_data["sn"]["$(interface.subnet_index)"]["bus"]))
            busmult=busmult+1
            new_bus_idx = Int(interface.bus + 900*10^(busmult+floor(log10(interface.bus))))
         end
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"] = Dict{String,Any}()
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["bus_i"] = new_bus_idx
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["index"] = new_bus_idx
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["bus_type"] = 1
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["vmax"] = interface.vmax
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["vmin"] = 0.8
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["zone"] = copy(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"]["zone"])
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["area"] = copy(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"]["area"])
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["va"] = copy(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"]["va"])
         mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["vm"] = copy(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"]["vm"])
         if "base_kv" in names(interfaces)
            mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["base_kv"] = interface.base_kv
         elseif "base_kv" in keys(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"])
            mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["base_kv"] = copy(mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(interface.bus)"]["base_kv"])
         end
         # Include current limit in bus data
         # mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["converter_imax"] = Dict{Int64,Float64}()
         # mn_data["sn"]["$(interface.subnet_index)"]["bus"]["$(new_bus_idx)"]["converter_imax"][interface.index] = interface.imax

         # Now define a branch at the terminal
         new_branch_idx = maximum(parse.(Int,keys(mn_data["sn"]["$(interface.subnet_index)"]["branch"]))) + 1
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"] = Dict{String,Any}()
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["index"] = new_branch_idx
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["f_bus"] = new_bus_idx
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["t_bus"] = interface.bus
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["angmax"] = pi/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["angmin"] = -pi/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["br_r"] = interface.R
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["br_x"] = interface.X
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["g_to"] = interface.G/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["g_fr"] = interface.G/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["b_to"] = interface.B/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["b_fr"] = interface.B/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["g_to"] = interface.G/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["g_fr"] = interface.G/2
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["rate_a"] = interface.smax
         mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["br_status"] = 1
         if ("transformer" in names(interfaces)) && !(ismissing(interface.transformer))
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["transformer"] = interface.transformer
         else
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["transformer"] = false
         end
         if ("tap" in names(interfaces)) && !(ismissing(interface.tap))
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["tap"] = interface.tap
         else
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["tap"] = 1
         end
         if ("shift" in names(interfaces)) && !(ismissing(interface.shift))
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["shift"] = interface.shift
         else
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["shift"] = 0
         end
         if ("rate_b" in names(interfaces)) && !(ismissing(interface.rate_b))
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["rate_b"] = interface.rate_b
         else
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["rate_b"] = interface.smax
         end
         if ("rate_c" in names(interfaces)) && !(ismissing(interface.rate_c))
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["rate_c"] = interface.rate_c
         else
            mn_data["sn"]["$(interface.subnet_index)"]["branch"]["$(new_branch_idx)"]["rate_c"] = interface.smax
         end
         interface_bus = (interface.subnet_index,new_bus_idx)
         # if !((interface.index,interface.subnet_index,new_bus_idx) in keys(mn_data["converter"]))
         #    mn_data["converter"][interface.index,interface.subnet_index,new_bus_idx] = Dict{String,Any}("converter_buses" => Array{Tuple{Int64,Int64},1}())
         # end
         # if !((interface.subnet_index,new_bus_idx) in mn_data["converter"][interface.index,interface.subnet_index,new_bus_idx]["converter_buses"])
         #    push!(mn_data["converter"][interface.index,interface.subnet_index,new_bus_idx]["converter_buses"], (interface.subnet_index,new_bus_idx))
         # end
      else
         interface_bus = (interface.subnet_index,interface.bus)
      end

      # println("interface_bus: $(interface_bus)")
      # Create array to store indices of the converters connected to this bus,
      # along with dictionaries for all the parameters, indexed by the converter index
      if !("converter_index" in keys(mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]))
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_index"] = Array{Tuple,1}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_imax"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_vmax"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c1"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c2"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c3"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw1"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw2"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw3"] = Dict{Int64,Float64}()
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_M"] = Dict{Int64,Float64}()
      end
      param_labels = [:imax,:vmax,:c1,:c2,:c3,:sw1,:sw2,:sw3,:M]
      converter_params = Dict{Symbol,Any}()
      for label in param_labels
         # println("Interface labels: $(names(interface))")
         # println("Label: $label")
         if string(label) in names(interface)
            converter_params[label] = interface[label]
         else
            println("The converter parameter $label is not specified for interface $(interface.index).")
            if label == :vmax
               if "v_max" in names(interface)
                  println("Found v_max. Applying this limit. Change the label to \"vmax\" in interfaces.csv to ensure future functionality.")
                  converter_params[label] = interface[!,:v_max]
               else
                  converter_params[label] = 1.2
                  println("Set to 1.2.")
               end
            elseif label == :imax
               converter_params[label] = Inf
               println("Set limit to Inf.")
            elseif label == :M
               converter_params[label] = 1
               println("Set to 1.")
            else
               converter_params[label] = 0
               println("Set to 0.")
            end
         end
      end

      if !(interface.index in mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_index"])
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_imax"][interface.index] = converter_params[:imax]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_vmax"][interface.index] = converter_params[:vmax]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c1"][interface.index]   = no_converter_loss ? 0.0 : converter_params[:c1]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c2"][interface.index]   = no_converter_loss ? 0.0 : converter_params[:c2]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_c3"][interface.index]   = no_converter_loss ? 0.0 : converter_params[:c3]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw1"][interface.index]  = no_converter_loss ? 0.0 : converter_params[:sw1]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw2"][interface.index]  = no_converter_loss ? 0.0 : converter_params[:sw2]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_sw3"][interface.index]  = no_converter_loss ? 0.0 : converter_params[:sw3]
         mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_M"][interface.index]    = converter_params[:M]
         if (interface.index,interface_bus[1],interface_bus[2]) ∉ mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_index"]
            push!(mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_index"], (interface.index,interface_bus[1],interface_bus[2]))
         end
      end
      # println("converter_imax: $(mn_data["sn"]["$(interface_bus[1])"]["bus"]["$(interface_bus[2])"]["converter_imax"])")

      if !(interface.index in keys(mn_data["converter"]))
         mn_data["converter"][interface.index] = Dict{String,Any}("converter_buses" => Array{Tuple{Int64,Int64},1}())
      end
      if !(interface_bus in mn_data["converter"][interface.index]["converter_buses"])
         push!(mn_data["converter"][interface.index]["converter_buses"], interface_bus)
      end

      if !("imax" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["imax"] = Dict(interface_bus => converter_params[:imax])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["imax"]))
         mn_data["converter"][interface.index]["imax"][interface_bus] = converter_params[:imax]
      end

      if !("vmax" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["vmax"] = Dict(interface_bus => converter_params[:vmax])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["vmax"]))
         mn_data["converter"][interface.index]["vmax"][interface_bus] = converter_params[:vmax]
      end

      if !("c1" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["c1"] = Dict(interface_bus => converter_params[:c1])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["c1"]))
         mn_data["converter"][interface.index]["c1"][interface_bus] = converter_params[:c1]
      end

      if !("c2" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["c2"] = Dict(interface_bus => converter_params[:c2])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["c2"]))
         mn_data["converter"][interface.index]["c2"][interface_bus] = converter_params[:c2]
      end

      if !("c3" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["c3"] = Dict(interface_bus => converter_params[:c3])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["c3"]))
         mn_data["converter"][interface.index]["c3"][interface_bus] = converter_params[:c3]
      end

      if !("sw1" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["sw1"] = Dict(interface_bus => converter_params[:sw1])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["sw1"]))
         mn_data["converter"][interface.index]["sw1"][interface_bus] = converter_params[:sw1]
      end

      if !("sw2" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["sw2"] = Dict(interface_bus => converter_params[:sw2])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["sw2"]))
         mn_data["converter"][interface.index]["sw2"][interface_bus] = converter_params[:sw2]
      end

      if !("sw3" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["sw3"] = Dict(interface_bus => converter_params[:sw3])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["sw3"]))
         mn_data["converter"][interface.index]["sw3"][interface_bus] = converter_params[:sw3]
      end

      if !("M" in keys(mn_data["converter"][interface.index]))
         mn_data["converter"][interface.index]["M"] = Dict(interface_bus => converter_params[:M])
      elseif !(interface_bus in keys(mn_data["converter"][interface.index]["M"]))
         mn_data["converter"][interface.index]["M"][interface_bus] = converter_params[:M]
      end

         # if "converter_buses" not in keys(mn_data["converter"][interface.index])
         #    mn_data["converter"][interface.index]["converter_buses"] = Array{Tuple(Int64,Int64)}()
         # end
      # print("parsed $(directory_path)/$(file_name)\n")
   end

   return mn_data
end

function to_powermodels_json(input_network, output_dir=nothing)
   filename = split(input_network, "/")[end]
   println("$filename")
   filetype = split(lowercase(filename), '.')[end]
   base_name = filename[1:(end-length(filetype)-1)]
   (base_name, filetype) = splitext(filename)
   filetype = Unicode.normalize(filetype, casefold=true)

   if output_dir == nothing
      output_dir = join(split(input_network, "/")[1:end-1],"/")
   end

   data = PowerModels.parse_file(input_network)

   stringnet = JSON.json(data)
   open("$output_dir/$(base_name).json", "w") do f
      write(f, stringnet)
   end
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demodercrv()

Demonstrates the construction of a general curve and determine of the
derivative.

# Examples:
```julia
julia> demodercrv()
```
"""
function demodercrv()
# make and draw nurbs test curve
crv = nrbtestcrv();
nrbplot(crv,48,
        title = "First derivatives along a test curve.",
        lw =1.0,
        c = :green,
        framestyle=:box,
        legend = false);

npts = 9;
tt = collect(range(0.0,1.0,length=npts));

# first derivative
dcrv = nrbderiv(crv);
p1, dp = nrbdeval(crv,dcrv,tt);

p2 = vecnorm_toolbox(dp);

plot!(p1[1,:],p1[2,:],
      m = (3, :transparent, stroke(1, :red)),
      st = :scatter,
      legend = false);
plot!(p1[1,:],p1[2,:],
      st=:quiver,
      quiver=(p2[1,:],p2[2,:]),
      c=:red);
end # demodercrv


text:
struct SNA_BINDICES
    j1 :: Int
    j2 :: Int
    j  :: Int
end
struct SNA_ZINDICES
    j1 :: Int
    j2 :: Int
    j  :: Int
    ma1min :: Int
    ma2max :: Int
    mb1min :: Int
    mb2max :: Int
    na     :: Int
    nb     :: Int
    jju    :: Int
end
struct PrebuiltArrays{T}
    # Array indices, sizes
    idxcg_block :: Array{Int}
    idxu_block  :: Array{Int}
    idxb_block  :: Array{Int}
    idxz_block  :: Array{Int}
    idxcg_max :: Int
    idxu_max  :: Int
    idxb_max  :: Int
    idxz_max  :: Int

    idxb      :: Vector{SNA_BINDICES}
    idxz      :: Vector{SNA_ZINDICES}
    ncoeff    :: Int
    num_elements :: Int

    # Clebsch Gordan
    rootpqarray :: Array{T}
    cglist      :: Vector{T}

    # bzero
    wself       :: T
    bzero       :: Vector{T}
end

# /* ---------------------------------------------------------------------- */
## Initialize

function initialize_prebuilt_arrays(twojmax::Int, n_elements::Int, chem_flag::Bool, bzero_flag::Bool, bnorm_flag::Bool)
    jdim = twojmax+1
    ncount = 0

    wself = 1.0
    bzero = zeros(jdim)
    if (bzero_flag)
        www = wself*wself*wself;
        for j = 0:twojmax
            if (bnorm_flag)
                bzero[j+1] = www;
            else
                bzero[j+1] = www*(j+1);
            end
        end
    end

    if (chem_flag)
        num_elements = n_elements
    else
        num_elements = 1
    end

    idxcg_block = zeros(Int, jdim, jdim, jdim)
    idxu_block = zeros(Int, jdim)
    idxb_block = zeros(Int, jdim, jdim, jdim)
    idxz_block = zeros(Int, jdim, jdim, jdim)

    idxcg_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1-j2
            while j <= minimum([twojmax, j1+j2])
                idxcg_block[j1+1,j2+1,j+1] = idxcg_count;
                for m1 = 0:j1
                    for m2 = 0:j2
                        idxcg_count+=1;
                    end
                end
                j+=2;
            end
        end
    end
    idxcg_max = idxcg_count
    #   // index list for uarray
    #   // need to include both halves

    idxu_count = 0;

    for j = 0:twojmax
        idxu_block[j+1] = idxu_count;
        for mb = 0:j
            for ma = 0:j
                idxu_count+=1;
            end
        end
    end

    idxu_max = idxu_count;

    # // index list for beta and B

    idxb_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1 - j2
            while j <= minimum([twojmax, j1+j2])
                if (j >= j1)
                    idxb_count+=1;
                end
                j+=2;
            end
        end
    end

    idxb_max = idxb_count;
    idxb = Vector{SNA_BINDICES}(undef, idxb_max)
    for k = 1:idxb_max
        idxb[k] = SNA_BINDICES(0, 0, 0)
    end
    idxb_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1-j2
            while j <= minimum([twojmax, j1+j2])
                if (j >= j1)
                    idxb[idxb_count+1] = SNA_BINDICES(j1, j2, j)
                    idxb_count+=1;
                end
                j +=2;
            end
        end
    end

    # // reverse index list for beta and b

    idxb_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1 - j2
            while j <= minimum([twojmax, j1+j2])
                if (j >= j1)
                    idxb_block[j1+1,j2+1,j+1] = idxb_count;
                    idxb_count+=1;
                end
                j +=2
            end
        end
    end

    # // index list for zlist
    idxz_count = 0;

    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1 - j2
            while j <= minimum([twojmax, j1+j2])
                mb = 0
                while 2*mb <= j
                    for ma = 0:j
                        idxz_count+=1;
                    end
                    mb += 1;
                end
                j+=2;
            end
        end
    end

    idxz_max = idxz_count;
    idxz = Vector{SNA_ZINDICES}(undef, idxz_max)
    for k = 1:idxz_max
        idxz[k] = SNA_ZINDICES(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    end

    idxz_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1-j2
            while j <= minimum([twojmax, j1+j2])
                idxz_block[j1+1,j2+1,j+1] = idxz_count;

                # // find right beta[jjb] entry
                # // multiply and divide by j+1 factors
                # // account for multiplicity of 1, 2, or 3
                mb = 0
                while 2*mb <= j
                    for ma = 0:j
                        ma1min = maximum([ 0, (2 * ma - j - j2 + j1) ÷ 2]);
                        ma2max = (2 * ma - j - (2 * ma1min - j1) + j2) ÷ 2;
                        na = minimum( [j1, (2 * ma - j + j2 + j1) ÷ 2] ) - ma1min + 1;
                        mb1min = maximum([0, (2 * mb - j - j2 + j1) ÷ 2]);
                        mb2max = (2 * mb - j - (2 * mb1min - j1) + j2) ÷ 2;
                        nb = minimum([ j1, (2 * mb - j + j2 + j1) ÷ 2] ) - mb1min + 1;
                        # // apply to z(j1,j2,j,ma,mb) to unique element of y(j)

                        jju = idxu_block[j+1] + (j+1)*mb + ma;

                        idxz[idxz_count+1] = SNA_ZINDICES(j1, j2, j,
                                        ma1min, ma2max, mb1min, mb2max,
                                        na, nb, jju)
                        idxz_count+=1;
                    end
                    mb +=1;
                end
                j+=2;
            end
        end
    end

    # Initialize Root array
    rootpqarray = zeros(twojmax+2, twojmax+2)        # rootpqarray
    rootpqarray = init_rootpqarray(twojmax, rootpqarray)

    # Initialize Clebsch Gordan
    cglist      = zeros(idxcg_max)                  # cglist
    cglist      = init_clebsch_gordan(twojmax, cglist)

    # Calculate number of coefficients
    ncoeff = ncount * num_elements * num_elements * num_elements

    return PrebuiltArrays(
    idxcg_block,
    idxu_block,
    idxb_block,
    idxz_block,
    idxcg_max,
    idxu_max,
    idxb_max,
    idxz_max,
    idxb,
    idxz,
    ncoeff,
    num_elements,
    rootpqarray,
    cglist,
    wself,
    bzero)
end

function deltacg(j1::Int, j2::Int, j::Int)
    sfaccg = factorial( (j1+j2+j) ÷ 2 + 1  )
    return sqrt(
            factorial( (j1+j2 -j) ÷ 2 ) *
            factorial( (j1-j2+j) ÷ 2  )  *
            factorial( (-j1+j2+j) ÷ 2 ) / sfaccg
    )
end

function init_clebsch_gordan(twojmax::Int, cglist::Vector{T}) where T<: AbstractFloat
    idxcg_count = 0;
    for j1 = 0:twojmax
        for j2 = 0:j1
            j = j1 - j2
            while j <= minimum([twojmax, j1+j2])
                for m1 = 0:j1
                    aa2 = 2 * m1 - j1;

                    for m2 = 0:j2
                        bb2 = 2 * m2 - j2;
                        m = (aa2 + bb2 + j) ÷ 2 ;

                        if (m < 0) || (m > j)
                            cglist[idxcg_count+1] = 0.0;
                            idxcg_count+=1;
                            continue
                        else
                            # println("j $j")
                            sum = 0.0;
                            min_z = maximum([0, maximum([-(j-j2+aa2) ÷ 2, -(j-j1-bb2) ÷ 2])])
                            max_z = minimum([ (j1 + j2 - j) ÷ 2, minimum([ (j1-aa2) ÷ 2, (j2+bb2) ÷ 2 ])])
                            for z = min_z:max_z
                                ifac = ( (z % 2)==0 ) ? 1.0 : -1.0;
                                sum += ifac /
                                    (factorial(Int(z)) *
                                    factorial( (j1 + j2 - j) ÷ 2 - z) *
                                    factorial( (j1 - aa2) ÷ 2 - z) *
                                    factorial( (j2 + bb2) ÷ 2 - z) *
                                    factorial( (j - j2 + aa2) ÷ 2 + z) *
                                    factorial( (j - j1 - bb2) ÷ 2 + z)) ;
                            end

                            cc2 = 2 * m - j;
                            dcg = deltacg(j1, j2, j);
                            sfaccg = sqrt(factorial( (j1 + aa2) ÷ 2)   *
                                        factorial( (j1 - aa2) ÷ 2)  *
                                        factorial( (j2 + bb2) ÷ 2)  *
                                        factorial( (j2 - bb2) ÷ 2)  *
                                        factorial( (j  + cc2) ÷ 2)  *
                                        factorial( (j  - cc2) ÷ 2) *
                                        (j + 1));

                            # println("j1 $j1, j2 $j2, m1 $m1, m2 $m2, j $j, m, $m")

                            #     println("dcg $dcg, sfaccg $sfaccg")
                            # cg = float(CG(j1, j2, j, m1, m2, m))
                            # println("cg $cg")
                            cglist[idxcg_count+1] = sum * dcg * sfaccg;
                            # cglist[idxcg_count+1] = cg
                            idxcg_count+=1;
                        end
                    end
                end
                j +=2;
            end
        end
    end
    return cglist
end

function init_rootpqarray(twojmax::Int, rootpqarray::Matrix{T}) where T<:AbstractFloat
    for p = 1:twojmax
        for q = 1:twojmax
            rootpqarray[p+1, q+1] = sqrt(p / q)
        end
    end
    return rootpqarray
end


text:
import RetroSignalModel as rs

conds = load_conditions()

const condkeys = (:Rtg1, :Rtg2, :Rtg3, :Mks, :s, :gfp, :Trans2Nuc)

@testset "Name availability (Conditions): $k" for cond in conds, k in condkeys
    @test k in keys(cond)
end

const paramkeys = (
    rs.n_S,
    rs.ksV,
    rs.ksD,
    rs.k2I,
    rs.k2M,
    rs.kn2M,
    rs.kBM,
    rs.knBM,
    rs.k13I,
    rs.k13IV,
    rs.k13ID,
    rs.k3A_c,
    rs.k3I_c,
    rs.k3I_n,
    rs.k13_c,
    rs.kn13_c,
    rs.k13_n,
    rs.kn13_n,
    rs.k1in,
    rs.k1out,
    rs.k3inA,
    rs.k3outA,
    rs.k3inI,
    rs.k3outI)

optimparams = load_parameters("solution_rtgMTK_optim.csv")

# Solutions of rtgM4
@testset "Name availability (Optim Parameters): $k" for param in optimparams, k in paramkeys
    @test k in keys(param)
end

randparams = load_parameters("solution_rtgM4.csv")

# Solutions of rtgM4
@testset "Name availability (Random Parameters): $k" for param in randparams, k in paramkeys
    @test k in keys(param)
end


text:
#=
runtime: in O(n^2.81)
accepts square matrices A,B (assumes nxn, mxm where n,m are powers of two)
returns matrix C
=#

function partition(A)
  #add try-catch
  #accepts square matrix A (assumes nxn, where n is power of two)
  n = length(A)[1]
  W = new matrix(  )
  X = new matrix(  )
  Y = new matrix(  )
  Z = new matrix(  )
  return W, X, Y, Z
end

function recursiveMultiply(A,B)
  n = length(A)[1]
  #matrix C
  C = new zeroesMatrix( length(A)[1,:],length(B)[:,1] )
  #trivial case
  if n == 1
     return A[1:1] * B[1:1]
  else
    partition(A)
    partition(B)
    #create 10 matrices (= Theta(n^2))
    S1 = B12 - B22
    S2 = A11 + A22
    S3 = A21 + A22
    S4 = B21 - B11
    S5 = A11 + A22
    S6 = B11 + B12
    S7 = A11 - A22
    S8 = B21 + B22
    S9 = A11 - A21
    S10 = B11 + B12
    #multiplications
    P1 = recursiveMultiply(A11, S1)
    P2 = recursiveMultiply(S2, B22)
    P3 = recursiveMultiply(S3, B11)
    P4 = recursiveMultiply(A22, S4)
    P5 = recursiveMultiply(S5, S6)
    P6 = recursiveMultiply(S7, S8)
    P7 = recursiveMultiply(S9, S10)
    #defining C
    C11 = P5 + P4 - P2 + P6
    C12 = P1 + P2
    C21 = P3 + P4
    C22 = P5 + P1 - P3 - P7
    return #build via blocks
  end
end

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing


text:
#!/usr/bin/env julia

p(n) = (reduce(+,map((x->(4/(1+(x*x)))),(0.5/n):(1/n):1))/n)

function picalc(numsteps)

  # Get number of processes.
  # np = nworkers()
  np = 1

  println("Calculating PI using:")
  println("  ", numsteps, " slices")
  println("  ", np, " worker(s)")
  start = time()

  mypi = p(numsteps)

  stop = time()

  elapsed = (stop - start)

  println("Obtained value of PI: ", mypi)
  println("Time taken: ", elapsed, " seconds")

end

numsteps=50000000

if length(ARGS) > 0
  numsteps = parse(Int, ARGS[1])
end

picalc(numsteps)


text:
struct Config
    beta_min::Float64
    beta_max::Float64
    delta_min::Float64
    delta_max::Float64
    mutate_add_rate::Float64
    mutate_delete_rate::Float64
    dist_coef_ids::Float64
    dist_coef_enh::Float64
    dist_coef_inh::Float64
    crossover::Symbol
    speciation_thresh::Float64
    init_species::Int64
    ga_population::Int64
    ga_crossover_rate::Float64
    ga_mutation_rate::Float64
    ga_num_generations::Int64
end

function Config(;beta_min=0.02, beta_max=0.5, delta_min=0.02, delta_max=0.5,
                mutate_add_rate=0.25, mutate_delete_rate=0.25,
                dist_coef_ids=0.75, dist_coef_enh=0.25, dist_coef_inh=0.25,
                crossover=:onepoint_crossover, speciation_thresh=0.15,
                init_species=20, ga_population=100, ga_crossover_rate=0.25,
                ga_mutation_rate=0.5, ga_num_generations=10)
    Config(beta_min, beta_max, delta_min, delta_max, mutate_add_rate, mutate_delete_rate,
           dist_coef_ids, dist_coef_enh, dist_coef_inh, crossover, speciation_thresh,
           init_species, ga_population, ga_crossover_rate, ga_mutation_rate,
           ga_num_generations)
end