text:
include("../utils.jl")
include("group_utils.jl")
using NLsolve

@testset "Special Linear group" begin
    @testset "basic properties" begin
        G = SpecialLinear(3)
        @test G === SpecialLinear(3, ℝ)
        @test repr(G) == "SpecialLinear(3, ℝ)"
        @test base_manifold(G) === SpecialLinear(3)
        @test decorated_manifold(G) == GeneralLinear(3)
        @test number_system(G) === ℝ
        @test manifold_dimension(G) == 8
        @test representation_size(G) == (3, 3)
        Gc = SpecialLinear(2, ℂ)
        @test decorated_manifold(Gc) == GeneralLinear(2, ℂ)
        @test repr(Gc) == "SpecialLinear(2, ℂ)"
        @test number_system(Gc) == ℂ
        @test manifold_dimension(Gc) == 6
        @test representation_size(Gc) == (2, 2)
        Gh = SpecialLinear(4, ℍ)
        @test repr(Gh) == "SpecialLinear(4, ℍ)"
        @test number_system(Gh) == ℍ
        @test manifold_dimension(Gh) == 4 * 15
        @test representation_size(Gh) == (4, 4)

        @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
        @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(false)
        @test is_default_metric(
            MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())),
        ) === true
        @test @inferred(Manifolds.default_metric_dispatch(G, EuclideanMetric())) ===
              Val(true)
        @test @inferred(
            Manifolds.default_metric_dispatch(
                G,
                InvariantMetric(EuclideanMetric(), LeftAction()),
            )
        ) === Val(true)
        @test @inferred(
            Manifolds.default_metric_dispatch(
                MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())),
            )
        ) === Val(true)
        @test Manifolds.allocation_promotion_function(Gc, exp!, (1,)) === complex
    end

    @testset "Real" begin
        G = SpecialLinear(3)

        @test_throws DomainError is_manifold_point(G, randn(2, 3), true)
        @test_throws DomainError is_manifold_point(G, Float64[2 1; 1 1], true)
        @test_throws DomainError is_manifold_point(G, [1 0 im; im 0 0; 0 -1 0], true)
        @test_throws DomainError is_manifold_point(G, zeros(3, 3), true)
        @test_throws DomainError is_manifold_point(G, Float64[1 3 3; 1 1 2; 1 2 3], true)
        @test_throws DomainError is_manifold_point(
            G,
            make_identity(SpecialLinear(2), ones(2, 2)),
            true,
        )
        @test is_manifold_point(G, Float64[1 1 1; 2 2 1; 2 3 3], true)
        @test is_manifold_point(G, make_identity(G, ones(3, 3)), true)
        @test_throws DomainError is_tangent_vector(
            G,
            Float64[2 3 2; 3 1 2; 1 1 1],
            randn(3, 3),
            true;
            atol=1e-6,
        )
        @test_throws DomainError is_tangent_vector(
            G,
            Float64[2 1 2; 3 2 2; 2 2 1],
            Float64[2 1 -1; 2 2 1; 1 1 -1],
            true;
            atol=1e-6,
        )
        @test is_tangent_vector(
            G,
            Float64[2 1 2; 3 2 2; 2 2 1],
            Float64[-1 -1 -1; 1 -1 2; -1 -1 2],
            true;
            atol=1e-6,
        )

        types = [Matrix{Float64}]
        pts =
            [[2 -1 -3; 4 -1 -6; -1 1 2], [0 2 1; 0 -3 -1; 1 0 2], [-2 0 -1; 1 0 0; -1 -1 2]]
        vpts = [[0 -1 -5; 1 2 0; 1 2 -2], [0 -2 1; -2 1 2; -4 2 -1]]

        retraction_methods = [
            Manifolds.GroupExponentialRetraction(LeftAction()),
            Manifolds.GroupExponentialRetraction(RightAction()),
        ]

        inverse_retraction_methods = [
            Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
            Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
        ]

        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true)
            test_manifold(
                G,
                gpts;
                test_reverse_diff=false,
                test_forward_diff=false,
                test_injectivity_radius=false,
                test_project_point=true,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                exp_log_atol_multiplier=1e10,
                retraction_atol_multiplier=1e8,
                is_tangent_atol_multiplier=1e10,
            )
        end

        @testset "project" begin
            p = randn(3, 3)
            @test !is_manifold_point(G, p)
            q = project(G, p)
            @test is_manifold_point(G, q)
            @test project(G, q) ≈ q

            X = randn(3, 3)
            @test !is_tangent_vector(G, q, X; atol=1e-6)
            Y = project(G, q, X)
            @test is_tangent_vector(G, q, Y; atol=1e-6)
            @test project(G, q, Y) ≈ Y
        end
    end

    @testset "Complex" begin
        G = SpecialLinear(2, ℂ)

        @test_throws DomainError is_manifold_point(G, randn(ComplexF64, 2, 3), true)
        @test_throws DomainError is_manifold_point(G, randn(2, 2), true)
        @test_throws DomainError is_manifold_point(
            G,
            ComplexF64[1 0 im; im 0 0; 0 -1 0],
            true,
        )
        @test_throws DomainError is_manifold_point(G, ComplexF64[1 im; im 1], true)
        @test_throws DomainError is_manifold_point(
            G,
            make_identity(SpecialLinear(2), ones(ComplexF64, 2, 2)),
            true,
        )
        @test is_manifold_point(G, ComplexF64[im 1; -2 im], true)
        @test is_manifold_point(G, make_identity(G, ones(3, 3)), true)
        @test_throws DomainError is_tangent_vector(
            G,
            ComplexF64[-1+im -1; -im 1],
            ComplexF64[1-im 1+im; 1 -1+im],
            true;
            atol=1e-6,
        )
        @test_throws DomainError is_tangent_vector(
            G,
            ComplexF64[1 1+im; -1+im -1],
            ComplexF64[1-im -1-im; -im im],
            true;
            atol=1e-6,
        )
        @test is_tangent_vector(
            G,
            ComplexF64[1 1+im; -1+im -1],
            ComplexF64[1-im 1+im; 1 -1+im],
            true;
            atol=1e-6,
        )

        types = [Matrix{ComplexF64}]
        pts = [
            [0+1im 0+2im; 0+1im 0+1im],
            [-2+0im 0+1im; 0-1im -1+0im],
            [0+0im 0-1im; 0-1im -1+3im],
        ]
        vpts = [[0-1im -1+1im; -1+0im 0+1im], [1+1im 0+0im; 0+1im -1-1im]]

        retraction_methods = [
            Manifolds.GroupExponentialRetraction(LeftAction()),
            Manifolds.GroupExponentialRetraction(RightAction()),
        ]

        inverse_retraction_methods = [
            Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
            Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
        ]

        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true)
            test_manifold(
                G,
                gpts;
                test_reverse_diff=false,
                test_forward_diff=false,
                test_injectivity_radius=false,
                test_project_point=true,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                exp_log_atol_multiplier=1e8,
                retraction_atol_multiplier=1e8,
                is_tangent_atol_multiplier=1e8,
            )
        end

        @testset "project" begin
            p = randn(ComplexF64, 2, 2)
            @test !is_manifold_point(G, p)
            q = project(G, p)
            @test is_manifold_point(G, q)
            @test project(G, q) ≈ q

            X = randn(ComplexF64, 2, 2)
            @test !is_tangent_vector(G, q, X; atol=1e-6)
            Y = project(G, q, X)
            @test is_tangent_vector(G, q, Y; atol=1e-6)
            @test project(G, q, Y) ≈ Y
        end
    end
end


text:
"""
    StrongWolfeLineSearch{C1, C2, S, T}

`StrongWolfeLineSearch` computes a step length that satifies both the Armijo sufficient
decrease and the curvature conditions.

# Fields
- `c₁::C1`: Armijo condition constant
- `c₂::C2`: curvature condition constant
- `scale::S`: scale factor to expand step length bracket
- `α_max::T`: maximum step length
"""
struct StrongWolfeLineSearch{C1, C2, S, T}
    c₁::C1
    c₂::C2
    scale::S
    αₘₐₓ::T

    function StrongWolfeLineSearch(c₁, c₂, scale, αₘₐₓ)
        0 < c₁ < c₂ < 1 || error("Constants must satisfy 0 < `c₁` < `c₂` < 1")
        scale > 1 || error("`scale` factor must be larger than 1")
        C1 = typeof(c₁)
        C2 = typeof(c₂)
        S = typeof(scale)
        T = typeof(αₘₐₓ)
        return new{C1, C2, S, T}(c₁, c₂, scale, αₘₐₓ)
    end
end

"""
    StrongWolfeLineSearch(; kwargs...)

Initiate `StrongWolfeLineSearch`.

# Keywords
- `c1=1e-4`: Armijo condition constant
- `c2=0.9: curvature condition constant
- `scale=2.0`: scale factor to expand step length bracket
- `alphamax=100_000`: maximum step length
"""
function StrongWolfeLineSearch(;
        c1=1e-4,
        c2=0.9,
        scale=2.0,
        alphamax=100_000)
    return StrongWolfeLineSearch(c1, c2, scale, alphamax)
end

function (sw::StrongWolfeLineSearch)(f, state, p, α₀::T) where {T}
    x = state.x
    ϕ₀ = state.f
    dϕ₀ = state.∇f ⋅ p
    c₁ = sw.c₁
    c₂ = sw.c₂
    scale = T(sw.scale)
    αₘₐₓ = T(sw.αₘₐₓ)
    αₖ = zero(T)
    ϕαₖ = ϕ₀
    αₖ₊₁ = α₀
    ϕ(α) = f(x + α*p)
    iter = 0
    while αₖ₊₁ < αₘₐₓ
        ϕαₖ₊₁ = ϕ(αₖ₊₁)
        if ϕαₖ₊₁ > ϕ₀ + c₁*αₖ₊₁*dϕ₀ || (ϕαₖ₊₁ ≥ ϕαₖ && iter > 1)
            return _zoom(ϕ, αₖ, αₖ₊₁, ϕ₀, dϕ₀, ϕαₖ, c₁, c₂)
        end
        dϕαₖ₊₁ = Zygote.gradient(ϕ, αₖ₊₁)[1]
        abs(dϕαₖ₊₁) ≤ -c₂*dϕ₀ && return αₖ₊₁
        if dϕαₖ₊₁ ≥ 0
            return _zoom(ϕ, αₖ₊₁, αₖ, ϕ₀, dϕ₀, ϕαₖ₊₁, c₁, c₂)
        end
        αₖ = αₖ₊₁
        ϕαₖ = ϕαₖ₊₁
        αₖ₊₁ *= scale
        iter += 1
    end
    return αₘₐₓ
end

function _zoom(ϕ, αₗₒ, αₕᵢ, ϕ₀, dϕ₀, ϕαₗₒ, c₁, c₂)
    while true
        α = (αₗₒ + αₕᵢ) / 2
        ϕα = ϕ(α)
        if ϕα > ϕ₀ + c₁*α*dϕ₀ || ϕα ≥ ϕαₗₒ
            αₕᵢ = α
        else
            dϕα = Zygote.gradient(ϕ, α)[1]
            abs(dϕα) ≤ -c₂*dϕ₀ && return α
            if dϕα*(αₕᵢ - αₗₒ) ≥ 0
                αₕᵢ = αₗₒ
            end
            αₗₒ = α
            ϕαₗₒ = ϕα
        end
    end
end

ext: jl


text:
module EOS
using GSW

export rho_s, rho_c

function rho_s(temp::Float64, salt::Float64, pres::Float64, lon= -75.0, lat=37.7)
    saltA = GSW.gsw_sa_from_sp(salt, pres, lon, lat)
    ctmp = GSW.gsw_ct_from_t(saltA, temp, pres)
    rho = GSW.gsw_rho(saltA, ctmp, pres)
    return rho
end

function rho_c(temp::Float64, cond::Float64, pres::Float64, lon= -75.0, lat=37.7)
    salt = GSW.gsw_sp_from_c(cond, temp, pres)
    saltA = GSW.gsw_sa_from_sp(salt, pres, lon, lat)
    ctmp = GSW.gsw_ct_from_t(saltA, temp, pres)
    rho = GSW.gsw_rho(saltA, ctmp, pres)
    return rho
end

end


text:
# Include everything common between our LLVM versions (That's a lot)
include("../llvm_common.jl")

# Build the tarballs, then upload to Yggdrasil releases

# Unfortunately, it looks like these are just too big for GitHub.  :(
#=
name = "LLVMBootstrapDebug"
version = v"8.0.1"
sources, script, products, dependencies = llvm_build_args(;version=version, llvm_build_type="RelWithDebInfo")
build_info = build_tarballs(ARGS, name, version, sources, script, [host_platform], products, dependencies; skip_audit=true)
upload_and_insert_shards("JuliaPackaging/Yggdrasil", name, version, build_info)
=#

name = "LLVMBootstrap"
version = v"8.0.1"
sources, script, products, dependencies = llvm_build_args(;version=version)
build_info = build_tarballs(ARGS, name, version, sources, script, [host_platform], products, dependencies; skip_audit=true)
upload_and_insert_shards("JuliaPackaging/Yggdrasil", name, version, build_info)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("./dce_models/tofts.jl")
include("./dce_models/uptake.jl")
include("./dce_models/exchange.jl")
include("./dce_models/filtration.jl")
include("./dce_models/reference_region.jl")

function fit_model(modelname, modelopt = :default; kwargs...)
    return model_dict[modelname][modelopt](; kwargs...)
end

const model_dict = Dict{Symbol,Dict{Symbol,Function}}(
    :tofts => Dict{Symbol,Function}(
        :lls => fit_tofts_lls,
        :nls => fit_tofts_nls,
        :default => fit_tofts_nls,
    ),
    :extendedtofts => Dict{Symbol,Function}(
        :lls => fit_extendedtofts_lls,
        :nls => fit_extendedtofts_nls,
        :default => fit_extendedtofts_nls,
    ),
    :uptake => Dict{Symbol,Function}(
        :lls => fit_uptake_lls,
        :nls => fit_uptake_nls,
        :default => fit_uptake_nls,
    ),
    :exchange => Dict{Symbol,Function}(
        :lls => fit_exchange_lls,
        :nls => fit_exchange_nls,
        :default => fit_exchange_nls,
    ),
    :filtration => Dict{Symbol,Function}(
        :lls => fit_filtration_lls,
        :nls => fit_filtration_nls,
        :default => fit_filtration_nls,
    ),
    :rrm => Dict{Symbol,Function}(
        :lls => fit_rrm_lls,
        :nls => fit_rrm_nls,
        :default => fit_rrm_nls,
    ),
    :crrm => Dict{Symbol,Function}(
        :lls => fit_crrm_lls,
        :nls => fit_crrm_nls,
        :default => fit_crrm_nls,
    ),
    :errm => Dict{Symbol,Function}(:lls => fit_errm_lls, :default => fit_errm_lls),
    :cerrm => Dict{Symbol,Function}(:lls => fit_cerrm_lls, :default => fit_cerrm_lls),
    :rrift => Dict{Symbol,Function}(
        :crrm => fit_rrift_with_crrm,
        :cerrm => fit_rrift_with_cerrm,
        :default => fit_rrift_with_crrm,
    ),
)


text:
module GravityWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractForceModel = @jimport org.orekit.forces.AbstractForceModel
const CelestialBody = @jimport org.orekit.bodies.CelestialBody
const Class = @jimport java.lang.Class
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeDerivativesEquations = @jimport org.orekit.propagation.numerical.FieldTimeDerivativesEquations
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const GravityFields = @jimport org.orekit.forces.gravity.potential.GravityFields
const HolmesFeatherstoneAttractionModel = @jimport org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel
const IERSConventions = @jimport org.orekit.utils.IERSConventions
const JString = @jimport java.lang.String
const NewtonianAttraction = @jimport org.orekit.forces.gravity.NewtonianAttraction
const NormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider
const NormalizedSphericalHarmonicsProvider_NormalizedSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider$NormalizedSphericalHarmonics
const Object = @jimport java.lang.Object
const OceanTides = @jimport org.orekit.forces.gravity.OceanTides
const OceanTidesField = @jimport org.orekit.forces.gravity.OceanTidesField
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Relativity = @jimport org.orekit.forces.gravity.Relativity
const SingleBodyAbsoluteAttraction = @jimport org.orekit.forces.gravity.SingleBodyAbsoluteAttraction
const SingleBodyRelativeAttraction = @jimport org.orekit.forces.gravity.SingleBodyRelativeAttraction
const SolidTides = @jimport org.orekit.forces.gravity.SolidTides
const SolidTidesField = @jimport org.orekit.forces.gravity.SolidTidesField
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const ThirdBodyAttraction = @jimport org.orekit.forces.gravity.ThirdBodyAttraction
const ThirdBodyAttractionEpoch = @jimport org.orekit.forces.gravity.ThirdBodyAttractionEpoch
const TideSystem = @jimport org.orekit.forces.gravity.potential.TideSystem
const TimeDerivativesEquations = @jimport org.orekit.propagation.numerical.TimeDerivativesEquations
const UT1Scale = @jimport org.orekit.time.UT1Scale
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export HolmesFeatherstoneAttractionModel
export NewtonianAttraction
export OceanTides
export OceanTidesField
export Relativity
export SingleBodyAbsoluteAttraction
export SingleBodyRelativeAttraction
export SolidTides
export SolidTidesField
export ThirdBodyAttraction
export ThirdBodyAttractionEpoch
export acceleration
export add_contribution
export depends_on_position_only
export equals
export get_ae
export get_class
export get_derivatives_to_epoch
export get_events_detectors
export get_field_events_detectors
export get_max_degree
export get_max_order
export get_mu
export get_offset
export get_parameter_driver
export get_parameters
export get_parameters_drivers
export get_reference_date
export get_tide_system
export gradient
export hash_code
export init
export is_supported
export non_central_part
export notify
export notify_all
export on_date
export to_string
export value
export wait

include("PotentialWrapper/PotentialWrapper.jl")
include("holmes_featherstone_attraction_model.jl")
include("newtonian_attraction.jl")
include("ocean_tides.jl")
include("ocean_tides_field.jl")
include("relativity.jl")
include("single_body_absolute_attraction.jl")
include("single_body_relative_attraction.jl")
include("solid_tides.jl")
include("solid_tides_field.jl")
include("third_body_attraction.jl")
include("third_body_attraction_epoch.jl")

end


text:
#!/usr/bin/julia

# Trizen
# Date: 11 June 2019
# https://github.com/trizen

# Generate the smallest extended Chernick-Carmichael number with n prime factors.

# OEIS sequence:
#   https://oeis.org/A318646 -- The least Chernick's "universal form" Carmichael number with n prime factors.

# See also:
#   https://oeis.org/wiki/Carmichael_numbers
#   http://www.ams.org/journals/bull/1939-45-04/S0002-9904-1939-06953-X/home.html

using Primes

function trial_pretest(k::UInt64)

    if ((k %  3)==0 || (k %  5)==0 || (k %  7)==0 || (k % 11)==0 ||
        (k % 13)==0 || (k % 17)==0 || (k % 19)==0 || (k % 23)==0)
        return (k <= 23)
    end

    return true
end

function gcd_pretest(k::UInt64)

    if (k <= 107)
        return true
    end

    gcd(29*31*37*41*43*47*53*59*61*67, k) == 1 &&
    gcd(71*73*79*83*89*97*101*103*107, k) == 1
end

function is_chernick(n::Int64, m::UInt64)

    t = 9*m

    if (!trial_pretest(6*m + 1))
        return false
    end

    if (!trial_pretest(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!trial_pretest((t << i) + 1))
            return false
        end
    end

    if (!gcd_pretest(6*m + 1))
        return false
    end

    if (!gcd_pretest(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!gcd_pretest((t << i) + 1))
            return false
        end
    end

    if (!isprime(6*m + 1))
        return false
    end

    if (!isprime(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!isprime((t << i) + 1))
            return false
        end
    end

    return true
end

function chernick_carmichael(n::Int64, m::UInt64)
    prod = big(1)

    prod *= 6*m + 1
    prod *= 12*m + 1

    for i in 1:n-2
        prod *= ((big(9)*m)<<i) + 1
    end

    prod
end

function cc_numbers(from, to)

    for n in from:to

        multiplier = 1

        if (n > 4) multiplier = 1 << (n-4) end
        if (n > 5) multiplier *= 5 end

        m = UInt64(multiplier)

        while true

            if (is_chernick(n, m))
                println("a(", n, ") = ", chernick_carmichael(n, m))
                break
            end

            m += multiplier
        end
    end
end

cc_numbers(3, 10)

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
module Arg

	export arg_create
	export arg_key
	export arg_has
	export arg_get

	type ArgPart
		start::Int
		sz::Int
		pid::Int

		ArgPart() = ( x = new(); x.start = 1; x.sz = 0; x.pid = -1; return x; )
	end

	type ArgKey
		keys::Vector{String}
		iter::Vector{Int}
		match::Vector{Int}

		ArgKey() = new()
	end

	arg_key() = ( x = ArgKey(); x.keys = Array(String, 0); x.iter = Array(Int, 0); x.match = Array(Int, 0); return x; )
	arg_key(k::String) = ( x = arg_key(); add_key(x, k); return x; )
	arg_key(prev_ak::ArgKey, k::String) = ( x = arg_key(); add_key(x, prev_ak, k); return x; )

	function add_key(ak::ArgKey, k::String, it = 0, mtch = -1)
		push!(ak.keys, k); push!(ak.iter, it); push!(ak.match, mtch);
	end

	function add_key(ak::ArgKey, prev_ak::ArgKey, k::String)
		for i=1:length(prev_ak.keys)
			add_key(ak, prev_ak.keys[i]. prev_ak.iter[i], prev_ak.match[i])
		end
		add_key(ak, k)
	end

	function key_length(ak::ArgKey)
		return length(ak.keys)
	end

	type Node
		pid::Int
		id::Int
		key_start::Int
		key_sz::Int
		key_first::Char
		val_is_node::Bool
		val_start::Int
		val_sz::Int
		depth::Int
	end

	type Tree
		str::Vector{Char}
		len::Int
		str_orig::String
		nodes::Vector{Node}

		Tree() = ( x = new(); x.str = []; x.str_orig = ""; x.len = 0; x.nodes = Node[]; return x; )
	end

	function _parse_recurse(args, depth::Int, pid::Int, str_i::Int, next_str_i::Array{Int, 1})
		const char0 = char(0)
		const strchar = char('\'')
		str = args.str
		while(isempty(str) == false && str[str_i] != char0)
			if (str[str_i] == '/')
				str[str_i] = char0
				id = length(args.nodes)
				key_start = str_i+1; i = key_start;
				while (str[i] != char0 && str[i] != ':' && str[i] != '/' && str[i] != '}') i = i+1; end;
				node; success = false;
				if (str[i] == char0 || str[i] == '/' || str[i] == '}')
					node2 = Node(pid, id, key_start, i-key_start, str[key_start], false, -1, 0, depth); node = node2; success = true;
					str_i = i
				elseif (str[i] == ':')
					str[i] = char0
					if (str[i+1] != '{')
						val_start = i+1
						vi = val_start;
						in_string = (str[i+1] == strchar)
						while (str[vi] != char0 && (in_string || (str[vi] != '/' && str[vi] != '}')))
							vi = vi+1
							if (str[vi] == strchar)
								in_string = !in_string
							end
						end
						val_sz = vi-val_start
						if (val_sz >= 2 && str[val_start] == strchar && str[val_start+val_sz-1] == strchar)
							val_start = val_start+1; val_sz = val_sz-2;
						end
						node2 = Node(pid, id, key_start, i-key_start, str[key_start], false, val_start, val_sz, depth); node = node2; success = true;
						str_i = vi
					elseif (str[i+1] == '{')
						cont_str_i = [0]
						val_start = length(args.nodes); if (!_parse_recurse(args, depth+1, id, i+2, cont_str_i)) return false; end; val_sz = length(args.nodes)-val_start;
						node2 = Node(pid, id, key_start, i-key_start, str[key_start], args.nodes, val_start, val_sz, depth); node = node2; success = true;
						str_i = cont_str_i[0]
					end
				end

				if (success)
					push!(args.nodes, node)
				end

				if (str[str_i] == '}')
					str[str_i] = char0
					if (length(next_str_i) > 0 ) next_str_i[1] = str_i + 1; end
					return true
				end
			elseif (length(next_str_i) > 0 && str[str_i] == '}')
				str[str_i] = char0
				next_str_i[1] = str_i + 1
				return true
			else
				return false
			end
		end
		return true
	end

	function _sub_str(str, start::Int, sz::Int)
		return join(str[start:start+sz-1])
	end

	function _matches(str, node::Node, key, depth::Int)
		if (node.depth == depth && node.key_first == key[1] && node.key_sz == length(key))
			for j = 1:node.key_sz
			 	if (key[j] != str[node.key_start+j-1]) return false; end;
			end
			return true
		end
		return false
	end

	function _get_empty_sub()
		return ArgPart()
	end

	function _get_root_sub(args::Tree)
		ret = ArgPart(); ret.start = 1; ret.sz = length(args.nodes); ret.pid = -1; return ret;
	end

	function _to_part(args::Tree, ni::Int)
		ret = _get_empty_sub()
		if (ni >= 1 && ni <= length(args.nodes))
			ret.start = args.nodes[ni].val_start; ret.sz = args.nodes[ni].val_sz; ret.pid = args.nodes[ni].id;
		end
		return ret
	end

	function _find_recurse(args::Tree, sub::ArgPart, key::ArgKey, ki::Int)
		if (ki > key_length(key)) return -1; end;
		if (key.match[ki] != -1)
			found_ni = key.match[ki]
			if (ki == key_length(key)) return found_ni; end;
			return _find_recurse(args, _to_part(args, found_ni), key, ki+1)
		end
		matchi = 0;
		for ni = sub.start:(sub.start+sub.sz-1)
			found_ni = -1
			if (_matches(args.str, args.nodes[ni], key.keys[ki], ki))
				if (matchi == key.iter[ki])
					if (ki == key_length(key))
						found_ni = ni
					else
						if (args.nodes[ni].val_is_node)
							found_ni = _find_recurse(args, _to_part(args, ni), key, ki+1)
						end
					end
				else
					found_ni = -1
					matchi = matchi+1
				end
			end
			if (found_ni != -1) key.match[ki] = found_ni; return found_ni; end;
		end
		return -1
	end

	function _find(args::Tree, key::ArgKey)
		return _find_recurse(args, _get_root_sub(args), key, 1)
	end

	function _iter_key(args::Tree, index::Int, key::ArgKey)
		found_ni = _find(args, key)
		if (found_ni < 0 || args.nodes[found_ni].val_is_node == false || args.nodes[found_ni].val_sz <= index) return -1; end;
		return args.nodes[found_ni].val_start + index
	end

	function arg_create(str::String)
		args = Tree()
		args.str_orig = str
		args.str = Array(Char, length(str)+1); args.str[end] = char(0);
		for i = 1:length(str) args.str[i] = str[i]; end
		_parse_recurse(args, 1, -1, 1, Array(Int, 0))
		return args
	end

	function root(args::Tree)
		return _get_root_sub(args)
	end

	function to_key(k::Any)
		return typeof(k) == ArgKey ? k : arg_key(k)
	end

	function get_part(args::Tree, key)
		key = to_key(key)
		ni = _find(args, key); val = _to_part(args, ni);
		return ((ni != -1 && !args.nodes[ni].val_is_node), val)
	end

	function arg_has(args::Tree, key)
		key = to_key(key)
		return get_part(args, key)[1]
	end

	function has_sub(args::Tree, key)
		key = to_key(key)
		ni = _find(args, key);
		return (ni != -1 && args.nodes[ni].val_is_node)
	end

	function iter(args::Tree, key)
		key = to_key(key)
		if (key.match[key_length(key)] == -1) return has_sub(key); end;
		if (key_length(key)) key.iter[end] = key.iter[end]+1; key.match[end] = -1; end;
		return has_sub(args, key)
	end

	function parse_s(args::Tree, key, dflt::String)
		key = to_key(key)
		found, part = get_part(args, key)
		if (found == false) return dflt; end;
		return _sub_str(args.str, part.start, part.sz)
	end

	function conv_get(v::Any, dflt::Any)
		ret = v
		if (typeof(v) != typeof(dflt))
			typ = typeof(dflt)
			try
				ret = convert(typ, eval(Base.parse(string(v))))
			catch
			end
		end
		return ret
	end

	function dict_get(dict::Dict, key, dflt::Any)
		val = get(dict, key, dflt)
		return conv_get(val, dflt)
	end

	function arg_get(args::Tree, key, dflt::Any)
		key = to_key(key)
		has, part = get_part(args, key)
		if (has)
			val = _sub_str(args.str, part.start, part.sz)
			return conv_get(val, dflt)
		end
		return dflt
	end

	function arg_get(args::Tree, keys::Dict{Any, Any})
		for k in keys
			keys[k[1]] = arg_get(args, k[1], k[2])
		end

		for i = 1:length(args.nodes)
			node = args.nodes[i]
			if (node.val_is_node == false)
				key = _sub_str(args.str, node.key_start, node.key_sz)
				if (haskey(keys, key) == false)
					val = parse_s(args, key, "")
					val = (length(val) == 0) ? true : val
					keys[key] = val
				end
			end
		end
	end

	#TODO. convert the rest of the c++ code.

end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
type VcovWhite <: AbstractVcovMethod end

type VcovWhiteData <: AbstractVcovMethodData end

VcovMethodData(::VcovWhite, ::AbstractDataFrame) = VcovWhiteData()

function vcov!(v::VcovWhiteData, x::VcovData)
    S = shat!(v, x)
    return sandwich(x.crossmatrix, S)
end

function shat!{T}(::VcovWhiteData, x::VcovData{T, 1})
    X = x.regressors
    res = x.residuals
    Xu = scale!(res, X)
    S = At_mul_B(Xu, Xu)
    scale!(S, size(X, 1) / x.df_residual)
    return S
end

# S_{(l-1) * K + k, (l'-1)*K + k'} = \sum_i X[i, k] res[i, l] X[i, k'] res[i, l']
function shat!{T}(::VcovWhiteData, x::VcovData{T, 2})
    X = x.regressors
    res = x.residuals
    nobs = size(X, 1)
    dim = size(X, 2) * size(res, 2)
    S = fill(zero(Float64), (dim, dim))
    temp = fill(zero(Float64), nobs, dim)
    index = zero(Int)
    @inbounds for k in 1:size(X, 2), l in 1:size(res, 2)
        index += 1
        @simd for i in 1:nobs
            temp[i, index] = X[i, k]* res[i, l]
        end
    end
    S = At_mul_B(temp, temp)
    scale!(S, nobs / x.df_residual)
    return S
end


text:
using RenyiEntropyOfMagic

N = parse(Int, ARGS[1])
prepare_states(N, 0.1:0.1:2.5)


text:
"""
See FIRBasis for an examples

    a BasisFunction should implement:
    kernel()
    collabel() [default "colname_basis"]
    colnames()
    times()
    name()
    shiftOnset() [default 0]
"""
abstract type BasisFunction end

"""
Defines a FIRBasisfunction which can be called for each event, defining the time-expanded basis kernel
$(TYPEDEF)
$(TYPEDSIGNATURES)
$(FIELDS)
# Examples
```julia-repl
julia>  b = BasisFunction(kernelfunction,"derivative",["f(x)"],range(0,(length(kernelfunction([0, 1]))-1)*TR,step=TR),"hrf_kernel","basis_A",0)
```
"""
struct FIRBasis <: BasisFunction
    "a design-matrix kernel function used to timeexpand, given a timepoint in 'sample' timeunits"
    kernel::Function

    " vector of times along rows of kernel-output (in seconds)"
    times::AbstractVector
    "name of the event, random 1:1000 if unspecified"
    name::String
    "by how many samples do we need to shift the event onsets? This number is determined by how many 'negative' timepoints the basisfunction defines"
    shiftOnset::Integer
end

collabel(basis::FIRBasis) = :time
colnames(basis::FIRBasis) = basis.times

struct SplineBasis <: BasisFunction
    "a design-matrix kernel function used to timeexpand, given a timepoint in 'sample' timeunits"
    kernel::Function
    "name of column dimension (e.g 'time' for FIR, 'derivative', for HRF etc.)"

    "vector of names along columns of kernel-output"
    colnames::AbstractVector
    " vector of times along rows of kernel-output (in seconds)"
    times::AbstractVector
    "name of the event, random 1:1000 if unspecified"
    name::String
    "by how many samples do we need to shift the event onsets? This number is determined by how many 'negative' timepoints the basisfunction defines"
    shiftOnset::Integer
end

struct HRFBasis <: BasisFunction
    "a design-matrix kernel function used to timeexpand, given a timepoint in 'sample' timeunits"
    kernel::Function
    "vector of names along columns of kernel-output"
    colnames::AbstractVector
    " vector of times along rows of kernel-output (in seconds)"
    times::AbstractVector
    "name of the event, random 1:1000 if unspecified"
    name::String
end

function Base.show(io::IO, obj::BasisFunction)
    println(io, "name: $(name(obj))")
    println(io, "collabel: $(collabel(obj))")
    println(io, "colnames: $(colnames(obj))")
    println(io, "kerneltype: $(typeof(obj))")
    println(io, "times: $(times(obj))")
    println(io, "shiftOnset: $(shiftOnset(obj))")
end

"""
$(SIGNATURES)
Generate a sparse FIR basis around the *τ* timevector at sampling rate *sfreq*. This is useful if you cannot make any assumptions on the shape of the event responses. If unrounded events are supplied, they are split between samples. E.g. event-latency = 1.2 will result in a "0.8" and a "0.2" entry.

# Examples
Generate a FIR basis function from -0.1s to 0.3s at 100Hz
```julia-repl
julia>  f = firbasis([-0.1,0.3],100)
```
Evaluate at an event occuring at sample 103.3
```julia-repl
julia>  f(103.3)
```

"""
function firbasis(τ, sfreq, name::String)
    τ = round_times(τ, sfreq)
    times = range(τ[1], stop = τ[2], step = 1 ./ sfreq)

    kernel = e -> firkernel(e, times)

    shiftOnset = Int64(floor(τ[1] * sfreq))

    return FIRBasis(kernel, times, name, shiftOnset)
end
# cant multiple dispatch on optional arguments
#firbasis(;τ,sfreq)           = firbasis(τ,sfreq)
firbasis(; τ, sfreq, name) = firbasis(τ, sfreq, name)
firbasis(τ, sfreq) = firbasis(τ, sfreq, "basis_" * string(rand(1:10000)))

"""
$(SIGNATURES)
Calculate a sparse firbasis
# Examples

```julia-repl
julia>  f = firkernel(103.3,range(-0.1,step=0.01,stop=0.31))
```
"""
function firkernel(e, times)
    @assert ndims(e) <= 1 #either single onset or a row vector where we will take the first one :)
    if size(e, 1) > 1
        # XXX we will soon assume that the second entry would be the duration
        e = Float64(e[1])
    end
    e = [1 .- (e .% 1) e .% 1]
    e[isapprox.(e, 0, atol = 1e-15)] .= 0
    ksize = length(times) # kernelsize

    kernel =
        spdiagm(ksize + 1, ksize, 0 => repeat([e[1]], ksize), -1 => repeat([e[2]], ksize))

    return (kernel)

end

splinebasis(; τ, sfreq, nsplines, name) = splinebasis(τ, sfreq, nsplines, name)
splinebasis(τ, sfreq, nsplines) =
    splinebasis(τ, sfreq, nsplines, "basis_" * string(rand(1:10000)))

shiftOnset(basis::Union{SplineBasis,FIRBasis}) = basis.shiftOnset

function splinebasis(τ, sfreq, nsplines, name::String)
    τ = Unfold.round_times(τ, sfreq)
    times = range(τ[1], stop = τ[2], step = 1 ./ sfreq)
    kernel = e -> splinekernel(e, times, nsplines - 2)

    shiftOnset = Int64(floor(τ[1] * sfreq))
    colnames = spl_breakpoints(times, nsplines)
    return SplineBasis(kernel, "spline", colnames, times, name, shiftOnset)
end

function spl_breakpoints(times, nsplines)
    # calculate the breakpoints, evenly spaced
    return collect(range(minimum(times), stop = maximum(times), length = nsplines))
end

function splinekernel(e, times, nsplines)
    breakpoints = spl_breakpoints(times, nsplines)
    basis = BSplineBasis(4, breakpoints)  # 4= cubic
    return sparse(Unfold.splFunction(times, basis))
end

"""
$(SIGNATURES)
Generate a Hemodynamic-Response-Functio (HRF) basis with inverse-samplingrate "TR" (=1/FS)

Optional Parameters p:
                                                           defaults
                                                          {seconds}
        p(1) - delay of response (relative to onset)          6
        p(2) - delay of undershoot (relative to onset)       16
        p(3) - dispersion of response                         1
        p(4) - dispersion of undershoot                       1
        p(5) - ratio of response to undershoot                6
        p(6) - onset {seconds}                                0
        p(7) - length of kernel {seconds}                    32

# Examples
Generate a HRF basis function object with Sampling rate 1/TR. And evaluate it at an event occuring at TR 103.3 with duration of 4.1 TRs
```julia-repl
julia>  f = hrfbasis(2.3)
julia>  f(103.3,4.1)

```

"""
function hrfbasis(
    TR::Float64;
    parameters = [6.0 16.0 1.0 1.0 6.0 0.0 32.0],
    name::String = "basis_" * string(rand(1:10000)),
)
    # Haemodynamic response function adapted from SPM12b "spm_hrf.m"
    # Parameters:
    #                                                           defaults
    #                                                          {seconds}
    #        p(1) - delay of response (relative to onset)          6
    #        p(2) - delay of undershoot (relative to onset)       16
    #        p(3) - dispersion of response                         1
    #        p(4) - dispersion of undershoot                       1
    #        p(5) - ratio of response to undershoot                6
    #        p(6) - onset {seconds}                                0
    #        p(7) - length of kernel {seconds}                    32
    kernel = e -> hrfkernel(e, TR, parameters)
    return HRFBasis(
        kernel,
        ["f(x)"],
        range(0, (length(kernel([0, 1])) - 1) * TR, step = TR),
        name,
    )
end

shiftOnset(basis::HRFBasis) = 0

collabel(basis::HRFBasis) = :derivative
collabel(basis::SplineBasis) = :splineTerm

collabel(uf::UnfoldModel) = collabel(formula(uf))
collabel(form::FormulaTerm) = collabel(form.rhs)
collabel(t::Tuple) = collabel(t[1]) # MixedModels has Fixef+ReEf
collabel(term::Array{<:AbstractTerm}) = collabel(term[1].rhs)  # in case of combined formulas
#collabel(form::MatrixTerm) = collabel(form[1].rhs)

# typical defaults
colnames(basis::BasisFunction) = basis.colnames
kernel(basis::BasisFunction) = basis.kernel

times(basis::BasisFunction) = basis.times
name(basis::BasisFunction) = basis.name

"""
$(SIGNATURES)
Calculate a HRF kernel. Input e can be [onset duration]
# Examples

```julia-repl
julia>  f = hrfkernel(103.3,2.3,[6. 16. 1. 1. 6. 0. 32.])
```
"""
function hrfkernel(e, TR, p)
    @assert ndims(e) <= 1 #either single onset or a row vector where we will take the first one :)

    # code adapted from SPM12b
    mt = 32
    dt = TR / mt
    firstElem = Int(round((1 - (e[1] % 1)) .* mt))

    if size(e, 1) == 2 && e[2] > 0.0
        duration = e[2]
    else
        duration = 1.0 / mt
    end
    duration_mt = Int(ceil(duration .* mt))
    box_mt = [zeros(firstElem)' ones(duration_mt)']'
    #box_mt = box_mt ./ (sum(box_mt))
    # u   = [0:ceil(p(7)/dt)] - p(6)/dt;
    u = range(0, stop = ceil((p[7] + duration) / dt)) .- p[6] / dt

    #hrf = spm_Gpdf(u,p(1)/p(3),dt/p(3))
    # Note the inverted scale parameter compared to SPM12.
    g1 = Gamma(p[1] ./ p[3], p[3] ./ dt)
    #spm_Gpdf(u,p(2)/p(4),dt/p(4))
    g2 = Gamma(p[2] ./ p[4], p[4] ./ dt)
    # g1 - g2/p(5);
    hrf = pdf.(g1, u) .- pdf.(g2, u) ./ p[5]

    # hrf = hrf([0:floor(p(7)/RT)]*fMRI_T + 1);
    hrf = hrf ./ sum(hrf)
    hrf = conv(box_mt, hrf)'
    #println(box_mt)
    hrf = hrf[((range(0, stop = Int(floor(p[7] ./ TR + duration)))*mt)).+1]

    return (SparseMatrixCSC(sparse(hrf)))

end


text:
@testset "Simple tests" begin
    B = FEDVR(range(0,stop=20,length=71), 10)
    C = FEDVR(range(0,stop=20,length=71), 10)

    @test B == C
    @test axes(B) == (Inclusion(0..20), 1:631)
    @test size(B) == (ℵ₁, 631)

    @test nel(B) == 70

    @test B[0.0,1] ≠ 0.0
    @test B[20.0,end] ≠ 0.0

    @testset "Restricted bases" begin
        B̃ = B[:, 2:end-1]
        @test axes(B̃) == (Inclusion(0..20), 1:629)
        @test size(B̃) == (ℵ₁, 629)

        @test leftendpoint(B̃) == 0
        @test rightendpoint(B̃) == 20

        @test B̃[0.0,1] == 0.0
        @test B̃[20.0,end] == 0.0

        @test B isa FEDVR
        @test B isa CompactBases.FEDVROrRestricted
        @test B' isa QuasiAdjoint{<:Any,<:FEDVR}
        @test B' isa CompactBases.AdjointFEDVROrRestricted
        @test !(B̃ isa FEDVR)
        @test B̃ isa CompactBases.RestrictedFEDVR
        @test B̃ isa CompactBases.FEDVROrRestricted
        @test B̃' isa CompactBases.AdjointRestrictedFEDVR
        @test B̃' isa CompactBases.AdjointFEDVROrRestricted

        @test unrestricted_basis(B) == B
        @test unrestricted_basis(B̃) == B

        @test restriction_extents(B) == (0,0)
        @test restriction_extents(B̃) == (1,1)

        @test locs(B̃) == locs(B)[2:end-1]
    end
end

@testset "Pretty printing" begin
    B = FEDVR(range(0,stop=20,length=71), 10)
    C = FEDVR(range(0,stop=20,length=71), 10, t₀=10.0, ϕ=π/3)
    @test occursin("FEDVR{Float64} basis with 70 elements on 0.0..20.0", string(B))
    Ts = string(ComplexF64)
    @test occursin("FEDVR{$(Ts)} basis with 70 elements on 0.0..20.0 with ECS @ 60.00° starting at 10.00", string(C))
end

@testset "Element access" begin
    B = FEDVR(range(0,stop=20,length=71), 10)
    @test first(@elem(B, x, 1)) == 0
    @test last(@elem(B, x, 1)) == first(@elem(B, x, 2))
    @test last(@elem(B, x, 69)) == first(@elem(B, x, 70))
    @test_throws BoundsError @elem(B, x, 71)
end


text:
module RoundFast

import Base:(+),(-),(*),(/),sqrt
export square

include("Adjacent.jl")
using .Adjacent

include("eftArith.jl")
include("eftRound.jl")
include("rounding.jl")

end # module


text:
include("generator_utils.jl")

# repl_latex_symbols
include("repl_latex_symbols.jl")

latex_name_table = reverse(repl_latex_symbols)

write_dict_to_file("latex_name_table.jl", "latex_name_table", latex_name_table)


text:
using Behavior
using Niancat
using Niancat.Games.NiancatGames
using Niancat.Formatters

@given("that the team defaultteam uses Slack formatting in Swedish") do context
    service = context[:service]

    updateformatting!(service, "defaultteam", "Slack", "sv")
end

@when("we format a Correct(\"PUSSGURKA\") response") do context
    service = context[:service]

    formatter = getformatter(service, "defaultteam")
    response = NiancatGames.Correct("PUSSGURKA")

    context[:text] = format(formatter, response)
end

@when("the notification CorrectNotification for user Erik is formatted") do context
    service = context[:service]

    user = getuser(service, "Erik", "defaultteam")
    notification = NiancatGames.CorrectNotification(user)

    formatter = getformatter(service, "defaultteam")
    context[:text] = format(formatter, notification)
end

@then("the formatted response is") do context
    expected = context[:block_text]

    response = context[:text]

    @expect response == expected
end


text:
"""
    WorldClim{Union{BioClim,Climate,Weather}} <: RasterDataSource

Data from WorldClim datasets, either [`BioClim`](@ref), [`Climate`](@ref) or
[`Weather`](@ref).

See: [www.worldclim.org](https://www.worldclim.org)
"""
struct WorldClim{X} <: RasterDataSource end

const WORLDCLIM_URI = URI(scheme="https", host="biogeo.ucdavis.edu", path="/data/worldclim/v2.1")

resolutions(::Type{<:WorldClim}) = ("30s", "2.5m", "5m", "10m")
defres(::Type{<:WorldClim}) = "10m"

rasterpath(::Type{WorldClim{T}}) where T = joinpath(rasterpath(), "WorldClim", string(nameof(T)))
rasterpath(T::Type{<:WorldClim}, layer; kw...) =
    joinpath(rasterpath(T), string(layer), rastername(T, layer; kw...))

_zipfile_to_read(raster_name, zf) = first(filter(f -> f.name == raster_name, zf.files))

max_forks_repo_name: cesaraustralia/RasterDataSources


text: println("hello word from julia!")

max_stars_repo_name: Mxrcon/helloword-nextflow

max_issues_repo_name: Mxrcon/helloword-nextflow

max_forks_repo_name: Mxrcon/helloword-nextflow


text:
# ------------------------------------------------------------------------------
# convenience utilities
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# compute symbol matrix
# ------------------------------------------------------------------------------

"""
```julia
computesymbolsoverrange(operator, θ; mass=nothing, θ_min=-π/2)
```

Compute the eigenvalues and eigenvectors of the symbol matrix for an operator over
  a range of θ from -π/2 to 3π/2

# Arguments:
- `operator`:       Finite element operator to compute symbol matrices for
- `numbersteps1d`:  Number of values of θ to sample in each dimension
                        Note: numbersteps1d^dimension symbol matrices will be computed

# Keyword Arguments:
- `mass`:   Mass operator to invert for comparison to analytic solution,
                default: nothing
- `θ_min`:  Bottom of range of θ, shifts range to [θ_min, θ_min + 2π],
                default: -π / 2

# Returns:
- Values of θ sampled
- Eigenvalues of symbol matrix at θ sampled
- Eigenvectors of symbol matrix at θ sampled

# Examples:
```jldoctest
using LFAToolkit;

numbersteps1d = 5

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);

    # compute symbols
    (_, eigenvalues, _) = computesymbolsoverrange(diffusion, numbersteps1d);

    # verify
    eigenvalues = real(eigenvalues);
    if dimension == 1
        @assert minimum(eigenvalues[4, :]) ≈ 1
        @assert maximum(eigenvalues[4, :]) ≈ 4/3
    elseif dimension == 2
        @assert minimum(eigenvalues[19, :]) ≈ 2/3
        @assert maximum(eigenvalues[19, :]) ≈ 64/45
    elseif dimension == 3
        @assert minimum(eigenvalues[94, :]) ≈ 1/3
        @assert maximum(eigenvalues[94, :]) ≈ 256/225
    end
end

# output

```

```jldoctest
using LFAToolkit;

# setup
numbersteps1d = 3
mesh = Mesh2D(0.5, 0.5)
p = 6

# operators
mass = GalleryOperator("mass", p + 1, p + 1, mesh)
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# compute symbols
(_, eigenvalues, _) = computesymbolsoverrange(
    diffusion,
    numbersteps1d;
    mass = mass,
    θ_min = -π
);
eigenvalues = real(eigenvalues);

@assert minimum(eigenvalues[2, :]) ≈ π^2

# output

```
"""
function computesymbolsoverrange(
    operator::Operator,
    numbersteps1d::Int;
    mass::Union{Operator,Nothing} = nothing,
    θ_min::Float64 = -π / 2,
)
    # setup range
    dimension = operator.dimension
    numbersteps = numbersteps1d^dimension
    θ_max = θ_min + 2π
    θ_range1d = LinRange(θ_min, θ_max, numbersteps1d)
    θ_range = zeros(numbersteps, dimension)

    # setup
    numbermodes = size(operator.rowmodemap)[1]
    eigenvalues = zeros(ComplexF64, numbersteps, numbermodes)
    eigenvectors = zeros(ComplexF64, numbersteps, numbermodes, numbermodes)
    rangeiterator = []
    if dimension == 1
        rangeiterator = θ_range1d
    elseif dimension == 2
        rangeiterator = Iterators.product(θ_range1d, θ_range1d)
    elseif dimension == 3
        rangeiterator = Iterators.product(θ_range1d, θ_range1d, θ_range1d)
    end

    # compute
    for (step, θ) in enumerate(rangeiterator)
        θ = collect(θ)
        A = computesymbols(operator, θ)
        if mass != nothing
            A = computesymbols(mass, θ) \ A
        end
        currenteigen = eigen(A)
        eigenvalues[step, :] = currenteigen.values
        eigenvectors[step, :, :] = currenteigen.vectors
        θ_range[step, :] = θ
    end

    # return
    return (θ_range, eigenvalues, eigenvectors)
end

"""
```julia
computesymbolsoverrange(preconditioner, ω, θ; mass=nothing, θ_min=-π/2)
```

Compute the eigenvalues and eigenvectors of the symbol matrix for a preconditioned
  operator over a range of θ from -π/2 to 3π/2

# Arguments:
- `preconditioner`:  Preconditioner to compute symbol matries for
- `ω`:               Smoothing parameter array
- `numbersteps1d`:   Number of values of θ to sample in each dimension
                         Note: numbersteps1d^dimension symbol matrices will be computed

# Keyword Arguments:
- `mass`:   Mass operator to invert for comparison to analytic solution,
                default: nothing
- `θ_min`:  Bottom of range of θ, shifts range to [θ_min, θ_min + 2π],
                default: -π / 2

# Returns:
- Values of θ sampled
- Eigenvalues of symbol matrix at θ sampled
- Eigenvectors of symbol matrix at θ sampled

# Example:
```jldoctest
using LFAToolkit;

numbersteps1d = 5

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);

    # preconditioner
    chebyshev = Chebyshev(diffusion);

    # compute symbols
    (_, eigenvalues, _) = computesymbolsoverrange(chebyshev, [1], numbersteps1d);

    # verify
    eigenvalues = real(eigenvalues);
    if dimension == 1
        @assert minimum(eigenvalues[4, :]) ≈ 0.15151515151515105
        @assert maximum(eigenvalues[4, :]) ≈ 0.27272727272727226
    elseif dimension == 2
        @assert minimum(eigenvalues[19, :]) ≈ -0.25495098334134725
        @assert maximum(eigenvalues[19, :]) ≈ -0.17128758445192374
    elseif dimension == 3
        @assert minimum(eigenvalues[94, :]) ≈ -0.8181818181818181
        @assert maximum(eigenvalues[94, :]) ≈ -0.357575757575757
    end
end

# output

```
"""
function computesymbolsoverrange(
    preconditioner::AbstractPreconditioner,
    ω::Array,
    numbersteps1d::Int;
    mass::Union{Operator,Nothing} = nothing,
    θ_min::Float64 = -π / 2,
)
    # setup range
    dimension = preconditioner.operator.dimension
    numbersteps = numbersteps1d^dimension
    θ_max = θ_min + 2π
    θ_range1d = LinRange(θ_min, θ_max, numbersteps1d)
    θ_range = zeros(numbersteps, dimension)

    # setup
    numbermodes = size(preconditioner.operator.rowmodemap)[1]
    eigenvalues = zeros(ComplexF64, numbersteps, numbermodes)
    eigenvectors = zeros(ComplexF64, numbersteps, numbermodes, numbermodes)
    rangeiterator = []
    if dimension == 1
        rangeiterator = θ_range1d
    elseif dimension == 2
        rangeiterator = Iterators.product(θ_range1d, θ_range1d)
    elseif dimension == 3
        rangeiterator = Iterators.product(θ_range1d, θ_range1d, θ_range1d)
    end

    # compute
    for (step, θ) in enumerate(rangeiterator)
        θ = collect(θ)
        A = computesymbols(preconditioner, ω, θ)
        if mass != nothing
            A = computesymbols(mass, θ) \ A
        end
        currenteigen = eigen(A)
        eigenvalues[step, :] = currenteigen.values
        eigenvectors[step, :, :] = currenteigen.vectors
        θ_range[step, :] = θ
    end

    # return
    return (θ_range, eigenvalues, eigenvectors)
end

"""
```julia
computesymbolsoverrange(multigrid, ω, v, θ; mass=nothing, θ_min=-π/2)
```

Compute the eigenvalues and eigenvectors of the symbol matrix for a multigrid
  preconditioned operator over a range of θ from -π/2 to 3π/2

# Arguments:
- `multigrid`:      Preconditioner to compute symbol matries for
- `p`:              Smoothing parameter array
- `v`:              Pre and post smooths iteration count array, 0 indicates no pre or post smoothing
- `numbersteps1d`:  Number of values of θ to sample in each dimension
                       Note: numbersteps1d^dimension symbol matrices will be computed

# Keyword Arguments:
- `mass`:   Mass operator to invert for comparison to analytic solution,
                default: nothing
- `θ_min`:  Bottom of range of θ, shifts range to [θ_min, θ_min + 2π],
                default: -π / 2

# Returns:
- Values of θ sampled
- Eigenvalues of symbol matrix at θ sampled
- Eigenvectors of symbol matrix at θ sampled

# Example:
```jldoctest
using LFAToolkit;

numbersteps1d = 5

using LinearAlgebra

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    ctofbasis = TensorH1LagrangeBasis(3, 5, 1, dimension; collocatedquadrature=true);

    # operators
    finediffusion = GalleryOperator("diffusion", 5, 5, mesh);
    coarsediffusion = GalleryOperator("diffusion", 3, 5, mesh);

    # smoother
    jacobi = Jacobi(finediffusion);

    # preconditioner
    multigrid = PMultigrid(finediffusion, coarsediffusion, jacobi, [ctofbasis]);

    # compute symbols
    (_, eigenvalues, _) = computesymbolsoverrange(multigrid, [1.0], [1, 1], numbersteps1d);

    # verify
    eigenvalues = real(eigenvalues)
    if dimension == 1
       @assert maximum(eigenvalues[4, :]) ≈ 0.64
    elseif dimension == 2
       @assert maximum(eigenvalues[19, :]) ≈ 0.9082562365654528
    elseif dimension == 3
       @assert maximum(eigenvalues[94, :]) ≈ 1.4359882222222669
    end
end

# output

```
"""
function computesymbolsoverrange(
    multigrid::Multigrid,
    p::Array,
    v::Array,
    numbersteps1d::Int;
    mass::Union{Operator,Nothing} = nothing,
    θ_min::Float64 = -π / 2,
)
    # setup range
    dimension = multigrid.fineoperator.dimension
    numbersteps = numbersteps1d^dimension
    θ_max = θ_min + 2π
    θ_range1d = LinRange(θ_min, θ_max, numbersteps1d)
    θ_range = zeros(numbersteps, dimension)

    # setup
    numbermodes = size(multigrid.fineoperator.rowmodemap)[1]
    eigenvalues = zeros(ComplexF64, numbersteps, numbermodes)
    eigenvectors = zeros(ComplexF64, numbersteps, numbermodes, numbermodes)
    rangeiterator = []
    if dimension == 1
        rangeiterator = θ_range1d
    elseif dimension == 2
        rangeiterator = Iterators.product(θ_range1d, θ_range1d)
    elseif dimension == 3
        rangeiterator = Iterators.product(θ_range1d, θ_range1d, θ_range1d)
    end

    # compute
    for (step, θ) in enumerate(rangeiterator)
        θ = [abs(θ_i) > 100 * eps() ? θ_i : 100 * eps() for θ_i in θ]
        A = computesymbols(multigrid, p, v, θ)
        if mass != nothing
            A = computesymbols(mass, θ) \ A
        end
        currenteigen = eigen(A)
        eigenvalues[step, :] = currenteigen.values
        eigenvectors[step, :, :] = currenteigen.vectors
        θ_range[step, :] = θ
    end

    # return
    return (θ_range, eigenvalues, eigenvectors)
end

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


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using CSV
using PyPlot
ioff() # hide
using HTTP
using StableRNGs

MLJ.color_off() # hide

req = HTTP.get("https://raw.githubusercontent.com/rupakc/UCI-Data-Analysis/master/Airfoil%20Dataset/airfoil_self_noise.dat");

df = CSV.read(req.body, DataFrames.DataFrame; header=[
                   "Frequency","Attack_Angle","Chord+Length",
                   "Free_Velocity","Suction_Side","Scaled_Sound"
                   ]
              );
df[1:5, :] |> pretty

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
mutable struct GenericModel
  is_directed::Bool
  vertices::Vector{Variable}           # An indexable container of vertices (variables)
  edges::Vector{Dependence}            # An indexable container of edges (dependencies)
  finclist::Vector{Vector{Dependence}} # Forward incidence list
  binclist::Vector{Vector{Dependence}} # Backward incidence list
  ofkey::Dict{Symbol, Integer}         # Dictionary storing index of vertex (variable) of corresponding key
end

getindex(m::GenericModel, k::Symbol) = m.vertices[m.ofkey[k]]

is_directed(m::GenericModel) = m.is_directed

num_vertices(m::GenericModel) = length(m.vertices)
vertices(m::GenericModel) = m.vertices

num_edges(m::GenericModel) = length(m.edges)
edges(m::GenericModel) = m.edges

vertex_index(v::Integer, m::GenericModel) = (v <= m.vertices[end] ? v : 0)
vertex_index(v::Variable, m::GenericModel) = vertex_index(v)

edge_index(d::Dependence, m::GenericModel) = edge_index(d)
source(d::Dependence, m::GenericModel) = d.source
target(d::Dependence, m::GenericModel) = d.target

keys(m::GenericModel) = keys(m.vertices)
indexes(m::GenericModel) = indexes(m.vertices)

out_edges(v::Variable, m::GenericModel) = m.finclist[vertex_index(v, m)]
out_degree(v::Variable, m::GenericModel) = length(out_edges(v, m))

in_edges(v::Variable, m::GenericModel) = m.binclist[vertex_index(v, m)]
in_degree(v::Variable, m::GenericModel) = length(in_edges(v, m))

function add_vertex!(m::GenericModel, v::Variable, n::Integer=num_vertices(m)+1)
    push!(m.vertices, v)

    push!(m.finclist, Integer[])
    push!(m.binclist, Integer[])

    m.ofkey[v.key] = n

    v
end

function add_vertex!(m::GenericModel, vs::VariableVector, n::Integer=num_vertices(m)+1)
  nvertices = n
  for v in vs
    add_vertex!(m, v, nvertices)
    nvertices += 1
  end
end

function set_vertex!(m::GenericModel, v::Variable)
  m.vertices[v.index] = v

  push!(m.finclist, Integer[])
  push!(m.binclist, Integer[])

  m.ofkey[v.key] = v.index

  v
end

function set_vertex!(m::GenericModel, vs::VariableVector)
  for v in vs
    set_vertex!(m, v)
  end
end

make_edge(m::GenericModel, s::Variable, t::Variable) = Dependence(num_edges(m)+1, s, t)

function add_edge!(m::GenericModel, u::Variable, v::Variable, d::Dependence)
    ui = vertex_index(u, m)
    vi = vertex_index(v, m)

    push!(m.edges, d)
    push!(m.finclist[ui], d)
    push!(m.binclist[vi], d)

    if !is_directed(m)
        rev_d = revedge(d)
        push!(m.finclist[vi], rev_d)
        push!(m.binclist[ui], rev_d)
    end

    d
end

add_edge!(m::GenericModel, d::Dependence) = add_edge!(m, source(d, m), target(d, m), d)
add_edge!(m::GenericModel, u::Variable, v::Variable) = add_edge!(m, u, v, make_edge(m, u, v))

function GenericModel(vs::VariableVector, ds::DependenceVector=Dependence[]; isdirected::Bool=true, isindexed::Bool=true)
  n = length(vs)

  m = GenericModel(
    isdirected,
    Variable[],
    Dependence[],
    multivecs(Dependence, n),
    multivecs(Dependence, n),
    Dict{Symbol, Integer}()
  )
  if isindexed
    add_vertex!(m, sort_by_index(vs), 1)
  else
    add_vertex!(m, vs, 1)
    for i in 1:n
      m.vertices[i].index = i
    end
  end

  for d in ds
    add_edge!(m, d)
  end

  return m
end

function GenericModel(vs::VariableVector, ds::Vector{Dependence}, isdirected::Bool)
  n = length(vs)

  m = GenericModel(
    isdirected,
    Array{typeof(vs)}(n),
    Dependence[],
    multivecs(Dependence, n),
    multivecs(Dependence, n),
    Dict{Symbol, Integer}()
  )

  set_vertex!(m, vs)

  for d in ds
    add_edge!(m, d)
  end

  return m
end

GenericModel(isdirected::Bool=true) = GenericModel(Variable[], Dependence[], isdirected)

GenericModel(vs::Vector{V}, ds::Matrix{V}; isdirected::Bool=true, isindexed::Bool=true) where {V<:Variable} =
  GenericModel(vs, [Dependence(i, ds[i, 1], ds[i, 2]) for i in 1:size(ds, 1)], isdirected=isdirected, isindexed=isindexed)

function GenericModel(vs::Dict{Symbol, DataType}, ds::Dict{Symbol, Symbol}, isdirected::Bool=true)
  i = 0
  m = GenericModel(Variable[vs[k](k, i+=1) for k in keys(vs)], Dependence[], isdirected, true)

  i = 0
  for k in keys(ds)
    add_edge!(m, Dependence(i+=1, m[k], m[ds[k]]))
  end

  return m
end

function show(io::IO, model::GenericModel)
  isdirected = is_directed(model) ? "directed": "undirected"
  print(io, "GenericModel: $(num_vertices(model)) variables, $(num_edges(model)) dependencies ($isdirected graph)")
end

function model2dot(stream::IOStream, model::GenericModel)
  graphkeyword, edgesign = is_directed(model) ? ("digraph", "->") : ("graph", "--")
  dotindentation, dotspacing = "  ", " "

  write(stream, "$graphkeyword GenericModel {\n")

  for v in vertices(model)
    write(stream, string(dotindentation, v.key, dotspacing, "[shape=", dotshape(v), "]\n"))
  end

  for d in edges(model)
    write(stream, string(dotindentation, d.source.key, dotspacing, edgesign, dotspacing, d.target.key, "\n"))
  end

  write(stream, "}\n")
end

function model2dot(filename::AbstractString, model::GenericModel, mode::AbstractString="w")
  stream = open(filename, mode)
  model2dot(stream, model)
  close(stream)
end


text:
function MLJModelInterface.fit(model::EvoTypes, verbosity::Int, X, y)
    Xmatrix = MLJModelInterface.matrix(X)
    fitresult, cache = init_evotree(model, Xmatrix, y, verbosity = verbosity)
    grow_evotree!(fitresult, cache, verbosity = verbosity)
    report = nothing
    return fitresult, cache, report
end

function okay_to_continue(new, old)
    new.nrounds - old.nrounds >= 0 &&
    new.loss == old.loss &&
    new.λ == old.λ &&
    new.γ == old.γ &&
    new.max_depth  == old.max_depth &&
    new.min_weight == old.min_weight &&
    new.rowsample ==  old.rowsample &&
    new.colsample ==  old.colsample &&
    new.nbins ==  old.nbins &&
    new.α ==  old.α &&
    new.metric ==  old.metric
end

function MLJModelInterface.update(model::EvoTypes, verbosity::Integer, fitresult, cache, X, y)

    if okay_to_continue(model, cache.params)
        grow_evotree!(fitresult, cache, verbosity = verbosity)
    else
        Xmatrix = MLJModelInterface.matrix(X)
        fitresult, cache = init_evotree(model, Xmatrix, y, verbosity = verbosity)
        grow_evotree!(fitresult, cache, verbosity = verbosity)
    end
    report = nothing
    return fitresult, cache, report
end

function predict(model::EvoTreeRegressor, fitresult, Xnew)
    Xnew = MLJModelInterface.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return pred
end

function predict(model::EvoTreeClassifier, fitresult, Xnew)
    Xnew = MLJModelInterface.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return [MLJModelInterface.UnivariateFinite(fitresult.levels, pred[i,:]) for i in 1:size(pred,1)]
end

function predict(model::EvoTreeCount, fitresult, Xnew)
    Xnew = MLJModelInterface.matrix(Xnew)
    λ = predict(fitresult, Xnew)
    return [Distributions.Poisson(λᵢ) for λᵢ ∈ λ]
end

function predict(model::EvoTreeGaussian, fitresult, Xnew)
    Xnew = MLJModelInterface.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return [Distributions.Normal(pred[i]...) for i in 1:size(pred,1)]
end

# Metadata
const EvoTreeRegressor_desc = "Regression models with various underlying methods: least square, quantile, logistic."
const EvoTreeClassifier_desc = "Multi-classification with softmax and cross-entropy loss."
const EvoTreeCount_desc = "Poisson regression fitting λ with max likelihood."
const EvoTreeGaussian_desc = "Gaussian maximum likelihood of μ and σ²."

MLJModelInterface.metadata_pkg.((EvoTreeRegressor, EvoTreeClassifier, EvoTreeCount, EvoTreeGaussian),
    name="EvoTrees",
    uuid="f6006082-12f8-11e9-0c9c-0d5d367ab1e5",
    url="https://github.com/Evovest/EvoTrees.jl",
    julia=true,
    is_wrapper=false)

MLJModelInterface.metadata_model(EvoTreeRegressor,
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{<:MLJModelInterface.Continuous},
    weights=false,
    descr=EvoTreeRegressor_desc)

MLJModelInterface.metadata_model(EvoTreeClassifier,
        input=MLJModelInterface.Table(MLJModelInterface.Continuous),
        target=AbstractVector{<:MLJModelInterface.Finite},
        weights=false,
        descr=EvoTreeClassifier_desc)

MLJModelInterface.metadata_model(EvoTreeCount,
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{<:MLJModelInterface.Count},
    weights=false,
    descr=EvoTreeCount_desc)

MLJModelInterface.metadata_model(EvoTreeGaussian,
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{<:MLJModelInterface.Continuous},
    weights=false,
    descr=EvoTreeGaussian_desc)

# shared metadata
# MLJModelInterface.package_name(::Type{<:EvoTypes}) = "EvoTrees"
# MLJModelInterface.package_uuid(::Type{<:EvoTypes}) = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
# MLJModelInterface.package_url(::Type{<:EvoTypes}) = "https://github.com/Evovest/EvoTrees.jl"
# MLJModelInterface.is_pure_julia(::Type{<:EvoTypes}) = true
#
# MLJModelInterface.load_path(::Type{<:EvoTreeRegressor}) = "EvoTrees.EvoTreeRegressor"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeRegressor}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeRegressor}) = AbstractVector{<:MLJModelInterface.Continuous}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeCount}) = "EvoTrees.EvoTreeCount"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeCount}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeCount}) = AbstractVector{<:MLJModelInterface.Count}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeClassifier}) = "EvoTrees.EvoTreeClassifier"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeClassifier}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeClassifier}) = AbstractVector{<:MLJModelInterface.Finite}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeGaussian}) = "EvoTrees.EvoTreeGaussian"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeGaussian}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeGaussian}) = AbstractVector{<:MLJModelInterface.Continuous}

# function MLJ.clean!(model::EvoTreeRegressor)
#     warning = ""
#     if model.nrounds < 1
#         warning *= "Need nrounds ≥ 1. Resetting nrounds=1. "
#         model.nrounds = 1
#     end
#     if model.λ < 0
#         warning *= "Need λ ≥ 0. Resetting λ=0. "
#         model.λ = 0.0
#     end
#     if model.γ < 0
#         warning *= "Need γ ≥ 0. Resetting γ=0. "
#         model.γ = 0.0
#     end
#     if model.η <= 0
#         warning *= "Need η > 0. Resetting η=0.001. "
#         model.η = 0.001
#     end
#     if model.max_depth < 1
#         warning *= "Need max_depth ≥ 0. Resetting max_depth=0. "
#         model.max_depth = 1
#     end
#     if model.min_weight < 0
#         warning *= "Need min_weight ≥ 0. Resetting min_weight=0. "
#         model.min_weight = 0.0
#     end
#     if model.rowsample < 0
#         warning *= "Need rowsample ≥ 0. Resetting rowsample=0. "
#         model.rowsample = 0.0
#     end
#     if model.rowsample > 1
#         warning *= "Need rowsample <= 1. Resetting rowsample=1. "
#         model.rowsample = 1.0
#     end
#     if model.colsample < 0
#         warning *= "Need colsample ≥ 0. Resetting colsample=0. "
#         model.colsample = 0.0
#     end
#     if model.colsample > 1
#         warning *= "Need colsample <= 1. Resetting colsample=1. "
#         model.colsample = 1.0
#     end
#     if model.nbins > 250
#         warning *= "Need nbins <= 250. Resetting nbins=250. "
#         model.nbins = 250
#     end
#     return warning
# end


text:
module StochasticCubicDescent

export optimize!, ObjectiveFunction, GradientDescent, GDState, SCD, SCDState

abstract type Solver end
abstract type State end

struct ObjectiveFunction
    val::Function     # evaluates the objective function
    grad!::Function   # evaluates the gradient
    hv!::Function     # hessian-vector product (may be left undefined)
    ξ!::Function      # the noise (may be undefined for deterministic functions)
    ρ::AbstractFloat  # a Lipschitz constant for the Hessian
end

noop() = None
ObjectiveFunction(val, grad!) = ObjectiveFunction(val, grad!, noop, noop, -1.0)

# Gradient Descent (GD)

struct GDState <: State
    iter::Ref{Int}
    x::Vector
    grad::Vector

    GDState(x) = new(0, x, similar(x))
end

struct GradientDescent <: Solver
    η::AbstractFloat
    max_iters::Int
    store_history::Bool
    history::Vector{GDState}
end

function GradientDescent(η, max_iters)
    GradientDescent(η, max_iters, false, GDState[])
end

function optimize!(solver::GradientDescent,
                   state::GDState,
                   f::ObjectiveFunction)
    while state.iter[] < solver.max_iters
        f.grad!(state.x, state.grad)
        solver.store_history && push!(solver.history, deepcopy(state))
        state.iter[] += 1
        state.x .-= solver.η * state.grad
    end
end

# Stochastic Cubic Descent (SCD)

struct SCDState <: State
    iter::Ref{Int}
    x::Vector
    grad::Vector
    hv::Vector
    ξ::Vector

    # TODO: allow ξ to have a different dimension than x
    SCDState(x0) = new(0, x0, similar(x0), similar(x0), similar(x0))
end

struct SCD <: Solver
    max_iters::Int
    store_history::Bool
    history::Vector{SCDState}
end

function SCD(max_iters)
    SCD(max_iters, false, SCDState[])
end

function optimize!(solver::SCD, state::SCDState, f::ObjectiveFunction)
    while state.iter[] < solver.max_iters
        f.ξ!(state.x, state.ξ)
        f.grad!(state.x, state.ξ, state.grad)
        solver.store_history && push!(solver.history, deepcopy(state))

        β = f.ρ  # Here f.ρ is an upper bound on the β from Carmon-Duchi 2017
        R = β / (2f.ρ) + sqrt((β / 2f.ρ)^2 + norm(state.grad) / f.ρ)
        η = (4(β + f.ρ * R))^(-1)
        gd_solver = GradientDescent(η, 1000)

        Rc_lower_bound = -β / f.ρ + R
        gd_x0 = -Rc_lower_bound / norm(state.grad) * state.grad
        gd_state = GDState(gd_x0)

        val(x) = begin
            f.hv!(state.x + x, state.ξ, x, state.hv)
            0.5(x' * state.hv) + x' * state.grad + (f.ρ / 3) * norm(x)^3
        end
        grad!(x, grad) = begin
            f.hv!(state.x + x, state.ξ, x, state.hv)
            grad[:] = state.hv + state.grad + f.ρ * norm(x) * x
        end
        gd_subproblem = ObjectiveFunction(val, grad!)

        f.ξ!(state.x, state.ξ)
        optimize!(gd_solver, gd_state, gd_subproblem)
        state.x .+= gd_state.x

        state.iter[] += 1
    end
end

end # module


text:
function structure(x)
  fields = fieldnames(typeof(x))
  if isempty(fields)
    isbits(x) ?
      Row(typeof(x), Text(" "), x) :
      Row(typeof(x), Text("()"))
  else
    LazyTree(typeof(x), () -> [SubTree(Text("$f → "), structure(getfield′(x, f, UNDEF))) for f in fields])
  end
end

structure(xs::Vector) =
  LazyTree(Row(eltype(xs), fade("[$(length(xs))]")),
       () -> [isassigned(xs, i) ? structure(xs[i]) : UNDEF for i = 1:length(xs)])

structure(s::Symbol) = s
structure(s::Ptr) = s
structure(s::String) = s
# TODO: do this recursively
structure(x::Array) = x


text:
module Routes

immutable Route
 params::Array
 parts::Array
 regex::Regex
end

abstract Param
immutable NamedParam
 partindex::Int
 name::Symbol
 regex::Regex
end
immutable PositionalParam
 partindex::Int
 regex::Regex
end

immutable Part
 paramindex::Int
end

export route
"""
r1 = route(\"../RESULT/\",
    :day=>r\"\\d\\d\\d\\d-\\d\\d-\\d\\d\", \".gz\")

Create a route. Parameters may be:

    string - will joined as is,

    Pair(Symbol,Regex) - parameter for creating concrete file name

"""
function route(args...)
 params_arr = params(args...)
 parts_arr = [part(i,x,params_arr) for (i,x) in enumerate(args)]
 pattern_re = pattern(args...)
 Route(params_arr, parts_arr, pattern_re)
end

params(args...) = enumerate(args)|>
    e->[param(i,a) for(i,a) in e] |>
    arr->filter(p->p!=nothing, arr) |> collect

param(i::Int, p::Pair{Symbol,Regex}) = NamedParam(i, p[1], p[2])
param(i::Int, re::Regex) = PositionalParam(i, re)
param(i::Int, x)=nothing

part(i::Int,x::AbstractString,params) = x
part(i::Int,x,params) = findfirst(p->p.partindex==i, params)|> i->Part(i)

pattern(args...) = Regex(join(map(string,args),""))
import Base.string
string(p::Pair{Symbol,Regex}) = "(?P<$(p[1])>$(p[2].pattern))"
string(r::Regex) = "($(r.pattern))"

export file
"""
r1 = route(\"../RESULT\", :day=>r\"\\d\\d\\d\\d-\\d\\d-\\d\\d\",
    \".gz\") # named param

r1 = route(\"../RESULT\", r\"\\d\\d\\d\\d-\\d\\d-\\d\\d\",
    \".gz\") # positional param

f1 = file(r1, \"2016-08-09\") # \"../RESULT/2016-08-09.gz\"

r1.regex # r\"../RESULT/(?P<day>\d\d\d\d-\d\d-\d\d).gz\"

match(r1.regex,f1) # RegexMatch(\"../RESULT/2016-08-09.gz\", day=\"2016-08-09\")

"""
function file(r::Route, pars...)
 if length(pars)!=length(r.params)
    error("Count of params not equal. $(r.params|>length) != $pars")
 end
 join([ yeld(r,p,pars) for p in r.parts ],"")

end

check(r::Regex,s::AbstractString) = ismatch(r,s)
check(r::Regex,x) = check(r,string(x))

yeld(r::Route, p::AbstractString, pars) = p
function yeld(r::Route, p::Part, pars)
    userparam = pars[p.paramindex]
    routeparam = r.params[p.paramindex]
    if !check(routeparam.regex, userparam)
	error("Param #$(p.paramindex) (:$(routeparam)) must match with $(routeparam.regex)")
    end
    userparam
end

end # module

max_stars_repo_name: closescreen/Routes

max_issues_repo_name: closescreen/Routes

max_forks_repo_name: closescreen/Routes


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ 7681f7f4-5830-4860-a67d-c2acd1ed8218
using Plots

# ╔═╡ b79ef21d-e232-4c10-a0ba-6b4c54ad71aa
begin
	using PlutoUI
	almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))
	keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]));
	alert(text=md"Alert!") = Markdown.MD(Markdown.Admonition("danger", "Alert!", [text]));
	exercise(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("correct", "Exercise", [text]));
	warning(text=md"Warning.") = Markdown.MD(Markdown.Admonition("warning", "Warning!", [text]));
	correct(text=md"Great! You got the right answer! Let's move on to the next section.") = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]));
	hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));
end;

# ╔═╡ 69b42df0-8ae2-11eb-0bba-1983308b31c6
md"""
# Preliminary Weight, CG and Static Margin Estimation

"""

# ╔═╡ e76379e2-4d75-4b74-800d-8672941dc9a5
plotlyjs()

# ╔═╡ 8da00018-3cac-40d0-8d49-14a147467387
begin
	S_W=27#wing reference area

	#canard
	AR_H=0.9*4.5
	lambda_H=0#LE sweep
	S_H=8.26#reference area
	#l_H=5.265967537251188 #cg to local  cg

	#vertical stabilizer_wingtip
	#AR_V1=1.38
	#lambda_V1=36
	#S_V1=1.38
	#l_V1=7.755

	#Vertical stabilizer_fin
	#AR_V2=0.58
	#lambda_V2=45
	#S_V2=0.58
	#l_V2=7.2

	#other
	Mach=0.235
	Eta=0.97
	Eta_H=0.9

	#Fuselage
	L_F=9
	w_F=1.8
	S_F=33.9
end;

# ╔═╡ 0b5a24d0-8b00-11eb-224a-8b6fb1aff81a
md"""
## Estimating Mean Aerodynamic Chord Location
![](https://godot-bloggy.xyz/post/diagrams/WingParams.svg)

From the trapezoidal geometry:
```math
x_{40\%~\text{MAC}} = x_{\text{LE},\ \text{MAC}} + 0.4\bar{c}, \quad \text{where} \quad x_{\text{LE},\ \text{MAC}} = x_{\text{LE},\ \text{root}} + \bar Y\tan\Lambda_{\text{LE}}
```

"""

# ╔═╡ e8dad1e6-7477-42cc-b8e3-c518af7d9b52
begin
	# Vertical stabilizer
	S_v = 604
	lambda_v = deg2rad(44.4)
end;

# ╔═╡ a3b5d8b0-89e2-11eb-1040-031e0142ca7c
struct LiftingSurface
	c_root 		:: Float64
	c_tip 		:: Float64
    x_RLE 		:: Float64
    b 			:: Float64
    sweep_LE 	:: Float64
    taper 		:: Float64
end

# ╔═╡ 780dc91e-e6cf-4013-a04a-b35fa0fa9755
struct VaryingTaperLiftingSurface
	yB:: Float64
	x_rootLE :: Float64
	c_root:: Float64
	c_B :: Float64
	b :: Float64
	c_tip :: Float64
	taper_inner :: Float64
	taper :: Float64
	sweep_LEinner :: Float64
	sweep_LEouter :: Float64
	sweep_TE :: Float64
end

# ╔═╡ c82ef9d0-89e5-11eb-02d1-db611f26e0b6
md"""
Define a function for MAC with linear taper
```math
\bar{c}=\frac{2}{3}\left(\frac{1+\lambda+\lambda^2}{1+\lambda}\right)c_\text{root},~ \bar{Y}=\frac{b}{6}\left(\frac{1+2\lambda}{1+\lambda}\right)
```
"""

# ╔═╡ 3d325ff0-2a5e-4267-86f1-41c2cc25d7b1
md"""
Defined function for MAC for wings with cranked leading and trailing edges
```math
c_{\text{MAC}}=\frac{2}{3} c_{r} \frac{(\lambda_{i}^2 + \lambda_{i}\lambda+\lambda^2) + \left [ \frac{y_{B}}{(b/2 - y_{B})} \right ] (1+\lambda_{i}+\lambda_{i}^2) }{(\lambda_{i}^2+\lambda) + \left [ \frac{y_{B}}{(b/2 - y_{B})} \right ](1+\lambda_{i}^2)}
```
```math
y_{\text{MAC}}=\frac{1}{3} \left ( \frac{b}{2} -y_{B}\right ) \frac{\lambda_{i} \left (1 + 2 \lambda_{i} \lambda \right) + \left [ \frac{y_{B}}{(b/2) -y_{B}} \right ] (2+\lambda_{i}\lambda ) + \left [ \frac{y_{B}}{(b/2) -y_{B}} \right ]^2 (1+2\lambda_{i})}{(\lambda_{i}+\lambda)+\left [ \frac{y_{B}}{(b/2) -y_{B}} \right ](1+\lambda_{i})}
```
The above equation is taken from https://www.sciencedirect.com/topics/engineering/mean-aerodynamic-chord

"""

# ╔═╡ 137918a0-8b06-11eb-0f23-2fd22e086911
md"The function defined for MAC along the chord for linear tapered wing"

# ╔═╡ 3053c570-89e3-11eb-1f51-8fb42324f651
mean_aerodynamic_chord(surf :: LiftingSurface) = 2/3 * (1 + surf.taper + surf.taper^2)/(1 + surf.taper) * surf.c_root

# ╔═╡ 11ddb3ae-db3a-41b0-a77d-387f99b0e9a0
md"The function defined for MAC along the chord for cranked LE & TE wing"

# ╔═╡ 09f81461-7713-41fc-a3b9-cc106018ccf5
c_MAC(surf :: VaryingTaperLiftingSurface)=(2/3)*surf.c_root*(((surf.taper_inner^2)+(surf.taper_inner*surf.taper)+(surf.taper^2)+((surf.yB/((surf.b/2)-surf.yB))*(1+surf.taper_inner+(surf.taper_inner^2))))/(surf.taper+(surf.taper_inner^2)+((surf.yB/((surf.b/2)-surf.yB))*(1+surf.taper_inner^2))))

# ╔═╡ 3a8b8e80-89e6-11eb-255d-df692399569f
md"Define a function for MAC for linear-tapered wing along the span-wise direction"

# ╔═╡ 68795e80-89e6-11eb-04da-c9e66888d483
Y_mac(surf :: LiftingSurface) = surf.b * 1/6 * (1 + 2 * surf.taper) / (1 + surf.taper)

# ╔═╡ 3a4723e4-ec74-4d59-b1cd-4d1e99f2434e
md"Define a function for MAC for cranked LE & TE wing along the span-wise direction"

# ╔═╡ 28ad10c1-678a-45e5-85b7-772a99108f0c
yMAC(surf :: VaryingTaperLiftingSurface)=(1/3)*((surf.b/2)-surf.yB)*(((surf.taper_inner*(1+(2*surf.taper_inner*surf.taper)))+((surf.yB/((surf.b/2)-surf.yB))*(2+(surf.taper_inner*surf.taper)))+((surf.yB/((surf.b/2)-surf.yB))^2)*(1+(2*surf.taper_inner)))/(surf.taper_inner+surf.taper+((surf.yB/((surf.b/2)-surf.yB))*(1+surf.taper_inner))))

# ╔═╡ 11a520d0-89e6-11eb-2337-a5c4040b8013
md"Define the values of the wing for its struct."

# ╔═╡ be1b261a-be10-4c60-adad-36614e7b5827
begin
	#Main wing with LERX
	yB = 2.24#y loc at where it changes sweep
	x_rootLE= 4.252#distance of wing to nose at root
	c_root = 4#root chord
	c_B = 2.1#chord at where it changes sweep
	b_W = 12.6#span of wing
	c_tip = 0.96#chord at tip
	taper_inner = 0.525
	taper_W = 0.24#overall taper
	sweep_LEinner = 56
	sweep_LEouter = 26
	sweep_TE = 12
end;

# ╔═╡ c3824140-28d3-4f0d-8d56-3dbf42dca2cf
function getTaperSweepofMAC(ymac)
	if ymac<yB
		relSweep=sweep_LEinner
	elseif ymac>yB
		relSweep=sweep_LEouter
	elseif ymac>b_W
		relSweep=0
	end
	relSweep
end

# ╔═╡ 12d99410-89e4-11eb-1bbb-6b20dcf1678a
begin

	#canards
	croot_h = 0.9
	ctip_h = 0.9
	xRLE_h = 0.675#root LE longitudinal distance from nose
	b_h = 5
	sweepLE_h = 0
	taper_h = ctip_h/croot_h
	S_c=croot_h*b_h

	#vertical stabilizer at wingtip #concerning wingtip and span
	croot_V1 = 1.65
	ctip_V1 = 0.45
	xRLE_V1 = x_rootLE+(yB*tan(sweep_LEinner))+((b_W-yB)*tan(sweep_LEouter))
	b_V1 = 2
	sweepLE_V1 = 36
	taper_V1 = ctip_V1/croot_V1
	S_wingtip=(((croot_V1+ctip_V1)*b_V1)/2)*2

	#vertical stabilizer at fin
	croot_V2=0.58
	ctip_V2=0.45
	xRLE_V2=7.2
	b_V2=0.42
	sweepLE_V2=45
	taper_V2=ctip_V2/croot_V2
	S_aftFin=((croot_V2+ctip_V2)*b_V2)/2
end;

# ╔═╡ 608bc750-8b06-11eb-2123-fb1b2a57a252
md"Structs of type ```LiftingSurface``` are created by applying the ```LiftingSurface``` type as a **function** with inputs for its fields:"

# ╔═╡ 30fe0cc4-7e57-49bc-89b2-ce3b00cc2a65
canard=LiftingSurface(croot_h,ctip_h,xRLE_h,b_h,sweepLE_h,taper_h);

# ╔═╡ 9e1c4b9c-4352-47e0-be98-65cc599474ac
vStab=LiftingSurface(croot_V1,ctip_V1,xRLE_V1,b_V1,sweepLE_V1,taper_V1);

# ╔═╡ cb11e4d7-233f-4779-b9ad-7059bff013c9
finStab=LiftingSurface(croot_V2,ctip_V2,xRLE_V2,b_V2,sweepLE_V2,taper_V2);

# ╔═╡ b356c994-2314-491d-a93d-b5625be04036
Mainwing = VaryingTaperLiftingSurface(yB, x_rootLE, c_root,c_B,b_W,c_tip,taper_inner,taper_W,sweep_LEinner,sweep_LEouter,sweep_TE);

# ╔═╡ 2a7232ae-70a5-4bee-9b63-6d9903701a40
mainWingMAC_c = c_MAC(Mainwing) #Main wing chord MAC

# ╔═╡ f90ac96f-200c-4c00-bb15-f24074865da4
canardMAC_c=mean_aerodynamic_chord(canard)#canard MAC

# ╔═╡ 463f1e3b-365a-4311-a4db-bdbefa66564d
vStabMAC_c=mean_aerodynamic_chord(vStab)#wingtip v. stabilizer MAC

# ╔═╡ 5072dbe8-a88e-4891-8cb6-6a765d055dae
finStabMAC_c=mean_aerodynamic_chord(finStab)#aft fin MAC

# ╔═╡ d76e65a8-9e8f-4665-8fa8-52382d8dfbb2
mainWingMAC_Y=yMAC(Mainwing) #Main wing MAC y.loc, within inboard wing

# ╔═╡ 45d4c80a-5798-4451-a591-c1811ba15fed
canardMAC_Y=Y_mac(canard)#canard MAC y.loc

# ╔═╡ fc3217ef-7fc1-433a-a9b2-84cfda8f04d8
vStabMAC_Y=Y_mac(vStab)+(b_W/2)#vertical stabilizer MAC y.loc

# ╔═╡ 2d4a5450-8b07-11eb-3b69-3dc27e505c50
md"Now you can obtain the CG location of the wing with respect to the reference mean aerodynamic chord location."

# ╔═╡ afaaa1f0-8a0f-11eb-00dc-47d7ecc444a7
wing_mac_origin(surf :: LiftingSurface, percent) = surf.x_RLE + Y_mac(surf) * tan(surf.sweep_LE) + percent / 100 * mean_aerodynamic_chord(surf)

# ╔═╡ 142c524a-3917-430d-aab8-1469dc8effa8
wingMAC_loc_origin(wing :: VaryingTaperLiftingSurface, perc) = wing.x_rootLE + yMAC(wing)*tand(wing.sweep_LEinner) + perc/100*c_MAC(wing)

# ╔═╡ aee47142-8ce8-4189-9bfd-3ee51fbb4cb1
Mainwing_cg40 = wingMAC_loc_origin(Mainwing,40)

# ╔═╡ 197916ee-fe63-4df2-869b-81c440d5ff84
Mainwing_cg25 = wingMAC_loc_origin(Mainwing,25)

# ╔═╡ 3e05470c-1161-4adc-9607-3638368a14d6
Canard_cg40=wing_mac_origin(canard,40)

# ╔═╡ 10ca0d20-b99e-4b6e-a6b3-ecb95d1bf1b0
Canard_cg25=wing_mac_origin(canard,25)

# ╔═╡ 2c683bb7-83d6-409d-850b-9f622cfeefe9
vStab_cg40=wing_mac_origin(vStab,40)

# ╔═╡ f1ab1278-4ab5-423c-a773-5d2e874d06f7
vStab_cg25=wing_mac_origin(vStab,25)

# ╔═╡ b9c58e91-11e8-47db-a331-14691993f212
finStab_cg40=wing_mac_origin(finStab,40)

# ╔═╡ 7a95efd5-9794-43cb-b6cc-970c0899f746
finStab_cg25=wing_mac_origin(finStab,25)

# ╔═╡ 9c13d810-8a13-11eb-07a5-3bcedbd77eac
md"""
## Center of Gravity Estimation

The position of the aircraft’s center of gravity (CG) is given by:
```math
\mathbf{x}_\text{cg} = \frac{\sum_i W_i \ (\mathbf{x}_{\text{cg}})_i}{\sum_i W_i}, \quad \mathbf{x} = \begin{bmatrix}
  x \\ y \\ z
\end{bmatrix}
```

where $W_i$ represents the weight for each component and $\mathbf x_{\text{cg}_i}$ is the longitudinal distance between the origin and the CG of the $i$th component. The product in the form of $W_i \mathbf x_{\text{cg}_i}$ is also referred to as the moment of the $i$th component.

Considering a takeoff gross weight of 2655 kg and using Raymer's reference data, the specific weight and position parameters of the structural components are shown in the following table:

Components | Weight Ratio | Reference Weight (kg) | Approximate Location
:-------- | :-----: | :----------:|----------:
Nose landing gear | 0.057 * 15% - W * 0.014 | 2655 | Centroid
Main landing gear | 0.057 * 85% - W * 0.014 | 2655 | Centroid
Installed engine     | 1.4 | 208.7     | Centroid
“All-else empty”    | 0.1  | 2655     | 40-50% Length

Note: We consider the **nose** as the reference point and **clockwise moments** as positive!
"""

# ╔═╡ 9b2da92a-d97c-485b-ab36-bc775a65a159
begin
	# Filled in quantities
	TOGW = 2858 	 # kg
	w_engine = 363*1.4 # kg
	reserveFuel=0.06
	FuelWF=0.31135388287326593
	TrappedFuel=TOGW*FuelWF*reserveFuel#kg
	UsableFuel=TOGW*FuelWF/(reserveFuel+1)#kg

	#weight of Components
	w_MainWing=278.55
	w_fuselage=275.8
	w_canard=82.7475#67.0255
	w_VtailTip=43.99
	w_VtailFin=13.529
	w_MainLG=(0.057*TOGW-(TOGW*0.014))*0.85
	w_NoseLG=(0.057*TOGW-(TOGW*0.014))*0.15
	w_SeatRow1=30*2
	w_SeatRow2=30*3
	w_Crew=90
	w_passenger1=90*1
	w_passenger2=90*3
	w_BRS=40

	#cg length
	cg_fuselage=4.05
	cg_MainLG=5.625
	cg_NoseLG=0.825
	cg_SeatRow1=3.825
	cg_SeatRow2=5.325
	cg_engine=8.040
	cg_fuel=5.79
	cg_BRS=6

end;

# ╔═╡ 7a540330-8a13-11eb-3e55-c15452e5004c
md"It's convenient to use a function for the moment calculations."

# ╔═╡ ea34f410-8a14-11eb-1556-3b8c3ea0307d
moment(dist, weight) = dist * weight

# ╔═╡ 00d09cb2-8a15-11eb-0638-dffe4c995856
md"For the previously generated wing, the total moment is:"

# ╔═╡ a9ab5c00-8a18-11eb-1992-93bb7c57c029
md"The weight and CG position of each component can hence be computed and included in a dictionary for convenience in calculations."

# ╔═╡ 89bf9247-00c9-4b30-b328-6b1777d474a7
weight_pos = Dict(	"MainWing"   => (w_MainWing,Mainwing_cg40),
						"fuselage"  => (w_fuselage,cg_fuselage),
        				"canard"  => (w_canard,Canard_cg40),
        				"vTailTip"   => (w_VtailTip,vStab_cg40),
        				"vTailFin" => (w_VtailFin,finStab_cg40),
        				"mainLG" => (w_MainLG,cg_MainLG),
        				"noseLG" => (w_NoseLG,cg_NoseLG),
						"SeatRow_1" =>(w_SeatRow1,cg_SeatRow1),
						"SeatRow_2" => (w_SeatRow2,cg_SeatRow2),
						"engine"=> (w_engine,cg_engine),
						"trappedFuel" =>(TrappedFuel,cg_fuel),
						"Crew" => (w_Crew,cg_SeatRow1),
						"usableFuel"=> (UsableFuel,cg_fuel),
						"Passenger1" => (w_passenger1,cg_SeatRow1),
						"Passenger2" => (w_passenger2,cg_SeatRow2),
						"BRS" => (w_BRS,cg_BRS),
						);

# ╔═╡ 4175e400-8a1a-11eb-305d-bf6a7ad893ed
md"Now we can calculate the total moments generated from all the component, i.e., $\sum_i W_i \mathbf x_{\text{cg}_i}$"

# ╔═╡ 5b9711fa-4950-4512-9bec-14d08f594324
moment_list = [ moment(weight, pos_x) for (weight, pos_x) in values(weight_pos) ];

# ╔═╡ b63575a7-d350-496e-a0a7-9653413bfb62
momentSum = sum(moment_list)

# ╔═╡ 84906350-8a1a-11eb-26b6-2324a98f390b
md"The same logic also applies to the total weight, i.e., $\sum_i W_i$"

# ╔═╡ 934c246e-b0b7-4421-8545-a9e7b3d16bd3
tot_weight = sum(weight for (weight, pos_x) in values(weight_pos))

# ╔═╡ a8e2beb4-817f-4de7-a883-f6777c2ab73c
xCG = momentSum/tot_weight

# ╔═╡ 3a64faf2-56cd-438a-a187-a213bcdfcd14
md"""## CG excursion calculation"""

# ╔═╡ edf063b2-54ea-46f5-9d13-854e5014657d
weight_Empty = Dict(	"MainWing"   => (w_MainWing,Mainwing_cg40),
						"fuselage"  => (w_fuselage,cg_fuselage),
        				"canard"  => (w_canard,Canard_cg40),
        				"vTailTip"   => (w_VtailTip,vStab_cg40),
        				"vTailFin" => (w_VtailFin,finStab_cg40),
        				"mainLG" => (w_MainLG,cg_MainLG),
        				"noseLG" => (w_NoseLG,cg_NoseLG),
						"SeatRow_1" =>(w_SeatRow1,cg_SeatRow1),
						"SeatRow_2" => (w_SeatRow2,cg_SeatRow2),
						"engine"=> (w_engine,cg_engine),
						);

# ╔═╡ 34c83825-cca7-430c-a70c-09cd85d99852
WeightComponent_fuel=Dict("usableFuel"=> (UsableFuel,cg_fuel));

# ╔═╡ e01584c3-1a50-4b0e-8f6c-008312e6e945
weightComponent_OPEmpty=Dict("trappedFuel" =>(TrappedFuel,cg_fuel),
						"Crew" => (w_Crew,cg_SeatRow1));

# ╔═╡ a5ba9d65-2ba9-4328-ad81-ac426f91320c
WeightComponent_payload=Dict("Passenger1" => (w_passenger1,cg_SeatRow1),
						"Passenger2" => (w_passenger2,cg_SeatRow2),
						"BRS" => (w_BRS,cg_BRS));

# ╔═╡ 1dbe47ad-952c-4074-aac3-80e643099531
emptyW_component = sum(weight for (weight, pos_x) in values(weight_Empty))

# ╔═╡ bf9d769c-1918-46a2-904c-199da2adb044
OpEmpt_component=sum(weight for (weight, pos_x) in values(weightComponent_OPEmpty))

# ╔═╡ d9cad966-9d5f-4c1f-a75a-56a66cfb4c98
fuel_component = UsableFuel

# ╔═╡ af74f497-4752-42fd-9a85-9618b787df96
payload_component = sum(weight for (weight, pos_x) in values(WeightComponent_payload))

# ╔═╡ 6df33c5b-fe54-480c-9d57-a2e806be631e
sum_moments_emptyComp = sum(moment(weight, pos_x) for (weight, pos_x) in values(weight_Empty));

# ╔═╡ e241a14b-8340-45dd-9bc7-04d2b426198c
sum_moments_OpEmpt_Comp = sum(moment(weight, pos_x) for (weight, pos_x) in values(weightComponent_OPEmpty))

# ╔═╡ df0597a3-492b-4735-9497-90179942a882
sum_moments_fuel_Comp = sum(moment(fuel_component, 4.8));

# ╔═╡ 7bab3e29-2558-4a20-b092-6e8921e427d5
sum_moments_payload = sum(moment(weight, pos_x) for (weight, pos_x) in values(WeightComponent_payload));

# ╔═╡ 3d3f2ec4-299e-4936-9fec-18d0542e273d
x_cg_Empt_Comp = sum_moments_emptyComp/emptyW_component

# ╔═╡ 6a8e4835-5643-4b3c-b9b5-28e596cd1170
x_cg_OpEmpt_Comp = sum_moments_OpEmpt_Comp/OpEmpt_component

# ╔═╡ e64398d7-6aeb-40b5-b774-e9a7ef50ab79
x_cg_fuel_Comp = sum_moments_fuel_Comp/fuel_component

# ╔═╡ 6df55772-ea54-48f0-a577-bb0a372dcffa
x_cg_Passenger = sum_moments_payload/payload_component

# ╔═╡ 33a141a5-3e02-419a-be64-5d654578982d
begin
	moment1=sum_moments_emptyComp
	moment2=sum_moments_OpEmpt_Comp
	moment3=sum_moments_fuel_Comp
	moment4=sum_moments_payload

	x_cg1=x_cg_Empt_Comp
	x_cg2=x_cg_OpEmpt_Comp
	x_cg3=x_cg_fuel_Comp
	x_cg4=x_cg_Passenger
end;

# ╔═╡ 61d42a8a-cb21-4f10-a310-66a06d3b2565
function excursion(m1,m2,m3,m4,x1,x2,x3,x4)
	excursion_xcg=[m1/(m1/x1),(m2+m1)/((m2/x2)+(m1/x1)),(m3+m2+m1)/((m3/x3)+(m2/x2)+(m1/x1)),(m4+m3+m2+m1)/((m4/x4)+(m3/x3)+(m2/x2)+(m1/x1))]
	excursion_weight=[m1/x1,(m2/x2)+(m1/x1),(m3/x3)+(m2/x2)+(m1/x1),(m4/x4)+(m3/x3)+(m2/x2)+(m1/x1)]
	excursion_xcg,excursion_weight
end

# ╔═╡ 30b2ad6c-88bc-48f0-a436-4e155cff23f7
#Case 1 - 1234
case1cg,case1weight=excursion(moment1,moment2,moment3,moment4,x_cg1,x_cg2,x_cg3,x_cg4);

# ╔═╡ 0ed8b8b0-0674-4eff-bd68-d05f9b8c2a33
#Case 2 - 1243
case2cg,case2weight=excursion(moment1,moment2,moment4,moment3,x_cg1,x_cg2,x_cg4,x_cg3);

# ╔═╡ f589ca20-8af5-4164-8fba-1124ff131d35
#Case 3 - 1324
case3cg,case3weight=excursion(moment1,moment3,moment2,moment4,x_cg1,x_cg3,x_cg2,x_cg4);

# ╔═╡ f000648b-a70e-4600-b2f7-0af98922da15
#Case 4 - 1342
case4cg,case4weight=excursion(moment1,moment3,moment4,moment2,x_cg1,x_cg3,x_cg4,x_cg2);

# ╔═╡ c4748573-4944-4000-a21e-e4254dac6f07
#Case 5 - 1423
case5cg,case5weight=excursion(moment1,moment4,moment2,moment3,x_cg1,x_cg4,x_cg2,x_cg3);

# ╔═╡ 872a4079-7ef5-4501-8326-7f542fa671ed
#Case 6 - 1432
case6cg,case6weight=excursion(moment1,moment4,moment3,moment2,x_cg1,x_cg4,x_cg3,x_cg2);

# ╔═╡ 1f3a5f6a-d99f-48b3-8d62-ef2b4cfb4cfb
begin
	cg_excursionPlot = plot(ylabel = "Weight, kg", xlabel = "Longitudinal distance from the tip of the nose, m",
		 			title = "CG Excursion Plot")
	plot!(case1cg, case1weight, label = "1234")
	plot!(case2cg, case2weight, label = "1243")
	plot!(case3cg, case3weight, label = "1324")
	plot!(case4cg, case4weight, label = "1342")
	plot!(case5cg, case5weight, label = "1423")
	plot!(case6cg,case6weight, label="1432")
end

# ╔═╡ 7eae14e0-8a1b-11eb-15f0-372c8c42804b
md"""## Static Margin Estimation

The forces depicted are the largest contributors to the longitudinal stability characteristics.

The static margin equation in a slightly modified form:
```math
\text{Static Margin} = \frac{l_h S_h}{\bar{c}S_w} \frac{C_{L_{\alpha_h}}}{C_{L_{\alpha_w}}} - \frac{\partial C_{m_{fus}}}{\partial C_L} - \frac{(x_{cg} - x_{25\%\ \text{MAC}})}{\bar{c}}
```
4 parameters are unknown after the sizing procedure:
1. The lift curve slope for the wing $C_{L_{\alpha_w}}$
2. The lift curve slope of the horizontal stabilizer $C_{L_{\alpha_h}}$
3. Derivative of pitching moment of fuselage (including other components) with respect to $C_L$ $\frac{\partial C_{m_{fus}}}{\partial C_L}$
4. Distance between wing aerodynamic center and the CG  $x_{cg} - x_{25\%\ \text{MAC}}$.

"""

# ╔═╡ d7449777-49c3-4acc-af6a-523fc6e67231
Canard_staticMargin(l_h,S_c,CLa_c,cbar,S_w,CLA_w,CmCL,x_cg,x_ac)=((l_h*S_c*CLa_c)/(cbar*S_w*CLA_w))+(CmCL)-((x_cg-x_ac)/cbar)

# ╔═╡ b89f03d2-0b78-45a9-9644-d5c1eb16f29d
md""" #### Downwash effect
"""

# ╔═╡ c77332d9-98eb-4cc5-a1c7-7fa60dd14dc3
begin
	downwashWing=0.41
	corrFactorWingDownwash=0.129
	CLwingIsolate=0.0688
	CLcanardIsolate=0.0693
	upwashCanard=0.05
end;

# ╔═╡ b620d985-d81b-466e-a489-6d6bf505f0c8
CLwing=CLwingIsolate*(1-(downwashWing*corrFactorWingDownwash))

# ╔═╡ 38fc9c2a-c84f-430e-aa2e-2b14e5cbe112
ClCanard_upwash=CLcanardIsolate*(1+upwashCanard)*S_c/S_W

# ╔═╡ 11dd1622-50c4-404e-88a2-c29c02760ca9
CLCanardWing=ClCanard_upwash+CLwing

# ╔═╡ e258ad9e-a39b-45c8-997b-d785977116fd
md"""### Fuselage Pitching Moment (Newton's Divided Differences Interpolation Polynomial) """

# ╔═╡ 55a66ce9-2c77-42dd-8aac-446cddb1274d
FusePitch(K_fus,w_fus,L_fus,S_w,cbar,CLA_w)=(K_fus*L_fus*(w_fus)^2)/(S_w*cbar*CLA_w)

# ╔═╡ fc926801-2060-4c9f-b103-466c30cb899a
getFuse(wingloc)=0.115+(0.57*(wingloc-0.1))+(5.75*(wingloc-0.1)*(wingloc-0.2))+((-24)*(wingloc-0.1)*(wingloc-0.2)*(wingloc-0.3))+(96.25*(wingloc-0.1)*(wingloc-0.2)*(wingloc-0.3)*(wingloc-0.4))+((-314.167)*(wingloc-0.1)*(wingloc-0.2)*(wingloc-0.3)*(wingloc-0.4)*(wingloc-0.5))+(890.2778*(wingloc-0.1)*(wingloc-0.2)*(wingloc-0.3)*(wingloc-0.4)*(wingloc-0.5)*(wingloc-0.6))

# ╔═╡ 5a22ad68-2241-4c8e-bfc5-bbaf9bf786e2
md"""### Tail volume parameters"""

# ╔═╡ cccce8d6-afa5-484f-8871-9106e71da354
#define moment arms
begin
	l_canard=abs(Mainwing_cg25-Canard_cg25)
	l_wingtip=abs(vStab_cg25-Mainwing_cg25)
	l_aftFin=abs(finStab_cg25-Mainwing_cg25)
end;

# ╔═╡ 84aedc8c-5a59-4e7a-8136-d209b16984bd
C_VT(L_vt,S_vt,b_w,S_w)=(L_vt*S_vt)/(b_w*S_w)

# ╔═╡ 87ad6166-a5f7-4fdd-9ca3-85ec1ac62fff
C_HT(L_ht,S_ht,cbar_w,S_w)=(L_ht*S_ht)/(cbar_w*S_w)

# ╔═╡ 217f7fb1-0cc4-4ead-8998-2fc5c13571a9
Volume_wingtip=C_VT(l_canard,S_wingtip,b_W,S_W)

# ╔═╡ 97c7bd4b-f860-4e66-bd5a-a901f9c014ee
Volume_aftFin=C_VT(l_aftFin,S_aftFin,b_W,S_W)

# ╔═╡ 8164635b-1ad6-4b74-8926-b90d6182cc69
Volume_Canard=C_HT(l_canard,S_c,mainWingMAC_c,S_W)

# ╔═╡ 995fb655-6e01-452b-91fd-93b5cb89a17e
md"""#### Value Checks"""

# ╔═╡ 67fc7cd4-afd1-4812-9e34-41527ca8e68a
testCmCL=FusePitch(getFuse(Mainwing_cg25/L_F),w_F,L_F,S_W,mainWingMAC_c,CLwing)

# ╔═╡ 8944cab7-3a9e-4573-93cc-3f605fcb545a
SM1=Canard_staticMargin(case1cg[1]-Canard_cg40,S_c,ClCanard_upwash,mainWingMAC_c,S_W,CLwing,testCmCL,case1cg[1],Mainwing_cg40)

# ╔═╡ d6562dfe-6656-494b-aae7-9807aec745f9
neutralPoint=(SM1/100)*mainWingMAC_c+case1cg[1]

# ╔═╡ Cell order:
# ╟─69b42df0-8ae2-11eb-0bba-1983308b31c6
# ╠═e76379e2-4d75-4b74-800d-8672941dc9a5
# ╠═7681f7f4-5830-4860-a67d-c2acd1ed8218
# ╠═8da00018-3cac-40d0-8d49-14a147467387
# ╟─0b5a24d0-8b00-11eb-224a-8b6fb1aff81a
# ╟─e8dad1e6-7477-42cc-b8e3-c518af7d9b52
# ╠═a3b5d8b0-89e2-11eb-1040-031e0142ca7c
# ╠═780dc91e-e6cf-4013-a04a-b35fa0fa9755
# ╠═c3824140-28d3-4f0d-8d56-3dbf42dca2cf
# ╟─c82ef9d0-89e5-11eb-02d1-db611f26e0b6
# ╟─3d325ff0-2a5e-4267-86f1-41c2cc25d7b1
# ╟─137918a0-8b06-11eb-0f23-2fd22e086911
# ╠═3053c570-89e3-11eb-1f51-8fb42324f651
# ╟─11ddb3ae-db3a-41b0-a77d-387f99b0e9a0
# ╠═09f81461-7713-41fc-a3b9-cc106018ccf5
# ╟─3a8b8e80-89e6-11eb-255d-df692399569f
# ╠═68795e80-89e6-11eb-04da-c9e66888d483
# ╟─3a4723e4-ec74-4d59-b1cd-4d1e99f2434e
# ╠═28ad10c1-678a-45e5-85b7-772a99108f0c
# ╟─11a520d0-89e6-11eb-2337-a5c4040b8013
# ╠═12d99410-89e4-11eb-1bbb-6b20dcf1678a
# ╠═be1b261a-be10-4c60-adad-36614e7b5827
# ╟─608bc750-8b06-11eb-2123-fb1b2a57a252
# ╠═30fe0cc4-7e57-49bc-89b2-ce3b00cc2a65
# ╠═9e1c4b9c-4352-47e0-be98-65cc599474ac
# ╠═cb11e4d7-233f-4779-b9ad-7059bff013c9
# ╠═b356c994-2314-491d-a93d-b5625be04036
# ╠═2a7232ae-70a5-4bee-9b63-6d9903701a40
# ╠═f90ac96f-200c-4c00-bb15-f24074865da4
# ╠═463f1e3b-365a-4311-a4db-bdbefa66564d
# ╠═5072dbe8-a88e-4891-8cb6-6a765d055dae
# ╠═d76e65a8-9e8f-4665-8fa8-52382d8dfbb2
# ╠═45d4c80a-5798-4451-a591-c1811ba15fed
# ╠═fc3217ef-7fc1-433a-a9b2-84cfda8f04d8
# ╟─2d4a5450-8b07-11eb-3b69-3dc27e505c50
# ╠═afaaa1f0-8a0f-11eb-00dc-47d7ecc444a7
# ╠═142c524a-3917-430d-aab8-1469dc8effa8
# ╠═aee47142-8ce8-4189-9bfd-3ee51fbb4cb1
# ╠═197916ee-fe63-4df2-869b-81c440d5ff84
# ╠═3e05470c-1161-4adc-9607-3638368a14d6
# ╠═10ca0d20-b99e-4b6e-a6b3-ecb95d1bf1b0
# ╠═2c683bb7-83d6-409d-850b-9f622cfeefe9
# ╠═f1ab1278-4ab5-423c-a773-5d2e874d06f7
# ╠═b9c58e91-11e8-47db-a331-14691993f212
# ╠═7a95efd5-9794-43cb-b6cc-970c0899f746
# ╟─9c13d810-8a13-11eb-07a5-3bcedbd77eac
# ╠═9b2da92a-d97c-485b-ab36-bc775a65a159
# ╟─7a540330-8a13-11eb-3e55-c15452e5004c
# ╠═ea34f410-8a14-11eb-1556-3b8c3ea0307d
# ╟─00d09cb2-8a15-11eb-0638-dffe4c995856
# ╟─a9ab5c00-8a18-11eb-1992-93bb7c57c029
# ╠═89bf9247-00c9-4b30-b328-6b1777d474a7
# ╟─4175e400-8a1a-11eb-305d-bf6a7ad893ed
# ╠═5b9711fa-4950-4512-9bec-14d08f594324
# ╠═b63575a7-d350-496e-a0a7-9653413bfb62
# ╟─84906350-8a1a-11eb-26b6-2324a98f390b
# ╠═934c246e-b0b7-4421-8545-a9e7b3d16bd3
# ╠═a8e2beb4-817f-4de7-a883-f6777c2ab73c
# ╟─3a64faf2-56cd-438a-a187-a213bcdfcd14
# ╠═edf063b2-54ea-46f5-9d13-854e5014657d
# ╠═34c83825-cca7-430c-a70c-09cd85d99852
# ╠═e01584c3-1a50-4b0e-8f6c-008312e6e945
# ╠═a5ba9d65-2ba9-4328-ad81-ac426f91320c
# ╠═1dbe47ad-952c-4074-aac3-80e643099531
# ╠═bf9d769c-1918-46a2-904c-199da2adb044
# ╠═d9cad966-9d5f-4c1f-a75a-56a66cfb4c98
# ╠═af74f497-4752-42fd-9a85-9618b787df96
# ╠═6df33c5b-fe54-480c-9d57-a2e806be631e
# ╠═e241a14b-8340-45dd-9bc7-04d2b426198c
# ╠═df0597a3-492b-4735-9497-90179942a882
# ╠═7bab3e29-2558-4a20-b092-6e8921e427d5
# ╠═3d3f2ec4-299e-4936-9fec-18d0542e273d
# ╠═6a8e4835-5643-4b3c-b9b5-28e596cd1170
# ╠═e64398d7-6aeb-40b5-b774-e9a7ef50ab79
# ╠═6df55772-ea54-48f0-a577-bb0a372dcffa
# ╠═33a141a5-3e02-419a-be64-5d654578982d
# ╠═61d42a8a-cb21-4f10-a310-66a06d3b2565
# ╠═30b2ad6c-88bc-48f0-a436-4e155cff23f7
# ╠═0ed8b8b0-0674-4eff-bd68-d05f9b8c2a33
# ╠═f589ca20-8af5-4164-8fba-1124ff131d35
# ╠═f000648b-a70e-4600-b2f7-0af98922da15
# ╠═c4748573-4944-4000-a21e-e4254dac6f07
# ╠═872a4079-7ef5-4501-8326-7f542fa671ed
# ╠═1f3a5f6a-d99f-48b3-8d62-ef2b4cfb4cfb
# ╟─7eae14e0-8a1b-11eb-15f0-372c8c42804b
# ╠═d7449777-49c3-4acc-af6a-523fc6e67231
# ╟─b89f03d2-0b78-45a9-9644-d5c1eb16f29d
# ╠═c77332d9-98eb-4cc5-a1c7-7fa60dd14dc3
# ╠═b620d985-d81b-466e-a489-6d6bf505f0c8
# ╠═38fc9c2a-c84f-430e-aa2e-2b14e5cbe112
# ╠═11dd1622-50c4-404e-88a2-c29c02760ca9
# ╟─e258ad9e-a39b-45c8-997b-d785977116fd
# ╠═55a66ce9-2c77-42dd-8aac-446cddb1274d
# ╠═fc926801-2060-4c9f-b103-466c30cb899a
# ╟─5a22ad68-2241-4c8e-bfc5-bbaf9bf786e2
# ╠═cccce8d6-afa5-484f-8871-9106e71da354
# ╠═84aedc8c-5a59-4e7a-8136-d209b16984bd
# ╠═87ad6166-a5f7-4fdd-9ca3-85ec1ac62fff
# ╠═217f7fb1-0cc4-4ead-8998-2fc5c13571a9
# ╠═97c7bd4b-f860-4e66-bd5a-a901f9c014ee
# ╠═8164635b-1ad6-4b74-8926-b90d6182cc69
# ╟─995fb655-6e01-452b-91fd-93b5cb89a17e
# ╠═67fc7cd4-afd1-4812-9e34-41527ca8e68a
# ╠═8944cab7-3a9e-4573-93cc-3f605fcb545a
# ╠═d6562dfe-6656-494b-aae7-9807aec745f9
# ╟─b79ef21d-e232-4c10-a0ba-6b4c54ad71aa

max_stars_repo_path: Preliminary WeightandBalance.jl

max_stars_repo_name: rahmanfahd/aircraft-design

max_issues_repo_path: Preliminary WeightandBalance.jl

max_issues_repo_name: rahmanfahd/aircraft-design

max_forks_repo_path: Preliminary WeightandBalance.jl

max_forks_repo_name: rahmanfahd/aircraft-design


text:
using ImageProjectiveGeometry
using Test

include("test_utilities.jl")
include("test_transforms.jl")
include("test_projective.jl")
include("test_ransac.jl")
include("test_cornerfeatures.jl")
include("test_geometry.jl")