text:
export GitCredential, PlainTextCred, SSHKeyCred, SSHAgentCred, DefaultCred

abstract GitCredential

type PlainTextCred <: GitCredential
    ptr::Ptr{Void}

    PlainTextCred(ptr::Ptr{Void}) = begin
        @assert ptr != C_NULL
        cred = new(ptr)
        return cred
    end
end

PlainTextCred(;username::MaybeString=nothing,
               password::MaybeString=nothing) = begin
    is(username, nothing) && throw(ArgumentError("username must be defined"))
    cred_ptr = Ptr{Void}[0]
    @check ccall((:git_cred_userpass_plaintext_new, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Uint8}, Ptr{Uint8}),
                 cred_ptr, username, password != nothing ? password : C_NULL)
    return PlainTextCred(cred_ptr[1])
end

type SSHKeyCred <: GitCredential
    ptr::Ptr{Void}

    SSHKeyCred(ptr::Ptr{Void}) = begin
        @assert ptr != C_NULL
        cred = new(ptr)
        return cred
    end
end

SSHKeyCred(;username::MaybeString=nothing,
            publickey::MaybeString=nothing,
            privatekey::MaybeString=nothing,
            passphrase::MaybeString=nothing) = begin
    is(username, nothing) && throw(ArgumentError("username must be defined"))
    is(privatekey, nothing) && throw(ArgumentError("privatekey must be defined"))
    cred_ptr = Ptr{Void}[0]
    @check ccall((:git_cred_ssh_key_new, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}),
                 cred,
                 username,
                 publickey != nothing ? publickey : C_NULL,
                 privatekey,
                 passphrase != nothing ? passphrase : C_NULL)
    return SSHKeyCred(cred_ptr[1])
end

type SSHAgentCred <: GitCredential
    ptr::Ptr{Void}

    SSHAgentCred(ptr::Ptr{Void}) = begin
        @assert ptr != C_NULL
    end
end

SSHAgentCred(username::String) = begin
    cred_ptr = Ptr{Void}[0]
    @check ccall((:git_cred_ssh_key_from_agent, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Uint8}), cred_ptr, username)
    return SSHAgentCred(cred_ptr[1])
end

type DefaultCred <: GitCredential
    ptr::Ptr{Void}

    DefaultCred(ptr::Ptr{Void}) = begin
        @assert ptr != C_NULL
        cred = new(ptr)
        return cred
    end
end

DefaultCred() = begin
    cred_ptr = Ptr{Void}[0]
    @check ccall((:git_cred_default_new, libgit2), Cint, (Ptr{Ptr{Void}},), cred_ptr)
    return DefaultCred(cred_ptr[1])
end


text:
##  SCRIPT NAME:   mUtilW_k.jl
##  DATE:          March 2015
#
#   This script calculates the wife's marginal utility with
#   respect to the family's public consumption.
#
#
##  ALEXANDROS THELOUDIS, UCL

##  ---------------------------------------------------------------------------------------------

function mUtilW_k(cW, K)

  marginal_utility = (1-alphaCW) * (K^(-gammaKW))

  return marginal_utility

end

ext: jl

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using DomainSets: ×

"Does `z` agree with `val` up to the given `threshold`?"
isaccurate(z, val, threshold = 1e-5) = abs(z-val) < threshold

function test_readme_examples()
    # Examples of README.md
    threshold = 1e-5
    @test isaccurate(integral(cos, 0..1.0), 0.8414709848078965)
    @test isaccurate(integral(exp(x) for x in 2..3), 12.69648082425702)
    @test isaccurate(integral(exp(x+y) for (x,y) in (0..1)^2), 2.9524924420120535)
    @test isaccurate(integrate(cos(x) for x in UnionDomain(0..1, 2..3))[1], 0.07329356604208204)
    @test isaccurate(integral( (sin(log(abs(t))) for t in  -1..1), LogSingPoint(0.0)), -1.0000000021051316)
    @test isaccurate(integral( ( exp(log(abs(x-y))) for (x,y) in (2..3) × (1..4) ), SingularDiagonal()),
2.333333333333333)
    @test isaccurate(integral(cos, ChebyshevTWeight()), 2.403939430634413)
    @test isaccurate(integral(cos(t)*1/sqrt(1-t^2) for t in  -1.0..1.0), 2.403939410869398)
    @test isaccurate(integrate(QuadAdaptive(atol=1e-3, rtol = 1e-3), t->cos(t^2), 0..10)[1], 0.6011251848111901, 1e-3)
    @test isaccurate(integral(Q_GaussLaguerre(10), cos), 0.5000005097999486)
    @test isaccurate(integral(cos(t)*exp(-t) for t in HalfLine()), 0.5)
end

function test_some_integrals()
    z1 = sin(1.0)
    @test abs(integral(cos, 0..1.0) - z1) < 1e-8
    I, E = integrate(cos, 0..1.0)
    @test abs(I-z1) < 1e-8
    @test E < 1e-8
    @test abs(integral(cos,uniondomain(0..1.0, 2..3.0)) - (sin(1)+sin(3)-sin(2))) < 1e-6

    f2 = x -> exp(x[1]+x[2])
    z2 = hcubature(f2, (0,0), (1,1))[1]
    @test abs(integral(f2, (0..1)^2) - z2) < 1e-6

    @test integral(exp, Point(0.5)) == 0.0

    # a mapped domain
    @test abs(integral(cos, DomainSets.MappedDomain(LinearMap(1/2), 0..1))-sin(2)) < 1e-6

    # diagonal singularity
    @test abs(integral((log(abs(x-y)) for (x,y) in UnitSquare()), SingularDiagonal())+1.5) < 1e-5

    # a mapped weight
    testmap_d = 2.0..4.0
    testmap_μ = lebesguemeasure(testmap_d)
    testmap_m = LinearMap(4.0)
    I1 = integral(cos, testmap_d, testmap_μ)
    @test abs(integral(t -> cos(inverse(testmap_m,t)), testmap_m.(2.0..4.0),
        DomainIntegrals.mappedmeasure(testmap_m, testmap_μ)) - I1)  < 1e-6
    @test abs(integral(t -> cos(inverse(testmap_m,t)), MappedDomain(inverse(testmap_m), 2.0..4.0),
        DomainIntegrals.mappedmeasure(testmap_m, testmap_μ)) - I1)  < 1e-6

    @test integrate(cos, EmptySpace()) === (0.0, 0.0)
    @test integrate(DomainIntegrals.Integrand{ComplexF64}(cos), EmptySpace()) === (0.0+0.0im, 0.0)
    @inferred integrate(cos, EmptySpace())

    # product weights
    μ1 = ChebyshevTWeight()
    @test abs(integral(x->cos(x[1]+x[2]), (-0.5..0.5)^2, productmeasure(μ1, μ1)) - 1.0049593511549983) < 1e-6
    @test abs(integral(x->cos(x[1]+x[2]), (-0.5..0.5)^3, productmeasure(μ1, μ1, μ1)) - 1.0523909717390105) < 1e-6

    # parametric domains
    @test abs(integral(x->cos(x[1]+x[2]), UnitCircle()) - 3.51314344095576) < 1e-6
    @test abs(integral(x->cos(x[1]+x[2]), 5*UnitCircle()) - 9.41394510) < 1e-6

    g(z) = 5/z
    @test abs(1/(2*pi*im)*integral(g, ComplexUnitCircle()) - 5) < 1e-6

    strategy = QuadAdaptive(maxevals=1e6)
    @test abs(integral(strategy, x->1, UnitDisk()) - pi) < 1e-5
end

function test_integrals()
    test_some_integrals()
    test_readme_examples()
end


text:
using Documenter, MRIReco

makedocs(
    modules = [MRIReco],
    format = Documenter.HTML(prettyurls = false),
    sitename = "Julia MRI Package",
    authors = "Tobias Knopp, Mirco Grosser",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "gettingStarted.md",
        "Overview" => "overview.md",
        "Acquisition Data" => "acquisitionData.md",
        "File Handling" => "filehandling.md",
        "Images" => "image.md",
        "Trajectory" => "trajectories.md",
        "Imaging Operators" => "operators.md",
        "Offresonance" => "offresonance.md",
        "Parallel Imaging" => "SENSE.md",
        "Compressed Sensing" => "compressedSensing.md",
        "Customize" => "custom.md",
        "API" => "API.md",
    ],
)

deploydocs(
    repo = "github.com/MagneticResonanceImaging/MRIReco.jl.git",
    target = "build",
)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the Direction Cosine Matrix (DCM).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export ddcm, orthonormalize

"""
    orthonormalize(dcm::DCM)

Perform the Gram-Schmidt orthonormalization process in the `dcm` and return the
new matrix.

!!! warning
    This function does not check if the columns of the input matrix span a
    three-dimensional space. If not, then the returned matrix should have `NaN`.
    Notice, however, that such input matrix is not a valid direction cosine
    matrix.

# Example

```julia-repl
julia> D = DCM(3I)

julia> orthonormalize(D)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
```
"""
function orthonormalize(dcm::DCM)
    e₁ = dcm[:, 1]
    e₂ = dcm[:, 2]
    e₃ = dcm[:, 3]

    en₁  = e₁ / norm(e₁)
    enj₂ = e₂ - (en₁ ⋅ e₂) * en₁
    en₂  = enj₂ / norm(enj₂)
    enj₃ = e₃ - (en₁ ⋅ e₃) * en₁
    enj₃ = enj₃ - (en₂ ⋅ enj₃) * en₂
    en₃  = enj₃ / norm(enj₃)

    return hcat(en₁, en₂, en₃)
end

################################################################################
#                                  Kinematics
################################################################################

"""
    ddcm(Dba::DCM, wba_b::AbstractArray)

Compute the time-derivative of the `dcm` that rotates a reference frame `a` into
alignment with the reference frame `b` in which the angular velocity of `b` with
respect to `a`, and represented in `b`, is `wba_b`.

# Example

```julia-repl
julia> D = DCM(1.0I);

julia> ddcm(D, [1, 0, 0])
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 0.0   0.0  0.0
 0.0   0.0  1.0
 0.0  -1.0  0.0
```
"""
function ddcm(Dba::DCM, wba_b::AbstractArray)
    # Auxiliary variable.
    w = wba_b

    # Check the dimensions.
    if length(wba_b) != 3
        throw(ArgumentError("The angular velocity vector must have three components."))
    end

    wx = SMatrix{3, 3}(
          0  , -w[3], +w[2],
        +w[3],   0  , -w[1],
        -w[2], +w[1],   0,
    )'

    # Return the time-derivative.
    return -wx * Dba
end


text:
module test_julia_intfuncs

using Test

@test ndigits( 5) == 1
@test ndigits(35) == 2

end # module test_julia_intfuncs


text:
"""
    module CplTestingBL
 Defines an equation set for testing coupling

 Defines kernels to evaluates RHS of

  d ϕ / dt = - div κ( ϕ_init ) grad ϕ

 subject to specified prescribed gradient or flux bc's

 - [`CplTestBL`](@ref)
     balance law struct created by this module

"""
module CplTestingBL

export CplTestBL
export PenaltyNumFluxDiffusive
export ExteriorBoundary
export CoupledPrimaryBoundary, CoupledSecondaryBoundary

const τ = 60 * 86400
using ClimateMachine.BalanceLaws: Auxiliary, BalanceLaw, Gradient, GradientFlux, Prognostic

import ClimateMachine.BalanceLaws:
    boundary_conditions,
    boundary_state!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    init_state_prognostic!,
    nodal_init_state_auxiliary!,
    source!,
    vars_state,
    wavespeed

using ClimateMachine.Mesh.Geometry: LocalGeometry
using ClimateMachine.MPIStateArrays

using ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux

import ClimateMachine.DGMethods.NumericalFluxes: numerical_boundary_flux_second_order!, numerical_flux_second_order!

using ClimateMachine.VariableTemplates

using LinearAlgebra
using StaticArrays

"""
    CplTestBL
Type that holds specification for a diffusion equation balance law instance.
"""
struct CplTestBL{BLP, BCS} <: BalanceLaw
    bl_prop::BLP
    boundaryconditions::BCS
end

l_type = CplTestBL

"""
  Extend the NumericalFluxSecondOrder to include a penalty term numerical
  flux formulation.
"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end

"""
  Set a default set of properties and their default values
  - init_aux_geom   :: function to initialize geometric terms stored in aux.
  - init_theta      :: function to set initial θ values.
  - source_theta    :: function to add a source term to θ.
  - calc_kappa_diff :: function to set diffusion coeffiecient(s).
  - get_wavespeed   :: function to return a wavespeed for Rusanov computations (there aren't any in this model)
  - get_penalty_tau :: function to set timescale on which to bring state+ and state- together
  - theta_shadow_boundary_flux :: function to set boundary flux into shadow variable for passing to coupler
"""
function prop_defaults()
    bl_prop = NamedTuple()

    function init_aux_geom(npt, elnum, x, y, z)
        return npt, elnum, x, y, z
    end
    # init_aux_geom(_...)=(return 0., 0., 0., 0., 0.)
    bl_prop = (bl_prop..., init_aux_geom = init_aux_geom)

    init_theta(_...) = (return 0.0)
    bl_prop = (bl_prop..., init_theta = init_theta)

    source_theta(_...) = (return 0.0)
    bl_prop = (bl_prop..., source_theta = source_theta)

    calc_kappa_diff(_...) = (return 0.0, 0.0, 0.0)
    bl_prop = (bl_prop..., calc_kappa_diff = calc_kappa_diff)

    get_wavespeed(_...) = (return 0.0)
    bl_prop = (bl_prop..., get_wavespeed = get_wavespeed)

    get_penalty_tau(_...) = (return 1.0)
    bl_prop = (bl_prop..., get_penalty_tau = get_penalty_tau)

    theta_shadow_boundary_flux(_...) = (return 0.0)
    bl_prop = (bl_prop..., theta_shadow_boundary_flux = theta_shadow_boundary_flux)

    coupling_lambda(_...) = (return 0.0)
    bl_prop = (bl_prop..., coupling_lambda = coupling_lambda)

    bl_prop = (bl_prop..., LAW = CplTestBL)
end

"""
  Declare single prognostic state variable, θ
  and a shadow variable, θ_boundary_export, for accumulating boundary flux over a single
  timestep.

  "Shadow" variables XX_boundary_export will be used to capture boundary fluxes that we
  want to export to coupling as time integrals. We use a shadow variable
  because we want to integrate over whatever timestepper is being used.
  Eventually we should have the ability to potentially use a 2d field here.
  The shadow variable needs to be zeroed at the start of each
  coupling cycle for a component.
"""
function vars_state(bl::l_type, ::Prognostic, FT)
    @vars begin
        θ::FT
        F_accum::FT # accumulated flux across boundary (atmosphere export)
    end
end

"""
  Auxiliary state variable symbols for array index and real world coordinates and for
  θ value at reference time used to compute κ.
"""
function vars_state(bl::l_type, ::Auxiliary, FT)
    @vars begin
        npt::Int
        elnum::Int
        xc::FT
        yc::FT
        zc::FT
        θⁱⁿⁱᵗ::FT
        θ_secondary::FT # stores opposite face for primary (atmospheric import)
        F_prescribed::FT # stores prescribed flux for secondary (ocean import)
    end
end

"""
  Gradient computation stage input (and output) variable symbols
"""
function vars_state(bl::l_type, ::Gradient, FT)
    @vars begin
        ∇θ::FT
        ∇θⁱⁿⁱᵗ::FT
    end
end

"""
  Flux due to gradient accumulation stage variable symbols
"""
function vars_state(bl::l_type, ::GradientFlux, FT)
    @vars begin
        κ∇θ::SVector{3, FT}
    end
end

"""
  Initialize prognostic state variables
"""
function init_state_prognostic!(bl::l_type, Q::Vars, A::Vars, geom::LocalGeometry, FT)
    npt = getproperty(geom, :n)
    elnum = getproperty(geom, :e)
    x = geom.coord[1]
    y = geom.coord[2]
    z = geom.coord[3]
    Q.θ = bl.bl_prop.init_theta(x, y, z, npt, elnum)
    Q.F_accum = 0
    nothing
end

"""
  Initialize auxiliary state variables
"""
function nodal_init_state_auxiliary!(bl::l_type, A::Vars, tmp::Vars, geom::LocalGeometry, _...)
    npt = getproperty(geom, :n)
    elnum = getproperty(geom, :e)
    x = geom.coord[1]
    y = geom.coord[2]
    z = geom.coord[3]
    A.npt, A.elnum, A.xc, A.yc, A.zc = bl.bl_prop.init_aux_geom(npt, elnum, x, y, z)
    A.θⁱⁿⁱᵗ = 0
    A.θ_secondary = 0
    A.F_prescribed = 0
    nothing
end

#====

Atmos

----

Land

====#

"""
  Set any source terms for prognostic state external sources.
  Also use to record boundary flux terms into shadow variables
  for export to coupler.
"""
function source!(bl::l_type, S::Vars, Q::Vars, G::Vars, A::Vars, _...)
    #S.θ=bl.bl_prop.source_theta(Q.θ,A.npt,A.elnum,A.xc,A.yc,A.zc,A.θ_secondary)
    # Record boundary condition fluxes as needed by adding to shadow
    # prognostic variable
    S.F_accum = (Q.θ - A.θ_secondary) * bl.bl_prop.coupling_lambda()
    nothing
end

"""
  No flux first order for diffusion equation, but we must define a stub.
"""
function flux_first_order!(bl::l_type, _...)
    nothing
end

"""
  Set values to have gradients computed.
"""
function compute_gradient_argument!(bl::l_type, G::Vars, Q::Vars, A::Vars, t)
    G.∇θ = Q.θ
    G.∇θⁱⁿⁱᵗ = A.θⁱⁿⁱᵗ
    nothing
end

"""
  Compute diffusivity tensor times computed gradient to give net gradient flux.
"""
function compute_gradient_flux!(bl::l_type, GF::Vars, G::Grad, Q::Vars, A::Vars, t)
    # "Non-linear" form (for time stepped)
    ### κ¹,κ²,κ³=bl.bl_prop.calc_kappa_diff(G.∇θ,A.npt,A.elnum,A.xc,A.yc,A.zc)
    # "Linear" form (for implicit)
    κ¹, κ², κ³ = bl.bl_prop.calc_kappa_diff(G.∇θⁱⁿⁱᵗ, A.npt, A.elnum, A.xc, A.yc, A.zc)
    # Maybe I should pass both G.∇θ and G.∇θⁱⁿⁱᵗ?
    GF.κ∇θ = Diagonal(@SVector([κ¹, κ², κ³])) * G.∇θ
    nothing
end

"""
  Pass flux components for second order term into update kernel.
"""
function flux_second_order!(bl::l_type, F::Grad, Q::Vars, GF::Vars, H::Vars, A::Vars, t)
    F.θ += GF.κ∇θ
    nothing
end

# Boundary conditions

"""
  Define boundary condition flags/types to iterate over, for now keep it simple.
"""
function boundary_conditions(bl::l_type, _...)
    bl.boundaryconditions
end

"""
  Set any first order numerical flux to null
  NumericalFluxFirstOrder is an abstract type that currently generalizes
  RusanovNumericalFlux, CentralNumericalFluxFirstOrder, RoeNumericalFlux, HLLCNumericalFlux.
"""
# No first order fluxes so numerical flux needed. NumericalFluxFirstOrder
function boundary_state!(
    nF::NumericalFluxFirstOrder,
    bc,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    nothing
end

## ExteriorBoundary
# flux is 0 across the boundary
struct ExteriorBoundary end

"""
  Zero normal gradient boundary condition.
"""
function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::ExteriorBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    nothing
end
function numerical_boundary_flux_second_order!(
    numerical_flux::NumericalFluxSecondOrder,
    bctype::ExteriorBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}
    fluxᵀn.θ = 0
end

## CoupledPrimaryBoundary
# compute flux based on opposite face
# also need to accumulate net flux across boundary
struct CoupledPrimaryBoundary end

function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::CoupledPrimaryBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ # Q⁺.θ=A⁺.θ_secondary
    nothing
end
function numerical_boundary_flux_second_order!(
    numerical_flux::NumericalFluxSecondOrder,
    bctype::CoupledPrimaryBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}

    fluxᵀn.θ = (state_prognostic⁻.θ - state_auxiliary⁺.θ_secondary) * balance_law.bl_prop.coupling_lambda() # W/m^2
end

# θ: J / m^3
# ∇θ: J / m^4
# κ: (m^2/s)
# κ∇θ: J/(m^2 s) = W/m^2 = J/m^3 * m/s
#

#  - clean up and write primer
#  - imports and exports
#  - add boundary auxiliary interface
#  - add boundary flux accumulation interface
#  - vector quantity
#  - run on sphere to check normal terms are all correct
#  - single stack integration
#  - Held-Suarez with wind stress and heat

## CoupledSecondaryBoundary
# use prescribed flux computed in primary
struct CoupledSecondaryBoundary end
function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::CoupledSecondaryBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    nothing
end
function numerical_boundary_flux_second_order!(
    numerical_flux::NumericalFluxSecondOrder,
    bctype::CoupledSecondaryBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}
    fluxᵀn.θ = -state_auxiliary⁺.F_prescribed
end

function wavespeed(bl::l_type, _...)
    # Used in Rusanov term.
    # Only active if there is a flux first order term?
    bl.bl_prop.get_wavespeed()
end

"""
  Penalty flux formulation of second order numerical flux. This formulation
  computes the CentralNumericalFluxSecondOrder term first (which is just the average
  of the + and - fluxes and an edge), and then adds a "penalty" flux that relaxes
  the edge state + and - toward each other.
"""
function numerical_flux_second_order!(
    ::PenaltyNumFluxDiffusive,
    bl::l_type,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    diff⁻::Vars{D},
    hyperdiff⁻::Vars{HD},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    diff⁺::Vars{D},
    hyperdiff⁺::Vars{HD},
    aux⁺::Vars{A},
    t,
) where {S, HD, D, A}

    numerical_flux_second_order!(
        CentralNumericalFluxSecondOrder(),
        bl,
        fluxᵀn,
        n,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        t,
    )

    Fᵀn = parent(fluxᵀn)
    FT = eltype(Fᵀn)
    tau = bl.bl_prop.get_penalty_tau()
    Fᵀn .+= tau * (parent(state⁻) - parent(state⁺))
end

end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using SubstitutionTilings
using SubstitutionTilings.Nilpotent

using Test
using LinearAlgebra

bhp = @time Dict(substitute(bhp_subst(), [a()], 3))

center_a = Dict([
    a(-2,-2,-2),
    a(0,-2,-2),
    a(2,-2,-2),
    a(-2,0,-2),
    a(0,0,-2),
    a(2,0,-2),
    a(-2,2,-2),
    a(0,2,-2),
    a(2,2,-2),
    b(-2,-2,0),
    b(0,-2,0),
    b(2,-2,0),
    b(-2,0,0),
    a(0,0,0),
    b(2,0,0),
    b(-2,2,0),
    b(0,2,0),
    b(2,2,0),
    a(-2,-2,2),
    a(0,-2,2),
    a(2,-2,2),
    a(-2,0,2),
    a(0,0,2),
    a(2,0,2),
    a(-2,2,2),
    a(0,2,2),
    a(2,2,2),
])
@time empirical_frequency(center_a, bhp)

patch = Dict([
    b(0,-2,0),
    a(0,0,0),
    b(0,2,0),
])
@time float(empirical_frequency(patch, bhp))

(collars, S) = total_collaring(bhp_subst(), center_a)
# 263 tiles

@time frequency(bhp_subst(), center_a, patch, 2)


text:
using Sparklines
using Base.Test

# Original tests from https://github.com/holman/spark/blob/master/spark-test.sh
# Copyright (c) Zach Holman, http://zachholman.com, MIT license

@test sprint(spark, [1,5,22,13,5])                == "▁▂█▅▂"
@test sprint(spark, [0,30,55,80,33,150])          == "▁▂▃▄▂█"
@test sprint(spark, [5.5,20])                     == "▁█"
@test sprint(spark, [1,2,3,4,100,5,10,20,50,300]) == "▁▁▁▁▃▁▁▁▂█"
@test sprint(spark, [1,50,100])                   == "▁▄█"
@test sprint(spark, [2,4,8])                      == "▁▃█"
@test sprint(spark, [1,2,3,4,5])                  == "▁▂▄▆█"

@test sprint(sparkln, [1,2,3,4,5])                == "▁▂▄▆█\n"


text:
# Model m4.3as.jl

using Pkg, DrWatson

begin
    using StanQuap
    using StatisticalRethinking
end

begin
    df = CSV.read(sr_datadir("Howell1.csv"), DataFrame; delim=';')
    df = filter(row -> row[:age] >= 18, df);
    mean_weight = mean(df.weight)
    #df.weight_c = df.weight .- mean_weight
end;

stan4_3 = "
data {
    int<lower=1> N;
    vector[N] weight;
    vector[N] height;
}
parameters {
    real a;
    real<lower=0> b;
    real<lower=0, upper=50> sigma;
}
model {
    // Define mu as a vector.
    vector[N] mu;

    // Priors for mu and sigma
    sigma ~ uniform(0 , 50);
    a ~ normal($(mean_weight), 20);
    b ~ lognormal(0, 1);

    // Observed heights
    for (i in 1:N) {
    mu[i] = a + b * (weight[i] - $(mean_weight));
    }
    height ~ normal(mu, sigma);
}
";

data = (N = length(df.height), height = df.height, weight = df.weight)
init = (a = 180.0, b = 1.0, sigma = 10.0)
q4_3s, m4_3s, o4_3s = stan_quap("m4.2s", stan4_3; data, init);

if q4_3s.converged
    quap4_3s_df = sample(q4_3s)          # DataFrame with samples
    precis(quap4_3s_df)
    post4_3s = read_samples(m4_3s)
    post4_3s |> display
end

# End of m4.3as.jl


text:
import Base.Pair
Pair(x::Modifier) = Pair{Symbol, Modifier}(nameof(typeof(x)), x)
Pair(x::NodeModifier) = Pair{Symbol, NodeModifier}(nameof(typeof(x)), x)
Pair(x::EdgeModifier) = Pair{Symbol, EdgeModifier}(nameof(typeof(x)), x)

function moddict(x::Modifier)
    Dict(Pair(x))
end

function moddict(x::Vector{M} where {M <: Modifier})
    Dict(Pair.(x)...)
end

moddict(x::Dict{Symbol, Modifier}) = x

modifiers(x) = x.modifiers

ModifiedEdge(edge::ModifiedEdge, modifier::EdgeModifier) = ModifiedEdge(edge.edge, merge(modifiers(edge), moddict(modifier)))
ModifyingNode(node::ModifyingNode, modifier::EdgeModifier) = ModifyingNode(node.node, merge(modifiers(node), moddict(modifier)))
ModifiedNode(node::ModifiedNode, modifier::NodeModifier) = ModifiedNode(node.node, merge(modifiers(node), moddict(modifier)))

ModifiedEdge(edge::ModifiedEdge, modifier::VecOrMat{M} where {M <: EdgeModifier}) = ModifiedEdge(edge.edge, merge(modifiers(edge), moddict(vec(modifier))))
ModifyingNode(node::ModifyingNode, modifier::VecOrMat{M} where {M <: EdgeModifier}) = ModifyingNode(node.node, merge(modifiers(node), moddict(vec(modifier))))
ModifiedNode(node::ModifyingNode, modifier::VecOrMat{M} where {M <: NodeModifier}) = ModifiedNode(node.node, merge(modifiers(node), moddict(vec(modifier))))

ModifiedEdge(edge::Edge, modifier::EdgeModifier) = ModifiedEdge(edge, moddict(modifier))
ModifyingNode(node::Node, modifier::EdgeModifier) = ModifyingNode(node, moddict(modifier))
ModifiedNode(node::Node, modifier::NodeModifier) = ModifiedNode(node, moddict(modifier))

ModifiedEdge(edge::Edge, modifier::VecOrMat{M} where {M <: EdgeModifier}) = ModifiedEdge(edge, moddict(vec(modifier)))
ModifyingNode(node::Node, modifier::VecOrMat{M} where {M <: EdgeModifier}) = ModifyingNode(node, moddict(vec(modifier)))
ModifiedNode(node::Node, modifier::VecOrMat{M} where {M <: NodeModifier}) = ModifiedNode(node, moddict(vec(modifier)))

ModifiedEdge(edge::Arrow, modifier) = Arrow(ModifiedEdge.(unarrow(edge), Ref(modifier)), edge.lhs, edge.rhs)
import Base.*

@communative function *(x::AbstractEdge, y)
    ModifiedEdge(x, y)
end

@communative function *(x::AbstractEdge, y::VecOrMat{T} where T)
    ModifiedEdge(x, vec(y))
end

@communative function *(y::AbstractNode, x::EdgeModifier)
    ModifyingNode(y, x)
end

@communative function *(y::AbstractNode, x::NodeModifier)
    ModifiedNode(y, x)
end

@communative function *(x::AbstractNode, y::VecOrMat{T} where T)
    ModifyingNode(x, vec(y))
end

function *(y::EdgeModifier, x::EdgeModifier)
    vec([y x])
end

@communative function *(y::EdgeModifier, x::VecOrMat{T} where {T <: EdgeModifier})
    vec([y x...])
end

function *(y::VecOrMat{M1} where {M1<:EdgeModifier}, x::VecOrMat{M2} where {M2<:EdgeModifier})
    vec([vec(y)... vec(x)...])
end

import Base.^

@communative function ^(x::AbstractNode, y::NodeModifier)
    ModifiedNode(x, y)
end

@communative function ^(x::AbstractNode, y::VecOrMat{T} where T)
    ModifiedNode(x, vec(y))
end

@communative function ^(y::NodeModifier, x::VecOrMat{T} where {T <: NodeModifier})
    vec([y x...])
end

function ^(y::NodeModifier, x::NodeModifier)
    vec([y x])
end

function ^(y::VecOrMat{M1} where {M1<:NodeModifier}, x::VecOrMat{M2} where {M2<:NodeModifier})
    vec([vec(y)... vec(x)...])
end

# Edge(Node(:a), Node(:b)) ≠ Edge(Node(:b), Node(:a))
# not communative

edges = (:DirectedEdge, :UndirectedEdge)
for e in edges
    @eval $e(y::AbstractNode, x::ModifyingNode) = ModifiedEdge($e(y, x.node), modifiers(x))
    @eval $e(y::ModifyingNode, x::AbstractNode) = ModifiedEdge($e(y.node, x), modifiers(y))
    @eval $e(y::ModifyingNode, x::ModifyingNode) = ModifiedEdge($e(y.node, x.node), merge(modifiers(y), modifiers(x)))
end


text:
"""
    function exponentiate(A, t::Number, x; kwargs...)
    function exponentiate(A, t::Number, x, algorithm)

Compute ``y = exp(t*A) x``, where `A` is a general linear map, i.e. a `AbstractMatrix` or
just a general function or callable object and `x` is of any Julia type with vector like
behavior.

### Arguments:
The linear map `A` can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object that implements the action of the linear map on a vector. If `A` is an
`AbstractMatrix`, `x` is expected to be an `AbstractVector`, otherwise `x` can be of any
type that behaves as a vector and supports the required methods (see KrylovKit docs).

The time parameter `t` can be real or complex, and it is better to choose `t` e.g. imaginary
and `A` hermitian, then to absorb the imaginary unit in an antihermitian `A`. For the
former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to
the exponential action of the linear map is obtained. The argument `x` can be of any type
and should be in the domain of `A`.

### Return values:
The return value is always of the form `y, info = exponentiate(...)` with
*   `y`: the result of the computation, i.e. `y = exp(t*A)*x`
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: 0 or 1 if the solution `y` was approximated up to the
        requested tolerance `tol`.
    -   `info.residual::Nothing`: value `nothing`, there is no concept of a residual in
        this case
    -   `info.normres::Real`: a (rough) estimate of the error between the approximate and
        exact solution
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
    By default (i.e. if `verbosity = 0`, see below), no warning is printed if the solution
    was not found with the requested precion, so be sure to check `info.converged == 1`.

### Keyword arguments:
Keyword arguments and their default values are given by:
*   `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
*   `krylovdim = 30`: the maximum dimension of the Krylov subspace that will be constructed.
    Note that the dimension of the vector space is not known or checked, e.g. `x₀` should
    not necessarily support the `Base.length` function. If you know the actual problem
    dimension is smaller than the default value, it is useful to reduce the value of
    `krylovdim`, though in principle this should be detected.
*   `tol = 1e-12`: the requested accuracy per unit time, i.e. if you want a certain
    precision `ϵ` on the final result, set `tol = ϵ/abs(t)`. If you work in e.g. single
    precision (`Float32`), you should definitely change the default value.
*   `maxiter::Int = 100`: the number of times the Krylov subspace can be rebuilt; see below
    for further details on the algorithms.
*   `issymmetric`: if the linear map is symmetric, only meaningful if `T<:Real`
*   `ishermitian`: if the linear map is hermitian
The default value for the last two depends on the method. If an `AbstractMatrix` is used,
`issymmetric` and `ishermitian` are checked for that matrix, ortherwise the default values
are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`.

### Algorithm
This is actually a simple wrapper over more general method [`expintegrator`](@ref) for
for integrating a linear non-homogeneous ODE.
"""
function exponentiate end

exponentiate(A, t::Number, v; kwargs...) = expintegrator(A, t, v; kwargs...)
exponentiate(A, t::Number, v, alg::Union{Lanczos,Arnoldi}) = expintegrator(A, t, (v,), alg)

# old implmentation: replaced with more general and accurate expintegrator
# function exponentiate(A, t::Number, v, alg::Lanczos)
#     # process initial vector and determine result type
#     β = norm(v)
#     Av = apply(A, v) # used to determine return type
#     numops = 1
#     T = promote_type(eltype(Av), typeof(β), typeof(t))
#     S = real(T)
#     w = mul!(similar(Av, T), v, 1/β)
#
#     # krylovdim and related allocations
#     krylovdim = alg.krylovdim
#     UU = Matrix{S}(undef, (krylovdim, krylovdim))
#     yy1 = Vector{T}(undef, krylovdim)
#     yy2 = Vector{T}(undef, krylovdim)
#
#     # initialize iterator
#     iter = LanczosIterator(A, w, alg.orth, true)
#     fact = initialize(iter; verbosity = alg.verbosity-2)
#     numops += 1
#     sizehint!(fact, krylovdim)
#
#     # time step parameters
#     sgn = sign(t)
#     τ::S = abs(t)
#     Δτ::S = τ
#
#     # tolerance
#     η::S = alg.tol / τ # tolerance per unit step
#     totalerr = zero(η)
#
#     δ::S = 0.9 # safety factor
#
#     # start outer iteration loop
#     maxiter = alg.maxiter
#     numiter = 0
#     while true
#         numiter += 1
#         Δτ = numiter == maxiter ? τ : min(Δτ, τ)
#
#         # Lanczos or Arnoldi factorization
#         while normres(fact) > η && length(fact) < krylovdim
#             fact = expand!(iter, fact; verbosity = alg.verbosity-2)
#             numops += 1
#         end
#         K = fact.k # current Krylov dimension
#         V = basis(fact)
#
#         # Small matrix exponential and error estimation
#         U = copyto!(view(UU, 1:K, 1:K), I)
#         H = rayleighquotient(fact) # tridiagonal
#         D, U = tridiageigh!(H, U)
#
#         # Estimate largest allowed time step
#         ϵ::S = zero(η)
#         while true
#             ϵ₁ = zero(eltype(H))
#             ϵ₂ = zero(eltype(H))
#             @inbounds for k = 1:K
#                 ϵ₁ += U[K,k] * exp(sgn * Δτ/2 * D[k]) * conj(U[1,k])
#                 ϵ₂ += U[K,k] * exp(sgn * Δτ * D[k]) * conj(U[1,k])
#             end
#             ϵ = normres(fact) * ( 2*abs(ϵ₁)/3 + abs(ϵ₂)/6 )
#             # error per unit time: see Lubich
#
#             if ϵ < δ * η || numiter == maxiter
#                 break
#             else # reduce time step
#                 Δτ = round(δ * (η / ϵ)^(1/krylovdim) * Δτ; sigdigits=2)
#             end
#         end
#
#         # Apply time step
#         totalerr += Δτ * ϵ
#         y1 = view(yy1, 1:K)
#         y2 = view(yy2, 1:K)
#         @inbounds for k = 1:K
#             y1[k] = exp(sgn*Δτ*D[k])*conj(U[1,k])
#         end
#         y2 = mul!(y2, U, y1)
#
#         # Finalize step
#         w = mul!(w, V, y2)
#         τ -= Δτ
#
#         if alg.verbosity > 1
#             msg = "Lanczos exponentiate in iteration $numiter: "
#             msg *= "reached time " * @sprintf("%.2e", abs(t)-τ)
#             msg *= ", total error = " * @sprintf("%.4e", totalerr)
#             @info msg
#         end
#
#         if iszero(τ) # should always be true if numiter == maxiter
#             w = rmul!(w, β)
#             converged = totalerr < alg.tol ? 1 : 0
#             if alg.verbosity > 0
#                 if converged == 0
#                     @warn """Lanczos exponentiate finished without convergence after $numiter iterations:
#                     total error = $totalerr"""
#                 else
#                     @info """Lanczos exponentiate finished after $numiter iterations: total error = $totalerr"""
#                 end
#             end
#             return w, ConvergenceInfo(converged, nothing, totalerr, numiter, numops)
#         else
#             normw = norm(w)
#             β *= normw
#             w = rmul!(w, inv(normw))
#             fact = initialize!(iter, fact)
#         end
#     end
# end


text:
module project_euler
    export prime_factors

    # It isn't the best algorithm, but is easy to read
    function prime_factors(limit::Int)
        factors = []
        p = 2
        n = limit

        while n >= p*p
            if 0 == n%p
                push!(factors, p)
                n = n/p
            else
                p += 1
            end
        end

        return push!(factors, n)
    end
end

max_stars_repo_name: Fazendaaa/project-euler

max_issues_repo_name: Fazendaaa/project-euler

max_forks_repo_name: Fazendaaa/project-euler


text:
using StabTrussTopOpt
sto = StabTrussTopOpt

import JSON
using Dates
using Printf
using LinearAlgebra: norm

function parse_TO_problem(ground_fp, load_supp_fp)::sto.TOProblemPar
    ndim, n, m, X, T, B = parse_TO_ground_json(ground_fp)
    F, Bound = parse_TO_support_load_json(load_supp_fp)

    dof_map = zeros(Int, m, 2*ndim) # e_id -> nodal dof
    id_minus = collect(ndim-1:-1:0)
    for i=1:m
        dof_map[i, 1:ndim] = ndim*T[i,1] * ones(ndim) - id_minus
        dof_map[i, ndim+1:2*ndim] = ndim*T[i,2] * ones(ndim) - id_minus
    end

    return sto.TOProblemPar(ndim, n, m, X, T, B, F, Bound, dof_map)
end

function parse_TO_ground_json(file_path::String)
    data = Dict()
    open(file_path, "r") do f
        data_txt = read(f, String)
        data = JSON.parse(data_txt)
    end
    ndim = data["dimension"]
    n = data["node_num"]
    m = data["element_num"]

    X = Matrix{Float64}(undef, n, ndim)
    T = Matrix{Int}(undef, m, 2)
    B = Matrix{Float64}(undef, m, ndim*2)

    # get node coord
    for i=1:n
        X[i,:] = data["node_list"][i]["point"]'
    end

    # get element node ids
    for i=1:m
        T[i,:] = data["element_list"][i]'
    end

    # build compatability matrix
    for i=1:m
        u_id = T[i,1]
        v_id = T[i,2]
        # compute length
        L = norm(X[v_id,:] - X[u_id,:])
        # compute direction cosine
        for k=1:ndim
            B[i, 1:ndim] = (X[v_id,:] - X[u_id,:]) ./ L
            B[i, ndim+1:2*ndim] = -B[i, 1:ndim]
        end
    end

    return ndim, n, m, X, T, B
end

function parse_TO_support_load_json(file_path::String)
    data = Dict()
    open(file_path, "r") do f
        data_txt = read(f, String)
        data = JSON.parse(data_txt)
    end

    ndim = data["dimension"]
    n = data["node_num"]

    n_load_nodes = length(data["point_load_list"])
    @assert(n_load_nodes > 0)
    F = zeros(Float64, n, ndim) # load
    for i=1:n_load_nodes
        load_v = data["point_load_list"][i]["node_id"]
        F[load_v, :] = data["point_load_list"][i]["force"]'
    end

    n_supp_nodes = length(data["support_node_list"])
    @assert(n_supp_nodes > 0)
    Bound = zeros(Int, n, ndim)
    for i=1:n_supp_nodes
        supp_v = data["support_node_list"][i]["node_id"]
        Bound[supp_v, :] = data["support_node_list"][i]["fixities"]'
    end

    # TODO: self-weight
    return F, Bound
end


text:
abstract type Problem end

struct TrajectoryOptimizationProblem <: Problem
    model  # model

    h      # time step
    T::Int # planning horizon

    obj    # objective

    con    # constraints
    xl     # state lower bound
    xu     # state upper bound
    ul     # control lower bound
    uu     # control upper bound

    num_var::Int      # number of decision variables
    num_con::Int      # number of constraints

    idx    # indices
end

function trajectory_optimization(model, obj, T;
        h = 0.0,
        xl = [-1.0 * Inf * ones(model.n) for t = 1:T],
        xu = [Inf * ones(model.n) for t = 1:T],
        ul = [-1.0 * Inf * ones(model.m) for t = 1:T-1],
        uu = [Inf * ones(model.m) for t = 1:T-1],
        con = EmptyConstraints(),
        dynamics = true)

    # indices
    idx = indices(model.n, model.m, T)

    # decision variables
    num_state = model.n * T
    num_control = model.m * (T - 1)
    num_var = num_state + num_control

    # constraints
    dynamics && (con = add_constraints(con, dynamics_constraints(model, T)))

    num_con = con.n

    prob = TrajectoryOptimizationProblem(
            model,
            h,
            T,
            obj,
            con, xl, xu, ul, uu,
            num_var, num_con,
            idx
            )

   return prob
end

function trajectory_optimization_problem(model, obj, T;
        xl = [-1.0 * Inf * ones(model.n) for t = 1:T],
        xu = [Inf * ones(model.n) for t = 1:T],
        ul = [-1.0 * Inf * ones(model.m) for t = 1:T-1],
        uu = [Inf * ones(model.m) for t = 1:T-1],
        h = 0.0,
        con = EmptyConstraints(),
        dynamics = true)

   prob = trajectory_optimization(model, obj, T,
                h = h,
                xl = xl,
                xu = xu,
                ul = ul,
                uu = uu,
                con = con,
                dynamics = dynamics)

   prob_moi = moi_problem(prob)

   return prob_moi
end

function pack(X, U, prob::TrajectoryOptimizationProblem)
    T = prob.T
    idx = prob.idx

    Z = zeros(prob.num_var)

    for t = 1:T
        Z[idx.x[t]] = X[t]

        t == T && continue

        Z[idx.u[t]] = U[t]
    end

    return Z
end

function unpack(Z, prob::TrajectoryOptimizationProblem)
    T = prob.T
    idx = prob.idx

    X = [Z[idx.x[t]] for t = 1:T]
    U = [Z[idx.u[t]] for t = 1:T-1]

    return X, U
end

function primal_bounds(prob::TrajectoryOptimizationProblem)
    T = prob.T
    idx = prob.idx

    Zl = -Inf * ones(prob.num_var)
    Zu = Inf * ones(prob.num_var)

    for t = 1:T
        Zl[idx.x[t]] = prob.xl[t]
        Zu[idx.x[t]] = prob.xu[t]

        t == T && continue

        Zl[idx.u[t]] = prob.ul[t]
        Zu[idx.u[t]] = prob.uu[t]
    end

    return Zl, Zu
end

function constraint_bounds(prob::TrajectoryOptimizationProblem)

    cl = zeros(prob.num_con)
    cu = zeros(prob.num_con)

    cu[prob.con.ineq] .= Inf

    return cl, cu
end

function eval_objective(prob::TrajectoryOptimizationProblem, Z)
    objective(Z, prob)
end

function eval_objective_gradient!(∇J, Z, prob::TrajectoryOptimizationProblem)
    ∇J .= 0.0
    objective_gradient!(∇J, Z, prob)
    return nothing
end

function eval_constraint!(c, Z, prob::TrajectoryOptimizationProblem)
    constraints!(c, Z, prob)
    return nothing
end

function eval_constraint_jacobian!(∇c, Z, prob::TrajectoryOptimizationProblem)
    constraints_jacobian!(∇c, Z, prob)
    return nothing
end

function sparsity_jacobian(prob::TrajectoryOptimizationProblem)
    constraints_sparsity(prob)
end

"""
    objective
"""
objective(Z, prob::TrajectoryOptimizationProblem) = objective(Z, prob.obj, prob.idx, prob.T)
objective_gradient!(∇J, Z, prob::TrajectoryOptimizationProblem) = objective_gradient!(∇J, Z, prob.obj, prob.idx, prob.T)

"""
    constraints
"""
constraints!(c, Z, prob::TrajectoryOptimizationProblem) = constraints!(c, Z, prob.con, prob.model, prob.idx, prob.h, prob.T)
constraints_jacobian!(∇c, Z, prob::TrajectoryOptimizationProblem) = constraints_jacobian!(∇c, Z, prob.con, prob.model, prob.idx, prob.h, prob.T)
constraints_sparsity(prob::TrajectoryOptimizationProblem; shift_row = 0, shift_col = 0) = constraints_sparsity(prob.con, prob.model, prob.idx, prob.T; shift_row = shift_row, shift_col = shift_col)

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# TODO: line and column give character offset not byte index

import Base: match, rsearch, length, showerror, isequal, show

export _Compound, LazyReference, unicode, IncompleteParseError, showerror, ParseError, Expression, Literal, Regex, Sequence, OneOf, Not, Optional, ZeroOrMore, OneOrMore, parse, Lookahead, isequal, @p_str, @p_mstr, show, as_rule, match

abstract Expression # -ism
abstract ParseException

type ParseError <: ParseException
    text::String
    expr::Expression
    pos::Int
end

# TODO: DRY
# A future Julia feature will be abstract types with fields
type IncompleteParseError <: ParseException
    text::String
    expr::Expression
    pos::Int
end

function showerror(io::IO, e::ParseError)
    print_escaped(io, "Rule '$(e.expr.name)' didn't match at '$(e.text[e.pos:min(end,e.pos+20)])' (line $(line(e)), column $(column(e))).", "")
end

function unicode(e::ParseException)
    io = IOString()
    showerror(io, e)
    seekstart(io)
    readline(io)
end

function showerror(io::IO, e::IncompleteParseError)
    print_escaped(io, "Rule '$(e.expr.name)' matched in its entirety, but it didn't consume all the text. The non-matching portion of the text begins with '$(e.text[e.pos+1:min(end,e.pos + 20)])' (line $(line(e)), column $(column(e))).", "")
end

# TODO: replace all pos[:] with SubString
line(e::ParseException) = 1 + count(SubString(e.text, 1, prevind(e.text, e.pos))) do c c == '\n' end

column(e::ParseException) = ind2chr(e.text, e.pos - rsearch(e.text, '\n', prevind(e.text, e.pos)))

function parse(expr::Expression, text::String, pos::Int=1)
    node = match(expr, text, pos)
    #@show node
    # workaround here
    if max(0, _end(node)) != endof(text)
        throw(IncompleteParseError(text, expr, _end(node)))
    end
    node
end

function match(expr::Expression, text::String, pos::Int=1)
    err::ParseError = ParseError(text, expr, pos)
    cache = Dict{Uint, AnyNode}()
    sizehint(cache, 10*length(text))
    node = _match(expr, text, pos, cache, err)
    if isempty(node)
        throw(err)
    end
    node
end

function hasfield(t, field)
    try
        getfield(t, field)
    catch
        return false
    end
    true
end

function _match(expr::Expression, text::String, pos::Int, cache::Dict, err::ParseError)
    # This is interesting. Tuple keys are slow so
    #   key = (object_id(expr), pos)
    # results in a much slower program
    # This key, however, is subject to collision, however unlikely so is a bug in the now much-faster program
    # When tuple keys are fixed, go back to the proper key.
    key::Uint = (uint(object_id(expr) << 30) | pos)
    local node
    if !haskey(cache, key)
        # TODO: hottest of all hot spots here
        #@show expr.name, typeof(expr), pos
        #println(text[pos:min(pos + 20, end)])
        node = cache[key] = _uncached_match(expr, text, pos, cache, err)

        #@show node
        #if isa(node, MatchNode)
        #    @show nodetext(node, text)
        #end
        # For fun, to see memory savings replace above line with:
        # node = _uncached_match(expr, text, pos, cache, err)
        # Surprisingly, memory usage went up to 183M without the cache from 139M with.
        # Either way that is a lot of fucking memory allocated to parse a few kB of text.
        # Hopefully @time reports it some wierd way and it's not really chewing up so
        # much memory.
    else
        node = cache[key]
    end

    if isempty(node) && pos >= err.pos && (
        !isempty(expr.name) || !hasfield(err.expr, :name))
        err.expr = expr
        err.pos = pos
    end

    node
end

function as_rule(expr::Expression)
    if !isempty(expr.name)
        return "$(expr.name) = $(_as_rhs(expr))"
    end
    return _as_rhs(expr)
end

function _as_rhs(expr::Expression)
    error("_as_rhs not implemented in "*string(typeof(expr))*" only subclasses")
end

type Literal <: Expression
    literal::String
    name::String

    Literal(literal, name) = new(literal, name)
end

# TODO: get rid of all keyword default constructors because slow
# Also, try not to use keyword constructors in Grammars.jl
Literal(literal::String; name::String="") = Literal(literal, name)

function _uncached_match(literal::Literal, text::String, pos::Int, cache::Dict, err::ParseError)
    if beginswith(SubString(text, pos, endof(text)), literal.literal)
        return Node(literal.name, text, pos, prevind(text, pos + endof(literal.literal)))
    end
    EmptyNode()  # "Empty" node == no match
end

function _as_rhs(literal::Literal)
    "\"" * literal.literal * "\""
end

# Stefan won't accept my anchor-regex patch. Fine
function MyBaseRegex(pattern::String, flags::String)
    options = Base.DEFAULT_OPTS
    PCRE = Base.PCRE
    for f in flags
        options |= f=='i' ? PCRE.CASELESS  :
                   f=='m' ? PCRE.MULTILINE :
                   f=='s' ? PCRE.DOTALL    :
                   f=='x' ? PCRE.EXTENDED  :
                   f=='a' ? PCRE.ANCHORED  :
                   error("unknown regex flag: $f")
    end
    Base.Regex(pattern, options)
end

type Regex <: Expression
    name::String
    re::Base.Regex
    options
    original_re

    function Regex(pattern, name="", options="")
        original_re = pattern
        @assert(length(pattern) > 0)

        new(name, MyBaseRegex(pattern, options * "a"), options, original_re)
    end
end

Regex(pattern::String; name="", options="") = Regex(pattern, name, options)

# This is really something magical. Say you want an unescaped string, just say:
#    s = p"\'"
# And you will have a 2 character string where the characters are '\' and '''.
# How does the macro do it? I don't know. Thanks `regex.jl`

function _uncached_match(regex::Regex, text::String, pos::Int, cache::Dict, err::ParseError)
# TODO: want but can't because bug
    m = Base.match(regex.re, text, pos)
    if m == nothing
        return EmptyNode()
    end
    Node(regex.name, text, pos, prevind(text, pos + endof(m.match))) # no longer special , m)
end

function _as_rhs(regex::Regex)
    "~\"$(regex.original_re)\"$(regex.options)"
end

abstract _Compound <: Expression
# has members

# A bit convoluted. You can efficiently create:
#   Sequence(name, members...)
#   Sequence(members...)  # no name
# And less-efficiently
#   Sequence(members..., name="name")
type Sequence <: _Compound
    name
    members

    Sequence(name::String, members::Expression...) = new(name,members)
    Sequence(members::Expression...) = new("", members)
end

# TODO: DRY: generic constructor for _Compound things?
Sequence(members::Expression...; name::String="") = Sequence(name, members...)

function isequal(a::Expression, b::Expression)
    typeof(a) == typeof(b) || return false

    for field in names(a)
# special handling for the re field
        if field == :re
            a.re.pattern == b.re.pattern || return false
        else
            getfield(a, field) == getfield(b, field) || return false
        end
    end
    true
end

using Base.Test
function Regex(pattern::Base.Regex; name="", options="")
# If someone uses a regex string, unpack it and use it anyway
    Regex(pattern.pattern, name, options)
end

# TODO: Why can't isequal just work for any two same-type objects?
# TODO: ask julia-users
function isequal(a::_Compound, b::_Compound)
    typeof(a) == typeof(b) || return false
    length(a.members) == length(b.members) || return false
    for (i, j) in zip(a.members, b.members)
        i == j || return false
    end
    for field in names(a)
        if field != :members
            getfield(a, field) == getfield(b, field) || return false
        end
    end
    true
end

function _uncached_match(sequence::Sequence, text::String, pos::Int, cache::Dict, err::ParseError)
    new_pos = pos
    length_of_sequence = 0
    first = true
    local children
    for m in sequence.members
        node = _match(m, text, new_pos, cache, err)
        if isempty(node)
            return node
        end
        if first
            children = MatchNode[node]
            first = false
        else
            push!(children, node)
        end
        new_pos += textlength(node)
        length_of_sequence += textlength(node)
    end
    # Hooray! We got through all the members!
    return Node(sequence.name, text, pos, prevind(text, pos + length_of_sequence), tuple(children...))
end

function _as_rhs(sequence::Sequence)
    return join([m.name for m in sequence.members], " ")
end

type OneOf <: _Compound
    name::String
    members # ::(Expression...) -- crashes resolve_refs

    OneOf(name::String, members::Expression...) = new(name,members)
    OneOf(members::Expression...) = new("", members)
end

OneOf(members::Expression...; name::String="") = OneOf(name, members...)

function _uncached_match(oneof::OneOf, text::String, pos::Int, cache::Dict, err::ParseError)
    local members = oneof.members
    for (i, m) in enumerate(members)
        node = _match(m, text, pos, cache, err)
        if !isempty(node)
            return Node(oneof.name, text, pos, _end(node), (node,))
        end
    end
    return EmptyNode()
end

function _as_rhs(e::OneOf)
    join(e.members, " / ")
end

type Lookahead <: _Compound
    name
    members

    Lookahead(name::String, members::Expression...) = new(name,members)
    Lookahead(members::Expression...) = new("", members)
end

Lookahead(members::Expression...; name::String="") = Lookahead(name, members...)

function _uncached_match(self::Lookahead, text::String, pos::Int, cache::Dict, err::ParseError)
    node = _match(self.members[1], text, pos, cache, err)
    if !isempty(node)
        return Node(self.name, text, pos, prevind(text, pos))
    end
    return EmptyNode()
end

function _as_rhs(e::Lookahead)
    "&" * join(e.members, " / ")
end

type Not <: _Compound
    name
    members

    Not(name::String, members::Expression...) = new(name,members)
    Not(members::Expression...) = new("", members)
end

Not(members::Expression...; name::String="") = Not(name, members...)

function _uncached_match(self::Not, text::String, pos::Int, cache::Dict, err::ParseError)
    node = _match(self.members[1], text, pos, cache, err)
    if isempty(node)
        return Node(self.name, text, pos, prevind(text, pos))
    end
    return EmptyNode()
end

function _as_rhs(e::Not)
    "!" * join(e.members, " ")
end

type Optional <: _Compound
    name
    members

    Optional(name::String, members::Expression...) = new(name,members)
    Optional(members::Expression...) = new("", members)
end

Optional(members::Expression...; name::String="") = Optional(name, members...)

function _uncached_match(self::Optional, text::String, pos::Int, cache::Dict, err::ParseError)
    node = _match(self.members[1], text, pos, cache, err)
    if isempty(node)
        return Node(self.name, text, pos, prevind(text, pos))
    end
    return Node(self.name, text, pos, _end(node), (node,))
end

function _as_rhs(e::Optional)
    join(e.members, " ") * "?"
end

type ZeroOrMore <: _Compound
    name
    members

    ZeroOrMore(name::String, members::Expression...) = new(name,members)
    ZeroOrMore(members::Expression...) = new("", members)
end

ZeroOrMore(members::Expression...; name::String="") = ZeroOrMore(name, members...)

function _uncached_match(self::ZeroOrMore, text::String, pos::Int, cache::Dict, err::ParseError)
    new_pos = pos
    first = true
    local children
    while true
        node = _match(self.members[1], text, new_pos, cache, err)
        if isempty(node) || textlength(node) == 0
            first && return Node(self.name, text, pos, prevind(text, new_pos))
            return Node(self.name, text, pos, _end(children[end]), tuple(children...))
        end
        if first
            children = MatchNode[node]
            first = false
        else
            push!(children, node)
        end
        new_pos += textlength(node)
    end
    return EmptyNode()
end

function _as_rhs(e::ZeroOrMore)
    join(e.members, " ") * "*"
end

type OneOrMore <: _Compound
    name::String
    _min::Int
    members

    # TODO: too many constructors -- members should always be first because it's the only thing that is always there
    OneOrMore(name::String="", _min::Int=1, members::Expression...=Expression[]) = new(name, _min, members)
    OneOrMore(members::Expression...=Expression[], name::String="", _min::Int=1) = new(name, _min, members)
end

OneOrMore(members::Expression...; _min::Int=1, name::String="") = OneOrMore(name, _min, members...)

function _uncached_match(self::OneOrMore, text::String, pos::Int, cache::Dict, err::ParseError)
    new_pos = pos
    first = true
    local children
    while true
        node = _match(self.members[1], text, new_pos, cache, err)
        if isempty(node)
            break
        end
        if first
            children = MatchNode[node]
            first = false
        else
            push!(children, node)
        end
        len = textlength(node)
        if len == 0
            break
        end
        new_pos += len
    end
    if !first && length(children) >= self._min
        return Node(self.name, text, pos, prevind(text, new_pos), tuple(children...))
    end
    return EmptyNode()
end

type LazyReference <: Expression
    name::String
    label::String
end

_as_rhs(lr::LazyReference) = "<LazyReference to $(lr.label)>"

unicode(lr::LazyReference) = lr.name

function _as_rhs(e::OneOrMore)
    join(e.members, " ") * "+"
end

function as_rule(expr::Expression)
    if expr.name == ""
        return _as_rhs(expr)
    end

    "$(expr.name) = $(_as_rhs(expr))"
end

function unicode(expr::Expression)
    class = string(typeof(expr))
    rulestr = as_rule(expr)
    id = hex(object_id(expr))
    "<$(class) $(rulestr) at 0x$(id)>"
end

indent(stack) = "  .  "^length(stack)
exprid(expr) = string(typeof(expr)) * " <" * expr.name * ">"
function show(io::IO, expr::Expression)
    function reshow(io::IO, expr::Expression, stack::Set)
        println(io, indent(stack), exprid(expr), ": ")
        for field in sort(names(expr))
            if(field == :name)
# already taken care of
            elseif(field == :members)
                push!(stack, exprid(expr))
                for m in expr.members
                    if in(exprid(m), stack)
                        println(io, indent(stack), "RECURSIVE " * exprid(expr) * " ...")
                    else
                        reshow(io, m, stack)
                    end
                end
                delete!(stack, exprid(expr))
            else
                println(io, indent(stack), ' ', string(field), '(', getfield(expr, field), ')')
            end
        end
    end
    reshow(io, expr, Set())
end


text:
"""
    function strength(
        team::Team,
        method::String
    )::Float64

Calculates the combined strength of a team.

Arguments
=========
- `team`:
    a team of players.
- `method`
    the method for calculating team strength
        current methods include:
            - "Simple Average"
                assigns strength to be average of individual player strengths.
            - "Weighted Average"
"""

function strength(team::Team, method::String)::Float64
    if method == "Simple Average"
        p1_strength = strength(team.player1)
        p2_strength = strength(team.player2)
        return (p1_strength + p2_strength) / 2.0
    else
        throw(DomainError(method, "this method isn't supported"))
    end
end

"""
    function probabilityOfWinning(
        team1::Team,
        team2::Team,
        method::String
    )::Float64

Calculates the probability of team1 beating team2.

Arguments
=========
- `team1`:
    the team for whose probability is calculated.
- `team2`
    the opposing team.
- `method`:
    method to determine the probability team 1 beats team 2.
    options:
        1) "BT" - Bradley Terry Model
        2) "BT_Exp" - Bradely Terry Exponential Model
"""

function probabilityOfWinning(
    team1::Team,
    team2::Team,
    method::String
)::Float64
    team1_strength = strength(team1, "Simple Average")
    team2_strength = strength(team2, "Simple Average")
    if method == "BT"
        return (team1_strength / (team1_strength + team2_strength))
    elseif method == "BT_Exp"
        return (exp(team1_strength) / (exp(team1_strength) + exp(team2_strength)))
    else
        throw(DomainError(method, "this method isn't supported"))
    end
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 2
coefs = Float64.([0.25; 0.125; 0.06])
dense_stencil = @def_stencil_expression(
            v*@sum(i,-2,2,
                @sum(j,-2,2,
                    @sum(k,-2,2, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))
st_def = CreateStencilDefinition(dense_stencil, coefs)
st_inst = NewStencilInstance(st_def)
## Input Data size Definition
radius = 2
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)
vsc = CreateData(dx,dy,dz, init=0.5, mid=1, rng=5)

t_steps = 1
## GPU Implementation
gpu_out = ApplyStencil(st_inst, data, t_steps, vsq=vsc)

## CPU result
cpu_coefs = fill(0.06, 5,5,5)
cpu_coefs[2:4, 2:4, 2:4] .= 0.125
cpu_coefs[3,3,3] = 0.25
cpu_out = cpu_dense_stencil_vsc(data, vsc, radius, cpu_coefs, t_steps)

## Check
mean_er = sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out)))
@assert(mean_er < 0.001)

## Visualize Error
using Plots

shift = 0
Int(round(dz/2))+shift
heatmap(gpu_out[:,:,Int(round(dz/2))+shift]-cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])


text:
#
println("Ci) Test of the DecayYield module with ASF from an internally generated multiplet.")
#
setDefaults("print summary: open", "zzz-DecayYield.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(12.); properties=[JAC.Yields],
                        grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600),
                        configs=[Configuration("1s 2s^2 2p^6")],
                        yieldSettings=DecayYield.Settings("SCA", true, false, Int64[]) )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
@testset "HO-SVD" begin
    r = 2
    T = rand_kruskal3(r, (10, 20, 30), false)

    @testset "no residuals calculation" begin
        @time factors = hosvd(T, r)
        @test length(factors.factors) == ndims(T)
        @test size.(factors.factors) == ntuple(i -> (size(T, i), r), ndims(T))
        @test size(factors.core) == (r, r, r)
        @test isnan(rel_residue(factors))
    end

    @testset "core reconstruction and residuals" begin
        @time factors = hosvd(T, r, compute_error=true)
        @test size(factors.core) == (r, r, r)
        @test rel_residue(factors) < 1e-5
    end

    @testset "core dimension equal to the original dimension" begin
        @time factors = hosvd(randn(10, 20, 30), (10, 15, 5), compute_error=true)
        @test size(factors.core) == (10, 15, 5)
    end

end


text:
"""
    RCMethod

Abstract type for nesting the various methods that can be used to perform a Reality Check. Subtypes are: \n
    RCBoot \n
The subtypes have entries in the help (?) menu.
"""
abstract type RCMethod end

"""
    RCBoot(data ; alpha::Float64=0.05, bootinput_kwargs...)

Method type for doing a Reality Check using a dependent bootstrap. A constructor
that takes the data and several keyword arguments is provided. The data refers
to the loss differences for each forecast relative to the basecase. Use ?rc for more detail.
Relevant keyword arguments are:
    alpha=0.05 <- The confidence level to use in the test \n
    kernelfunction=KernelEpanechnikov() <- The kernel function to use with HAC variance estimator. See ?hacvariance for more detail. \n
    bandwidth=-1 <- The bandwidth for HAC variance estimator. If bandwidth <= -1 then bandwidth is estimated using Politis (2003) "Adaptive Bandwidth Choice" \n
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
See the DependentBootstrap package docs for more info on bootstrap input keyword arguments.
"""
struct RCBoot <: RCMethod
    alpha::Float64
	bootinput::BootInput
    function RCBoot(alpha::Float64, bootinput::BootInput)
        !(0.0 < alpha < 0.5) && error("Confidence level set to $(alpha) which is not on the (0, 0.5) interval")
        new(alpha, bootinput)
    end
end
RCBoot(data ; alpha::Float64=0.05, kwargs...)::RCBoot = RCBoot(alpha, BootInput(data ; kwargs...))

"""
    RCTest(rejH0::Int, pvalue::Float64)

Output type from a Reality Check test.
The fields of this type follow: \n
    rejH0 <- true if the null is rejected, false otherwise
    pvalue <- p-value from the test
"""
struct RCTest
    rejH0::Bool
    pvalue::Float64
end
function Base.show(io::IO, a::RCBoot)
    println("Reality Check test results:")
    println("    Reject null hypothesis: $(a.rejH0)")
    println("    p-value: $(a.pvalue)")
end

"""
	rc(lD::Matrix{T}, method ; kwargs)
    rc(lD::Matrix{T} ; kwargs)

This function implements the test proposed in White (2000) "A Reality Check for Data Snooping" following the methodology in Hansen (2005). \n
Let x_0 denote a base-case forecast, x_k, k = 1, ..., K, denote K alternative forecasts, and y denote the forecast target.
Let L(., .) denote a loss function. The first argument of rc is a matrix where the kth column of the matrix is created by the operation: \n
L(x_k, y) - L(x_0, y) \n
Note that the forecast loss comes first and the base case loss comes second. This is the opposite to what is described in White's paper. \n
The second method argument determines which methodology to use. Currently, only RCBoot is available and if
this input type is provided, the keyword arguments are not needed.
Alternatively, the user can omit the second argument, and then any keyword arguments will be passed to the
RCBoot constructor. See ?RCBoot for more detail. The most relevant keyword arguments are:
    alpha=0.05 <- The confidence level to use in the test \n
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
The output of a Reality Check test is of type RCTest. Use ?RCTest for more information.
"""
function rc(lD::Matrix{<:Number}, method::RCBoot)::RCTest
    lD *= -1 #White's loss differentials have base case first
	numObs = size(lD, 1)
	numModel = size(lD, 2)
	numResample = method.bootinput.numresample
	numObs < 2 && error("Number of observations = $(numObs) which is not enough to perform a reality check")
	numModel < 1 && error("Input dataset is empty")
	inds = dbootinds(Float64[], method.bootinput) #Bootstrap indices
    #Get mean loss differentials and bootstrapped mean loss differentials
	mld = Float64[ mean(view(lD, 1:numObs, k)) for k = 1:numModel ]
	mldBoot = Array{Float64}(undef, numModel, numResample)
	for j = 1:numResample
		for k = 1:numModel
			mldBoot[k, j] = mean(view(lD, 1:numObs, k)[inds[j]])
		end
	end
	#Get RC test statistic and bootstrapped density under the null
	v = maximum(sqrt(numObs) * mld)
	vBoot = maximum(sqrt(numObs) * (mldBoot .- mld), dims=1)
	#Calculate p-value and return (as vector)
	pVal = sum(vBoot .> v) / numResample
    pVal < method.alpha ? (rejH0 = true) : (rejH0 = false)
	return RCTest(rejH0, pVal)
end
#Keyword method
rc(lD::Matrix{<:Number} ; kwargs...)::RCTest = rc(lD, RCBoot(lD ; kwargs...))
function rc(lD::Matrix{<:Number}, method ; kwargs...)::RCTest
    method == :boot && return rc(lD, RCBoot(lD ; kwargs...))
    error("Invalid method: $(method)")
end


text:
__precompile__()

"""
# Bootstrap package

Statistical bootstrapping for Julia
"""
module Bootstrap

using Compat

using StatsBase
using Distributions
using DataFrames

import Base: start, next, done, eltype, length

export
    Model,
    nvar,
    nrun,
    draw!,
    rademacher,
    mammen,
    iquantile,
    bootstrap,
    BasicSampling,
    AntitheticSampling,
    BalancedSampling,
    ExactSampling,
    ResidualSampling,
    WildSampling,
    BootstrapSample,
    NonParametricBootstrapSample,
    ParametricBootstrapSample,
    original,
    straps,
    data,
    model,
    statistic,
    sampling,
    noise,
    bias,
    se,
    ci,
    BasicConfInt,
    PercentileConfInt,
    NormalConfInt,
    BCaConfInt,
    StudentConfInt,
    level

include("stats.jl")
include("nobs.jl")
include("draw.jl")
include("bootsampling.jl")
include("get.jl")
include("show.jl")
include("confint.jl")
include("datasets/Datasets.jl")

end


text:
using NLopt, ForwardDiff, DiffResults, Optim , Sobol, MappedArrays

#here i will procede to define an VTD solver:

#tangent plane distance functions of a potential
#temperature is constant, and the volume is reduced in the case of VTD to
#density = n/v, so the potencial is _helmholtzd(d)

#in TPD,
_TPD(fn,x,y) = sum(ForwardDiff.gradient(fn,x) .* (x - y)) - fn(x) + fn(y)
_dTPDdy(fn,x,y) =   ForwardDiff.gradient(fn,y) .- ForwardDiff.gradient(fn,x)
_dTPDdy!(fn,F,dx,y) = begin F .=  ForwardDiff.gradient(fn,y) .- dx end
_st_twophases(fn,x,y,β) = β*fn(y) + (1-β)*fn((x .- β .* y)./ (1-β)) - fn(x)
_a_twophases(fn,x,y,β) = β*fn(y) + (1-β)*fn((x .- β .* y) ./ (1-β))
_alpha_phase_2(α,x0,x1) = (x0-α*x1)/(1-α)
#auxiliary funtion: nlopt form to input functions with gradient
function nlopt_form(f,x,g,diffresult_cache)
    if length(g) > 0
        ForwardDiff.gradient!(diffresult_cache,f,x)
        g .= DiffResults.gradient(diffresult_cache)
        return DiffResults.value(diffresult_cache)
    else
        return f(x)
    end
end

function autodiff_optim(f,x0)
    g! = (g,x) -> begin
        ForwardDiff.gradient!(g,f,x)
    end

    fg! = (g,x) -> begin
        ForwardDiff.gradient!(g,f,x)
        f(x)
    end
    h! = (H,x) -> begin
        ForwardDiff.hessian!(H,f,x)
    end
    return TwiceDifferentiable(f, g!, fg!, h!,x0)
end

#molar volume and mol_number
function stvtd(model::AbstractHelmholtzModel,v,T,x0;verbose=false,test_phase=nothing,local_iters=nothing)

    x = mol_fraction(x0,model)
    len = length(model)
    c0 =  x ./ v
    min_c = fill(1e-03,len)

    bi = covolumes(model)
    max_c = x./bi
    fn = c-> _helmholtzd(model,T,c)
    tpd = z -> _TPD(fn,c0,z)

    function _pressure(c) #pressure from molar concentrations
        v1 = inv(sum(c))
        x1 = c*v1
        return core_pressure(model,v1,T0,x1)
    end
    #NLopt part. the strategy is the following:
    # a relatively low ammount of local optimizations is made, starting from random points.
    #if those optimizations fail to find a feasible point, a global optimization is performed.
    #if the global optimization fail to find a feasible point, there is not phase split.
    diffresults_cache_st = DiffResults.GradientResult(x)
    opt_local = Opt(:LD_SLSQP, len) #sequencial quadratic programming
    opt_local.lower_bounds = min_c
    opt_local.upper_bounds = max_c
    opt_local.xtol_rel = 1e-4
    opt_local.min_objective = (x,grad) -> nlopt_form(tpd,x,grad,diffresults_cache_st)

    isnothing(local_iters) && (local_iters = 10*len)

    s = SobolSeq(min_c, max_c) #random distributed numbers with sobol sequences
    skip(s,local_iters) #recomended by the package
    crand = zeros(len)
    one_phase = true
    minf = Inf
    new_phase = zeros(len)

    for i in 1:local_iters
        next!(s,crand) #new random point
        minf,new_phase,ret = NLopt.optimize(opt_local, crand)

        length(intersect(new_phase,min_c,max_c)) == 0 && (one_phase = false)
        minf <0.0                                     && (one_phase = false)
        (abs(minf)>sqrt(eps(Float64)))                && (one_phase = false)
        _pressure(new_phase) > 0.0                    && (one_phase = false)

        if one_phase == false
            if verbose == true
                println("Candidate phase found after $i local optimizations")
                println("Phase(moles/m3):: $new_phase")
                println("Potential Energy:: ",fn(new_phase))
            end
            break
        end
    end

    if one_phase == true
        if verbose == true
            println("Local optimizations failed to find a candidate phase,trying global optimization.")
        end
        opt_global = Opt(:GN_DIRECT_L, len) #sequencial quadratic programming
        opt_global.lower_bounds = min_c
        opt_global.upper_bounds = max_c
        opt_global.xtol_rel = 5e-5
        opt_global.min_objective = (x,grad) -> nlopt_form(tpd,x,grad,diffresults_cache_st)
        next!(s,crand)
        minf,new_phase,ret = NLopt.optimize(opt_local, crand)
        if length(intersect(new_phase,min_c,max_c)) == 0 &&
             minf <0.0 &&
             (abs(minf)>sqrt(eps(Float64))) &&
             _pressure(new_phase) > 0.0
            one_phase = false
        end
    end

    if one_phase == false
        if verbose == true
            println("Candidate phase found in global optimization")
            println("Phase(moles/m3):: $new_phase")
        end
    else
        if verbose == true
            println("Global optimization failed to find a candidate phase,returning original phase.")
        end
        return helmholtz_phase(model,v,T,x) # the material is definitely one phase
    end

    βmin = [0.0]
    βmax = [min(1.0,minimum(c0./new_phase))]

    opt_eq0 = Opt(:GN_DIRECT_L, 1) #Global minimization in alpha
    opt_eq0.lower_bounds = βmin
    opt_eq0.upper_bounds = βmax
    opt_eq0.xtol_rel = 5e-2
    β0 = 0.5*(βmin + βmax)/2
    diffresults_cache_eq0 = DiffResults.GradientResult([1])
    fn_eq0 = β -> _st_twophases(fn,c0,new_phase,β[1])
    opt_eq0.min_objective = (x,grad) -> nlopt_form(fn_eq0,x,grad,diffresults_cache_eq0)

    minf_eq0,minβ,ret = NLopt.optimize(opt_eq0, β0)

    if !((minf_eq0<0.0) &&
        (!(minβ in [βmin,βmax,β0])) &&
        (ret != :FORCED_STOP) &&
        (abs(minf_eq0)>sqrt(eps(Float64))))

        if verbose == true
            println("Candidate phase does not meet the criteria to be an equilibrium phase")
            println("Phase(moles/m3):: $new_phase")
        end
        return helmholtz_phase(model,v,T,x)
    else
        if verbose == true
            println("Candidate phase meets the criteria to be an equilibrium phase")
            println("Proceeding to phase equilibria")
            println("Phase(moles/m3):: $new_phase")
            println("fraction:: ",minβ[1])
        end
    end

    #println(c0)
    #println(minβ[1]*new_phase)
    #println((mol_fraction(n2),mol_fraction(n1),minβ[1]))
    #println(_helmholtzn(model,minβ[1]v,T,minβ[1]*new_phase*v)+_helmholtzn(model,(1-minβ[1])*v,T,x0.-minβ[1]*new_phase*v))
    #return (minf,minβ[1],minβ[1]*new_phase)

    #here ends the stability test and the search of a candidade phase. now we proceed to calculate
    #the real phase split
    equilibrium_phase0 = [new_phase...,minβ[1]]

    function fn_eq(z)
        α = z[end];
        c1 = view(z,1:(length(z)-1))
        fc2(x0,x1) = _alpha_phase_2(α,x0,x1)
        c2 = mappedarray(fc2, c0, c1)
        return α*fn(c1) + (1-α)*fn(c2)
    end

    function alpha_max(z)
        α = z[end];
        c1 = view(z,1:(length(z)-1))
        return α - minimum(c0./c1)
    end

    function covolume_max(z)
        #c1 is restricted via box bounds
        α = z[end];
        c1 = view(z,1:(length(z)-1))
        fc2(x0,x1) = _alpha_phase_2(α,x0,x1)
        c2 = mappedarray(fc2, c0, c1)
        c2[1]*bi[1]+c2[2]*bi[2] - 1

    end

    min_eq = vcat(fill(0.0,len),sqrt(eps(Float64)))
    max_eq = vcat(fill(2.0^256,len),1.0-sqrt(eps(Float64)))
    diffresults_cache_eq = DiffResults.GradientResult(min_eq)
    diffresults_cache_const_eq = DiffResults.GradientResult(min_eq)

    opt_eq = Opt(:LD_SLSQP, len+1) #sequencial quadratic programming again, but with stricter tolerance
    opt_eq.lower_bounds = min_eq
    opt_eq.upper_bounds = max_eq
    opt_eq.ftol_abs = 1e-16
    objective_eq = (x,grad) -> nlopt_form(fn_eq,x,grad,diffresults_cache_eq)

    grad = zeros(len+1)

    exe = objective_eq(equilibrium_phase0,grad)
    println("c0:: ",equilibrium_phase0)
    println("function:: grad: ",grad,", f0: ",exe)
    println(minβ[1])

    opt_eq.min_objective = objective_eq
    inequality_constraint!(opt_eq,
        (x,grad) -> nlopt_form(alpha_max,x,grad,diffresults_cache_const_eq)
    )

    #phase equilibria
    minf_eq,equilibrium_phase,ret = NLopt.optimize(opt_eq,equilibrium_phase0)
    phase_success = true

    function phase_to_helmholtz(c1,α)
        v1 = inv(sum(c1))
        n1 = mol_number(c1*v1*α)
        return helmholtz_phase(model,v1,T0,n1)
    end

    αfinal = equilibrium_phase[end]
    c1 = @view equilibrium_phase[1:end-1]
    fc2(x0,x1) = _alpha_phase_2(αfinal,x0,x1)
    c2 = mappedarray(fc2, c0, c1)

    v1 = inv(sum(c1))
    v2 = inv(sum(c2))
    isapprox(v1,v2;rtol = 1e-6) && (phase_success = false)
    ret != :FTOL_REACHED && (phase_success = false)
    #if the phase equilibria failed, return the original phase.
    phase_success = false && return helmholtz_phase(model,v,T,x)

    phase1  = phase_to_helmholtz(c1,αfinal)
    phase2  = phase_to_helmholtz(c2,1-αfinal)
    return [phase1,phase2]
end

function autodiff(f)
    res =(x,grad=[]) ->begin
        if length(grad) > 0
            diffresult_cache = DiffResults.GradientResult(x)
            ForwardDiff.gradient!(diffresult_cache,f,x)
            grad .= DiffResults.gradient(diffresult_cache)
            return DiffResults.value(diffresult_cache)
        else
            return f(x)
        end
    end
    return res
    end

function _pressure(res)
    x1 = res[1:2]
    v1 = inv(sum(res))
    x1 = x1*v1
    return pressure(minigerg,v1,T0,x1)
end


text:
using RobustAndOptimalControl, ControlSystems

G = tf(200, [10, 1])*tf(1, [0.05, 1])^2     |> ss
Gd = tf(100, [10, 1])                       |> ss
W1 = tf([1, 2], [1, 1e-6])                  |> ss
Gs = G*W1
Ks, γ, info = glover_mcfarlane(Gs, 1.1)
@test info.γmin ≈ 2.34 atol=0.005

if isinteractive()
    bodeplot([G, Gs, Gs*Ks]) |> display

    plot( step(Gd*feedback(1, G*W1), 3))
    plot!(step(Gd*feedback(1, G*W1*Ks), 3)) |> display

    nyquistplot([G*W1, G*W1*Ks], ylims=(-2,1), xlims=(-2, 1), Ms_circles=1.5) |> display
end

W1h = hanus(W1)
@test W1h.nu == 2W1.nu

##

P = tf([1, 5], [1, 2, 10])
W1 = tf(1,[1, 0])

K,γ,info = glover_mcfarlane(ss(P); W1)

Gcl0 = let
    tempA = [-2.0 -2.5 2.0 0.0 0.0 0.0; 4.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -0.7189 0.0453 -2.9175; 0.0493 0.0616 0.0 -2.0493 -2.5616 2.0; 0.399 0.4988 0.0 3.601 -0.4988 0.0; 0.5 0.625 0.0 -1.2189 -0.5797 -2.9175]
    tempB = [0.0 2.0; 0.0 0.0; 0.0 0.0; 0.0985 0.0; 0.7981 0.0; 1.0 0.0]
    tempC = [0.5 0.625 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0]
    tempD = [1.0 0.0; 0.0 0.0]
    ss(tempA, tempB, tempC, tempD)
end

K0 = let
    tempA = [0.0 -0.7189 0.0453 -2.9175; 0.0 -2.0493 -2.5616 2.0; 0.0 3.601 -0.4988 0.0; 0.0 -1.2189 -0.5797 -2.9175]
    tempB = [0.0; 0.0985; 0.7981; 1.0;;]
    tempC = [1.0 0.0 0.0 0.0]
    tempD = [0.0;;]
    ss(tempA, tempB, tempC, tempD)
end

Gcl = extended_gangoffour(P, K)
@test hinfnorm2(Gcl-Gcl0)[1] < 1e-4
@test hinfnorm2(info.Gcl-Gcl0)[1] < 1e-4
@test nugap(K, -K0)[1] < 1e-4
@test info.margin ≈ 0.6325 atol=1e-3

@test ncfmargin(P, K)[1] ≈ 0.4472 atol=1e-3

p = ss(tf(4, [1, -0.001]))
cL = 1
cH = 10
@test ncfmargin(p,cL)[1] ≈ 0.7069 atol=1e-3

P = ss(tf(1, [1,1]))
K = ss(1)
S, PS, CS, T = RobustAndOptimalControl.gangoffour2(P,K)
# gangoffourplot(P, K)
# bodeplot!(extended_gangoffour(P, K), plotphase=false)

isstable.((S, PS, CS, T))
gof = extended_gangoffour(P, K)
isstable(gof)
@test nugap(S, gof[1,1])[1] < 1e-6
@test nugap(PS, gof[1,2])[1] < 1e-6
@test nugap(CS, -gof[2,1])[1] < 1e-6 # NOTE: slightly disturbing to have - here
@test nugap(T, -gof[2,2])[1] < 1e-6


text:
function update!(cur_reactor::AbstractReactor)

  if !cur_reactor.is_good
    _NanizeReactor!(cur_reactor)
    return cur_reactor
  end

  cur_reactor.sigma_v = calc_sigma_v(cur_reactor)

  secondary_constraint = getfield(Fussy, Symbol("$(cur_reactor.constraint)_equation"))(cur_reactor)

  cur_equation_set = EquationSet(cur_reactor, secondary_constraint)

  try

    cur_reactor.B_0 = convert(Real, cur_equation_set.B_0)
    cur_reactor.R_0 = convert(Real, cur_equation_set.R_0)
    cur_reactor.n_bar = convert(Real, calc_n_bar(cur_reactor))
    cur_reactor.tau_E = tau_E(cur_reactor)

  catch cur_error

    is_caught_error = isa(cur_error, InexactError)
    is_caught_error |= isa(cur_error, DomainError)
    is_caught_error || rethrow(cur_error)

    cur_reactor.is_good = false
    _NanizeReactor!(cur_reactor)

    return cur_reactor

  end

  cur_reactor.p_bar = p_bar(cur_reactor)
  cur_reactor.P_F = P_F(cur_reactor)

  cur_reactor.f_BS = f_BS(cur_reactor)
  cur_reactor.f_CD = f_CD(cur_reactor)
  cur_reactor.f_IN = f_IN(cur_reactor)

  cur_reactor.beta_N = beta_N(cur_reactor)
  cur_reactor.q_95 = q_95(cur_reactor)
  cur_reactor.q_DV = q_DV(cur_reactor)
  cur_reactor.P_E = P_E(cur_reactor)
  cur_reactor.P_W = P_W(cur_reactor)

  cur_reactor.norm_beta_N = beta_N(cur_reactor)
  cur_reactor.norm_q_95 = q_95(cur_reactor)
  cur_reactor.norm_q_DV = q_DV(cur_reactor)
  cur_reactor.norm_P_E = P_E(cur_reactor)
  cur_reactor.norm_P_W = P_W(cur_reactor)

  cur_reactor.norm_beta_N /= cur_reactor.max_beta_N
  cur_reactor.norm_q_95 /= cur_reactor.max_q_95
  cur_reactor.norm_q_DV /= cur_reactor.max_q_DV
  cur_reactor.norm_P_E /= cur_reactor.max_P_E
  cur_reactor.norm_P_W /= cur_reactor.max_P_W

  cur_reactor.norm_q_95 ^= -1

  cur_reactor.is_valid = true

  for (cur_limit, cur_variable) in secondary_params
    in(cur_limit, cur_reactor.ignored_limits) && continue

    cur_norm_symbol = Symbol("norm_$(cur_variable)")

    cur_norm_value = getfield(cur_reactor, cur_norm_symbol)

    cur_norm_value <= 1.001 && continue

    cur_reactor.is_valid = false
    break
  end

  cur_reactor.volume = volume(cur_reactor)
  cur_reactor.W_M = W_M(cur_reactor)

  cur_reactor.cost = gray_cost(cur_reactor)

  cur_reactor.resistance = R_P(cur_reactor)
  cur_reactor.voltage = V_loop(cur_reactor)
  cur_reactor.inductance = L_P(cur_reactor)

  cur_reactor.a = a(cur_reactor)
  cur_reactor.b = b(cur_reactor)

  try

    @assert !isnan(cur_reactor.f_BS)

    cur_reactor.c = c(cur_reactor)
    cur_reactor.d = d(cur_reactor)

    cur_reactor.R_CS = R_CS(cur_reactor)
    cur_reactor.h_CS = h_CS(cur_reactor)

    @assert cur_reactor.c >= 0
    @assert cur_reactor.d >= 0

    @assert cur_reactor.R_CS >= 0
    @assert cur_reactor.h_CS >= 0

  catch cur_error

    is_caught_error = isa(cur_error, InexactError)
    is_caught_error |= isa(cur_error, DomainError)
    is_caught_error |= isa(cur_error, AssertionError)

    is_caught_error || rethrow(cur_error)

    cur_reactor.is_valid = false
    return cur_reactor

  end

  cur_reactor

end


text:
module TheDataMustFlow

using MacroTools

abstract type DataRemote end

abstract type AbstractMetaData end
abstract type AbstractDataSet end
abstract type AbstractInput end
abstract type AbstractProblem end

abstract type ProjectTag end
abstract type Tag end
abstract type TableTag <: Tag end

include("dataset.jl")

defaultconfig!(::ProjectTag, d::AbstractDict) = d

function config(::Type{Dict}, ::ProjectTag; kwargs...)
    d = isempty(kwargs) ? Dict{Symbol,Any}() : convert(Dict{Symbol,Any}, Dict(kwargs))
    defaultconfig!(tag, d)
end
config(tag::ProjectTag; kwargs...) = config(Dict, tag; kwargs...)

export ProjectTag, Tag, TableTag
export directory!
export @tag, @tabletag
export alltags, iscomplete, initialize!, initialize
export load, load!, save

end # module


text:
using Documenter, SCEDC

makedocs(;
    modules=[SCEDC],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tclements/SCEDC.jl/blob/{commit}{path}#L{line}",
    sitename="SCEDC.jl",
    authors="Tim Clements <thclements@g.harvard.edu>",
    assets=String[],
)

deploydocs(;
    repo="github.com/tclements/SCEDC.jl",
)


text:
using DetectionIoTools
using Test

@testset "DetectionIoTools.jl" begin
    @testset "Save interesting" begin
        @info "Testing Save interesting"
        allsounds = [randn(100) for _ in 1:10]

        sounddir = save_interesting(allsounds, [1,3])
        @test readdir(sounddir) |> length == 3

    end
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

a = 1.4; b = 0.3
he = Systems.henon(; a, b)
fig = figure()

x0, y0 = columns(trajectory(he, 5000, Ttr = 100))
x1 = x0
y1 = @. 1 - a*x0^2 + y0

x2 = b .* x1
y2 = y1
y3 = x2
x3 = y2

ax1 = subplot(141)
ax1.scatter(x0, y0, color = COLORS[1], s = 1, label = "initial")
ax1.text(0.1, 0.05, "initial", transform=ax1.transAxes)

αprev = 0.75
ax2 = subplot(142)
ax2.scatter(x0, y0, color = COLORS[1], s = 0.25, alpha = αprev)
ax2.scatter(x1, y1, color = COLORS[2], s = 1, label = "stretch/bend")
ax2.text(0.1, 0.05, "bend/stretch", transform=ax2.transAxes)
ax2.arrow(0,-0.19, 0, 0.8, head_width=0.1, head_length=0.1, linewidth=4, color=COLORS[5], length_includes_head=true)
ax2.arrow(-1.5, 0.0,  0, -0.5, head_width=0.1, head_length=0.1, linewidth=4, color=COLORS[5], length_includes_head=true)
ax2.arrow(1.5, 0.0,  0, -0.5, head_width=0.1, head_length=0.1, linewidth=4, color=COLORS[5], length_includes_head=true)

ax3 = subplot(143)
ax3.scatter(x1, y1, color = COLORS[2], s = 0.25, alpha = αprev)
ax3.scatter(x2, y2, color = COLORS[3], s = 1)
ax3.text(0.1, 0.05, "squeeze", transform=ax3.transAxes)
ax3.arrow(-1.5, 0.0,  1, 0, head_width=0.1, head_length=0.1, linewidth=4, color=COLORS[5], length_includes_head=true)
ax3.arrow(1.5, 0.0,  -1, 0, head_width=0.1, head_length=0.1, linewidth=4, color=COLORS[5], length_includes_head=true)

ax4 = subplot(144)
ax4.scatter(x2, y2, color = COLORS[3], s = 0.25, alpha = αprev)
ax4.scatter(x3, y3, color = COLORS[4], s = 1, label = "reflect")
ax4.text(0.1, 0.05, "reflect", transform=ax4.transAxes)

ax4.plot([-2,2], [-2,2], ls = "dashed", lw = 2, color = COLORS[5])

# ticks and stuff
ax1.set_ylabel("\$y\$")
for ax in (ax1, ax2, ax3, ax4)
    ax.set_xlim(-2, 2)
    ax.set_ylim(-1.5, 1.5)
    ax.set_yticks(-1.5:0.5:1.5)
    ax.set_xticks(-1.5:1:1.5)
    ax.grid("on")
    ax.set_yticklabels([])
    ax.set_xticklabels([])
    ax.set_xlabel("\$x\$", labelpad = -10)
    # ax.legend(markerscale = 4, scatterpoints=1)
end
tight_layout(pad = 0.33)
add_identifiers!(fig; xloc = 0)
wsave(plotsdir("3", "stretchinghenon"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module TestNamefinders

    using Test
    using NCBITaxonomy

    dipl = descendants(ncbi"Diplectanidae")

    lam = taxon(namefilter(dipl), "Lamellodiscus elegans")

    @test typeof(lam) == NCBITaxon
    @test lam.name == "Lamellodiscus elegans"

    dipl_2 = namefilter(descendants(ncbi"Diplectanidae"))

    @test !isnothing(taxon(dipl_2, "Lamellodiscus elegans"))
    @test_throws NameHasNoDirectMatch taxon(dipl_2, "Gallus gallus")

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# helper function for table view
function viewtable(table, rows, cols)
  t = Tables.columns(table)
  v = map(cols) do c
    col = Tables.getcolumn(t, c)
    c => view(col, rows)
  end
  (; v...)
end

"""
    AbstractData

An abstract type to aid with custom spatial data types. Users can
subtype their domains from this type, and implement the methods in
`geotraits/data.jl`.
"""
abstract type AbstractData end

"""
    sdata₁ == sdata₂

Tells whether or not the spatial `sdata₁` and `sdata₂` are equal.
"""
==(sdata₁::AbstractData, sdata₂::AbstractData) =
  domain(sdata₁) == domain(sdata₂) &&
  values(sdata₁) == values(sdata₂)

# -----------
# TABLES API
# -----------

Tables.istable(::Type{<:AbstractData}) = true
Tables.materializer(sdata::AbstractData) = Tables.materializer(values(sdata))
Tables.columnaccess(sdata::AbstractData) = Tables.columnaccess(values(sdata))
Tables.rowaccess(sdata::AbstractData) = Tables.rowaccess(values(sdata))
Tables.schema(sdata::AbstractData) = Tables.schema(values(sdata))
Tables.columns(sdata::AbstractData) = Tables.columns(values(sdata))
Tables.columnnames(sdata::AbstractData) = Tables.columnnames(values(sdata))
Tables.getcolumn(sdata::AbstractData, c::Symbol) = Tables.getcolumn(values(sdata), c)
Tables.rows(sdata::AbstractData) = Tables.rows(values(sdata))

# -------------
# VARIABLE API
# -------------

function variables(sdata::AbstractData)
  s = Tables.schema(sdata)
  ns, ts = s.names, s.types
  @. Variable(ns, nonmissing(ts))
end

Base.getindex(sdata::AbstractData, var::Symbol) =
  Tables.getcolumn(sdata, var)

# ---------
# VIEW API
# ---------

Base.view(sdata::AbstractData, inds::AbstractVector{Int}) =
  DataView(sdata, inds, collect(name.(variables(sdata))))
Base.view(sdata::AbstractData, vars::AbstractVector{Symbol}) =
  DataView(sdata, 1:nelms(sdata), vars)
Base.view(sdata::AbstractData, inds, vars) =
  DataView(sdata, inds, vars)

# ------------
# IO methods
# ------------
function Base.show(io::IO, sdata::AbstractData)
  N = ncoords(sdata)
  T = coordtype(sdata)
  n = nelms(sdata)
  print(io, "$n SpatialData{$T,$N}")
end

MIMES = [MIME"text/plain", MIME"text/html"]

for MIME in MIMES
  @eval function Base.show(io::IO, ::$MIME, sdata::AbstractData)
    𝒟 = domain(sdata)
    𝒯 = values(sdata)
    s = Tables.schema(𝒯)
    vars = zip(s.names, s.types)
    println(io, 𝒟)
    println(io, "  variables")
    varlines = ["    └─$var ($V)" for (var,V) in vars]
    print(io, join(sort(varlines), "\n"))
  end
end

#------------------
# IMPLEMENTATIONS
#------------------
"""
    SpatialData(domain, data)

Tabular `data` georeferenced in a spatial `domain`.
"""
struct SpatialData{𝒟,𝒯} <: AbstractData
  domain::𝒟
  table::𝒯

  function SpatialData{𝒟,𝒯}(domain, table) where {𝒟,𝒯}
    ne = nelms(domain)
    nr = length(Tables.rows(table))
    @assert ne == nr "number of table rows ≠ number of mesh elements"
    new(domain, table)
  end
end

SpatialData(domain::𝒟, table::𝒯) where {𝒟,𝒯} =
  SpatialData{𝒟,𝒯}(domain, table)

domain(sdata::SpatialData) = sdata.domain
values(sdata::SpatialData) = sdata.table

"""
    DataView(sdata, inds, vars)

Return a view of spatial data `sdata` at `inds` and `vars`.
"""
struct DataView{𝒮,I,V} <: AbstractData
  data::𝒮
  inds::I
  vars::V
end

domain(dv::DataView) = view(domain(dv.data), dv.inds)
values(dv::DataView) = viewtable(values(dv.data), dv.inds, dv.vars)

coordinates!(buff::AbstractVector, dv::DataView, ind::Int) =
  coordinates!(buff, dv.data, dv.inds[ind])

==(dv₁::DataView, dv₂::DataView) =
  dv₁.data == dv₂.data &&
  dv₁.inds == dv₂.inds &&
  dv₁.vars == dv₂.vars

Base.view(dv::DataView, inds::AbstractVector{Int}) =
  DataView(dv.data, dv.inds[inds], dv.vars)
Base.view(dv::DataView, vars::AbstractVector{Symbol}) =
  DataView(dv.data, dv.inds, vars)
Base.view(dv::DataView, inds, vars) =
  DataView(dv.data, dv.inds[inds], vars)

"""
    collect(dataview)

Materialize spatial `dataview` into a new block of memory.
"""
Base.collect(dv::DataView) = georef(values(dv), coordinates(dv))

function Base.show(io::IO, dv::DataView)
  N = ncoords(dv)
  T = coordtype(dv)
  n = nelms(dv)
  print(io, "$n DataView{$T,$N}")
end