text:
# You should set up a Git repo with an example Julia package for the integration tests.
# You should create a bot user GitHub account, and give the bot user GitHub account push access to the testing repo.
# You should generate a personal access token (PAT) for the bot user GitHub account.
# You should add the PAT as an encrypted Travis CI environment variable named `BCBI_TEST_USER_GITHUB_TOKEN`.
# You should generate an SSH deploy key. Upload the public key as a deploy key for the testing repo.
# You should make the private key an encrypted Travis CI environment variable named `COMPATHELPER_PRIV`. Probably you will want to Base64-encode it.
const GLOBAL_PR_TITLE_PREFIX = Random.randstring(8)
const GITHUB = "GITHUB"
const GITLAB = "GITLAB"

function run_integration_tests(
    url::AbstractString, env::AbstractDict, ci_cfg::CompatHelper.CIService
)
    @testset "master_1" begin
        with_master_branch(templates("master_1"), url, "master") do master_1
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-1c] ",
                    master_branch=master_1,
                    entry_type=KeepEntry(),
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_2" begin
        with_master_branch(templates("master_2"), url, "master") do master_2
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-2c] ",
                    master_branch=master_2,
                    entry_type=KeepEntry(),
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_3" begin
        with_master_branch(templates("master_3"), url, "master") do master_3
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-3a] ",
                    master_branch=master_3,
                    entry_type=DropEntry(),
                )

                sleep(1)  # Prevent hitting the GH Secondary Rate Limits

                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-3b] ",
                    master_branch=master_3,
                    entry_type=KeepEntry(),
                )

                sleep(1)  # Prevent hitting the GH Secondary Rate Limits

                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-3c] ",
                    master_branch=master_3,
                    entry_type=KeepEntry(),
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_4" begin
        with_master_branch(templates("master_4"), url, "master") do master_4
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-4c] ",
                    master_branch=master_4,
                    entry_type=KeepEntry(),
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_5" begin
        with_master_branch(templates("master_5"), url, "master") do master_5
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-5a] ",
                    master_branch=master_5,
                    entry_type=DropEntry(),
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_6" begin
        with_master_branch(templates("master_6"), url, "master") do master_6
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-6c] ",
                    master_branch=master_6,
                    entry_type=KeepEntry(),
                    subdirs=["subdir_1", "subdir_2"],
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_7" begin
        with_master_branch(templates("master_7"), url, "master") do master_7
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-7a] ",
                    master_branch=master_7,
                    entry_type=DropEntry(),
                    bump_compat_containing_equality_specifier=false,
                )

                sleep(1)  # Prevent hitting the GH Secondary Rate Limits

                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-7b] ",
                    master_branch=master_7,
                    entry_type=DropEntry(),
                    bump_compat_containing_equality_specifier=true,
                )
            end
        end
    end

    sleep(1)  # Prevent hitting the GH Secondary Rate Limits

    @testset "master_8" begin
        with_master_branch(templates("master_8"), url, "master") do master_8
            withenv(env...) do
                CompatHelper.main(
                    ENV,
                    ci_cfg;
                    pr_title_prefix="$(GLOBAL_PR_TITLE_PREFIX) [test-8a] ",
                    master_branch=master_8,
                    entry_type=DropEntry(),
                    use_existing_registries=true,
                )
            end
        end
    end

    return _cleanup_old_branches(url)
end

@testset "$(service)" for service in [GITHUB, GITLAB]
    personal_access_token = ENV["INTEGRATION_PAT_$(service)"]
    test_repo = ENV["INTEGRATION_TEST_REPO_$(service)"]

    env = Dict(
        "$(service)_REPOSITORY" => test_repo,
        "$(service)_TOKEN" => personal_access_token,
        "CI_PROJECT_PATH" => ENV["INTEGRATION_TEST_REPO_GITLAB"],
    )

    # Otherwise auto_detect_ci_service() will think we're testing on GitHub
    if service == GITLAB
        delete!(env, "GITHUB_REPOSITORY")
        env["GITLAB_CI"] = "true"
    end

    ci_cfg = CompatHelper.auto_detect_ci_service(; env=env)
    api, repo = CompatHelper.get_api_and_repo(ci_cfg; env=env)
    url = CompatHelper.get_url_with_auth(api, ci_cfg, repo)

    run_integration_tests(url, env, ci_cfg)
end


text:
## OUTPUT FUNCTION
"""
# Duty Cycle
`set_duty(obj::T;ch="1",func="SQUare",duty=0.01) where {T<:DG1000Z}`
- func = {"SQUare", "PULSe"})
"""
set_duty(obj::T;ch="1",func="SQUare",duty=50) where {T<:DG1000Z} = write(obj,"SOURce$ch:FUNC:$func:DCYC $duty")

## BURST
"""
# Burst on
`set_burst_stat(obj::T;ch="1",st="off") where {T<:DG1000Z}`
- st = {"on", "off"}
"""
set_burst_stat(obj::T;ch="1",st="off") where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:STAT $st")

"""
# Burst mode
`set_burst_mode(obj::T;ch="1",md="triggered") where {T<:DG1000Z}`
- md = {"TRIGgered", "INFinity", "GATed"}
"""
set_burst_mode(obj::T;ch="1",md="triggered") where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:MODE $md")

"""
# Burst polarity
`set_burst_pol(obj::T;ch="1",pol="normal") where {T<:DG1000Z}`
- pol = {"NORMal", "INVerted"}
"""
set_burst_pol(obj::T;ch="1",pol="normal") where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:GATE:POL $pol")

"""
# Burst ncycles
`set_burst_ncyc(obj::T;ch="1",nc=1) where {T<:DG1000Z}`
- nc = 1 to 1'000'000 (external or manual trigger) 1 to 500'000 (internal trigger)
"""
set_burst_ncyc(obj::T;ch="1",nc=1) where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:NCYC $nc")

"""
# Burst period
`set_burst_per(obj::T;ch="1",per=0.01) where {T<:DG1000Z}`
- per = 2.0166us to 500s (in s)
"""
set_burst_per(obj::T;ch="1",per=0.01) where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:INT:PER $per")

"""
# Burst phase
`set_burst_phase(obj::T;ch="1",unit="deg",phase=0) where {T<:DG1000Z}``
- unit: *unused* on DG1000Z
- phase = 0 to 360
"""
set_burst_phase(obj::T;ch="1",unit="deg",phase=0) where {T<:DG1000Z} = write(obj,"SOURce$ch:BURS:PHAS $phase")

ext: jl


text:
using Pathnames
using Test
using Pathnames: PathnameWrapper as p

@testset "Pathnames.jl" begin
    home = expanduser(p"~")
    @test home isa Pathname
    @test home == pathname(ENV["HOME"])

    root = p"/"
    usr = root/"usr"

    @test root isa Pathname
    @test usr isa Pathname
    @test root/usr isa Pathname
    @test root/usr == usr
    @test root/usr/"local" isa Pathname
    @test pathname(root, usr, "local") == p"/usr/local"
    @test pathname(root, usr, p"local", "bin") == p"/usr/local/bin"

    testpath = p.mktempdir()
    cd(testpath)

    @test testpath isa Pathname
    @test realpath(testpath) isa Pathname
    @test p.pwd() == realpath(testpath)

    testpath = realpath(testpath)
    foofile = testpath/"foo"

    @test dirname(foofile) == testpath
    @test abspath(p"foo") == foofile
    @test touch(foofile) isa Pathname
    @test isfile(foofile)

    barfile = testpath/"bar"
    symlink(foofile, testpath/"bar")
    @test islink(barfile)
    @test readlink(barfile) == foofile
    rm(barfile)

    data = "Hello, world!"
    open(foofile, "w") do io
        write(io, data)
    end

    cp(foofile, barfile)
    @test isfile(barfile)
    open(barfile, "r") do io
        @test read(io, String) == data
    end

    bazfile = testpath/"baz"
    mv(barfile, bazfile)
    @test !isfile(barfile)
    @test isfile(bazfile)
    open(bazfile, "r") do io
        @test read(io, String) == data
    end
end


text:
using DiskJockey.EnsembleSampler
using Distributed
using LinearAlgebra: logdet

# The EnsembleSampler was written to *only* work in parallel, so in order to test it properly, we
# need to also define all of our test lnprobs to work using Distributed,
# even though they are quick enough to work on their own.

# Rosenbrock function
@everywhere function lnprob(p)
    x, y = p

    return -((1.0 - x)^2 + 100 * (y - x^2)^2)
end

@testset "sample Rosenbrock" begin
    nwalkers = 40
    ndim = 2

    sampler = Sampler(nwalkers, ndim, lnprob)

    # pos0 is the starting position, it needs to be a (ndim, nwalkers array)
    pos0 = Array{Float64}(undef, ndim, nwalkers)
    for i = 1:nwalkers
        pos0[:,i] = randn(ndim)
    end

    run_schedule(sampler, pos0, 10, 1, outputdir)
end

# Now, let's try sampling a function with more parameters than walkers (just for testing purposes)
# A 4D Gaussian to check against
@everywhere function lnprob(p)

    mu = [1.0, -2.0, 0.0, 3.4]
    Sigma = [[0.5^2,  0.00, 0.0, 0.0] [0.00, 0.3^2, 0.0, 0.0] [0.0, 0.0, 0.7^2, 0.0] [0.0, 0.0, 0.0, 1.0^2] ]

    R = p - mu # Residual vector
    ld = logdet(Sigma)

    lnp = -0.5 * ((R' * (Sigma \ R))[1] + ld + 2 * log(2 * pi))
    return lnp

end

@testset "sample Gaussian" begin
    nwalkers = 4
    ndim = 4

    sampler = Sampler(nwalkers, ndim, lnprob, true)

    # pos0 is the starting position, it needs to be a (ndim, nwalkers array)
    pos0 = Array{Float64}(undef, ndim, nwalkers)
    for i = 1:nwalkers
        pos0[:,i] = randn(ndim)
    end

    run_schedule(sampler, pos0, 1, 1, outputdir)
end


text:
#=
  ClickGraph example from PSI I paper, page 4

  Where the objective is to determine how similar the
  clicks of A and B are.

  From the PSI model
    // Orig: E[simAll]: 14475/23569 (~0.6141)
    clicksA := [1,1,1,0,0];
    clicksB := [1,1,1,0,0];

    // Variant: E[simAll]: 10581/38342 (~0.2759636)
    clicksA := [1,1,1,0,1,1,1,0];
    clicksB := [0,0,0,1,0,0,0,1];

   This Turing.jl model:
   Variant 1 ("very similar")  : simAll    0.6159
   Variant 2 ("very different"): simAll    0.2389

  cf ~/psi/clickGraph_psi1.psi
     ~/blog/clickGraph_psi1.blog
     ~/webppl/clickGraph.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function clickGraph(clickAs,clickBs)
    n = length(clickAs)

    simAll ~ Uniform(0,1)

    sim ~ filldist(flip(simAll),n)
    pA ~ filldist(Uniform(0,1),n)
    pB = tzeros(n)
    for i in 1:n
        if sim[i] == true
            pB[i] ~ Dirac(pA[i])
        else
            pB[i] ~ Uniform(0,1)
        end
    end

    for i in 1:n
        clickAs[i] ~ flip(pA[i])
        clickBs[i] ~ flip(pB[i])
    end

end

# Very similar
# clickAs = [true,true,true,false,false]
# clickBs = [true,true,true,false,false]

# Variant: completely different
clickAs = [true  ,true,true ,false,true ,true ,true ,false]
clickBs = [false,false,false,true ,false,false,false,true]
println("clickAs: ", clickAs)
println("clickBs: ", clickBs)
model = clickGraph(clickAs,clickBs)
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
import Base.convert

abstract type DynamicSystem end

show(io::IO, ds::DynamicSystem) = print(
    io,
    """$(typeof(ds))
    x     = $(get_x(ds))
    dx/dt = $(get_x_dot(ds))
    """
)

get_x(ds::DynamicSystem) = ds.x
get_x(ds::DynamicSystem, state::State) = get_x(convert(typeof(ds), state))
get_x_dot(ds::DynamicSystem) = ds.x_dot
get_x_dot(ds::DynamicSystem, state::State) = get_x(convert(typeof(ds), state))
get_state_equation(ds::DynamicSystem) = ds.state_equation
get_state_equation_ode_wrapper(ds::DynamicSystem) = ds.state_equation_ode_wrapper!

get_x_count(ds::DynamicSystem) = length(get_x(ds))
get_state(ds::DynamicSystem) = convert(State, ds)

set_x!(ds::DynamicSystem, x) = ds.x[:] = x
set_x_dot!(ds::DynamicSystem, x_dot) = ds.x_dot[:] = x_dot

# ---------------------- SixDOFEulerFixedMass ----------------------
struct SixDOFEulerFixedMass <: DynamicSystem
    x :: Array{T, 1} where T<:Number
    x_dot :: Array{T, 1} where T<:Number
    state_equation :: Function
    state_equation_ode_wrapper! :: Function
end

function six_dof_euler_fixed_mass_ode_wrapper!(dx, x, p, t)
    mass = p[1]
    inertia = p[2]
    forces = p[3]
    moments = p[4]
    h = p[5]
    dx[:] = six_dof_euler_fixed_mass(x, mass, inertia, forces, moments, h)
end

SixDOFEulerFixedMass(x, x_dot) = SixDOFEulerFixedMass(
    x,
    x_dot,
    six_dof_euler_fixed_mass,
    six_dof_euler_fixed_mass_ode_wrapper!
    )
SixDOFEulerFixedMass(x) = SixDOFEulerFixedMass(x, fill(NaN, 12))
SixDOFEulerFixedMass() = SixDOFEulerFixedMass(fill(NaN, 12))

function convert(::Type{SixDOFEulerFixedMass}, state::State)
    x = [
        get_body_velocity(state)...,
        get_body_ang_velocity(state)...,
        get_euler_angles(state)...,
        get_xyz_earth(state)...,
    ]

    x_dot = [
        get_body_accel(state)...,
        get_body_ang_accel(state)...,
        get_euler_angles_rates(state)...,
        get_horizon_velocity(state)...,
    ]

    SixDOFEulerFixedMass(x, x_dot)
end

function convert(state::State, ds::SixDOFEulerFixedMass)
    x = get_x(ds)
    x_dot = get_x_dot(ds)
    State(
        EarthPosition(x[10:12]..., get_position(state).ref_point),
        Attitude(x[7:9]...),
        x[1:3],
        x[4:6],
        x_dot[1:3],
        x_dot[4:6]
        )
end

# ---------------------- SixDOFQuaternionFixedMass ----------------------
struct SixDOFQuaternionFixedMass <: DynamicSystem
    x :: Array{T, 1} where T<:Number
    x_dot::Array{T, 1} where T<:Number
    state_equation :: Function
    state_equation_ode_wrapper! :: Function
    # TODO: must include k
    # k :: Number  # orthonormality error factor (see six_dof_quaternion_fixed_mass)
end

function six_dof_quaternion_fixed_mass_ode_wrapper!(dx, x, p, t, k=0.0)
    mass = p[1]
    inertia = p[2]
    forces = p[3]
    moments = p[4]
    h = p[5]
    dx[:] = six_dof_quaternion_fixed_mass(x, mass, inertia, forces, moments, h, k)
end

SixDOFQuaternionFixedMass(x, x_dot) = SixDOFQuaternionFixedMass(
    x,
    x_dot,
    six_dof_quaternion_fixed_mass,
    six_dof_quaternion_fixed_mass_ode_wrapper!
    )
SixDOFQuaternionFixedMass(x) = SixDOFQuaternionFixedMass(x, fill(NaN, 13))
SixDOFQuaternionFixedMass() = SixDOFQuaternionFixedMass(fill(NaN, 13))

function convert(state::State, ds::SixDOFQuaternionFixedMass)
    x = get_x(ds)
    x_dot = get_x_dot(ds)
    State(
        EarthPosition(x[11:13]..., get_position(state).ref_point),
        Attitude(x[7:10]...),
        x[1:3],
        x[4:6],
        x_dot[1:3],
        x_dot[4:6]
        )
end

function convert(::Type{SixDOFQuaternionFixedMass}, state::State)
    x = [
        get_body_velocity(state)...,
        get_body_ang_velocity(state)...,
        get_quaternions(state)...,
        get_xyz_earth(state)...,
    ]

    x_dot = [
        get_body_accel(state)...,
        get_body_ang_accel(state)...,
        get_quaternions_rate(state)...,
        get_horizon_velocity(state)...,
    ]

    SixDOFQuaternionFixedMass(x, x_dot)
end
# ---------------------- SixDOFECEFQuaternionFixedMass ----------------------
struct SixDOFECEFQuaternionFixedMass <: DynamicSystem
    x :: Array{T, 1} where T<:Number
    x_dot::Array{T, 1} where T<:Number
    state_equation :: Function
    state_equation_ode_wrapper! :: Function
    # TODO: must include k
    # k :: Number  # orthonormality error factor (see six_dof_ecef_quaternion_fixed_mass)
end

function six_dof_ecef_quaternion_fixed_mass_ode_wrapper!(dx, x, p, t, k=0.0)
    mass = p[1]
    inertia = p[2]
    forces = p[3]
    moments = p[4]
    h = p[5]
    dx[:] = six_dof_ecef_quaternion_fixed_mass(x, mass, inertia, forces, moments, h, k)
end

SixDOFECEFQuaternionFixedMass(x, x_dot) = SixDOFECEFQuaternionFixedMass(
    x,
    x_dot,
    six_dof_ecef_quaternion_fixed_mass,
    six_dof_ecef_quaternion_fixed_mass_ode_wrapper!
    )
SixDOFECEFQuaternionFixedMass(x) = SixDOFECEFQuaternionFixedMass(x, fill(NaN, 13))
SixDOFECEFQuaternionFixedMass() = SixDOFECEFQuaternionFixedMass(fill(NaN, 13))

function convert(state::State, ds::SixDOFECEFQuaternionFixedMass)
    x = get_x(ds)
    x_dot = get_x_dot(ds)
    State(
        ECEFPosition(x[11:13]...),
        Attitude(x[7:10]...),
        x[1:3],
        x[4:6],
        x_dot[1:3],
        x_dot[4:6]
        )
end

function convert(::Type{SixDOFECEFQuaternionFixedMass}, state::State)
    x = [
        get_body_velocity(state)...,
        get_body_ang_velocity(state)...,
        get_quaternions(state)...,
        get_xyz_ecef(state)...,
    ]

    x_dot = [
        get_body_accel(state)...,
        get_body_ang_accel(state)...,
        get_quaternions_rate(state)...,
        get_ecef_velocity(state)...,
    ]

    SixDOFECEFQuaternionFixedMass(x, x_dot)
end

# ---------------------- SixDOFAeroEulerFixedMass ----------------------
struct SixDOFAeroEulerFixedMass <: DynamicSystem
    x :: Array{T, 1} where T<:Number
    x_dot :: Array{T, 1} where T<:Number
    state_equation :: Function
    state_equation_ode_wrapper! :: Function
end

function six_dof_aero_euler_fixed_mass_ode_wrapper!(dx, x, p, t)
    mass = p[1]
    inertia = p[2]
    forces = p[3]
    moments = p[4]
    h = p[5]
    dx[:] = six_dof_aero_euler_fixed_mass(x, mass, inertia, forces, moments, h)
end

SixDOFAeroEulerFixedMass(x, x_dot) = SixDOFAeroEulerFixedMass(
    x,
    x_dot,
    six_dof_aero_euler_fixed_mass,
    six_dof_aero_euler_fixed_mass_ode_wrapper!
    )
SixDOFAeroEulerFixedMass(x) = SixDOFAeroEulerFixedMass(x, fill(NaN, 12))
SixDOFAeroEulerFixedMass() = SixDOFAeroEulerFixedMass(fill(NaN, 12))

function convert(::Type{SixDOFAeroEulerFixedMass}, state::State)
    x = [
        uvw_to_tasαβ(get_body_velocity(state)...)...,
        get_body_ang_velocity(state)...,
        get_euler_angles(state)...,
        get_xyz_earth(state)...,
    ]

    x_dot = [
        tas_α_β_dot_from_uvw_dot(
            get_body_velocity(state)...,
            get_body_accel(state)...
        )...,
        get_body_ang_accel(state)...,
        get_euler_angles_rates(state)...,
        get_horizon_velocity(state)...,
    ]

    SixDOFAeroEulerFixedMass(x, x_dot)
end

function convert(state::State, ds::SixDOFAeroEulerFixedMass)
    x = get_x(ds)
    x_dot = get_x_dot(ds)

    tas = x[1]
    α = x[2]
    β = x[3]

    tas_dot = x_dot[1]
    α_dot = x_dot[2]
    β_dot = x_dot[3]

    State(
        EarthPosition(
            x[10:12]..., get_position(state).ref_point
            ),
        Attitude(x[7:9]...),
        wind2body(tas, 0.0, 0.0, α, β),
        x[4:6],
        uvw_dot_from_tas_α_β(tas, α, β, tas_dot, α_dot, β_dot),
        x_dot[4:6]
        )
end


text:
# This code shows the functions to calculate the transfer matrices in high speed PM machines
# author: Jaime Renedo Anglada (renedo.jaime@ieee.org)

function TransferMatrix_comp( R_1, R_2, mu, sigma, omega, q)
# Calculation of the transfer matrix of a region. Complete methodology.
#   Detailed explanation goes here

    mu_0=4.0*pi*10.0^-7; # [m kg s^-2 A^-2]
    k_p=sqrt(im*omega*mu*sigma);

    F_2=mu_0*im*k_p*q/(mu*R_1)*(besselk(q,k_p*R_1)*besseli_d(q,k_p*R_1)-besseli(q,k_p*R_1)*besselk_d(q,k_p*R_1));
    M_2=[-mu_0*k_p/mu*besselk_d(q,k_p*R_1) -im*q/R_1*besselk(q,k_p*R_1); mu_0*k_p/mu*besseli_d(q,k_p*R_1) im*q/R_1*besseli(q,k_p*R_1)];
    N_2=[im*q/R_2*besseli(q,k_p*R_2) im*q/R_2*besselk(q,k_p*R_2); -mu_0*k_p/mu*besseli_d(q,k_p*R_2) -mu_0*k_p/mu*besselk_d(q,k_p*R_2)];
    T_mat=[1.0 0.0; 0.0 -1.0]*N_2*M_2/F_2*[1.0 0.0; 0.0 -1.0];

return T_mat
end

function besseli_d(q,x)
# Derivative of the modified Bessel function.
    I_prima=besseli(q-1,x)-q/x*besseli(q,x)
return I_prima
end

function besselk_d(q,x)
# Derivative of the modified Bessel function.
    K_prima=-besselk(q-1,x)-q/x*besselk(q,x)
return K_prima
end

function TransferMatrix_medium_simp( R_1, R_2, mu, sigma, omega, q, N_div)
# Calculation of the transfer matrix of a region. Medium simplified methodology.
mu_0=4.0*pi*10.0^-7; # [m kg s^-2 A^-2]

S_2l=R_2-R_1;
S_2d=S_2l/N_div;

T_2=[1.0 0.0; 0.0 1.0];
for count2 in 0:(N_div)
    delta=S_2d;
    r_2a=R_1+delta*count2+delta/2;
    k_2=q/r_2a;
    d_2=1/sqrt(sigma*mu*omega);
    gamma_2h=sqrt(k_2^2+im/d_2^2);
    beta_2h=gamma_2h/(im*mu*k_2);
#    T_2h=[1.0 gamma_2h*S_2d/beta_2h/mu_0; mu_0*beta_2h*gamma_2h*S_2d 1.0];
    T_2h=[cosh(gamma_2h*S_2d) sinh(gamma_2h*S_2d)/beta_2h/mu_0; mu_0*beta_2h*sinh(gamma_2h*S_2d) cosh(gamma_2h*S_2d)];

    T_2=T_2*T_2h;
end
return T_2
end

function TransferMatrix_simp( R_1, R_2, mu, sigma, omega, q, N_div)
# Calculation of the transfer matrix of a region. Medium simplified methodology.
#   Detailed explanation goes here

mu_0=4.0*pi*10.0^-7; # [m kg s^-2 A^-2]

S_2l=R_2-R_1;
S_2d=S_2l/N_div;

T_2=[1.0 0.0; 0.0 1.0];
for count2 in 0:(N_div)
    delta=S_2d;
    r_2a=R_1+delta*count2+delta/2;
    k_2=q/r_2a;
    d_2=1/sqrt(sigma*mu*omega);
    gamma_2h=sqrt(k_2^2+im/d_2^2);
    beta_2h=gamma_2h/(im*mu*k_2);
    T_2h=[1.0 gamma_2h*S_2d/beta_2h/mu_0; mu_0*beta_2h*gamma_2h*S_2d 1.0];
#     T_2h=[cosh(gamma_2h*S_2d) sinh(gamma_2h*S_2d)/beta_2h/mu_0; mu_0*beta_2h*sinh(gamma_2h*S_2d) cosh(gamma_2h*S_2d)];

    T_2=T_2*T_2h;

end
return T_2
end

function TransferMatrix_super_simp( R_1, R_2, mu, sigma, omega, q, N_div)
#Calculation of the transfer matrix of a region. Medium simplified methodology.
#   Detailed explanation goes here

mu_0=4.0*pi*10.0^-7; # [m kg s^-2 A^-2]

S_2l=R_2-R_1;
S_2d=S_2l/N_div;

delta=S_2d;
r_2a=(R_1+R_2)/2;
k_2=q/r_2a;
d_2=1/sqrt(sigma*mu*omega);
gamma_2h=sqrt(k_2^2+im/d_2^2);
beta_2h=gamma_2h/(im*mu*k_2);

T_2h=[1.0 gamma_2h*S_2d/beta_2h/mu_0; mu_0*beta_2h*gamma_2h*S_2d 1.0];
#     T_2h=[cosh(gamma_2h*S_2d) sinh(gamma_2h*S_2d)/beta_2h/mu_0; mu_0*beta_2h*sinh(gamma_2h*S_2d) cosh(gamma_2h*S_2d)];

T_2=T_2h^N_div;

return T_2
end

function TransferMatrix_nc( R_1, R_2, mu, sigma, omega, q)
# Calculation of the transfer matrix of a region. Complete methodology.
#   Detailed explanation goes here

    mu_0=4.0*pi*10.0^-7; # [m kg s^-2 A^-2]
    k_p=sqrt(im*omega*mu*sigma);

    F_2=mu_0*im*q/(mu*R_1*R_2);
    M_2=[-mu_0*R_1^(-q)/mu -im*q*R_1^(-q-1); mu_0*R_1^(q)/mu im*q*R_1^(q-1)];
    N_2=[im*q*R_1^(q-1) im*q*R_1^(-q-1); -mu_0*R_1^(q)/mu -mu_0*R_1^(-q)/mu];
    T_mat=[1.0 0.0; 0.0 -1.0]*N_2*M_2/F_2*[1.0 0.0; 0.0 -1.0];

return T_mat
end


text:
# code from the example notebook for automatic testing

using PyPlot

using Coils
using Coils.CoilsPlot

g, vertex_positions = cuboid_system([1, 1, 1], [3, 3, 3])

figure(figsize = (6, 6))
plot_vertices(vertex_positions)
plot_edges(g, vertex_positions, standalone = false)

poi = cuboid_poi([0.75, 0.75, 0.75], [0.0, 0.0, 0.0], [10, 10, 10], filled = false)

figure(figsize = (6, 6))
plot_system(g, vertex_positions, poi)

Bgoal(x) = [0, 100e-6, 0]

cells = find_cells(g)

figure(figsize = (8, 8))
plot_vertices(vertex_positions, labels = false)
plot_cells(cells, vertex_positions)

M = system_matrix(poi, vertex_positions, cells)

Bpoi = vcat(Bgoal.(poi)...)

optI = M \ Bpoi

figure(figsize = (10, 10))
plot_vertices(vertex_positions, labels = false)
current_norm = 1000 / maximum(abs, optI)
plot_cells(cells, vertex_positions, optI * current_norm)

edgecurrents = find_edgecurrents(g, cells, optI)

figure(figsize = (8, 8))
plot_vertices(vertex_positions, labels = false)
plot_edge_currents(g, edgecurrents, vertex_positions)

simpleloops, simpleloopscurrents = find_all_simpleloops(g, edgecurrents)

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[1], vertex_positions, color = "orange")
plot_loop(simpleloops[2], vertex_positions, color = "orange")
title("Current: $(simpleloopscurrents[1] * current_norm)")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[3], vertex_positions, color = "green")
plot_loop(simpleloops[4], vertex_positions, color = "green")
title("Current: $(simpleloopscurrents[3] * current_norm)")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[5], vertex_positions, color = "blue")
plot_loop(simpleloops[6], vertex_positions, color = "blue")
title("Current: $(simpleloopscurrents[5] * current_norm)")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[7], vertex_positions, color = "indigo")
plot_loop(simpleloops[8], vertex_positions, color = "indigo")
title("Current: $(simpleloopscurrents[7] * current_norm)")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[9], vertex_positions, color = "grey")
plot_loop(simpleloops[10], vertex_positions, color = "grey")
title("Current: $(simpleloopscurrents[9] * current_norm)")

figure()
plot_deviation_histogram(poi, simpleloops, simpleloopscurrents, vertex_positions, Bgoal)

figure()
plot_loops_field(cells, optI, vertex_positions, [1, 2, 3], 0;
        n = 50, spanA = [-0.5, 0.5], spanB = [-0.5, 0.5], Bref = Bgoal, levels = -10:10)

g, vertex_positions = cuboid_system([1, 2, 1], [3, 6, 3],
    skipfaces = [false, false, true, true, false, false])
poi = cuboid_poi([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [10, 10, 10], filled = true)

# the wires used to construct the coil
elemcurrents = [10.0, 1.0, 0.1]

figure(figsize = (6, 6))
plot_system(g, vertex_positions, poi)

# Find out the location of the openings
extrema(getindex.(vertex_positions, 2))

# prepare the large cells
front_face_unordered = [ i for i in eachindex(vertex_positions) if vertex_positions[i][2] < -0.9 ]
front_face = order_cell(g, front_face_unordered)

back_face_unordered = [ i for i in eachindex(vertex_positions) if vertex_positions[i][2] > 0.9 ]
back_face = order_cell(g, back_face_unordered)

initialcells = [front_face, back_face]

Bgoal(x) = [√2 * 100e-6 * x[2], √2 * 100e-6 * x[1], 0]

simpleloops, simpleloopscurrents = solve_system(g, vertex_positions, poi, Bgoal,
    verbose = true, initialcells = [front_face, back_face],
    mincurrent = minimum(elemcurrents) / 2)

simpleloopscurrents_decomp = decompose_currents(simpleloopscurrents, elemcurrents)
simpleloopscurrents_real = real_decomposed_currents(simpleloopscurrents_decomp, elemcurrents)

figure()
plot_deviation_histogram(poi, simpleloops, simpleloopscurrents_real, vertex_positions, Bgoal)

figure()
plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions, [1, 2, 3], 0;
        n = 50, spanA = [-0.5, 0.5], spanB = [-1, 1], Bref = Bgoal, levels = -20:4:20)

figure()
plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions, [1, 2, 3], 0;
        n = 50, spanA = [-0.5, 0.5], spanB = [-1, 1], levels = -50:5:50)

for i in eachindex(simpleloopscurrents)
    print("#$(lpad(i,2)): ")
    print("$(lpad(simpleloopscurrents[i], 20)) = ")
    print("$(lpad(simpleloopscurrents_decomp[i][1], 2)) * 10A  +  ")
    print("$(simpleloopscurrents_decomp[i][2]) * 1A  +  ")
    println("$(simpleloopscurrents_decomp[i][3]) * 0.1A")
end

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 1:8
    plot_loop(simpleloops[i], vertex_positions, color = "C0")
end
title("Current: $(simpleloopscurrents[1])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 9:16
    plot_loop(simpleloops[i], vertex_positions, color = "C1")
end
title("Current: $(simpleloopscurrents[9])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 17:20
    plot_loop(simpleloops[i], vertex_positions, color = "C2")
end
title("Current: $(simpleloopscurrents[17])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 21:28
    plot_loop(simpleloops[i], vertex_positions, color = "C3")
end
title("Current: $(simpleloopscurrents[21])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 29:30
    plot_loop(simpleloops[i], vertex_positions, color = "C4")
end
title("Current: $(simpleloopscurrents[29])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 31:38
    plot_loop(simpleloops[i], vertex_positions, color = "C5")
end
title("Current: $(simpleloopscurrents[31])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 39:42
    plot_loop(simpleloops[i], vertex_positions, color = "C6")
end
title("Current: $(simpleloopscurrents[39])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 43:46
    plot_loop(simpleloops[i], vertex_positions, color = "C7")
end
title("Current: $(simpleloopscurrents[43])")

figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 47:48
    plot_loop(simpleloops[i], vertex_positions, color = "C8")
end
title("Current: $(simpleloopscurrents[47])")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletpolyline

"""
    ''_leafletpolyline(;kwargs...)

A LeafletPolyline component.

Keyword arguments:
- `id` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `pathOptions` (optional)
- `positions` (optional): . positions has the following type: Array of lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)s | Array of Array of lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)ss
"""
function ''_leafletpolyline(; kwargs...)
        available_props = Symbol[:id, :loading_state, :pathOptions, :positions]
        wild_props = Symbol[]
        return Component("''_leafletpolyline", "LeafletPolyline", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
using PolynomialRoots: roots
using Roots: Order2, Newton, newton, find_zeros, find_zero

using .FiniteStrains: FromEulerianStrain, FromNaturalStrain

export vsolve

function vsolve(
    eos::PressureEquation{<:Murnaghan1st},
    p;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
)
    @unpack v0, b0, b′0 = getparam(eos)
    soln = [v0 * (1 + b′0 / b0 * p)^(-1 / b′0)]
    return _clamp(soln, bounds)
end
function vsolve(
    eos::PressureEquation{<:Murnaghan2nd},
    p;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    h = sqrt(2b0 * b″0 - b′0^2)
    k = b″0 * p + b′0
    numerator = exp(-2 / h * atan(p * h / (2b0 + p * b′0))) * v0
    denominator = (abs((k - h) / (k + h) * (b′0 + h) / (b′0 - h)))^(1 / h)
    soln = [numerator / denominator]
    return _clamp(soln, bounds)
end
function vsolve(
    eos::EnergyEquation{<:BirchMurnaghan2nd},
    e;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
)
    @unpack v0, b0, e0 = getparam(eos)
    Δ = (e - e0) / v0 / b0
    if Δ >= 0
        f = sqrt(2 / 9 * Δ)
        soln = map(FromEulerianStrain(v0), [f, -f])
        return _clamp(soln, bounds)
    elseif Δ < 0
        return typeof(v0)[]  # Complex strains
    else
        @assert false "Δ == (e - e0) / v0 / b0 == $Δ. this should never happen!"
    end
end
function vsolve(
    eos::PressureEquation{<:BirchMurnaghan2nd},
    p;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    stopping_criterion = 1e-20,  # Unitless
    chop = eps(),
    rtol = sqrt(eps()),
)
    @unpack v0, b0 = getparam(eos)
    # Solve f for (3 B0 f (2f + 1))^2 == p^2
    fs = roots(
        [-(p / 3b0)^2, 0, 1, 10, 40, 80, 80, 32];
        polish = true,
        epsilon = stopping_criterion,
    )
    soln = _strain2volume(eos, v0, fs, p, chop, rtol)
    return _clamp(soln, bounds)
end
function vsolve(
    eos::BulkModulusEquation{<:BirchMurnaghan2nd},
    b;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    stopping_criterion = 1e-20,  # Unitless
    chop = eps(),
    rtol = sqrt(eps()),
)
    @unpack v0, b0 = getparam(eos)
    # Solve f for ((7f + 1) * (2f + 1)^(5/2))^2 == (b/b0)^2
    fs = roots(
        [1 - (b / b0)^2, 24, 229, 1130, 3160, 5072, 4368, 1568];
        polish = true,
        epsilon = stopping_criterion,
    )
    soln = _strain2volume(eos, v0, fs, b, chop, rtol)
    return _clamp(soln, bounds)
end
function vsolve(
    eos::EnergyEquation{<:BirchMurnaghan3rd},
    e;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    chop = eps(),
    rtol = sqrt(eps()),
)
    @unpack v0, b0, b′0, e0 = getparam(eos)
    # Constrcut ax^3 + bx^2 + d = 0, see https://zh.wikipedia.org/wiki/%E4%B8%89%E6%AC%A1%E6%96%B9%E7%A8%8B#%E6%B1%82%E6%A0%B9%E5%85%AC%E5%BC%8F%E6%B3%95
    if b′0 == 4
        @warn "`b′0 == 4` for a `BirchMurnaghan3rd` is just a `BirchMurnaghan2nd`!"
        return vsolve(EnergyEquation(BirchMurnaghan2nd(v0, b0, e0)), e; bounds = bounds)
    else
        a = b′0 - 4
        r = 1 / 3a  # b = 1
        d = (e0 - e) / (9 / 2 * b0 * v0)
        p, q = -r^3 - d / 2a, -r^2
        Δ = p^2 + q^3
        fs = -r .+ if Δ > 0
            [cbrt(p + √Δ) + cbrt(p - √Δ)]  # Only 1 real solution
        elseif Δ < 0
            SIN, COS = sincos(acos(p / abs(r)^3) / 3)
            [2COS, -COS - √3 * SIN, -COS + √3 * SIN] .* abs(r)  # Verified
        elseif Δ == 0
            if p == q == 0
                [0]  # 3 reals are equal
            else  # p == -q != 0
                [2cbrt(p), -cbrt(p)]  # 2 real roots are equal, leaving 2 solutions
            end
        else
            @assert false "Δ == p^2 + q^3 == $Δ. this should never happen!"
        end  # solutions are strains
        soln = _strain2volume(eos, v0, fs, e, chop, rtol)
        return _clamp(soln, bounds)
    end
end
function vsolve(
    eos::PressureEquation{<:BirchMurnaghan3rd},
    p;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    stopping_criterion = 1e-20,  # Unitless
    chop = eps(),
    rtol = sqrt(eps()),
)
    @unpack v0, b0, b′0 = getparam(eos)
    # Solve f for (f (2f + 1)^(5/2) [2 + 3f (b′0 - 4)])^2 - (p / (3b0/2))^2 = 0
    fs = roots(
        [
            -(2p / 3 / b0)^2,
            0,
            4,
            12b′0 - 8,
            9b′0^2 + 48b′0 - 176,
            90b′0^2 - 80(2 + 3b′0),
            360b′0^2 + 320(7 - 6b′0),
            720b′0^2 + 64(122 - 75b′0),
            720b′0^2 + 768(13 - 7b′0),
            288b′0^2 + 2304(2 - b′0),
        ];
        polish = true,
        epsilon = stopping_criterion,
    )
    soln = _strain2volume(eos, v0, fs, p, chop, rtol)
    return _clamp(soln, bounds)
end
function vsolve(
    eos::EnergyEquation{<:BirchMurnaghan4th},
    e;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    stopping_criterion = 1e-20,  # Unitless
    chop = eps(),
    rtol = sqrt(eps()),
)
    @unpack v0, b0, b′0, b″0, e0 = getparam(eos)
    h = b0 * b″0 + b′0^2
    fs = roots(
        [(e0 - e) / (3 / 8 * v0 * b0), 0, 12, 12(b′0 - 4), 143 - 63b′0 + 9h];
        polish = true,
        epsilon = stopping_criterion,
    )
    soln = _strain2volume(eos, v0, fs, e, chop, rtol)
    return _clamp(soln, bounds)
end
function vsolve(
    eos::EnergyEquation{<:PoirierTarantola2nd},
    e;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
)
    @unpack v0, b0, e0 = getparam(eos)
    Δ = (e - e0) / v0 / b0
    if Δ >= 0
        f = sqrt(2 / 9 * Δ)
        soln = map(FromNaturalStrain(v0), [f, -f])
        return _clamp(soln, bounds)
    elseif Δ < 0
        return typeof(v0)[]  # Complex strains
    else
        @assert false "Δ == (e - e0) / v0 / b0 == $Δ. this should never happen!"
    end
end
function vsolve(
    eos::EquationOfStateOfSolids,
    y;
    bounds = (zero(eos.param.v0), 4 * eos.param.v0),
    xrtol = eps(),
    rtol = 4eps(),
)
    # Bisection method
    try
        return find_zeros(v -> eos(v) - y, bounds; xrtol = xrtol, rtol = rtol)
    catch e
        @error "cannot find solution! Come across `$e`!"
        return typeof(eos.param.v0)[]
    end
end
function vsolve(
    eos::EquationOfStateOfSolids,
    y,
    vᵢ,
    method = Order2();
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    maxiter = 40,
    verbose = false,
    xrtol = eps(),
    rtol = 4eps(),
)
    try
        vᵣ = find_zero(
            v -> eos(v) - y,
            vᵢ,
            method;
            maxevals = maxiter,
            verbose = verbose,
            xrtol = xrtol,
            rtol = rtol,
        )
        return _clamp([vᵣ], bounds)
    catch e
        @error "cannot find solution! Come across `$e`!"
        return typeof(vᵢ)[]
    end
end
function vsolve(
    eos::EnergyEquation,
    e,
    vᵢ,
    ::Newton;
    bounds = (zero(eos.param.v0), Inf * eos.param.v0),
    kwargs...,
)
    try
        vᵣ = newton(v -> eos(v) - e, v -> -PressureEquation(eos)(v), vᵢ; kwargs...)
        return _clamp([vᵣ], bounds)
    catch e
        @error "cannot find solution! Come across `$e`!"
        return typeof(vᵢ)[]
    end
end

function _strain2volume(
    eos::EquationOfStateOfSolids{<:BirchMurnaghan},
    v0,
    fs,
    y,
    chop = eps(),
    rtol = sqrt(eps()),
)
    @assert _ispositive(chop) && _ispositive(rtol) "either `chop` or `rtol` is less than 0!"
    return fs |>
           Base.Fix1(map, FromEulerianStrain(v0)) |>
           Base.Fix1(filter, x -> abs(imag(x)) < chop * oneunit(imag(x))) .|>  # If `x` has unit
           real |>
           Base.Fix1(filter, x -> isapprox(eos(x), y; rtol = rtol))
end

function _clamp(soln, bounds)
    low, high = extrema(bounds)
    return filter(v -> low <= v <= high, soln)
end


text:
module MyOptim

using LinearAlgebra

function conjgrad(fn, grad, x0; maxstep=2.0, niters=50, conv_thr = 1e-2)

    println("CG START")
    x = deepcopy(x0)
    storage = zeros(size(x))

#    println("x0 ", x)

    println("CG RUN GRAD1")
    #iter 1
    f, g = grad(storage, x)

#    println("f ", f, " x   ",  x)

    dx = deepcopy(-g)

    dx1 = deepcopy(dx)

    s = deepcopy(dx)

    step_size = maxstep/10.0

    println("CG RUN LS1")

    x, step_size = linesearch(x, dx, fn, f, step_size)
    step_size = min(step_size, maxstep)

    for i = 1:niters

        f, g = grad(storage, x)

        if sqrt(sum(g.^2)) < conv_thr && f < 0.1
            println("yes conv")
            return x, f, g
        end

        println("MY CONJGRAD $i ", f, " sg ", sqrt(sum(g.^2)), " step_size $step_size --------------------")

        dx1[:] = dx[:]
        dx[:] = -g[:]

#        println("dx ", dx)
#        println("dx1 ", dx1)

        beta = ( dx' * (dx - dx1)) / (dx1' * dx1)  #PR
#        println("beta $beta")
        beta = max(0, beta)

        s = dx + beta * s

        x, step_size, good = linesearch(x, s, fn, f, step_size)
        if good
            step_size =min(step_size, maxstep)
        else
            beta = 0.0
        end
        println("MY CG $i ", f, " sg ", sqrt(sum(g.^2)), " step_size $step_size good $good ")

    end

    println("conv not achieved")
    return x, f, g

end

function linesearch(x, dx, fn, f0, step_size)
#    println("CG FN0 $f0")
#    println("CG dx $dx")
#    println("ls $x $dx $f0 $step_size")
#    println("CG FN1 ss $step_size")
    f0r, flag0r = fn( x + dx * step_size * 0.0)
#    println("CG real $f0r $flag0r")

    f1, flag1 = fn( x + dx * step_size * 0.5)
#    println("CG FN1 $f1 $flag1")
    if f1 > f0 || flag1 == false
#        println("MY LS uphill $f1 > $f0 flag $flag1, reduce step")
        return x, step_size/3.1, false
    end
#    println("CG FN2")
    f2, flag2 = fn( x + dx * step_size)
#    println("CG FN2 $f2 $flag2")

    if flag2 == false
#        println("MY LS flag2 $flag2")
        return x + dx * step_size * 0.5, step_size/1.75, false
    end

    c = f0
    a = 2 * (f2-2*f1+f0)
    b = f2 - a - c

#    println("a b c $a $b $c ", -b/(2*a))

    if a < 0
#        println("MY LS a negative $a, take max step increase stepsize")
        return x + dx * step_size * 1.0, step_size * 2.1, true
    else
        step = min(-b/(2*a), 1.0)
#        println("MY LS NORMAL $step")
        return x + dx * step * step_size, step * 2.0, true
    end

end

end #end module

#=
function myf(x)

    return -0.5*x'*A*x + 0.1*sum(x.^4)

end

using ForwardDiff

function grad(storage, x)
    sleep(0.1)

    return myf(x), ForwardDiff.gradient(myf, x)

end

function grad_only(storage, x)
    storage[:] = ForwardDiff.gradient(myf, x)
    sleep(0.1)
    return storage
end
=#


text:
# numerial integration trapezoidal method

  function integration(a,b,n::Int64, func::Function)
    sum = 0.
    step = (b-a) / n # b and a are, respectivelly, upper and lower limits. n is the number of steps
    fa = func(a) / 2
    fb = func(b) / 2
    for i in 1:n-1
      x = i*step + a
      sum += func(x)
    end
    sum = (sum + fb + fa)*step

    return sum
  end

# Numerical integration

  function width(array)
    dt = ((maximum(array) - minimum(array))/(length(array)-1))
    return dt
  end

# NUMERICAL INTEGRATION - TRAPEZIUM METHOD

  function integration(x::Vector, y::Vector)
    dx = width(x)

# I = (h/2) * (y(1) + 2*y(2) + ... + 2y(n) + y(n+1)
    I = 0.
    n = length(y) - 1
    for i in 2:n
      I = I + y[i]
    end
    I = (dx/2)*(y[1]+2*I+y[n+1])

    return I

  end

# Simpson's method

  function integration_simp(x::Vector, y::Vector)
    dx = width(x)
    n = length(y) - 1

    if n % 2 != 0
      n += 1
    end

    I = y[1]
    for i in 2:n
      if i % 2 == 0
        I = I + 4*y[i]
      elseif i % 2 == 1
        I = I + 2*y[i]
      end
    end

    I = (dx/3) * (I + y[n+1])
    return I
  end


text:
using Nonconvex, LinearAlgebra, Test

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

options = NLoptOptions(xtol_rel = 1e-4)

@testset "Simple constraints" begin
    m = Model(f)
    addvar!(m, [0.0, 0.0], [10.0, 10.0])
    add_ineq_constraint!(m, x -> g(x, 2, 0))
    add_ineq_constraint!(m, x -> g(x, -1, 1))

    alg = NLoptAlg(:LD_MMA)
    r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
    @test abs(r.minimum - sqrt(8/27)) < 1e-6
    @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
end

@testset "Block constraints" begin
    m = Model(f)
    addvar!(m, [0.0, 0.0], [10.0, 10.0])
    add_ineq_constraint!(m, FunctionWrapper(x -> [g(x, 2, 0), g(x, -1, 1)], 2))

    alg = NLoptAlg(:LD_MMA)
    r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
    @test abs(r.minimum - sqrt(8/27)) < 1e-6
    @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
end

@testset "Infinite bounds" begin
    @testset "Infinite upper bound" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [Inf, Inf])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = NLoptAlg(:LD_MMA)
        r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-6
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
    end
    @testset "Infinite lower bound" begin
        m = Model(f)
        addvar!(m, [-Inf, -Inf], [10, 10])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = NLoptAlg(:LD_MMA)
        r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-6
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
    end
    @testset "Infinite upper and lower bound" begin
        m = Model(f)
        addvar!(m, [-Inf, -Inf], [Inf, Inf])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = NLoptAlg(:LD_MMA)
        r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-6
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
    end
end


text:
export Transition

struct Transition end

@node Transition Stochastic [ out, in, a ]

conjugate_type(::Type{ <: Transition }, ::Type{ Val{ :out } }) = Categorical
conjugate_type(::Type{ <: Transition }, ::Type{ Val{ :in} })   = Categorical
conjugate_type(::Type{ <: Transition }, ::Type{ Val{ :a } })   = MatrixDirichlet

@average_energy Transition (q_out::Any, q_in::Any, q_a::MatrixDirichlet) = begin
    return -probvec(q_out)' * logmean(q_a) * probvec(q_in)
end

@average_energy Transition (q_out_in::Contingency, q_a::MatrixDirichlet) = begin
    return -tr(contingency_matrix(q_out_in)' * logmean(q_a))
end


text:
#  print_invoice.jl

# Example code to print invoices

# print header

using Printf # package for formatting numbers

include("./infrastructure/infrastructure.jl")

function show_invoice(invoice::UnpaidInvoice)
    basic(invoice)
    println("="^60)
    println()
end

#==
function show_invoice(invoice::PaidInvoice)
    println("P A I D  $symbol$(@sprintf("%.2f", invoice.entry.amount))") : println("U N P A I D")
    println("="^60)
    println()
end
==#

function basic(invoice)
    println(invoice.header.name)
    println("Attn. $(invoice.header.name_contact)")
    println("$(invoice.header.address) ")
    print("$(invoice.header.city), ")
    println("$(invoice.header.zip)")
    println("$(invoice.header.country)")
    println()
    #print("INVOICE")
    #println("Date: $(Date(inv.date))")
    println("Invoice number: $(invoice.header.invoice_nbr)")
    println("="^40)
    println("Reference: $(invoice.header.order_ref) - $(string(Date(invoice.body.date)))")
    println()
    println("Training: $(invoice.body.name_training) - $(string(Date(invoice.body.date)))")
    println()
    nbr_students = length(invoice.body.students)
    unitprice = round(invoice.body.price_per_student * invoice.meta.currency_ratio; digits=2)
    price = round(nbr_students * unitprice; digits=2)
    vat = round(price*invoice.body.vat_perc; digits=2)
    totalprice = round(price + vat; digits=2)
    symbol = invoice.meta.currency

    println("#     Price                        Total")
    println("="^40)
    println("$nbr_students     $(@sprintf("%.2f", unitprice))                   $symbol$(@sprintf("%.2f", price))")
    println()
    println("STUDENT(S)")
    println("-"^20)
    for name in invoice.body.students
        println(name)
    end
    println()
    println("VAT                              $symbol$(@sprintf("%.2f", vat))")
    println("="^40)
    println("Total price                     $symbol$(@sprintf("%.2f", totalprice))")
    println()
end


text:
# # Crosswalk

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/crosswalk.ipynb)

# In this notebook we demonstrate how to define a crosswalk with pedestrians using `AutomotiveDrivingModels`.
# To do this, we define a crosswalk area as well as a pedestrian agent type.

# ## Generate a crosswalk environment

# We define a new concrete type that will contain the roadway (where cars drive)
# and the crosswalk definition which is just a regular lane.

using AutomotiveDrivingModels
using AutoViz
AutoViz.colortheme["background"] = colorant"white"; # hide
using Random

struct CrosswalkEnv
    roadway::Roadway{Float64}
    crosswalk::Lane{Float64}
end

# The crosswalk lane consists of a straight road segment perpendicular to the road.
# We will define the roadway just as a straight road.

# geometry parameters
roadway_length = 50.
crosswalk_length = 20.
crosswalk_width = 6.0
crosswalk_pos = roadway_length/2

# Generate a straight 2-lane roadway and a crosswalk lane
roadway = gen_straight_roadway(2, roadway_length)
crosswalk_start = VecE2(crosswalk_pos, -crosswalk_length/2)
crosswalk_end = VecE2(crosswalk_pos, crosswalk_length/2)
crosswalk_lane = gen_straight_curve(crosswalk_start, crosswalk_end, 2)
crosswalk = Lane(LaneTag(2,1), crosswalk_lane, width = crosswalk_width)
cw_segment = RoadSegment(2, [crosswalk])
push!(roadway.segments, cw_segment) # append it to the roadway

# initialize crosswalk environment
env = CrosswalkEnv(roadway, crosswalk)
#md nothing # hide

# **Render the crosswalk**

# We will define a new method to render this new environment.
# The roadway part is just rendered regularly, we add specific instuction
# for the crosswalk part that will display the white stripes.

using Cairo

function AutoViz.add_renderable!(rendermodel::RenderModel, env::CrosswalkEnv)

    ## render the road without the crosswalk
    roadway = gen_straight_roadway(2, roadway_length)
    add_renderable!(rendermodel, roadway)

    ## render crosswalk
    curve = env.crosswalk.curve
    n = length(curve)
    pts = Array{Float64}(undef, 2, n)
    for (i,pt) in enumerate(curve)
        pts[1,i] = pt.pos.x
        pts[2,i] = pt.pos.y
    end

    add_instruction!(
        rendermodel, render_dashed_line,
        (pts, colorant"white", env.crosswalk.width, 1.0, 1.0, 0.0, Cairo.CAIRO_LINE_CAP_BUTT)
    )
    return rendermodel
end

snapshot = render([env])
#md write("crosswalk.svg", snapshot) # hide
#md # ![crosswalk](crosswalk.svg)

# ### Navigate the crosswalk example

# Cars will be navigating in the roadway just as before.
# For the pedestrian we can define a new vehicle definition where we specify
# the size of the bounding box represented by the pedestrian.

# field of the VehicleDef type
fieldnames(VehicleDef)
## Agent.Class is from AutomotiveDrivingModels
const PEDESTRIAN_DEF = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
#md nothing # hide

# Car definition
car_initial_state = VehicleState(VecSE2(5.0, 0., 0.), roadway.segments[1].lanes[1],roadway, 8.0)
car = Entity(car_initial_state, VehicleDef(), :car)

# Pedestrian definition using our new Vehicle definition
ped_initial_state = VehicleState(VecSE2(+24.5,-7.0,π/2), env.crosswalk, roadway, 0.5)
ped = Entity(ped_initial_state, PEDESTRIAN_DEF, :pedestrian)

scene = Scene([car, ped])

# visualize the initial state
snapshot = render([env, scene])
#md write("crosswalk_initial.svg", snapshot) # hide
#md # ![initial state of crosswalk](crosswalk_initial.svg)

# ### Simulate the scenario

# As before, associate a driver model to each vehicle in the scene.
# We will use the model defined in the intersection example for both agents.

mutable struct LinearDriver <: DriverModel{LaneFollowingAccel}
    a::LaneFollowingAccel
    p::Float64 # confidence on the pedestrian intention
    k::Float64 # gain
end

function AutomotiveDrivingModels.observe!(model::LinearDriver, scene::Scene, roadway::Roadway, egoid)
    model.a = LaneFollowingAccel(model.k*model.p)
    ## change the confidence based on some policy
    ## you can get the position of the pedestrian from the scene
    model.p = 100.0
end
Base.rand(rng::AbstractRNG, model::LinearDriver) = model.a

timestep = 0.1
nticks = 50

## define a model for each entities present in the scene
models = Dict{Symbol, DriverModel}()

## Constant speed model
models[:car] = LinearDriver(LaneFollowingAccel(0.0), 20.0, -0.02)
models[:pedestrian] = IntelligentDriverModel(v_des=1.0)

## execute the simulation
scenes = simulate(scene, roadway, models, nticks, timestep)
#md nothing # hide

# ## Generate a video with Reel.jl

using Reel

function animate_record(scenes::Vector{Scene{E}},dt::Float64, env::CrosswalkEnv) where {E<:Entity}
    duration = length(scenes)*dt::Float64
    fps = Int(1/dt)
    function render_rec(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        return render([env, scenes[frame_index]])
    end
    return duration, fps, render_rec
end

duration, fps, render_hist = animate_record(scenes, timestep, env)
film = roll(render_hist, fps=fps, duration=duration)
#md write("crosswalk_animated.gif", film) # hide
#md nothing # hide
#md # ![animated crosswalk](crosswalk_animated.gif)


text:
module POMDPXFiles

using POMDPs
using POMDPModelTools
import POMDPs: action, value

# import o avoid naming conflict in POMDPs.jl (value is overloaded in LightXML)
import LightXML: parse_file, root, get_elements_by_tagname, attribute, content

export
    AbstractPOMDPXFile,
    POMDPXFile,
    MOMDPXFile,
    Alphas,
    POMDPAlphas,

    read_pomdp,
    action,
    value

include("writer.jl")
include("policy.jl")
include("read.jl")

end # module


text:
export plotmetrics

"""
Plot the metrics after some training. This function will plot all the metrics
in a single graph.

In order to avoid Photon being dependend on Plots, the calling code will have
to provide that module as the first parameter.

# Usage

```julia
fit!(workout, mydata, epochs=10)

import Plots
plotmetrics(Plots, workout)
```
"""
function plotmetrics(Plots::Module, workout::Workout, metrics=[:loss, :val_loss])
      p = nothing

      for (idx, metric) in enumerate(metrics)
            h = history(workout, metric)

            if idx == 1
                  p = Plots.plot(h..., xlabel = "steps", ylabel="values", label=metric)
            else
                  Plots.plot!(h..., label=metric)
            end
      end
      return p # otherwise no plotting in Juno
end


text:
using DocumenterGitLabCI
using Documenter

DocMeta.setdocmeta!(DocumenterGitLabCI, :DocTestSetup, :(using DocumenterGitLabCI); recursive=true)

makedocs(;
    modules=[DocumenterGitLabCI],
    authors="tester",
    repo="https://github.com/tester/DocumenterGitLabCI.jl/blob/{commit}{path}#{line}",
    sitename="DocumenterGitLabCI.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tester.gitlab.io/DocumenterGitLabCI.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
using Logging
using GasSteadySim
using DelimitedFiles

logger = SimpleLogger(stdout, Logging.Error)

function get_pressure(x; eos=:ideal)
    data = GasSteadySim._parse_data("./data/2-node/")
    data["pipes"]["1"]["length"] = x
    ss = initialize_simulator(data, eos=eos)
    run_simulator!(ss)
    p_fr = control(ss, :node, 1)[2]
    slack_pressure = p_fr * nominal_values(ss, :pressure) * 1e-6
    return slack_pressure, ss.sol["nodal_pressure"][2] * 1e-6
end

function get_density(x; eos=:ideal)
    data = GasSteadySim._parse_data("./data/2-node/")
    data["pipes"]["1"]["length"] = x
    ss = initialize_simulator(data, eos=eos)
    run_simulator!(ss)
    p_fr = control(ss, :node, 1)[2]
    rho_fr = (eos == :ideal) ? GasSteadySim._pressure_to_density_ideal(p_fr, ss.nominal_values, ss.params) :
        GasSteadySim._pressure_to_density_full_cnga(p_fr, ss.nominal_values, ss.params)
    rho_fr *= nominal_values(ss, :density)
    p = ss.sol["nodal_pressure"][2] / nominal_values(ss, :pressure)
    rho = (eos == :ideal) ? GasSteadySim._pressure_to_density_ideal(p, ss.nominal_values, ss.params) :
        GasSteadySim._pressure_to_density_full_cnga(p, ss.nominal_values, ss.params)
    rho *= nominal_values(ss, :density)
    return rho_fr, rho
end

function write_csv(results::Dict{Symbol,Any}, output_path, lengths)
    to_write_header = ["distance", "ideal", "cnga"]
    to_write = Array{Any,2}(undef, length(lengths)+1, 3)
    x = [0.0]
    append!(x, lengths)
    for i in 1:length(x)
        d = x[i]/1000.0
        to_write[i, :] = [d, results[:ideal][i], results[:full_cnga][i]]
    end
    open(output_path, "w") do io
        writedlm(io, [permutedims(to_write_header); to_write], ',')
    end
end

function run!(eos_list, lengths, results; type=:pressure)
    for x in lengths
        for eos in eos_list
            p = (type == :pressure) ? get_pressure(x, eos=eos) : get_density(x, eos=eos)
            if length(results[eos]) == 0
                push!(results[eos], p[1])
            end
            push!(results[eos], p[2])
        end
    end
end

function run()
    eos_list = [:ideal, :full_cnga]
    lengths = collect(range(1.0, 70.0, step=1.0)) * 1000.0
    results = Dict{Symbol,Any}(:ideal => [], :full_cnga => [])
    run!(eos_list, lengths,  results, type=:pressure)
    write_csv(results, "./output/ideal-vs-cnga-pressure.csv", lengths)
    results = Dict{Symbol,Any}(:ideal => [], :full_cnga => [])
    run!(eos_list, lengths,  results, type=:density)
    write_csv(results, "./output/ideal-vs-cnga-density.csv", lengths)
end

run()

max_stars_repo_name: kaarthiksundar/GasFlowRuns

max_issues_repo_name: kaarthiksundar/GasFlowRuns

max_forks_repo_name: kaarthiksundar/GasFlowRuns


text:
using OffsetArrays: OffsetArray

# replace first entries in a by b
merge_tups(a::Tuple, b) = merge_tups(a, to_tuple(b)::Tuple)
merge_tups(a::Tuple, ::Tuple{}) = a
merge_tups(a::Tuple, b::Tuple) = (first(b), merge_tups(Base.tail(a), Base.tail(b))...)

# only take first few entries of a so that result is as long as b
cut_tup(a::Tuple, b) = cut_tup(a, to_tuple(b)::Tuple)
cut_tup(a::Tuple, ::Tuple{}) = ()
cut_tup(a::Tuple, b::Tuple) = (first(a), cut_tup(Base.tail(a), Base.tail(b))...)

to_indexarray(t::Tuple{AbstractArray}) = t[1]
to_indexarray(t::Tuple{AbstractArray, Vararg{AbstractArray}}) = CartesianIndices(t)

_view(a::AbstractArray{<:Any, M}, b::AbstractArray{<:Any, N}) where {M, N} = view(a, b, ntuple(_ -> :, M-N)...)
_view(a::AbstractArray{<:Any, N}, b::AbstractArray{<:Any, N}) where {N} = view(a, b)

offsetrange(v, offset, range=()) = offsetrange(v, to_tuple(offset)::Tuple, to_tuple(range)::Tuple)

function offsetrange(v, offset::Tuple, range::Tuple=())
    short_axes = cut_tup(axes(v), offset)
    padded_range = merge_tups(short_axes, range)
    rel_offset = map(short_axes, offset, padded_range) do ax, off, r
        - off + first(r) - first(ax)
    end
    OffsetArray(to_indexarray(padded_range), rel_offset)
end

aroundindex(v, args...) = _view(v, offsetrange(v, args...))

function initstats(stat, ranges)
    ranges = to_tuple(ranges)
    itr = Iterators.product(ranges...)
    vec = [copy(stat) for _ in itr]
    return OffsetArray(vec, ranges)
end

# TODO combine fititer! and fitvec!
function fitvec!(m, val, shared_indices = eachindex(m, val))
    for cart in shared_indices
        @inbounds fit!(m[cart], val[cart])
    end
    return m
end

function fitvecmany!(m, iter)
    for el in iter
        am, ael = axes(m), axes(el)
        shared = (am == ael) ? eachindex(m, el) : CartesianIndices(map(intersect, am, ael))
        fitvec!(m, el, shared)
    end
    return m
end


text:
mutable struct _GP{T} <: Abstract_GP{T}
    dim::Int
    μ::Vector{T}
    Σ::Matrix{T}
    kernel::Kernel
    σ_k::Float64
    μ₀::PriorMean{T}
    K::PDMat{T,Matrix{T}}
    opt_ρ::Union{Optimizer,Nothing}
    opt_σ::Union{Optimizer,Nothing}
end

function _GP{T}(dim::Int,kernel::Kernel,mean::PriorMean,σ_k::Real,opt_ρ::Union{Optimizer,Nothing}=Adam(α=0.01),opt_σ=deepcopy(opt_ρ)) where {T<:Real}
    _GP{T}(dim,
            zeros(T,dim),
            Matrix{T}(I,dim,dim),
            deepcopy(kernel),
            σ_k,
            deepcopy(mean),
            PDMat(Matrix{T}(I,dim,dim)),
            deepcopy(opt_σ),
            deepcopy(opt_ρ))
end

mutable struct _VGP{T} <: Abstract_GP{T}
    dim::Int
    μ::Vector{T}
    Σ::Matrix{T}
    η₁::Vector{T}
    η₂::Symmetric{T,Matrix{T}}
    kernel::Kernel
    σ_k::Float64
    μ₀::PriorMean{T}
    K::PDMat{T,Matrix{T}}
    opt_ρ::Union{Optimizer,Nothing}
    opt_σ::Union{Optimizer,Nothing}
end

function _VGP{T}(dim::Int,kernel::Kernel,mean::PriorMean,σ_k::Real,opt_ρ::Union{Optimizer,Nothing}=Adam(α=0.01),opt_σ=deepcopy(opt_ρ)) where {T<:Real}
    _VGP{T}(dim,
            zeros(T,dim),
            Matrix{T}(I,dim,dim),
            zeros(T,dim),
            Symmetric(Matrix{T}(-0.5*I,dim,dim)),
            deepcopy(kernel),
            σ_k,
            deepcopy(mean),
            PDMat(Matrix{T}(I,dim,dim)),
            deepcopy(opt_σ),
            deepcopy(opt_ρ))
end

mutable struct _SVGP{T} <: Abstract_GP{T}
    dim::Int
    μ::Vector{T}
    Σ::Matrix{T}
    η₁::Vector{T}
    η₂::Symmetric{T,Matrix{T}}
    kernel::Kernel
    σ_k::Float64
    μ₀::PriorMean{T}
    Z::InducingPoints
    K::PDMat{T,Matrix{T}}
    Knm::Matrix{T}
    κ::Matrix{T}
    K̃::Vector{T}
    opt_ρ::Union{Optimizer,Nothing}
    opt_σ::Union{Optimizer,Nothing}
end

function _SVGP{T}(  dim::Int,nSamplesUsed::Int,
                    Z::Union{AbstractMatrix,InducingPoints},
                    kernel::Kernel,mean::PriorMean,σ_k::Real,
                    opt_ρ::Union{Optimizer,Nothing}=Adam(α=0.01),opt_σ::Union{Optimizer,Nothing}=opt_ρ
                 ) where {T<:Real}
    _SVGP{T}(dim,
            zeros(T,dim),
            Matrix{T}(I,dim,dim),
            zeros(T,dim),
            Symmetric(Matrix{T}(-0.5*I,dim,dim)),
            deepcopy(kernel),
            σ_k,
            deepcopy(mean),
            deepcopy(Z),
            PDMat(Matrix{T}(I,dim,dim)),
            Matrix{T}(undef,nSamplesUsed,dim),
            Matrix{T}(undef,nSamplesUsed,dim),
            Vector{T}(undef,nSamplesUsed),
            deepcopy(opt_ρ),
            deepcopy(opt_σ))
end

mutable struct _MCGP{T} <: Abstract_GP{T}
    dim::Int
    f::Vector{T}
    kernel::Kernel
    σ_k::Float64
    μ₀::PriorMean{T}
    K::PDMat{T,Matrix{T}}
end

function _MCGP{T}(dim::Int,kernel::Kernel,mean::PriorMean,σ_k::Real) where {T<:Real}
    _MCGP{T}(dim,
            zeros(T,dim),
            deepcopy(kernel),
            σ_k,
            deepcopy(mean),
            PDMat(Matrix{T}(I,dim,dim)))
end

mean_f(model::AbstractGP) = mean_f.(model.f)

mean_f(gp::_GP) = gp.μ
mean_f(gp::_VGP) = gp.μ
mean_f(gp::_SVGP) = gp.κ*gp.μ

diag_cov_f(model::AbstractGP) = diag_cov_f.(model.f)
diag_cov_f(gp::_GP{T}) where {T} = zeros(T,gp.dim)
diag_cov_f(gp::_VGP) = diag(gp.Σ)
diag_cov_f(gp::_SVGP) = opt_diag(gp.κ*gp.Σ,gp.κ) + gp.K̃

compute_K!(gp::Union{_GP,_VGP,_MCGP},X::AbstractMatrix,jitter) = gp.K = PDMat(gp.σ_k*(kernelmatrix(gp.kernel,X,obsdim=1)+jitter*I))
compute_K!(gp::_SVGP,jitter) = gp.K = PDMat(gp.σ_k*(kernelmatrix(gp.kernel,gp.Z.Z,obsdim=1)+jitter*I))

function compute_κ!(gp::_SVGP,X,jitter)
    gp.Knm .= gp.σ_k * kernelmatrix(gp.kernel,X,gp.Z.Z,obsdim=1)
    gp.κ .= gp.Knm / gp.K.mat
    gp.K̃ .= gp.σ_k * (kerneldiagmatrix(gp.kernel,X,obsdim=1) .+ jitter) - opt_diag(gp.κ,gp.Knm)
    @assert all(gp.K̃ .> 0) "K̃ has negative values"
end


text:
module TruncatedGaussQuadrature

#=
To import module...
push!(LOAD_PATH, "/Path/To/Src/")
... then load the module
using TruncatedGaussQuadrature
.. Test the module with:
truncatedGaussQuad(0.0, 0.1, -Inf, 10.0, 5)
=#

import Distributions
import Combinatorics
using LinearAlgebra # Hermitian, chol, permute, eigen

# include Distributions # loads other packages that subsequent code depends ons
include("normalDistMoment.jl")
include("momentGaussQuad.jl")
include("spd.jl")

export truncatedGaussQuad, normalDistMoment, momentGaussQuad, buildHankel, spdMatrix

function truncatedGaussQuad(μ::Real, σ::Real, lb::Real, ub::Real, n::Int;matrixFix::Bool=true)

  # 1) Calculate moments
  mm = normalDistMoment(μ, σ, 2n, lb, ub)
  # 2) Build Hankel Matrix
  M = buildHankel(mm,matrixFix)
  # 3) Calculate nodes:
  x, w = momentGaussQuad(M)
  return x, w
end

function truncatedGaussQuad(lb::Real, ub::Real, n::Int;matrixFix::Bool=true)
  # N(0,1)
  x, w = truncatedGaussQuad(0.0,1.0,lb,ub,n,matrixFix=matrixFix)
  return x, w
end

function buildHankel(mm::Array{T,1},matrixFix::Bool) where T<:Real
  n = div((length(mm)-1),2)
  M = zeros(eltype(mm),n+1,n+1)
  for i = 1:n+1
    M[:,i] = mm[i:i+n]
  end
  if  ~isposdef(M) && matrixFix
      M = spdMatrix(M)
  end
  return M
end

end


text:
using Test

@class ConsNative begin
    field1::Int
    field2::Int
    ConsNative(f1) = begin
        self = new(f1)
        self.field2 = f1
        self
    end
end

cNative = ConsNative(2)
@test cNative.field1 == cNative.field2

@class ConInit1 begin
    field1::Int
    field2::Int
    __init__(self, f1) = begin
        self.field1 = f1
        self.field2 = f1
    end

end

cInit1 = ConInit1(2)
@test cInit1.field1 == cInit1.field2

@class ConInit2 begin
    field1::Int
    field2::Int
    __init__(self, f1, f2=1) = begin
        self.field1 = f1
        if f1>f2
            self.field2 = f1
        else
            self.field2 = f2
        end
    end
end

cInit2a = ConInit2(2)
@test cInit2a.field1 == 2
@test cInit2a.field2 == 2

cInit2b = ConInit2(2,3)
@test cInit2b.field1 == 2
@test cInit2b.field2 == 3

@class ConInit3 begin
    field1::Int
    field2::Int
    __init__(self, f1=1; f2=1) = begin
        self.field1 = f1
        if f1>f2
            self.field2 = f1
        else
            self.field2 = f2
        end
    end
end

cInit3a = ConInit3(2)
@test cInit3a.field1 == 2
@test cInit3a.field2 == 2

cInit3b = ConInit3(2,f2=3)
@test cInit3b.field1 == 2
@test cInit3b.field2 == 3


text:
"""
```julia
calculate_tgrad(sys::AbstractTimeDependentSystem)
```

Calculate the time gradient of a system.

Returns a vector of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_tgrad end

"""
```julia
calculate_gradient(sys::AbstractSystem)
```

Calculate the gradient of a scalar system.

Returns a vector of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_gradient end

"""
```julia
calculate_jacobian(sys::AbstractSystem)
```

Calculate the jacobian matrix of a system.

Returns a matrix of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_jacobian end

"""
```julia
calculate_control_jacobian(sys::AbstractSystem)
```

Calculate the jacobian matrix of a system with respect to the system's controls.

Returns a matrix of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_control_jacobian end

"""
```julia
calculate_factorized_W(sys::AbstractSystem)
```

Calculate the factorized W-matrix of a system.

Returns a matrix of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_factorized_W end

"""
```julia
calculate_hessian(sys::AbstractSystem)
```

Calculate the hessian matrix of a scalar system.

Returns a matrix of [`Num`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_hessian end

"""
```julia
generate_tgrad(sys::AbstractTimeDependentSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generates a function for the time gradient of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_tgrad end

"""
```julia
generate_gradient(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generates a function for the gradient of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_gradient end

"""
```julia
generate_jacobian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the jacobian matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_jacobian end

"""
```julia
generate_factorized_W(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the factorized W-matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_factorized_W end

"""
```julia
generate_hessian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the hessian matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_hessian end

"""
```julia
generate_function(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generate a function to evaluate the system's equations.
"""
function generate_function end

Base.nameof(sys::AbstractSystem) = getfield(sys, :name)

#Deprecated
function independent_variable(sys::AbstractSystem)
    Base.depwarn("`independent_variable` is deprecated. Use `get_iv` or `independent_variables` instead.",:independent_variable)
    isdefined(sys, :iv) ? getfield(sys, :iv) : nothing
end

#Treat the result as a vector of symbols always
function independent_variables(sys::AbstractSystem)
    systype = typeof(sys)
    @warn "Please declare ($systype) as a subtype of `AbstractTimeDependentSystem`, `AbstractTimeIndependentSystem` or `AbstractMultivariateSystem`."
    if isdefined(sys, :iv)
        return [getfield(sys, :iv)]
    elseif isdefined(sys, :ivs)
        return getfield(sys,:ivs)
    else
        return []
    end
end

independent_variables(sys::AbstractTimeDependentSystem) = [getfield(sys, :iv)]
independent_variables(sys::AbstractTimeIndependentSystem) = []
independent_variables(sys::AbstractMultivariateSystem) = getfield(sys, :ivs)

function structure(sys::AbstractSystem)
    s = get_structure(sys)
    s isa SystemStructure || throw(ArgumentError("SystemStructure is not yet initialized, please run `sys = initialize_system_structure(sys)` or `sys = alias_elimination(sys)`."))
    return s
end
for prop in [
             :eqs
             :noiseeqs
             :iv
             :states
             :ps
             :var_to_name
             :ctrls
             :defaults
             :observed
             :tgrad
             :jac
             :ctrl_jac
             :Wfact
             :Wfact_t
             :systems
             :structure
             :op
             :equality_constraints
             :inequality_constraints
             :controls
             :loss
             :bcs
             :domain
             :ivs
             :dvs
             :connection_type
             :preface
            ]
    fname1 = Symbol(:get_, prop)
    fname2 = Symbol(:has_, prop)
    @eval begin
        $fname1(sys::AbstractSystem) = getfield(sys, $(QuoteNode(prop)))
        $fname2(sys::AbstractSystem) = isdefined(sys, $(QuoteNode(prop)))
    end
end

Setfield.get(obj::AbstractSystem, ::Setfield.PropertyLens{field}) where {field} = getfield(obj, field)
@generated function ConstructionBase.setproperties(obj::AbstractSystem, patch::NamedTuple)
    if issubset(fieldnames(patch), fieldnames(obj))
        args = map(fieldnames(obj)) do fn
            if fn in fieldnames(patch)
                :(patch.$fn)
            else
                :(getfield(obj, $(Meta.quot(fn))))
            end
        end
        kwarg = :($(Expr(:kw, :checks, false))) # Inputs should already be checked
        return Expr(:block,
            Expr(:meta, :inline),
            Expr(:call, :(constructorof($obj)), args..., kwarg)
        )
    else
        error("This should never happen. Trying to set $(typeof(obj)) with $patch.")
    end
end

rename(x::AbstractSystem, name) = @set x.name = name

function Base.propertynames(sys::AbstractSystem; private=false)
    if private
        return fieldnames(typeof(sys))
    else
        names = Symbol[]
        for s in get_systems(sys)
            push!(names, getname(s))
        end
        has_states(sys) && for s in get_states(sys)
            push!(names, getname(s))
        end
        has_ps(sys) && for s in get_ps(sys)
            push!(names, getname(s))
        end
        has_observed(sys) && for s in get_observed(sys)
            push!(names, getname(s.lhs))
        end
        return names
    end
end

Base.getproperty(sys::AbstractSystem, name::Symbol; namespace=true) = wrap(getvar(sys, name; namespace=namespace))
function getvar(sys::AbstractSystem, name::Symbol; namespace=false)
    systems = get_systems(sys)
    if isdefined(sys, name)
        Base.depwarn("`sys.name` like `sys.$name` is deprecated. Use getters like `get_$name` instead.", "sys.$name")
        return getfield(sys, name)
    elseif !isempty(systems)
        i = findfirst(x->nameof(x)==name, systems)
        if i !== nothing
            return namespace ? rename(systems[i], renamespace(sys, name)) : systems[i]
        end
    end

    if has_var_to_name(sys)
        avs = get_var_to_name(sys)
        v = get(avs, name, nothing)
        v === nothing || return namespace ? renamespace(sys, v) : v
    else
        sts = get_states(sys)
        i = findfirst(x->getname(x) == name, sts)
        if i !== nothing
            return namespace ? renamespace(sys, sts[i]) : sts[i]
        end

        if has_ps(sys)
            ps = get_ps(sys)
            i = findfirst(x->getname(x) == name,ps)
            if i !== nothing
                return namespace ? renamespace(sys, ps[i]) : ps[i]
            end
        end
    end

    sts = get_states(sys)
    i = findfirst(x->getname(x) == name, sts)

    if has_observed(sys)
        obs = get_observed(sys)
        i = findfirst(x->getname(x.lhs)==name,obs)
        if i !== nothing
            return namespace ? renamespace(sys, obs[i]) : obs[i]
        end
    end

    throw(ArgumentError("System $(nameof(sys)): variable $name does not exist"))
end

function Base.setproperty!(sys::AbstractSystem, prop::Symbol, val)
    # We use this weird syntax because `parameters` and `states` calls are
    # potentially expensive.
    if (
        params = parameters(sys);
        idx = findfirst(s->getname(s) == prop, params);
        idx !== nothing;
       )
        get_defaults(sys)[params[idx]] = value(val)
    elseif (
            sts = states(sys);
            idx = findfirst(s->getname(s) == prop, sts);
            idx !== nothing;
           )
        get_defaults(sys)[sts[idx]] = value(val)
    else
        setfield!(sys, prop, val)
    end
end

abstract type SymScope end

struct LocalScope <: SymScope end
LocalScope(sym::Union{Num, Symbolic}) = setmetadata(sym, SymScope, LocalScope())

struct ParentScope <: SymScope
    parent::SymScope
end
ParentScope(sym::Union{Num, Symbolic}) = setmetadata(sym, SymScope, ParentScope(getmetadata(value(sym), SymScope, LocalScope())))

struct GlobalScope <: SymScope end
GlobalScope(sym::Union{Num, Symbolic}) = setmetadata(sym, SymScope, GlobalScope())

renamespace(sys, eq::Equation) = namespace_equation(eq, sys)

function renamespace(sys, x)
    x = unwrap(x)
    if x isa Symbolic
        let scope = getmetadata(x, SymScope, LocalScope())
            if scope isa LocalScope
                rename(x, renamespace(getname(sys), getname(x)))
            elseif scope isa ParentScope
                setmetadata(x, SymScope, scope.parent)
            else # GlobalScope
                x
            end
        end
    else
        Symbol(getname(sys), :₊, x)
    end
end

namespace_variables(sys::AbstractSystem) = states(sys, states(sys))
namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
namespace_controls(sys::AbstractSystem) = controls(sys, controls(sys))

function namespace_defaults(sys)
    defs = defaults(sys)
    Dict((isparameter(k) ? parameters(sys, k) : states(sys, k)) => namespace_expr(defs[k], sys) for k in keys(defs))
end

function namespace_equations(sys::AbstractSystem)
    eqs = equations(sys)
    isempty(eqs) && return Equation[]
    map(eq->namespace_equation(eq, sys), eqs)
end

function namespace_equation(eq::Equation, sys)
    _lhs = namespace_expr(eq.lhs, sys)
    _rhs = namespace_expr(eq.rhs, sys)
    _lhs ~ _rhs
end

function namespace_assignment(eq::Assignment, sys)
    _lhs = namespace_expr(eq.lhs, sys)
    _rhs = namespace_expr(eq.rhs, sys)
    Assignment(_lhs, _rhs)
end

function namespace_expr(O, sys) where {T}
    ivs = independent_variables(sys)
    O = unwrap(O)
    if any(isequal(O), ivs)
        return O
    elseif isvariable(O)
        renamespace(sys, O)
    elseif istree(O)
        renamed = map(a->namespace_expr(a, sys), arguments(O))
        if symtype(operation(O)) <: FnType
            renamespace(sys, O)
        else
            similarterm(O, operation(O), renamed)
        end
    else
        O
    end
end

function states(sys::AbstractSystem)
    sts = get_states(sys)
    systems = get_systems(sys)
    unique(isempty(systems) ?
           sts :
           [sts; reduce(vcat,namespace_variables.(systems))])
end

function parameters(sys::AbstractSystem)
    ps = get_ps(sys)
    systems = get_systems(sys)
    isempty(systems) ? ps : [ps; reduce(vcat,namespace_parameters.(systems))]
end

function controls(sys::AbstractSystem)
    ctrls = get_ctrls(sys)
    systems = get_systems(sys)
    isempty(systems) ? ctrls : [ctrls; reduce(vcat,namespace_controls.(systems))]
end

function observed(sys::AbstractSystem)
    obs = get_observed(sys)
    systems = get_systems(sys)
    [obs;
     reduce(vcat,
            (map(o->namespace_equation(o, s), observed(s)) for s in systems),
            init=Equation[])]
end

Base.@deprecate default_u0(x) defaults(x) false
Base.@deprecate default_p(x) defaults(x) false
function defaults(sys::AbstractSystem)
    systems = get_systems(sys)
    defs = get_defaults(sys)
    isempty(systems) ? defs : mapreduce(namespace_defaults, merge, systems; init=defs)
end

states(sys::AbstractSystem, v) = renamespace(sys, v)
parameters(sys::AbstractSystem, v) = toparam(states(sys, v))
for f in [:states, :parameters]
    @eval $f(sys::AbstractSystem, vs::AbstractArray) = map(v->$f(sys, v), vs)
end

flatten(sys::AbstractSystem) = sys

function equations(sys::ModelingToolkit.AbstractSystem)
    eqs = get_eqs(sys)
    systems = get_systems(sys)
    if isempty(systems)
        return eqs
    else
        eqs = Equation[eqs;
               reduce(vcat,
                      namespace_equations.(get_systems(sys));
                      init=Equation[])]
        return eqs
    end
end

function preface(sys::ModelingToolkit.AbstractSystem)
    has_preface(sys) || return nothing
    pre = get_preface(sys)
    systems = get_systems(sys)
    if isempty(systems)
        return pre
    else
        pres = pre === nothing ? [] : pre
        for sys in systems
            pre = get_preface(sys)
            pre === nothing && continue
            for eq in pre
                push!(pres, namespace_assignment(eq, sys))
            end
        end
        return isempty(pres) ? nothing : pres
    end
end

function islinear(sys::AbstractSystem)
    rhs = [eq.rhs for eq ∈ equations(sys)]

    all(islinear(r, states(sys)) for r in rhs)
end

struct AbstractSysToExpr
    sys::AbstractSystem
    states::Vector
end
AbstractSysToExpr(sys) = AbstractSysToExpr(sys,states(sys))
function (f::AbstractSysToExpr)(O)
    !istree(O) && return toexpr(O)
    any(isequal(O), f.states) && return nameof(operation(O))  # variables
    if isa(operation(O), Sym)
        return build_expr(:call, Any[nameof(operation(O)); f.(arguments(O))])
    end
    return build_expr(:call, Any[operation(O); f.(arguments(O))])
end

###
### System utils
###
function push_vars!(stmt, name, typ, vars)
    isempty(vars) && return
    vars_expr = Expr(:macrocall, typ, nothing)
    for s in vars
        if istree(s)
            f = nameof(operation(s))
            args = arguments(s)
            ex = :($f($(args...)))
        else
            ex = nameof(s)
        end
        push!(vars_expr.args, ex)
    end
    push!(stmt, :($name = $collect($vars_expr)))
    return
end

function round_trip_expr(t, var2name)
    name = get(var2name, t, nothing)
    name !== nothing && return name
    t isa Sym && return nameof(t)
    istree(t) || return t
    f = round_trip_expr(operation(t), var2name)
    args = map(Base.Fix2(round_trip_expr, var2name), arguments(t))
    return :($f($(args...)))
end
round_trip_eq(eq, var2name) = Expr(:call, :~, round_trip_expr(eq.lhs, var2name), round_trip_expr(eq.rhs, var2name))

function push_eqs!(stmt, eqs, var2name)
    eqs_name = gensym(:eqs)
    eqs_expr = Expr(:vcat)
    eqs_blk = Expr(:(=), eqs_name, eqs_expr)
    for eq in eqs
        push!(eqs_expr.args, round_trip_eq(eq, var2name))
    end

    push!(stmt, eqs_blk)
    return eqs_name
end

function push_defaults!(stmt, defs, var2name)
    defs_name = gensym(:defs)
    defs_expr = Expr(:call, Dict)
    defs_blk = Expr(:(=), defs_name, defs_expr)
    for d in defs
        n = round_trip_expr(d.first, var2name)
        v = round_trip_expr(d.second, var2name)
        push!(defs_expr.args, :($(=>)($n, $v)))
    end

    push!(stmt, defs_blk)
    return defs_name
end

###
### System I/O
###
function toexpr(sys::AbstractSystem)
    sys = flatten(sys)
    expr = Expr(:block)
    stmt = expr.args

    name = Meta.quot(nameof(sys))
    ivs = independent_variables(sys)
    ivname = gensym(:iv)
    for iv in ivs
        ivname = gensym(:iv)
        push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
    end

    stsname = gensym(:sts)
    sts = states(sys)
    push_vars!(stmt, stsname, Symbol("@variables"), sts)
    psname = gensym(:ps)
    ps = parameters(sys)
    push_vars!(stmt, psname, Symbol("@parameters"), ps)

    var2name = Dict{Any,Symbol}()
    for v in Iterators.flatten((sts, ps))
        var2name[v] = getname(v)
    end

    eqs_name = push_eqs!(stmt, equations(sys), var2name)
    defs_name = push_defaults!(stmt, defaults(sys), var2name)

    if sys isa ODESystem
        iv = get_iv(sys)
        ivname = gensym(:iv)
        push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
        push!(stmt, :($ODESystem($eqs_name, $ivname, $stsname, $psname; defaults = $defs_name, name = $name, checks = false)))
    elseif sys isa NonlinearSystem
        push!(stmt, :($NonlinearSystem($eqs_name, $stsname, $psname; defaults = $defs_name, name = $name, checks = false)))
    end

    striplines(expr) # keeping the line numbers is never helpful
end

Base.write(io::IO, sys::AbstractSystem) = write(io, readable_code(toexpr(sys)))

function Base.show(io::IO, ::MIME"text/plain", sys::AbstractSystem)
    eqs = equations(sys)
    if eqs isa AbstractArray
        Base.printstyled(io, "Model $(nameof(sys)) with $(length(eqs)) equations\n"; bold=true)
    else
        Base.printstyled(io, "Model $(nameof(sys))\n"; bold=true)
    end
    # The reduced equations are usually very long. It's not that useful to print
    # them.
    #Base.print_matrix(io, eqs)
    #println(io)

    rows = first(displaysize(io)) ÷ 5
    limit = get(io, :limit, false)

    vars = states(sys); nvars = length(vars)
    Base.printstyled(io, "States ($nvars):"; bold=true)
    nrows = min(nvars, limit ? rows : nvars)
    limited = nrows < length(vars)
    defs = has_defaults(sys) ? defaults(sys) : nothing
    for i in 1:nrows
        s = vars[i]
        print(io, "\n  ", s)

        if defs !== nothing
            val = get(defs, s, nothing)
            if val !== nothing
                print(io, " [defaults to ")
                show(IOContext(io, :compact=>true, :limit=>true, :displaysize=>(1,displaysize(io)[2])), val)
                print(io, "]")
            end
        end
    end
    limited && print(io, "\n⋮")
    println(io)

    vars = parameters(sys); nvars = length(vars)
    Base.printstyled(io, "Parameters ($nvars):"; bold=true)
    nrows = min(nvars, limit ? rows : nvars)
    limited = nrows < length(vars)
    for i in 1:nrows
        s = vars[i]
        print(io, "\n  ", s)

        if defs !== nothing
            val = get(defs, s, nothing)
            if val !== nothing
                print(io, " [defaults to ")
                show(IOContext(io, :compact=>true, :limit=>true, :displaysize=>(1,displaysize(io)[2])), val)
                print(io, "]")
            end
        end
    end
    limited && print(io, "\n⋮")

    if has_structure(sys)
        s = get_structure(sys)
        if s !== nothing
            Base.printstyled(io, "\nIncidence matrix:"; color=:magenta)
            show(io, incidence_matrix(s.graph, Num(Sym{Real}(:×))))
        end
    end
    return nothing
end

function split_assign(expr)
    if !(expr isa Expr && expr.head === :(=) && expr.args[2].head === :call)
        throw(ArgumentError("expression should be of the form `sys = foo(a, b)`"))
    end
    name, call = expr.args
end

function _named(name, call, runtime=false)
    has_kw = false
    call isa Expr || throw(Meta.ParseError("The rhs must be an Expr. Got $call."))
    if length(call.args) >= 2 && call.args[2] isa Expr
        # canonicalize to use `:parameters`
        if call.args[2].head === :kw
            call.args[2] = Expr(:parameters, Expr(:kw, call.args[2].args...))
            has_kw = true
        elseif call.args[2].head === :parameters
            has_kw = true
        end
    end

    if !has_kw
        param = Expr(:parameters)
        if length(call.args) == 1
            push!(call.args, param)
        else
            insert!(call.args, 2, param)
        end
    end

    kws = call.args[2].args

    if !any(kw->(kw isa Symbol ? kw : kw.args[1]) == :name, kws) # don't overwrite `name` kwarg
        pushfirst!(kws, Expr(:kw, :name, runtime ? name : Meta.quot(name)))
    end
    call
end

function _named_idxs(name::Symbol, idxs, call)
    if call.head !== :->
        throw(ArgumentError("Not an anonymous function"))
    end
    if !isa(call.args[1], Symbol)
        throw(ArgumentError("not a single-argument anonymous function"))
    end
    sym, ex = call.args
    ex = Base.Cartesian.poplinenum(ex)
    ex = _named(:(Symbol($(Meta.quot(name)), :_, $sym)), ex, true)
    ex = Base.Cartesian.poplinenum(ex)
    :($name = $map($sym->$ex, $idxs))
end

check_name(name) = name isa Symbol || throw(Meta.ParseError("The lhs must be a symbol (a) or a ref (a[1:10]). Got $name."))

"""
    @named y = foo(x)
    @named y[1:10] = foo(x)
    @named y 1:10 i -> foo(x*i)

Rewrite `@named y = foo(x)` to `y = foo(x; name=:y)`.

Rewrite `@named y[1:10] = foo(x)` to `y = map(i′->foo(x; name=Symbol(:y_, i′)), 1:10)`.

Rewrite `@named y 1:10 i -> foo(x*i)` to `y = map(i->foo(x*i; name=Symbol(:y_, i)), 1:10)`.

Examples:
```julia
julia> using ModelingToolkit

julia> foo(i; name) = i, name
foo (generic function with 1 method)

julia> x = 41
41

julia> @named y = foo(x)
(41, :y)

julia> @named y[1:3] = foo(x)
3-element Vector{Tuple{Int64, Symbol}}:
 (41, :y_1)
 (41, :y_2)
 (41, :y_3)

julia> @named y 1:3 i -> foo(x*i)
3-element Vector{Tuple{Int64, Symbol}}:
 (41, :y_1)
 (82, :y_2)
 (123, :y_3)
```
"""
macro named(expr)
    name, call = split_assign(expr)
    if Meta.isexpr(name, :ref)
        name, idxs = name.args
        check_name(name)
        esc(_named_idxs(name, idxs, :($(gensym()) -> $call)))
    else
        check_name(name)
        esc(:($name = $(_named(name, call))))
    end
end

macro named(name::Symbol, idxs, call)
    esc(_named_idxs(name, idxs, call))
end

function _config(expr, namespace)
    cn = Base.Fix2(_config, namespace)
    if Meta.isexpr(expr, :.)
        return :($getvar($(map(cn, expr.args)...); namespace=$namespace))
    elseif Meta.isexpr(expr, :function)
        def = splitdef(expr)
        def[:args] = map(cn, def[:args])
        def[:body] = cn(def[:body])
        combinedef(def)
    elseif expr isa Expr && !isempty(expr.args)
        return Expr(expr.head, map(cn, expr.args)...)
    elseif Meta.isexpr(expr, :(=))
        return Expr(:(=), map(cn, expr.args)...)
    else
        expr
    end
end

"""
$(SIGNATURES)

Rewrite `@nonamespace a.b.c` to
`getvar(getvar(a, :b; namespace = false), :c; namespace = false)`.

This is the default behavior of `getvar`. This should be used when inheriting states from a model.
"""
macro nonamespace(expr)
    esc(_config(expr, false))
end

"""
$(SIGNATURES)

Rewrite `@namespace a.b.c` to
`getvar(getvar(a, :b; namespace = true), :c; namespace = true)`.
"""
macro namespace(expr)
    esc(_config(expr, true))
end

"""
$(SIGNATURES)

Structurally simplify algebraic equations in a system and compute the
topological sort of the observed equations. When `simplify=true`, the `simplify`
function will be applied during the tearing process.
"""
function structural_simplify(sys::AbstractSystem; simplify=false)
    sys = initialize_system_structure(alias_elimination(sys))
    check_consistency(sys)
    if sys isa ODESystem
        sys = dae_index_lowering(sys)
    end
    sys = tearing(sys, simplify=simplify)
    fullstates = [map(eq->eq.lhs, observed(sys)); states(sys)]
    @set! sys.observed = topsort_equations(observed(sys), fullstates)
    return sys
end

@latexrecipe function f(sys::AbstractSystem)
    return latexify(equations(sys))
end

Base.show(io::IO, ::MIME"text/latex", x::AbstractSystem) = print(io, latexify(x))

struct InvalidSystemException <: Exception
    msg::String
end
Base.showerror(io::IO, e::InvalidSystemException) = print(io, "InvalidSystemException: ", e.msg)

struct ExtraVariablesSystemException <: Exception
    msg::String
end
Base.showerror(io::IO, e::ExtraVariablesSystemException) = print(io, "ExtraVariablesSystemException: ", e.msg)

struct ExtraEquationsSystemException <: Exception
    msg::String
end
Base.showerror(io::IO, e::ExtraEquationsSystemException) = print(io, "ExtraEquationsSystemException: ", e.msg)

AbstractTrees.children(sys::ModelingToolkit.AbstractSystem) = ModelingToolkit.get_systems(sys)
AbstractTrees.printnode(io::IO, sys::ModelingToolkit.AbstractSystem) = print(io, nameof(sys))
AbstractTrees.nodetype(::ModelingToolkit.AbstractSystem) = ModelingToolkit.AbstractSystem

function check_eqs_u0(eqs, dvs, u0; check_length=true, kwargs...)
    if u0 !== nothing
        if check_length
            if !(length(eqs) == length(dvs) == length(u0))
                throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))), and initial conditions ($(length(u0))) are of different lengths. To allow a different number of equations than states use kwarg check_length=false."))
            end
        elseif length(dvs) != length(u0)
            throw(ArgumentError("States ($(length(dvs))) and initial conditions ($(length(u0))) are of different lengths."))
        end
    elseif check_length && (length(eqs) != length(dvs))
        throw(ArgumentError("Equations ($(length(eqs))) and states ($(length(dvs))) are of different lengths. To allow these to differ use kwarg check_length=false."))
    end
    return nothing
end

###
### Connectors
###

function with_connection_type(expr)
    @assert expr isa Expr && (expr.head == :function || (expr.head == :(=) &&
                                       expr.args[1] isa Expr &&
                                       expr.args[1].head == :call))

    sig = expr.args[1]
    body = expr.args[2]

    fname = sig.args[1]
    args = sig.args[2:end]

    quote
        struct $fname
            $(gensym()) -> 1 # this removes the default constructor
        end
        function $fname($(args...))
            function f()
                $body
            end
            res = f()
            $isdefined(res, :connection_type) ? $Setfield.@set!(res.connection_type = $fname) : res
        end
    end
end

macro connector(expr)
    esc(with_connection_type(expr))
end

promote_connect_rule(::Type{T}, ::Type{S}) where {T, S} = Union{}
promote_connect_rule(::Type{T}, ::Type{T}) where {T} = T
promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_type(promote_connect_rule(t1, t2), ts...)
@inline function promote_connect_type(::Type{T}, ::Type{S}) where {T,S}
    promote_connect_result(
        T,
        S,
        promote_connect_rule(T,S),
        promote_connect_rule(S,T)
    )
end

promote_connect_result(::Type, ::Type, ::Type{T}, ::Type{Union{}}) where {T} = T
promote_connect_result(::Type, ::Type, ::Type{Union{}}, ::Type{S}) where {S} = S
promote_connect_result(::Type, ::Type, ::Type{T}, ::Type{T}) where {T} = T
function promote_connect_result(::Type{T}, ::Type{S}, ::Type{P1}, ::Type{P2}) where {T,S,P1,P2}
    throw(ArgumentError("connection promotion for $T and $S resulted in $P1 and $P2. " *
                        "Define promotion only in one direction."))
end

throw_connector_promotion(T, S) = throw(ArgumentError("Don't know how to connect systems of type $S and $T"))
promote_connect_result(::Type{T},::Type{S},::Type{Union{}},::Type{Union{}}) where {T,S} = throw_connector_promotion(T,S)

promote_connect_type(::Type{T}, ::Type{T}) where {T} = T
function promote_connect_type(T, S)
    error("Don't know how to connect systems of type $S and $T")
end

function connect(syss...)
    connect(promote_connect_type(map(get_connection_type, syss)...), syss...)
end

###
### Inheritance & composition
###
function Base.hash(sys::AbstractSystem, s::UInt)
    s = hash(nameof(sys), s)
    s = foldr(hash, get_systems(sys), init=s)
    s = foldr(hash, get_states(sys), init=s)
    s = foldr(hash, get_ps(sys), init=s)
    if sys isa OptimizationSystem
        s = hash(get_op(sys), s)
    else
        s = foldr(hash, get_eqs(sys), init=s)
    end
    s = foldr(hash, get_observed(sys), init=s)
    s = hash(independent_variables(sys), s)
    return s
end

"""
    $(TYPEDSIGNATURES)

entend the `basesys` with `sys`, the resulting system would inherit `sys`'s name
by default.
"""
function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
    T = SciMLBase.parameterless_type(basesys)
    ivs = independent_variables(basesys)
    if !(typeof(sys) <: T)
        if length(ivs) == 0
            sys = convert_system(T, sys)
        elseif length(ivs) == 1
            sys = convert_system(T, sys, ivs[1])
        else
            throw("Extending multivariate systems is not supported")
        end
    end

    eqs = union(equations(basesys), equations(sys))
    sts = union(states(basesys), states(sys))
    ps = union(parameters(basesys), parameters(sys))
    obs = union(observed(basesys), observed(sys))
    defs = merge(defaults(basesys), defaults(sys)) # prefer `sys`
    syss = union(get_systems(basesys), get_systems(sys))

    if length(ivs) == 0
        T(eqs, sts, ps, observed = obs, defaults = defs, name=name, systems = syss)
    elseif length(ivs) == 1
        T(eqs, ivs[1], sts, ps, observed = obs, defaults = defs, name = name, systems = syss)
    end
end

Base.:(&)(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys)) = extend(sys, basesys; name=name)

"""
    $(SIGNATURES)

compose multiple systems together. The resulting system would inherit the first
system's name.
"""
function compose(sys::AbstractSystem, systems::AbstractArray{<:AbstractSystem}; name=nameof(sys))
    nsys = length(systems)
    nsys >= 1 || throw(ArgumentError("There must be at least 1 subsystem. Got $nsys subsystems."))
    @set! sys.name = name
    @set! sys.systems = systems
    return sys
end
compose(syss::AbstractSystem...; name=nameof(first(syss))) = compose(first(syss), collect(syss[2:end]); name=name)
Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2)

UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where p = getproperty(sys, p; namespace=false)