text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset "test_check_TS" begin

    network = PSCOPF.Networks.Network()
    PSCOPF.add_new_bus!(network, "bus_1")
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "limitable_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100., #pmin, pmax : Not concerned ? min is always 0, max is the limitation
                                                0., 10., #start_cost, prop_cost : start cost is always 0 ?
                                                Dates.Second(0), Dates.Second(0)) #dmo, dp : always 0. ?
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "pilotable_1", PSCOPF.Networks.PILOTABLE,
                                                10., 200., #pmin, pmax
                                                45000., 30., #start_cost, prop_cost
                                                Dates.Second(3600), Dates.Second(0)) #dmo, dp

    TS = [DateTime("2015-01-01T14:00:00")]
    ECH = [DateTime("2015-01-01T07:00:00"), DateTime("2015-01-01T11:00:00")]

    @testset "definition_example" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                "S2" => 20.
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 0.,
                                                                "S2" => 16.
                                                                )
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                "S2" => 20.
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 10.,
                                                                "S2" => 16.
                                                                )
                                ) ,
                            ),
                        )

        @test PSCOPF.check_uncertainties(uncertainties, network)
    end

    @testset "all_entries_need_to_have_the_same_scenarios" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                "S2" => 20.
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 10.,
                                                                # "S2" => 16. #Missing Scenario
                                                                )
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                "S2" => 20.
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 10.,
                                                                "S2" => 16.
                                                                )
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_same_scenarios(uncertainties)
    end

    @testset "all_entries_need_to_have_the_same_timesteps" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                #different/missing timestep
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_same_timesteps(uncertainties)
    end

    # in the execution we allow having more TS than what we need !
    @testset "all_specified_ts_are_listed" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )
        @test PSCOPF.check_uncertainties_contain_ts(uncertainties, [DateTime("2015-01-01T14:00:00")])

        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                #DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                #DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                #DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ),
                            "limitable_1" =>SortedDict(
                                #DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.),
                                DateTime("2015-01-01T14:15:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )
        @test !PSCOPF.check_uncertainties_contain_ts(uncertainties,
                                                    [DateTime("2015-01-01T14:00:00"),DateTime("2015-01-01T14:15:00")])
    end

    @testset "all_values_must_be_>=0" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => -10., # negative value
                                                                )
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 15.,
                                                                )
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict(
                                                                "S1" => 10.,
                                                                )
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_values(uncertainties, network)
    end

    @testset "limitable_injection_uncertainties<=pmax" begin
        # Allow extra ECHs
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 115.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 110.) # exceeds pmax (100)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 115.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 90.)
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_values(uncertainties, network)
    end

    @testset "all_listed_ids_are_either_a_bus_or_a_limitable_of_the_network" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "bus_2" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "bus_2" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_values(uncertainties, network)
    end

    @testset "all_limitable_generators_are_listed_for_all_ech" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            # missing :
                            # "limitable_1" =>SortedDict(
                            #     DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                            #     ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )

        @test !PSCOPF.check_uncertainties_limitables(uncertainties, network)
    end

    @testset "all_buses_are_listed_for_all_ech" begin
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            # missing :
                            # "bus_1" => SortedDict(
                            #     DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                            #     ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )
        @test !PSCOPF.check_uncertainties_buses(uncertainties, network)
    end

    @testset "at_least_all_ech_are_described" begin
        # Allow extra ECHs
        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        # Extra ECH :
                        DateTime("2015-01-01T10:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        DateTime("2015-01-01T11:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        )

        @test PSCOPF.check_uncertainties_contains_ech(uncertainties, ECH)

        uncertainties = PSCOPF.Uncertainties(
                        DateTime("2015-01-01T07:00:00") => SortedDict(
                            "bus_1" => SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                                ),
                            "limitable_1" =>SortedDict(
                                DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                                ) ,
                            ),
                        # missing :
                        # DateTime("2015-01-01T11:00:00") => SortedDict(
                        #     "bus_1" => SortedDict(
                        #         DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 15.)
                        #         ),
                        #     "limitable_1" =>SortedDict(
                        #         DateTime("2015-01-01T14:00:00") => SortedDict("S1" => 10.)
                        #         ) ,
                        #     ),
                        )

        @test !PSCOPF.check_uncertainties_contains_ech(uncertainties, ECH)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using SpecialFunctions: gamma
using Interpolations: LinearInterpolation
using ProgressMeter

struct SedovData

     t::Float64
     m::Float64
     r::Vector{Float32}
     vr::Vector{Float32}
     rho::Vector{Float32}
     U::Vector{Float32}
     P::Vector{Float32}
     Pth::Vector{Float32}
     Pcr::Vector{Float32}
     E::Float64
     hsml::Vector{Float32}
     mach::Vector{Float32}
     r_shock_rho::Float64
     r_shock_P::Float64
     rho_out::Float64
     rho_s::Float64
     cs_out::Float64
     gamma::Float64
     ndim::Int64

     function SedovData(t::Float64,             m::Float64,
                        r::Vector{Float32},     vr::Vector{Float32},
                        rho::Vector{Float32},   U::Vector{Float32},
                        Pcr::Vector{Float32},   E::Float64,
                        hsml::Vector{Float32},  mach::Vector{Float32};
                        γ::Float64=5.0/3.0,     ndim::Int64=3)

        P = @. (γ - 1.0 ) * rho * U
        Pth = P

        if Pcr != zeros(Float32, length(P))
            P += Pcr
        end

        firstbin = Int(floor(length(rho) * 0.99)-1)

        rho_out = mean( rho[ firstbin:end ])
        P_out   = mean( P[ firstbin:end ] )

        cs_out  = sqrt( 5.0/3.0 * P_out/rho_out)

        rho_s = rho_out * ( γ + 1.0 )/( γ - 1.0 )

        k = findmax(rho)[2]
        r_shock_rho = r[k]

        k = findmax(P)[2]
        r_shock_P   = r[k]

        new(t, m, r, vr, rho, U, P, Pth, Pcr, E, hsml, mach, r_shock_rho, r_shock_P,
            rho_out, rho_s, cs_out, γ, ndim)
    end
end

struct SedovFit

    alpha::Float64
    w::Int64
    D
    P
    V

end

struct SedovSolution

    r::Vector{Float32}
    vr::Vector{Float64}
    rho::Vector{Float64}
    P::Vector{Float64}
    U::Vector{Float64}
    r_shock::Float64
    v_shock::Float64
    xs::Float64
    mach::Float64

end

function R_s_analytic(i::Int64, data::SedovData, fit::SedovFit)
    return (data.E * data.t^2 / ( fit.alpha * data.rho[i] ))^(1.0/ ( 2.0 + data.ndim ))
end

function R_s_dot_analytic(i::Int64, data::SedovData, fit::SedovFit)
    t = maximum([data.t, 1.e-30])
    return 2.0 * R_s_analytic(i, data, fit) / ( (data.ndim + 2.0 - fit.w) * t)
end

function v_s_analytic(i::Int64, data::SedovData, fit::SedovFit)
    return 2.0 * R_s_dot_analytic(i, data, fit) /
            ( data.gamma + 1.0 )
end

function P_s_analytic(i::Int64, data::SedovData, fit::SedovFit)
    return 2.0 * data.rho_out * R_s_dot_analytic(i, data, fit)^2 / ( data.gamma + 1.0 )
end

function rho_analytic(i::Int64, data::SedovData, fit::SedovFit)
    R = R_s_analytic(i, data, fit)
    xi = data.r[i]/R
    if xi < 1.0
        return data.rho_s * fit.D(xi)
    else
        return data.rho_out
    end
end

function P_analytic(i::Int64, data::SedovData, fit::SedovFit)
    R = R_s_analytic(i, data, fit)
    xi = data.r[i]/R
    if xi < 1.0
        return P_s_analytic(i, data, fit) * fit.P(xi)
    else
        return 0.0
    end
end

function v_analytic(i::Int64, data::SedovData, fit::SedovFit)
    R = R_s_analytic(i, data, fit)
    xi = data.r[i]/R
    if xi < 1.0
        return v_s_analytic(i, data, fit) * fit.V(xi)
    else
        return 0.0
    end
end

struct fit_parameters

    ndim::Int64
    γ::Float64

    w::Float64
    w1::Float64
    w2::Float64
    w3::Float64

    b0::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    b4::Float64
    b5::Float64
    b6::Float64
    b7::Float64
    b8::Float64

    C0::Float64
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64

end

struct w1Data
    ndim::Int64
    γ::Float64

    w::Float64
    w1::Float64
    w2::Float64
    w3::Float64

    b0::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    b4::Float64
    b5::Float64
    b6::Float64
    b7::Float64
    b8::Float64

    C0::Float64
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
end

struct w2Data
    ndim::Int64
    γ::Float64

    w::Float64
    w1::Float64
    w2::Float64
    w3::Float64

    b0::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    b4::Float64
    b5::Float64
    b6::Float64
    b7::Float64
    b8::Float64

    C0::Float64
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
end

struct w3Data
    ndim::Int64
    γ::Float64

    w::Float64
    w1::Float64
    w2::Float64
    w3::Float64

    b0::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    b4::Float64
    b5::Float64
    b6::Float64
    b7::Float64
    b8::Float64

    C0::Float64
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
end

struct elseData
    ndim::Int64
    γ::Float64

    w::Float64
    w1::Float64
    w2::Float64
    w3::Float64

    b0::Float64
    b1::Float64
    b2::Float64
    b3::Float64
    b4::Float64
    b5::Float64
    b6::Float64
    b7::Float64
    b8::Float64

    C0::Float64
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
end

"""
        w1 functions
"""
function xifunc(F::Float64, c::w1Data)
    return F^(-c.b6) * (c.C1 * (c.F - c.C2))^c.b2 * (c.C3 * (c.C4 - F))^(-c.b1)
end
function Dfunc(F::Float64, c::w1Data)
    return xifunc(F, c)^(c.ndim-2)
end
function Vfunc(F::Float64, c::w1Data)
    return xifunc(F, c)
end
function Pfunc(F::Float64, c::w1Data)
    return xifunc(F, c)^c.ndim
end

"""
        w2 functions
"""
function xifunc(F::Float64, c::w2Data)
    return  F^(-c.b6) * (c.C1 * (F - c.C2))^((c.γ - 1.0 )*c.b0) *
            exp(((c.γ + 1.0)*c.b0*(1.0 - F)) / (F - c.C2))
end
function Dfunc(F::Float64, c::w2Data)
    return  F^c.b7 * (c.C1 * (F - c.C2))^((4.0 - c.ndim - 2.0*c.γ)*c.b0) *
            (c.C5 * (c.C6 - F))^(-c.b5) *
            exp((-2.0*(c.γ + 1.0)*c.b0*(1.0 - F)) / (F - c.C2))
end
function Vfunc(F::Float64, c::w2Data)
    return xifunc(F, c) * F
end
function Pfunc(F::Float64, c::w2Data)
    return F^c.b8 * (c.C1 * (F - c.C2))^(-c.ndim * c.γ * c.b0) *
           (c.C5 * (c.C6 - F))^(1.0 - c.b5)
end

"""
        w3 functions
"""
function xifunc(F::Float64, c::w3Data)
    return F^(-c.b6) * (c.C1 * (F - c.C2))^c.b2 * (c.C5 * (c.C6 - F))^(-c.b1)
end
function Dfunc(F::Float64, c::w3Data)
    return F^c.b7 * (c.C1 * (F - c.C2))^(c.b3-c.w*c.b2) *
           (c.C5 * (c.C6 - F))^(1.0 - 4.0*c.b0) *
           exp((-c.ndim*c.γ*(c.γ + 1.0)*c.b0*(1.0 - F)) /
           (c.C6 - F))
end
function Vfunc(F::Float64, c::w3Data)
    return xifunc(F, c) * F
end
function Pfunc(F::Float64, c::w3Data)
    return F^c.b8 * (c.C5 * (c.C6 - F))^(2.0*(c.ndim*c.γ-c.ndim-c.γ)*c.b0) *
           exp((-c.ndim*c.γ*(c.γ+1.0)*c.b0*(1.0 - F)) / (c.C6 - F))
end

"""
    else functions
"""
function xifunc(F::Float64, c::elseData)
    return F^(-c.b6) * (c.C1 * (F - c.C2))^c.b2 * (c.C3 * (c.C4 - F))^(-c.b1)
end
function Dfunc(F::Float64, c::elseData)
    return F^c.b7 * (c.C1 * (F - c.C2))^(c.b3 - c.w*c.b2) *
           (c.C3 * (c.C4 - F))^(c.b4 + c.w*c.b1) *
           (c.C5 * (c.C6 - F))^(-c.b5)
end
function Vfunc(F::Float64, c::elseData)
    return xifunc(F, c) * F
end
function Pfunc(F::Float64, c::elseData)
    return F^(c.b8) * (c.C3 * (c.C4 - F))^(c.b4 + (c.w - 2.0)*c.b1) *
           (c.C5 * (c.C6 - F))^(1.0 - c.b5)
end

function get_ideal_sedov_fit(rho_s::Float64; ndim::Int64=3, γ::Float64=5.0/3.0)

    w = 0 # w parametrizes the background density as rho = A * r^(-w)
    # The following is copied from Book (1991)
    # The omega parameters
    w1 = (3*ndim - 2 + γ*(2-ndim)) / (γ + 1)
    w2 = (2*(γ-1) + ndim) / γ
    w3 = ndim * (2-γ)
    # The beta_n parameters:
    b0 = 1 / (ndim * γ - ndim + 2)
    b2 = (γ - 1) / (γ*(w2-w))
    b3 = (ndim - w) / (γ*(w2-w))
    b5 = (2*ndim - w*(γ+1)) / (w3 - w)
    b6 = 2 / (ndim + 2 -w)
    b1 = b2 + (γ+1)*b0 -b6
    b4 = b1 * ((ndim - w) * (ndim + 2 - w)) / (w3 - w)
    b7 = w * b6
    b8 = ndim * b6
    # The C parameters
    C0 = 2^ndim * π^((ndim-1)/2) * gamma((ndim+1)/2) / gamma(ndim)
    C5 = 2 / (γ-1)
    C6 = (γ+1) / 2
    C1 = γ * C5
    C2 = C6 / γ
    C3 = (ndim*γ - ndim + 2) / ((w1-w)*C6)
    C4 = (ndim + 2 - w) * b0 * C6
    # Eq. (8) to (11) in Book (1991) or their later alternatives

    # println("w = $w\nw1 = $w1\nw2 = $w2")
    if w == w1
        # println("w1!")
        c = w1Data(ndim, γ,
                   w, w1, w2, w3,
                   b0, b1, b2, b3, b4, b5, b6, b7, b8,
                   C0, C1, C2, C3, C4, C5, C6)
    elseif w == w2
        # println("w2!")
        c = w2Data(ndim, γ,
                   w, w1, w2, w3,
                   b0, b1, b2, b3, b4, b5, b6, b7, b8,
                   C0, C1, C2, C3, C4, C5, C6)
    elseif w == w3
        # println("w3!")
        c = w3Data(ndim, γ,
                   w, w1, w2, w3,
                   b0, b1, b2, b3, b4, b5, b6, b7, b8,
                   C0, C1, C2, C3, C4, C5, C6)
    else
        # println("none!")
        c = elseData(ndim, γ,
                   w, w1, w2, w3,
                   b0, b1, b2, b3, b4, b5, b6, b7, b8,
                   C0, C1, C2, C3, C4, C5, C6)
    end

    if w1 > w
        Fmin = C2
    else
        Fmin = C6
    end

    # We use F as a proxy to do the integration over xi
    F = range(Fmin, stop=1.0, length=Int(1e5))
    xi = zeros(length(F))
    for i = 1:length(F)
        xi[i] = xifunc(F[i], c)
    end

    k = sortperm(F)
    F = F[k] # sort F according to the xi-values it procudes
    xi = xi[k]

    D = zeros(length(F))
    V = zeros(length(F))
    P = zeros(length(F))
    for i = 1:length(F)
        D[i] = Dfunc(F[i], c)
        V[i] = Vfunc(F[i], c)
        P[i] = Pfunc(F[i], c)
    end
    D_interpolated = LinearInterpolation(xi, D)
    P_interpolated = LinearInterpolation(xi, P)
    V_interpolated = LinearInterpolation(xi, V)
    I = @. xi^(ndim-1) * (D*V^2 + P) # Integrand
    dxi = zeros(length(I))
    # Use middle points
    for i = 1:(length(I)-1)
        I[i] = 0.5 * (I[i] + I[i+1])
        dxi[i] = xi[i+1] - xi[i]
    end
    integral = sum( @. I * dxi )
    # Finally, get alpha
    alpha = ((8 * C0) / ((γ^2 - 1)*(ndim + 2 - w)^2)) * integral

    return SedovFit(alpha, w, D_interpolated, P_interpolated, V_interpolated)
end

function get_ideal_sedov_data(data::SedovData, fit::SedovFit)

    N   = length(data.r)
    rho = zeros(N)
    P   = zeros(N)
    v   = zeros(N)

    for i = 1:N
        rho[i] = rho_analytic(i, data, fit)
        P[i]   = P_analytic(i, data, fit)
        v[i]   = v_analytic(i, data, fit)
    end

    U = @. P / ( (data.gamma - 1.0 ) * rho )

    k = findmax(rho)[2]

    r_shock = data.r[k]
    xs = rho[k]/rho[end]
    v_shock = v[k]/(1.0 - 1.0/xs)
#    v_shock = v[k] + data.cs_out * sqrt( 1 + (data.gamma + 1.0)/data.gamma * (P[k]/data.P[end] - 1.0))
    mach = v_shock/data.cs_out

    return SedovSolution(data.r, v, rho, P, U, r_shock, v_shock, xs, mach)
end

function get_sedov_data_from_snapshot(fi::String, blast_center::Vector{Float64}=[3.0, 3.0, 3.0];
                                      CRs=false, Nbins::Int64=0)

    h = head_to_obj(fi)

    t = h.time
    m = h.massarr[1]

    x = read_block_by_name(fi, "POS", parttype=0)

    r = @. sqrt( (x[:,1] - blast_center[1])^2 +
                 (x[:,2] - blast_center[2])^2 +
                 (x[:,3] - blast_center[3])^2 )

    k = sortperm(r)

    r = Float32.(r[k])

    v = read_block_by_name(fi, "VEL", parttype=0)[k,:]

    vr = @. sqrt( v[:,1]^2 + v[:,2]^2 + v[:,3]^2 )

    rho = read_block_by_name(fi, "RHO", parttype=0)[k, 1]

    U = read_block_by_name(fi, "U", parttype=0)[k, 1]

    hsml = read_block_by_name(fi, "HSML", parttype=0)[k, 1]

    try
        mach = read_block_by_name(fi, "MACH", parttype=0)[k, 1]
    catch
        mach = zeros(Float32, length(U))
    end

    if CRs
        CRpP = read_block_by_name(fi, "CRpP", parttype=0)[k, 1]
        Ecr = @. CRpP/(1.0/3.0 * rho)
        γ = 7.0/5.0
    else
        CRpP = zeros(Float32, length(U))
        Ecr = zeros(Float32, length(U))
        γ = 5.0/3.0
    end

    E = sum( @. (Ecr * m + U * m + 0.5 * m * vr^2) )

    if Nbins != 0
        N = length(r)

        rbinwidth = maximum(r)/Nbins
        r_raw    = r
        v_raw    = vr
        rho_raw  = rho
        U_raw    = U
        CRpP_raw = CRpP
        hsml_raw  = hsml
        mach_raw  = mach

        r     = zeros(Float32, Nbins+1)
        vr    = zeros(Float32, Nbins+1)
        rho   = zeros(Float32, Nbins+1)
        U     = zeros(Float32, Nbins+1)
        CRpP  = zeros(Float32, Nbins+1)
        hsml  = zeros(Float32, Nbins+1)
        mach  = zeros(Float32, Nbins+1)
        Npart = zeros(Float32, Nbins+1)

        @showprogress "Binning..." for i = 1:N
            bin = floor(Int64, r_raw[i]/rbinwidth) + 1

            vr[bin]     += v_raw[i]
            rho[bin]    += rho_raw[i]
            U[bin]      += U_raw[i]
            CRpP[bin]   += CRpP_raw[i]
            hsml[bin]  += hsml_raw[i]
            mach[bin]  += mach_raw[i]
            Npart[bin]  += 1
        end

        for i = 1:(Nbins+1)
            r[i] = (i-1) * rbinwidth
            vr[i]   /= Npart[i]
            rho[i]  /= Npart[i]
            U[i]    /= Npart[i]
            CRpP[i] /= Npart[i]
            mach[i] /= Npart[i]
            hsml[i] /= Npart[i]
        end

    end

    k = findall(isnan.(rho))

    deleteat!(r, k)
    deleteat!(vr, k)
    deleteat!(rho, k)
    deleteat!(U, k)
    deleteat!(CRpP, k)
    deleteat!(hsml, k)
    deleteat!(mach, k)

    return SedovData(t, m, r, vr, rho, U, CRpP, E, hsml, mach, γ=γ)

end

function get_sedov_solution(filename::String, blast_center::Vector{Float64}=[3.0, 3.0, 3.0];
                            CRs::Bool=false, Nbins::Int64=0, Ndim::Int64=3)

    sedov_data = get_sedov_data_from_snapshot(filename, blast_center, CRs=CRs, Nbins=Nbins)

    if CRs
        γ=7.0/5.0
    else
        γ=5.0/3.0
    end

    sedov_fit  = get_ideal_sedov_fit(sedov_data.rho_s, ndim=Ndim, γ=γ)

    sedov_ideal = get_ideal_sedov_data(sedov_data, sedov_fit)

    return sedov_data, sedov_ideal

end

ext: jl

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Umatrix

makedocs(
    sitename = "Umatrix.jl",
    authors = "Tobias Frilling",
    format = Documenter.HTML(),
    modules = [Umatrix],
    pages = [
             "Home" => "index.md",
             "ESOM Training" => "esom.md",
             "Matrices" => "matrices.md",
             "Plotting" => "plotting.md",
             "Settings" => "settings.md",
             "API" => "api.md"
            ]
)

deploydocs(
    repo = "github.com/ckafi/Umatrix.jl"
)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function JuliaFEM.get_model(::Type{Val{Symbol("1x1 plane stress quad4 block")}})

    X = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])

    body = Problem(Elasticity, "body", 2)
    body.properties.formulation = :plane_stress
    body.elements = [Element(Quad4, [1, 2, 3, 4])]
    update!(body.elements, "geometry", X)
    update!(body.elements, "youngs modulus", 288.0)
    update!(body.elements, "poissons ratio", 1/3)

    # boundary conditions
    bc_13 = Problem(Dirichlet, "symmetry 13", 2, "displacement")
    bc_13.properties.dual_basis = true
    bc_13.elements = [Element(Seg2, [1, 2])]
    update!(bc_13.elements, "geometry", X)
    update!(bc_13.elements, "displacement 2", 0.0)

    bc_23 = Problem(Dirichlet, "symmetry 23", 2, "displacement")
    bc_23.properties.dual_basis = true
    bc_23.elements = [Element(Seg2, [4, 1])]
    update!(bc_23.elements, "geometry", X)
    update!(bc_23, "displacement 1", 0.0)

    solver = Solver(Nonlinear, "1x1 plane stress quad4 block")
    push!(solver, body, bc_13, bc_23)

    return solver
end

@testset "test dirichlet spc in point" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    solver = get_model("1x1 plane stress quad4 block")
    update!(solver["symmetry 13"], "displacement 1", 0.0)
    update!(solver["symmetry 23"], "displacement 2", 0.0)
    nodal_bc = Problem(Dirichlet, "dx=0.5", 2, "displacement")
    nodal_bc.elements = [Element(Poi1, [3])]
    update!(nodal_bc, "geometry", X)
    update!(nodal_bc, "displacement 1", 0.5)
    update!(nodal_bc, "displacement 2", 0.0)
    push!(solver, nodal_bc)
    solver()
    pel = nodal_bc.elements[1]
    la = pel("reaction force", [0.0], 0.0)
    info("reaction force: $la")
    info(solver["body"].assembly.u)
    @test isapprox(pel("displacement", [], 0.0), [0.5, 0.0])
end

@testset "test nodal point force" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    solver = get_model("1x1 plane stress quad4 block")
    update!(solver["symmetry 13"], "displacement 1", 0.0)
    update!(solver["symmetry 23"], "displacement 2", 0.0)
    point_load = Element(Poi1, [3])
    update!(point_load, "geometry", X)
    update!(point_load, "displacement traction force 1", 72.0)
    update!(point_load, "displacement traction force 2", 27.0)
    push!(solver["body"], point_load)
    solver()
    @test isapprox(point_load("displacement", [], 0.0), [0.5, 0.0])
end


text:
module Simulations

using CombinatorialSpaces
using Catlab
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Theories
using PreallocationTools
using LinearAlgebra

export gen_sim,
       BoxFunc, MatrixFunc, ElementwiseFunc, ArbitraryFunc, InPlaceFunc,
       ConstantFunc, TDInPlaceFunc

abstract type BoxFunc end
struct MatrixFunc <: BoxFunc end
struct ElementwiseFunc <: BoxFunc end
struct ArbitraryFunc <: BoxFunc end
struct ConstantFunc <: BoxFunc end
struct InPlaceFunc <: BoxFunc end
struct TDInPlaceFunc <: BoxFunc end

form2dim = Dict(:Scalar => x->1,
                :Form0 => nv,
                :Form1 => ne,
                :Form2 => ntriangles,
                :DualForm2 => nv,
                :DualForm1 => ne,
                :DualForm0 => ntriangles)
dims(x) = begin
  k = filter(i -> x[:type] isa ObExpr{i}, keys(form2dim))
  length(k) == 1 || error("Object $x has multiple keys in `form2dim`")
  form2dim[first(k)]
end

function gen_sim(dwd::WiringDiagram, name2func, s; form2dim=form2dim, params=[], autodiff=false)
  check_consistency(dwd)
  d = dwd.diagram

  # Generate cached memory. These variables are indexed by their respective
  # output ports
  mem = [zeros(Float64, dims(f)(s)) for f in d[:out_port_type]]
  if autodiff
    mem = [dualcache(zeros(Float64, dims(f)(s))) for f in d[:out_port_type]]
  end
  tgt2src = Dict{Int64, Int64}()
  for w in 1:nparts(d, :Wire)
    d[w, :tgt] ∈ keys(tgt2src) && error("Two wires input to port $(d[w, :src])")
    tgt2src[d[w, :tgt]] = d[w, :src]
  end

  # Get the indices in input/output arrays for specific input/output forms
  input_inds = Vector{Tuple{Int64, Int64, Bool}}()
  output_inds = Vector{Tuple{Int64, Int64}}()

  cur_ind = 0
  cur_param = 0
  for ft in d[:outer_in_port_type]
    mag = dims(ft)(s)
    if ft[:name] ∈ params
      push!(input_inds, (cur_param+1, mag + cur_param, true))
      cur_param += mag
    else
      push!(input_inds, (cur_ind+1, mag + cur_ind, false))
      cur_ind += mag
    end
  end
  cur_ind = 0
  for ft in d[:outer_out_port_type]
    mag = dims(ft)(s)
    push!(output_inds, (cur_ind+1, mag + cur_ind))
    cur_ind += mag
  end

  # Make local variables for any MatrixFunc/ArbitraryFunc functions to be
  # included in the function closure
  n2f = deepcopy(name2func)
  matrices = Vector{Any}()
  funcs = Vector{Function}()

  # First develop a list of needed functions/matrices
  used_funcs = unique(d[:value])

  for k in used_funcs
    ops = split("$k", "⋅")
    if length(ops) > 1
      mat_val = foldl(*, map(o->name2func[Symbol(o)][:operator], ops))
      push!(matrices, mat_val)
      n2f[k] = Dict(:operator => :(matrices[$(length(matrices))]), :type => MatrixFunc())
    else
      v = name2func[k]
      if(v[:type] isa MatrixFunc)
        push!(matrices, v[:operator])
        n2f[k] = Dict(:operator => :(matrices[$(length(matrices))]), :type => v[:type])
      elseif(v[:type] isa Union{ElementwiseFunc, ArbitraryFunc, InPlaceFunc, TDInPlaceFunc})
        push!(funcs, v[:operator])
        n2f[k] = Dict(:operator => :(funcs[$(length(funcs))]),:type => v[:type])
      end
    end
  end

  # Compute composed matrices
  for v in d[:value]

  end

  # Topological sort of DWD for scheduling
  execution_order = topological_sort(internal_graph(dwd))

  # Generate function expression
  body = quote
  end

  # Assign variables for each dual memory location
  # Allows for get_tmp to be called for each cached memory location
  if autodiff
    dual_init = map(1:length(mem)) do i
      :($(Symbol("d_$i")) = get_tmp(mem[$i], u))
    end
    append!(body.args, dual_init)
  end

  exec_dwd = map(execution_order) do b
    iw = in_wires(dwd, b)

    input_args = map(iw) do w
      p = w.source
      if p.box == -2
        if input_inds[p.port][3]
          :(view(p, $(input_inds[p.port][1]):$(input_inds[p.port][2])))
        else
          :(view(u, $(input_inds[p.port][1]):$(input_inds[p.port][2])))
        end
      else
        if autodiff
          :($(Symbol("d_$(incident(d, p.box, :out_port_box)[p.port])")))
        else
          :(mem[$(incident(d, p.box, :out_port_box)[p.port])])
        end
      end
    end
    input_args[[w.target.port for w in iw]] .= input_args

    output_args = map(incident(d, b, :out_port_box)) do p
      if autodiff
        :($(Symbol("d_$p")))
      else
        :(mem[$p])
      end
    end

    gen_func(input_args, output_args, n2f[d[b, :value]])
  end
  append!(body.args, exec_dwd)

  # set du values
  du_set = map(1:nparts(d, :OuterOutPort)) do i
    iw = incident(d, i, :out_tgt)
    length(iw) == 1 || error(Output)
    w = first(iw)
    pt = d[w, :out_tgt]
    ps = d[w, :out_src]
    if autodiff
      :(setindex!(du, Symbol("d_$ps"), $(collect(output_inds[pt][1]:output_inds[pt][2]))))
    else
      :(setindex!(du, mem[$ps], $(collect(output_inds[pt][1]:output_inds[pt][2]))))
    end
  end
  append!(body.args, du_set)

  # Currently this generates a function at the module (can be accessed from the module)
  # Is this alright, or is this going to just lead to memory problems?
  fname = gensym()
  res = :($fname(du, u, p, t, mem, matrices, funcs) = $body)
  eval(res)
  f_local = deepcopy(eval(fname))
  function sim_func(du, u, p, t)
    f_local(du, u, p, t, mem, matrices, funcs)
  end, res
end

function gen_func(input_args, output_args, func_info)
  _gen_func(func_info[:type], input_args, output_args, func_info[:operator])
end

function _gen_func(::MatrixFunc, input_args, output_args, func)
  length(input_args) == 1 || error("Length of input for Matrix function is $(length(input_args)) ≠ 1")
  length(output_args) == 1 || error("Length of output for Matrix function is $(length(output_args)) ≠ 1")
  :(mul!($(output_args[1]), $(func), $(input_args[1])))
end

function _gen_func(::ArbitraryFunc, input_args, output_args, func)

  length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))")
  :($(Meta.parse(join(output_args, ","))) = $(Expr(:call, func, input_args...)))
end

function _gen_func(::InPlaceFunc, input_args, output_args, func)
  length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))")
  :($(Expr(:call, func, vcat(output_args, input_args)...)))
end

function _gen_func(::TDInPlaceFunc, input_args, output_args, func)
  _gen_func(InPlaceFunc(), vcat(input_args, [:t]), output_args, func)
end

function _gen_func(::ConstantFunc, input_args, output_args, func)
  length(output_args) == 0 && error("Length of output for Constant function is $(length(output_args))")
  :($(Meta.parse(join(output_args, ","))) .= $(func))
end

function _gen_func(::ElementwiseFunc, input_args, output_args, func)

  length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))")
  body = :($(Meta.parse(join([Expr(:ref, o, :i) for o in output_args], ","))) = $(Expr(:call, func, [Expr(:ref, n, :i) for n in input_args]...)))
  quote
    for i in eachindex($(input_args[1]))
      $body
    end
  end
end

function check_consistency(dwd::WiringDiagram)
  d = dwd.diagram
  for w in 1:nparts(d, :Wire)
    sp = d[w, :src]
    tp = d[w, :tgt]
    stype = d[sp, :out_port_type]
    ttype = d[tp, :in_port_type]
    sb = d[sp, :out_port_box]
    tb = d[tp, :in_port_box]

    stype == ttype ||
      error("""Wire $w maps a port of type $stype to a port of type $ttype
            This wire is between box $(sb)($(d[sb, :value])) and box $(tb)($(d[tb, :value]))""")
  end
end

end


text:
function prod_dim_gen_impl{T, N}(x::Array{T, N})
ex = :(1)
for i = 1:N
ex = :(size(x, $i) * $ex)
end
return ex
end

max_stars_repo_path: Module 2/Chapter04/@prod_dim_gen.jl

max_issues_repo_path: Module 2/Chapter04/@prod_dim_gen.jl

max_forks_repo_path: Module 2/Chapter04/@prod_dim_gen.jl


text:
"""
Conjoin two SDDs
"""
@inline conjoin(::SddFalseNode, ::SddTrueNode) = false_sdd
@inline conjoin(::SddTrueNode, ::SddFalseNode) = false_sdd
@inline conjoin(s::Sdd, ::SddTrueNode) = s
@inline conjoin(::Sdd, ::SddFalseNode) = false_sdd
@inline conjoin(::SddTrueNode, s::Sdd) = s
@inline conjoin(::SddFalseNode, ::Sdd) = false_sdd
@inline conjoin(::SddTrueNode, ::SddTrueNode) = true_sdd
@inline conjoin(::SddFalseNode, ::SddFalseNode) = false_sdd

# const stats = Dict{Tuple{Int,Int},Int}()

function conjoin(s::SddLiteralNode, t::SddLiteralNode)::Sdd
    if mgr(s) === mgr(t)
        (s === t) ? s : false_sdd
    else
        conjoin_indep(s,t)
    end
end

# Note: attempts to make a special cache for conjunctions with literals have not yielded speedups

function conjoin(s::Sdd, t::Sdd)::Sdd
    if mgr(s) === mgr(t)
        conjoin_cartesian(t,s)
    elseif varsubset(s,t)
        conjoin_descendent(s,t)
    elseif varsubset(t,s)
        conjoin_descendent(t,s)
    else
        conjoin_indep(s,t)
    end
end

"""
Conjoin two SDDs when they respect the same vtree node
"""
function conjoin_cartesian(n1::Sdd⋁Node, n2::Sdd⋁Node)::Sdd
    if n1 === n2
        return n1
    elseif n1 === !n2
        return false_sdd
    end
    get!(mgr(n1).conjoin_cache, ApplyArgs(n1,n2)) do
        conjoin_cartesian_general(n1,n2)
    end::Sdd
end

function conjoin_cartesian_general(n1::Sdd⋁Node, n2::Sdd⋁Node)
    # vast majority of cases are 2x2 and 2x3 applies, yet specializing for those cases does not appear to speed things up

    out = XYPartition()
    sizehint!(out, num_children(n1) * num_children(n2))

    elems1 = children(n1)
    elems2 = children(n2)

    if length(elems1) > length(elems2)
        elems1, elems2 = elems2, elems1
    end

    maski = falses(length(elems1))
    maskj = falses(length(elems2))

    conjoin_cartesian_cheap(out, elems1, elems2, maski, maskj)
    conjoin_cartesian_expensive(out, elems1, elems2, maski, maskj)

    canonicalize(out, mgr(n1))::Sdd
end

function conjoin_cartesian_cheap(out, elems1, elems2, maski, maskj)
    # first do a pass checking pointer equality, no actual applies unless strictly needed
    for i in 1:length(elems1)
        if !maski[i]
            e1 = elems1[i]
            for j in 1:length(elems2)
                if !maski[i] && !maskj[j]
                    e2 = elems2[j]
                    if prime(e1) === prime(e2)
                        push!(out, Element(prime(e1), conjoin(sub(e1),sub(e2))))
                        maski[i] = true
                        maskj[j] = true
                    elseif prime(e1) === !prime(e2)
                        # p1 |= prime(e2) for all other primes p1 in e1, save the unnecessary conjunctions
                        for k in 1:length(elems1)
                            if i!=k && !maski[k]
                                e3 = elems1[k]
                                push!(out, Element(prime(e3), conjoin(sub(e2),sub(e3))))
                            end
                        end
                        # p2 |= prime(e1) for all other primes p2 in e2, save the unnecessary conjunctions
                        for k in 1:length(elems2)
                            if j!=k && !maskj[k]
                                e3 = elems2[k]
                                push!(out, Element(prime(e3), conjoin(sub(e1),sub(e3))))
                            end
                        end
                        maski[i] = true
                        maskj[j] = true
                        break # nothing more to do for j
                    end
                end
            end
        end
    end
end

function conjoin_cartesian_expensive(out, elems1, elems2, maski, maskj)
    # now do the actual applies as needed
    for i in 1:length(elems1)
        if !maski[i]
            e1 = elems1[i]
            for j in 1:length(elems2)
                if !maski[i] && !maskj[j]
                    e2 = elems2[j]
                    newprime = conjoin(prime(e1),prime(e2))
                    if newprime !== false_sdd
                        newsub = conjoin(sub(e1),sub(e2))
                        push!(out, Element(newprime, newsub))
                    end
                    if newprime === prime(e2)
                        # p2 |= p1 and therefore p2 will be mutex with all other p1-primes
                        maskj[j] = true
                    end
                    if newprime === prime(e1)
                        # p1 |= p2 and therefore p1 will be mutex with all other p2-primes
                        maski[i] = true
                        break # all other j will have UNSAT primes
                    end
                end
            end
        end
    end
end

"""
Conjoin two SDDs when one descends from the other
"""
function conjoin_descendent(d::Sdd, n::Sdd)
    get!(mgr(n).conjoin_cache, ApplyArgs(d,n)) do
        elems = children(n)
        if varsubset_left(d, n)
            out = XYPartition()
            sizehint!(out, length(elems)+1)
            i = findfirst(c -> prime(c) === d, elems)
            if issomething(i)
                # there is a prime equal to d, all other primes will conjoin to false
                if sub(elems[i]) === false_sdd
                    return false_sdd
                elseif sub(elems[i]) === true_sdd
                    return d
                else
                    push!(out, Element(d, sub(elems[i])))
                    push!(out, Element(!d, false_sdd))
                    # since d is not a constant, must be trimmed and compressed
                    return unique⋁(out, mgr(n))
                end
            end
            i = findfirst(c -> prime(c) === !d, elems)
            if issomething(i)
                # there is a prime equal to !d, all other primes will conjoin to themselves
                for j in eachindex(elems)
                    j!=i && push!(out,Element(prime(elems[j]), sub(elems[j])))
                end
            else
                for e in elems
                    newprime = conjoin(prime(e),d)
                    if (newprime !== false_sdd)
                        push!(out, Element(newprime, sub(e)))
                    elseif newprime === d
                        # all future conjunctions will yield false
                        break
                    end
                end
            end
            push!(out, Element(!d, false_sdd))
        else
            # @assert varsubset_right(d, n)
            # TODO: build vector and compress all at once...
            out = [Element(prime(e),conjoin(sub(e),d)) for e in elems]
        end
        canonicalize(out, mgr(n))
    end
end

"""
Conjoin two SDDs in separate parts of the vtree
"""
function conjoin_indep(s::Sdd, t::Sdd)
    # @assert GateType(s)!=ConstantGate() && GateType(t)!=ConstantGate()
    lca_mgr = lca(mgr(s),mgr(t))
    get!(lca_mgr.conjoin_cache, ApplyArgs(s,t)) do
        if varsubset_left(mgr(s), lca_mgr)
            @assert varsubset_right(mgr(t), lca_mgr)
            elements = Element[Element(s,t),Element(!s,false_sdd)]
        else
            @assert varsubset_left(mgr(t), lca_mgr)
            @assert varsubset_right(mgr(s), lca_mgr)
            elements = Element[Element(t,s),Element(!t,false_sdd)]
        end
        # TODO: the XY partition must already be compressed and trimmed
        unique⋁(elements, lca_mgr)
    end::Sdd⋁Node
end

"""
Disjoin two SDDs
"""

disjoin(s::Sdd, t::Sdd) = !conjoin(!s,!t)


text:
##
#   File:      sdo1.jl
#
#   Purpose:   Demonstrates how to solve a small mixed semidefinite and conic quadratic
#              optimization problem using the MOSEK Julia API.
##

using Mosek

printstream(msg::String) = print(msg)

# Bound keys for constraints
bkc = [MSK_BK_FX,
       MSK_BK_FX]

# Bound values for constraints
blc = [1.0, 0.5]
buc = [1.0, 0.5]

A = sparse( [1,2,2],[1,2,3],[1.0, 1.0, 1.0])
conesub = [1, 2, 3]

barci = [1, 2, 2, 3, 3]
barcj = [1, 1, 2, 2, 3]
barcval = [2.0, 1.0, 2.0, 1.0, 2.0]

barai   = Any[ [1, 2, 3],
               [1, 2, 3, 2, 3, 3] ]
baraj   = Any[ [1, 2, 3],
               [1, 1, 1, 2, 2, 3] ]
baraval = Any[ [1.0, 1.0, 1.0],
               [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] ]

numvar = 3
numcon = length(bkc)
barvardim = [3]

# Create a task object and attach log stream printer
maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,printstream)

    # Append 'numvar' variables.
    # The variables will initially be fixed at zero (x=0).
    appendvars(task,numvar)

    # Append 'numcon' empty constraints.
    # The constraints will initially have no bounds.
    appendcons(task,numcon)

    # Append matrix variables of sizes in 'BARVARDIM'.
    # The variables will initially be fixed at zero.
    appendbarvars(task,barvardim)

    # Set the linear term c_0 in the objective.
    putcj(task, 1, 1.0)

    # Set the bounds on variable j
    # blx[j] <= x_j <= bux[j]
    putvarboundslice(task,1,numvar+1,
                     [ MSK_BK_FR::Int32 for i in 1:numvar ],
                     [ -Inf             for i in 1:numvar ],
                     [ +Inf             for i in 1:numvar ])

    # Set the bounds on constraints.
    # blc[i] <= constraint_i <= buc[i]
    putconboundslice(task,1,numcon+1, bkc,blc,buc)

    # Input row i of A
    putacolslice(task,1,numvar+1,
                 A.colptr[1:numvar], A.colptr[2:numvar+1],
                 A.rowval,A.nzval)

    appendcone(task,MSK_CT_QUAD, 0.0, conesub)

    symc  = appendsparsesymmat(task,barvardim[1],
                               barci,
                               barcj,
                               barcval)

    syma0 = appendsparsesymmat(task,barvardim[1],
                               barai[1],
                               baraj[1],
                               baraval[1])

    syma1 = appendsparsesymmat(task,barvardim[1],
                               barai[2],
                               baraj[2],
                               baraval[2])

    putbarcj(task,1, [symc], [1.0])

    putbaraij(task,1, 1, [syma0], [1.0])
    putbaraij(task,2, 1, [syma1], [1.0])

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Solve the problem and print summary
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    # Get status information about the solution
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta == MSK_SOL_STA_OPTIMAL || solsta == MSK_SOL_STA_NEAR_OPTIMAL
        # Output a solution
        xx = getxx(task,MSK_SOL_ITR)
        barx = getbarxj(task,MSK_SOL_ITR, 1)

        @printf("Optimal solution: \n  xx = %s\n  barx = %s\n", xx',barx')
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif  solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif  solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
# This must be run from  CGP.jl.
# To run:  julia -p 4 -L src/CGP.jl runtests.jl # choose nprocs according to your computer
# As of 2/8/21 all files run, but many have no tests.
# Note that for test_random_walk.jl one test written as an assertion fails (now commented out)

println("running test/test_Entropy.jl")
@everywhere include("test/test_Entropy.jl")

println("running test/test_Entropy.jl")
@everywhere include("test/test_mutate.jl")

println("running test/test_degen_complexity.jl")
@everywhere include("test/test_degen_complexity.jl")

println("running test/test_inf_alleles.jl")
@everywhere include("test/test_inf_alleles.jl")

println("running test/test_mut_evolve.jl")
@everywhere include("test/test_mut_evolve.jl")

println("running test/test_pop_evolve.jl")
@everywhere include("test/test_pop_evolve.jl")

println("running test/test_geno_robust_evolvability.jl")
@everywhere include("test/test_geno_robust_evolvability.jl")  # no test included

println("running test/test_generate_random_functions.jl")
@everywhere include("test/test_generate_random_functions.jl")

println("running test/test_run_geno_pheno.jl")
@everywhere include("test/test_run_geno_pheno.jl")

println("running test/test_build_chromosome.jl")
@everywhere include("test/test_build_chromosome.jl")  # no test included

println("running test/test_mutrobust_evolvability.jl")
@everywhere include("test/test_mutrobust_evolvability.jl")    # runs functions with no tests

println("running test/test_next_chromosome.jl")
@everywhere include("test/test_next_chromosome.jl")    # runs functions with no tests

println("running test/test_random_walk.jl")
@everywhere include("test/test_random_walk.jl")  # runs but assertion (commented out) fails

println("running test/test_selection.jl")
@everywhere include("test/test_selection.jl")  # no test included


text:
export
        DEFAULT_GTK_CANVAS_WIDTH,
        DEFAULT_GTK_CANVAS_HEIGHT,
        DEFAULT_GTK_WINDOW_IS_RESIZABLE,
        GTK_RENDERMODEL,
        GTK_CANVAS,
        GTK_WINDOW,
        GTK_SCENE,
        GTK_ROADWAY,

        show_gtk_window,
        hide_gtk_window,
        set_gtk_scene,
        set_gtk_roadway,
        render_gtk_video_1d

const DEFAULT_GTK_CANVAS_WIDTH = 800 # [pix]
const DEFAULT_GTK_CANVAS_HEIGHT = 600 # [pix]
const DEFAULT_GTK_WINDOW_IS_RESIZABLE = false

const GTK_RENDERMODEL = RenderModel()

GTK_CANVAS = Gtk.@GtkCanvas()
GTK_WINDOW = Gtk.@GtkWindow(GTK_CANVAS, "Automotive Driving Models Simulation",
    DEFAULT_GTK_CANVAS_WIDTH, DEFAULT_GTK_CANVAS_HEIGHT, DEFAULT_GTK_WINDOW_IS_RESIZABLE)
Gtk.show(GTK_CANVAS)
Gtk.visible(GTK_WINDOW, false)

GTK_SCENE = Scene()
GTK_ROADWAY = Roadway()

Gtk.@guarded Gtk.draw(GTK_CANVAS) do widget
    ctx = getgc(GTK_CANVAS)
    canvas_height = height(GTK_CANVAS)
    canvas_width = width(GTK_CANVAS)

    clear_setup!(GTK_RENDERMODEL)

    render!(GTK_RENDERMODEL, GTK_ROADWAY)
    render!(GTK_RENDERMODEL, GTK_SCENE)
    render!(GTK_RENDERMODEL, CarFollowingStatsOverlay(2), GTK_SCENE, GTK_ROADWAY)

    camera_fit_to_content!(GTK_RENDERMODEL, canvas_width, canvas_height)
    render(GTK_RENDERMODEL, ctx, canvas_width, canvas_height)
end

function show_gtk_window()
    Gtk.visible(GTK_WINDOW, true)
end
function hide_gtk_window()
    Gtk.visible(GTK_WINDOW, false)
end

function set_gtk_scene(scene::Scene)
    global GTK_SCENE = scene
end
function set_gtk_roadway(roadway::Roadway)
    global GTK_ROADWAY = roadway
end

function render_gtk_video_1d(
    veh2_actions::Vector{Float64}, # actions are for model2, just longitudinal
    scene::Scene, # initial scene, is overwritten
    model1::DriverModel{LatLonAccel, IntegratedContinuous}, # IDM
    model2::DriverModel{LatLonAccel, IntegratedContinuous}, # StaticLognitudinal
    roadway::Roadway;
    framerate::Int=10,
    )

    global GTK_ROADWAY = roadway
    actions = Array{LatLonAccel}(length(scene))
    models = Dict{Int, DriverModel}()
    models[1] = model1
    models[2] = model2

    frame_index = length(veh2_actions)
    function drawsim(t::Float64)

        frame_index += 1
        if frame_index > length(veh2_actions)
            copy!(GTK_SCENE, scene)
            frame_index = 1
        end

        actions[1] = rand(observe!(model1, GTK_SCENE, roadway, 1))

        observe!(model2, GTK_SCENE, roadway, 1)
        model2.mlon.a = veh2_actions[frame_index]
        actions[2] = rand(model2)

        tick!(GTK_SCENE, GTK_ROADWAY, actions, models)
        Gtk.draw(GTK_CANVAS)
    end

    drawsim(1.0)

    show_gtk_window()

    ticks = fps(framerate)
    timestamps = map(_ -> time(), ticks)
    map(drawsim, timestamps)
end


text:
r = PH.ScalarPenaltyParameter(25.0)
atol = 1e-8
rtol = 1e-12
max_iter = 500
obj_val = 178.3537498401004
var_vals = Dict([
    "x[1]_{0,1,2,3}" => 7.5625,
    "x[2]_{0,1,2,3}" => 0.0,
    "x[3]_{0,1,2,3}" => 1.0,
    "y_{0,1}" => 1.75,
    "y_{2,3}" => 0.0,
    "z[1]_{0}" => -0.65625,
    "z[2]_{0}" => -0.65625,
    "z[1]_{1}" => 2.84375,
    "z[2]_{1}" => 2.84375,
    "z[1]_{2}" => -1.78125,
    "z[2]_{2}" => -1.78125,
    "z[1]_{3}" => 4.71875,
    "z[2]_{3}" => 4.71875,
])

@testset "Solve Extensive" begin
    efm = PH.solve_extensive(build_scen_tree(),
                             create_model,
                             ()->Ipopt.Optimizer();
                             opt_args=(print_level=0,tol=1e-12)
                             )

    @test JuMP.num_variables(efm) == length(keys(var_vals))

    cnum = 0
    for (f,s) in JuMP.list_of_constraint_types(efm)
        cnum += JuMP.num_constraints(efm, f, s)
    end
    @test cnum == 12

    @test JuMP.termination_status(efm) == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(efm), obj_val)

    for var in JuMP.all_variables(efm)
        @test isapprox(JuMP.value(var), var_vals[JuMP.name(var)])
    end

    struct FakeSubproblem <: AbstractSubproblem end
    fake_constructor(scen::Int) = FakeSubproblem()
    @test_throws(ProgressiveHedging.UnimplementedError,
                 PH.solve_extensive(build_scen_tree(),
                                    fake_constructor,
                                    ()->Ipopt.Optimizer(),
                                    subproblem_type=FakeSubproblem)
                 )
end

@testset "Solve (Scalar)" begin
    (n, err, rerr, obj, soln, phd) = PH.solve(build_scen_tree(),
                                              create_model,
                                              r,
                                              opt=Ipopt.Optimizer,
                                              opt_args=(print_level=0,tol=1e-12),
                                              atol=atol,
                                              rtol=rtol,
                                              max_iter=max_iter,
                                              report=0,
                                              timing=false,
                                              warm_start=false)

    @test err < atol
    @test isapprox(obj, obj_val)
    @test n < max_iter

    for row in eachrow(soln)
        var = row[:variable] * "_{" * row[:scenarios] * "}"
        @test isapprox(row[:value], var_vals[var], atol=1e-7)
    end
end

@testset "Solve (Proportional)" begin
    prop_max_iter = 500
    prop_atol = 1e-8
    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ProportionalPenaltyParameter(25.0),
                                              atol=prop_atol,
                                              rtol=1e-12,
                                              max_iter=prop_max_iter,
                                              report=0,
                                              timing=false,
                                              warm_start=false
                                              )

    @test err < prop_atol
    @test isapprox(obj, 8.25, atol=1e-6)
    @test n < prop_max_iter

    for row in eachrow(soln)
        if row[:variable] == "x"
            @test isapprox(row[:value], 3.0)
        elseif row[:variable] == "u"
            @test isapprox(row[:value], 1.0)
        elseif row[:variable] == "y"
            if row[:scenarios] == "0"
                @test isapprox(row[:value], 7.0)
            else
                @test isapprox(row[:value], 0.0, atol=1e-8)
            end
        end
    end

end

@testset "Solve (SEP)" begin
    prop_max_iter = 500
    prop_atol = 1e-8
    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.SEPPenaltyParameter(),
                                              atol=prop_atol,
                                              rtol=1e-12,
                                              max_iter=prop_max_iter,
                                              report=0,
                                              timing=false,
                                              warm_start=false
                                              )

    @test err < prop_atol
    @test isapprox(obj, 8.25, atol=1e-6)
    @test n < prop_max_iter

    for row in eachrow(soln)
        if row[:variable] == "x"
            @test isapprox(row[:value], 3.0)
        elseif row[:variable] == "u"
            @test isapprox(row[:value], 1.0)
        elseif row[:variable] == "y"
            if row[:scenarios] == "0"
                @test isapprox(row[:value], 7.0)
            else
                @test isapprox(row[:value], 0.0, atol=1e-8)
            end
        end
    end

end

@testset "Warm-start" begin
    (n, err, rerr, obj, soln, phd) = PH.solve(build_scen_tree(),
                                        create_model,
                                        r,
                                        opt=Ipopt.Optimizer,
                                        opt_args=(print_level=0,tol=1e-12),
                                        atol=atol,
                                        rtol=rtol,
                                        max_iter=max_iter,
                                        report=-5,
                                        timing=false,
                                        warm_start=true)

    @test err < atol
    @test isapprox(obj, obj_val)
    @test n < max_iter

    for row in eachrow(soln)
        var = row[:variable] * "_{" * row[:scenarios] * "}"
        @test isapprox(row[:value], var_vals[var], atol=1e-7)
    end

end

@testset "Max iteration termination" begin
    max_iter = 4
    regex = r".*"
    @info "Ignore the following warning."
    (n, err, rerr, obj, soln, phd) = @test_warn(regex,
                                                PH.solve(build_scen_tree(),
                                                         create_model,
                                                         r,
                                                         opt=Ipopt.Optimizer,
                                                         opt_args=(print_level=0,tol=1e-12),
                                                         atol=atol,
                                                         rtol=rtol,
                                                         max_iter=max_iter,
                                                         report=0,
                                                         timing=false)
                                                )
    @test err > atol
    @test rerr > rtol
    @test n == max_iter
end

@testset "Relative tolerance termination" begin
    rtol = 1e-6
    (n, err, rerr, obj, soln, phd) = PH.solve(build_scen_tree(),
                                        create_model,
                                        r,
                                        opt=Ipopt.Optimizer,
                                        opt_args=(print_level=0,tol=1e-12),
                                        atol=atol,
                                        rtol=rtol,
                                        max_iter=max_iter,
                                        report=0,
                                        timing=false,
                                        warm_start=true)

    # xmax = maximum(soln[soln[:,:stage] .!= 3,:value])
    @test err > atol
    @test n < max_iter
    #@test err < rtol * xmax
    @test rerr < rtol
end

@testset "Lower Bound Algorithm" begin

    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ScalarPenaltyParameter(2.0),
                                              atol=atol,
                                              rtol=rtol,
                                              max_iter=max_iter,
                                              report=0,
                                              lower_bound=5,
                                              timing=false,
                                              warm_start=false
                                              )

    @test err < atol
    @test isapprox(obj, 8.25, atol=1e-6)
    @test n < max_iter

    lb_df = PH.lower_bounds(phd)
    for row in eachrow(lb_df)
        @test row[:bound] <= 8.25
    end

    @test size(lb_df,1) == 12
    @test isapprox(lb_df[size(lb_df,1), "absolute gap"], 0.0, atol=1e-7)
    @test isapprox(lb_df[size(lb_df,1), "relative gap"], 0.0, atol=(1e-7/8.25))

end

@testset "Lower Bound Termination" begin

    gap_tol = 1e-3

    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ScalarPenaltyParameter(2.0),
                                              gap_tol=gap_tol,
                                              atol=atol,
                                              rtol=rtol,
                                              max_iter=max_iter,
                                              report=0,
                                              lower_bound=5,
                                              timing=false,
                                              warm_start=false
                                              )

    lb_df = PH.lower_bounds(phd)

    @test err > atol
    @test rerr > rtol
    @test n < max_iter
    @test lb_df[size(lb_df,1), "relative gap"] < gap_tol

end


text:
immutable FASTA <: Format end

# FASTA Parser
# ============

function _pre_readfasta(io::AbstractString)
    seqs = split(io, '>')
    N = length(seqs) - 1

    IDS  = Array(ASCIIString, N)
    SEQS = Array(ASCIIString, N)

    for i in 1:N
        fields = split(seqs[i+1], '\n')
        IDS[i] = fields[1]
        SEQS[i] = replace(fields[2], r"\s+", "")
    end

    (IDS, SEQS)
end

function _pre_readfasta(io::IO)
    IDS  = ASCIIString[]
    SEQS = ASCIIString[]
    for (name, seq) in FastaReader{ASCIIString}(io) # FastaIO
        push!(IDS,  name)
        push!(SEQS, seq )
    end
    (IDS, SEQS)
end

function parse(io::Union{IO,AbstractString}, format::Type{FASTA}, output::Type{AnnotatedMultipleSequenceAlignment};
               generatemapping::Bool=false, useidcoordinates::Bool=false,
               deletefullgaps::Bool=true, checkalphabet::Bool=false, keepinserts::Bool=false)
    IDS, SEQS = _pre_readfasta(io)
    annot = Annotations()
    if keepinserts
        _keepinserts!(SEQS, annot)
    end
    if generatemapping
        MSA, MAP = useidcoordinates  && hascoordinates(IDS[1]) ? _to_msa_mapping(SEQS, IDS) : _to_msa_mapping(SEQS)
        setannotfile!(annot, "NCol", string(size(MSA,2)))
        setannotfile!(annot, "ColMap", join(vcat(1:size(MSA,2)), ','))
        for i in 1:length(IDS)
            setannotsequence!(annot, IDS[i], "SeqMap", MAP[i])
        end
    else
        MSA = convert(Matrix{Residue}, SEQS)
    end
    msa = AnnotatedMultipleSequenceAlignment(IndexedArray(IDS), MSA, annot)
    if checkalphabet
        deletenotalphabetsequences!(msa, SEQS)
    end
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

function parse(io::Union{IO,AbstractString}, format::Type{FASTA}, output::Type{MultipleSequenceAlignment};
               deletefullgaps::Bool=true, checkalphabet::Bool=false)
    IDS, SEQS = _pre_readfasta(io)
    msa = MultipleSequenceAlignment(IndexedArray(IDS), convert(Matrix{Residue}, SEQS))
    if checkalphabet
        deletenotalphabetsequences!(msa, SEQS)
    end
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

function parse(io::Union{IO,AbstractString}, format::Type{FASTA}, output::Type{Matrix{Residue}};
               deletefullgaps::Bool=true, checkalphabet::Bool=false)
    IDS, SEQS = _pre_readfasta(io)
    _strings_to_msa(SEQS, deletefullgaps, checkalphabet)
end

parse(io::Union{IO,AbstractString}, format::Type{FASTA}; generatemapping::Bool=false,
      useidcoordinates::Bool=false, deletefullgaps::Bool=true,
      checkalphabet::Bool=false, keepinserts::Bool=false) = parse(io, FASTA,
                                                                  AnnotatedMultipleSequenceAlignment;
                                                                  generatemapping=generatemapping,
                                                                  useidcoordinates=useidcoordinates,
                                                                  deletefullgaps=deletefullgaps,
                                                                  checkalphabet=checkalphabet,
                                                                  keepinserts=keepinserts)

# Print FASTA
# ===========

function print(io::IO, msa::AbstractMultipleSequenceAlignment, format::Type{FASTA})
    for i in 1:nsequences(msa)
        id = msa.id[i]
        seq = asciisequence(msa, i)
        println(io, string(">", id, "\n", seq))
    end
end

print(msa::MultipleSequenceAlignment, format::Type{FASTA}) = print(STDOUT, msa, FASTA)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AzureAsyncOperationResult <: SwaggerModel
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    error::Any # spec type: Union{ Nothing, Error } # spec name: error

    function AzureAsyncOperationResult(;status=nothing, error=nothing)
        o = new()
        validate_property(AzureAsyncOperationResult, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(AzureAsyncOperationResult, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type AzureAsyncOperationResult

const _property_map_AzureAsyncOperationResult = Dict{Symbol,Symbol}(Symbol("status")=>Symbol("status"), Symbol("error")=>Symbol("error"))
const _property_types_AzureAsyncOperationResult = Dict{Symbol,String}(Symbol("status")=>"String", Symbol("error")=>"Error")
Base.propertynames(::Type{ AzureAsyncOperationResult }) = collect(keys(_property_map_AzureAsyncOperationResult))
Swagger.property_type(::Type{ AzureAsyncOperationResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AzureAsyncOperationResult[name]))}
Swagger.field_name(::Type{ AzureAsyncOperationResult }, property_name::Symbol) =  _property_map_AzureAsyncOperationResult[property_name]

const _allowed_AzureAsyncOperationResult_status = ["InProgress", "Succeeded", "Failed"]

function check_required(o::AzureAsyncOperationResult)
    true
end

function validate_property(::Type{ AzureAsyncOperationResult }, name::Symbol, val)
    if name === Symbol("status")
        Swagger.validate_param(name, "AzureAsyncOperationResult", :enum, val, _allowed_AzureAsyncOperationResult_status)
    end
end


text:
# Affine scaling layer
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
#

export AffineLayer

"""
    AL = AffineLayer(nx, ny, nc; logdet=false)

 Create a layer for an affine transformation.

 *Input*:

 - `nx`, `ny, `nc`: input dimensions and number of channels

 - `logdet`: bool to indicate whether to compute the logdet

 *Output*:

 - `AL`: Network layer for affine transformation.

 *Usage:*

 - Forward mode: `Y, logdet = AL.forward(X)`

 - Inverse mode: `X = AL.inverse(Y)`

 - Backward mode: `ΔX, X = AL.backward(ΔY, Y)`

 *Trainable parameters:*

 - Scaling factor `AL.s`

 - Bias `AL.b`

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct AffineLayer <: NeuralNetLayer
    s::Parameter
    b::Parameter
    logdet::Bool
end

@Flux.functor AffineLayer

# Constructor: Initialize with nothing
function AffineLayer(nx::Int64, ny::Int64, nc::Int64; logdet=false)
    s = Parameter(glorot_uniform(nx, ny, nc))
    b = Parameter(zeros(Float32, nx, ny, nc))
    return AffineLayer(s, b, logdet)
end

# Foward pass: Input X, Output Y
function forward(X, AL::AffineLayer)

    Y = X .* AL.s.data .+ AL.b.data

    # If logdet true, return as second ouput argument
    AL.logdet == true ? (return Y, logdet_forward(AL.s)) : (return Y)
end

# Inverse pass: Input Y, Output X
function inverse(Y, AL::AffineLayer; eps::Float32=0f0)
    X = (Y .- AL.b.data) ./ (AL.s.data .+ eps)   # avoid division by 0
    return X
end

# Backward pass: Input (ΔY, Y), Output (ΔY, Y)
function backward(ΔY, Y, AL::AffineLayer; set_grad::Bool=true)
    nx, ny, n_in, batchsize = size(Y)
    X = inverse(Y, AL)
    ΔX = ΔY .* AL.s.data
    Δs = sum(ΔY .* X, dims=4)[:,:,:,1]
    if AL.logdet
        set_grad ? (Δs -= logdet_backward(AL.s)) : (Δs_ = logdet_backward(AL.s))
    end
    Δb = sum(ΔY, dims=4)[:,:,:,1]
    if set_grad
        AL.s.grad = Δs
        AL.b.grad = Δb
    else
        Δθ = [Parameter(Δs), Parameter(Δb)]
    end
    if set_grad
        return ΔX, X
    else
        AL.logdet ? (return ΔX, Δθ, X, [Parameter(Δs_), Parameter(0f0.*Δb)]) : (return ΔX, Δθ, X)
    end
end

## Jacobian-related utils

function jacobian(ΔX, Δθ, X, AL::AffineLayer)
    Y = X .* AL.s.data .+ AL.b.data
    ΔY = ΔX .* AL.s.data + X .* Δθ[1].data .+ Δθ[2].data
    if AL.logdet
        return ΔY, Y, logdet_forward(AL.s), [Parameter(logdet_hessian(AL.s).*Δθ[1].data), 0f0*Δθ[2]]
    else
        return ΔY, Y
    end
end

adjointJacobian(ΔY, Y, AL::AffineLayer) = backward(ΔY, Y, AL; set_grad=false)

## Other utils

# Clear gradients
function clear_grad!(AL::AffineLayer)
    AL.s.grad = nothing
    AL.b.grad = nothing
end

# Get parameters
get_params(AL::AffineLayer) = [AL.s, AL.b]

# Logdet
logdet_forward(s) = sum(log.(abs.(s.data)))
logdet_backward(s) = 1f0 ./ s.data
logdet_hessian(s) = -1f0 ./ s.data.^2f0


text:
@universal function abs!(x::Ubound)
  if (x.lower.flags & UNUM_SIGN_MASK) != (x.upper.flags & UNUM_SIGN_MASK)
    #number straddles zero.
    x.upper = is_inward(x.lower, x.upper) ? x.upper : x.lower
    #ensure that the upper branch is positive.
    abs!(x.upper)
    #set the lower branch to zero.
    x.lower = zero(U)
  elseif (x.upper.flags & UNUM_SIGN_MASK) != 0  #we might be negative-definite
    (x.lower, x.upper) = (x.upper, x.lower)     #swap 'em
    x.lower.flags $= UNUM_SIGN_MASK
    x.upper.flags $= UNUM_SIGN_MASK
  end  #otherwise do nothing
  x
end

@universal Base.abs(x::Ubound) = abs!(deepcopy(x))

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
" Generic basis "
abstract GenBasis{T<:Number}

" Generic 1d basis "
abstract GenBasis1d{T<:Number} <: GenBasis{T}

"""
Structure that stores 1d Basis function information.

In FEM, integrals and derivatives are calculated in standard elements,
which in the case of 1D domains is the interval (-1,1).

This type stores the information obtained from the combination
of a set of basis functions and a set of quadrature rules.
"""
type Basis1d{T<:Number,B<:BasisFun1d} <: GenBasis1d{T}
    "Number of modes"
    M::Int

    "Number of quadrature nodes"
    Q::Int

    "Quadrature nodes"
    ξ::Vector{T}

    "Quadrature weights"
    w::Vector{T}

    "Derivative matrix for functions known at quadrature nodes"
    D::Array{T,2}

    "Basis functions at quadrature nodes"
    ϕ::Array{T,2}

    "Derivative of basis functions at quadrature nodes"
    dϕ::Array{T,2}

    "Weight X Basis"
    wϕ::Array{T,2}

    "Inverse of mass matrix"
    imass::Cholesky{T}

    "Actual basis function. Used to compute the basis function"
    bas::B
end

"""
Construct a Basis from a set of basis functions and quadrature rules.
"""
function Basis1d{T<:Number, B<:BasisFun1d}(b::B, q::QuadType, ::Type{T}=Float64)
    # Create the basis function
    m = nmodes(b)
    # Obter as informações da quadratura
    Q = q.Q
    ξ = q.z
    w = q.w
    D = q.D
    ϕ = zeros(T, Q, m)
    wϕ = zeros(T, Q, m)

    # Preencher as funções de base:

    for k = 1:m
        for i=1:Q
            ϕ[i,k] = basis1d(b, ξ[i], k)
            wϕ[i,k] = w[i] * ϕ[i,k]
        end
    end

    # Calcular as derivadas:
    dϕ = D * ϕ

    # calcular a matrix de massa
    mass = ϕ' * wϕ #std_mass_matrix1d(ϕ, w)

    imass = cholfact(Hermitian(mass))
    Basis1d{T,B}(m, Q, ξ, w, D, ϕ, dϕ, wϕ, imass, b)
end

"Return quadrature nodes"
qnodes(b::Basis1d) = b.ξ

"Return number of modes"
nmodes(b::Basis1d) = b.M

"Return number of quadrature nodes"
nquad(b::Basis1d) = b.Q

"Quadrature weights"
qweights(b::Basis1d) = b.w

"Return basis function"
basis(b::Basis1d) = b.bas

"Basis functions at quadrature nodes"
qbasis(b::Basis1d) = b.ϕ

"Derivative matrix at quadrature nodes"
diffmat(b::Basis1d) = b.D

"Derivative of basis functions at quadrature nodes"
dqbasis(b::Basis1d) = b.dϕ

locmap(b::Basis1d) = b.bas.lnum

basis1d(b::GenBasis1d, x, p) = basis1d(basis(b), x, p)

basis1d!(b::GenBasis1d, x::AbstractArray, y::AbstractArray, p) = basis1d!(basis(b), x, y, p)
basis1d(b::GenBasis1d, x::AbstractArray, p) = basis1d!(b, x, similar(x), p)

#call(b::GenBasis1d, x, p) = basis1d(basis(b), x, p)

Basis1d(m::Int, q::Int) = Basis1d(ModalC01d(m), QuadType(q))
Basis1d(m::Int) = Basis1d(m, m+1)

nbndry(b::GenBasis1d) = nbndry(basis(b))
ninterior(b::GenBasis1d) = ninterior(basis(b))

bndidx(b::GenBasis1d) = bndidx(basis(b))
intidx(b::GenBasis1d) = intidx(basis(b))

seq2bi!{T}(b::GenBasis1d, x::AbstractVector{T}, y::AbstractVector{T}) = seq2bi!(basis(b).lnum, x, y)
seq2bi{T}(b::GenBasis1d, x::AbstractVector{T}) = seq2bi(basis(b).lnum, x)
bi2seq!{T}(b::GenBasis1d, x::AbstractVector{T}, y::AbstractVector{T}) = bi2seq!(basis(b).lnum, x, y)
bi2seq{T}(b::GenBasis1d, x::AbstractVector{T}) = bi2seq(basis(b).lnum, x)
bi2seq!{T}(b::GenBasis1d, xb::AbstractVector{T}, xi::AbstractVector{T}, y::AbstractVector{T}) =
    bi2seq!(basis(b).lnum, xb, xi, y)
bi2seq{T}(b::GenBasis1d, xb::AbstractVector{T}, xi::AbstractVector{T}) = bi2seq(basis(b).lnum, xb, xi)
seq2b!{T}(b::GenBasis1d, x::AbstractVector{T}, y::AbstractVector{T}) = seq2b!(basis(b).lnum, x, y)
seq2i!{T}(b::GenBasis1d, x::AbstractVector{T}, y::AbstractVector{T}) = seq2i!(basis(b).lnum, x, y)
seq2b{T}(b::GenBasis1d, x::AbstractVector{T}) = seq2b(basis(b).lnum, x)
seq2i{T}(b::GenBasis1d, x::AbstractVector{T}) = seq2i(basis(b).lnum, x)

function std_mass_matrix1d{T<:Number}(ϕ::AbstractArray{T,2}, w::AbstractVector{T})
    m = size(ϕ,2)
    Q = length(w)
    mass = zeros(T,m,m)
    for k = 1:m
        for i = k:m
            mm = 0.0
            for j = 1:Q
                mm += ϕ[j,i] * ϕ[j,k] * w[j]
            end
            mass[i,k] = mm
            mass[k,i] = mm
        end
    end

    return mass
end

""" Calculate the mass matrix

"""
mass_matrix(b::GenBasis1d) =  std_mass_matrix1d(qbasis(b), qweights(b))

"""
Project a function know at quadrature nodes to a basis
"""
function project!(b::Basis1d, f::AbstractVector, u::AbstractVector)

  ϕ = qbasis(b)
  w = qweights(b)
  Q = nquad(b)
  M = nmodes(b)
  iM = b.imass

  for k = 1:M
    F = zero(eltype(f))
    for q = 1:Q
      F += f[q] * ϕ[q,k] * w[q]
    end
    u[k] = F
  end

  A_ldiv_B!(iM, u)

  return u

end

project(b::Basis1d, f::AbstractVector) = project!(b, f, similar(f, nmodes(b)))
"""
Project a function  to a basis

"""
project(b::Basis1d, f::Function) = project(b, f(qnodes(b)))

type SEM1d{T<:Number} <: GenBasis1d
    "Number of quadrature nodes"
    Q::Int

    "Quadrature nodes"
    ξ::Vector{T}

    "Quadrature weights"
    w::Vector{T}

    "Derivative matrix for functions known at quadrature nodes"
    D::Array{T,2}

    "Actual basis function. Used to compute the basis function"
    bas::Lagrange1d{T}

end

function SEM1d{T<:Number}(n::Integer, ::Type{T}=Float64)
    q = QuadType(n, Jacobi.GLJ, T)
    SEM1d{T}(q.Q, q.z, q.w, q.D, Lagrange1d(q.z))
end

qnodes(b::SEM1d) = b.ξ
nmodes(b::SEM1d) = b.Q
nquad(b::SEM1d) = b.Q
qweights(b::SEM1d) = b.w
basis(b::SEM1d) = b.bas
qbasis{T<:Number}(b::SEM1d{T}) = eye(T,nmodes(b))
dqbasis(b::SEM1d) = b.D
diffmat(b::SEM1d) = b.D

locmap(b::SEM1d) = b.bas.lnum

project!(b::SEM1d, f::AbstractVector, u::AbstractVector) = copy!(u, f)
project(b::SEM1d, f::AbstractVector) = project!(b, f, similar(f))
project(b::SEM1d, f::Function) = project(b, f(qnodes(b)))


text:
###############################################################################
#
# Make density plot
#
###############################################################################
"""
    plot_density!(
                ax::PyObject,
                xs::Array,
                ys::Array;
                cmap::String = "viridis",
                markersize::Number = 8,
                dmax::Number = NaN
    )
    plot_density!(
                ax::PyObject,
                df::DataFrame;
                cmap::String = "viridis",
                markersize::Number = 8,
                dmax::Number = NaN
    )

Plot density plot on axis, given
- `ax` Axis to plot on
- `xs` Array of X
- `ys` Array of Y
- `cmap` Optional. Color map scheme
- `markersize` Optional. Marker size dimension, scatter size is markersize^2
- `dmax` Maximal density. If `dmax` is not NaN, use dmax as maximum density
- `df` A dataframe with column names of ("X","Y","C")
"""
function plot_density!(
            ax::PyObject,
            xs::Array,
            ys::Array;
            cmap::String = "viridis",
            markersize::Number = 8,
            dmax::Number = NaN
)
    df = calculate_density(xs, ys);
    plot_density!(ax, df; cmap=cmap, markersize=markersize, dmax=dmax);

    return nothing
end

function plot_density!(
            ax::PyObject,
            df::DataFrame;
            cmap::String = "viridis",
            markersize::Number = 8,
            dmax::Number = NaN
)
    # plot the scatter on axis, if cmax is not NaN, use cmax
    cc = (isnan(dmax) ? (df).C : min.((df).C, dmax));
    ax.scatter((df).X, (df).Y, s=markersize^2, c=cc, cmap=cmap);

    return nothing
end


text:
"""
	text(cmd0::String="", arg1=nothing; kwargs...)

Plots text strings of variable size, font type, and orientation. Various map projections are
provided, with the option to draw and annotate the map boundaries.

Full option list at [`pstext`]($(GMTdoc)pstext.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **A** | **azimuths** :: [Type => Bool]

    Angles are given as azimuths; convert them to directions using the current projection.
    ($(GMTdoc)text.html#a)
- **C** | **clearance** :: [Type => Str]

    Sets the clearance between the text and the surrounding box [15%].
    ($(GMTdoc)text.html#c)
- **D** | **offset** :: [Type => Str]

    Offsets the text from the projected (x,y) point by dx,dy [0/0].
    ($(GMTdoc)text.html#d)
- **F** | **attrib** :: [Type => Str | Tuple]

    Specify up to three text attributes (font, angle, and justification).
    ($(GMTdoc)text.html#f)
- **G** | **fill** :: [Type => Str | Number]

    Sets the shade or color used for filling the text box [Default is no fill].
    ($(GMTdoc)text.html#g)
- $(GMT.opt_Jz)
- **L** | **list** :: [Type => Bool]

    Lists the font-numbers and font-names available, then exits.
    ($(GMTdoc)text.html#l)
- **M** | **paragraph** :: [Type => Str | []]

    Paragraph mode.
    ($(GMTdoc)text.html#m)
- **N** | **no_clip** | **noclip** :: [Type => Str | []]

    Do NOT clip text at map boundaries.
    ($(GMTdoc)text.html#n)
- $(GMT.opt_P)
- **Q** | **change_case** :: [Type => Str]

    Change all text to either lower or upper case.
    ($(GMTdoc)text.html#q)
- **T** | **text_box** :: [Type => Str]

    Specify the shape of the textbox when using G and/or W.
    ($(GMTdoc)text.html#t)
- **W** | **pen** :: [Type => Str]

    Sets the pen used to draw a rectangle around the text string.
    ($(GMTdoc)text.html#w)
- **Z** | **threeD** :: [Type => Str]

    For 3-D projections: expect each item to have its own level given in the 3rd column.
    ($(GMTdoc)text.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function text(cmd0::String="", arg1=nothing; first=true, kwargs...)

    gmt_proggy = (IamModern[1]) ? "text "  : "pstext "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	N_args = (arg1 === nothing) ? 0 : 1

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd, = parse_common_opts(d, cmd, [:c :e :f :p :t :JZ :UVXY :params], first)
	cmd  = parse_these_opts(cmd, d, [[:A :azimuths], [:L :list], [:M :paragraph],
	                                 [:N :no_clip :noclip], [:Q :change_case], [:T :text_box], [:Z :threeD]])
	cmd  = add_opt(d, cmd, 'C', [:C :clearance], (margin="#", round="_+tO", concave="_+tc", convex="_+tC"))

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)
	if (isa(arg1, Array{<:Number}))
		arg1 = [GMTdataset(arg1, Array{String,1}(), "", Array{String,1}(), "", "")]
	end

	cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, [:C :color], 'C', N_args, arg1)

	cmd = add_opt(d, cmd, 'D', [:D :offset], (away=("j", nothing, 1), corners=("J", nothing, 1), shift="", line=("+v",add_opt_pen)), true)
	cmd = add_opt(d, cmd, 'F', [:F :attrib],
		(angle="+a", Angle="+A", font=("+f", font), justify="+j", region_justify="+c", header="_+h", label="_+l", rec_number="_+r", text="+t", zvalues="_+z"), true, true)
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	cmd *= add_opt_pen(d, [:W :pen], "W", true)     # TRUE to also seek (lw,lc,ls)

	if (!occursin(" -F", cmd))		# Test if the GMTdataset has text or if a numeric column is to be used as such
		if ((isa(arg1, GMTdataset) && isempty(arg1.text)) || (isa(arg1, Array{<:GMTdataset,1}) && isempty(arg1[1].text)) )
			if (isa(arg1, GMTdataset))  arg1 = [arg1]  end
			for n = 1:length(arg1)
				nr, nc = size(arg1[n].data)
				(nc < 3) && error("TEXT: input file must have at least three columns")
				arg1[n].text = Array{String,1}(undef, nr)
				for k = 1:nr
					arg1[n].text[k] = @sprintf("%.16g", arg1[n].data[k,3])
				end
			end
		end
	end

	r = finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1, arg2)
	gmt("destroy")
	return r
end

# ---------------------------------------------------------------------------------------------------
#function text(arg::DataFrame, cols...; first=true, kw...)
#	arg1 = text_record([arg.cols[1] arg.col[2]], [arg.col[3]])
#	text("", arg1; first=first, kw...)
#end
text!(cmd0::String="", arg1=nothing; kw...) = text(cmd0, arg1; first=false, kw...)
text(arg1;  kw...) = text("", arg1; first=true, kw...)
text!(arg1; kw...) = text("", arg1; first=false, kw...)

const pstext  = text			# Alias
const pstext! = text!			# Alias


text:
infPixel(step::Int) = (iseven(step)) ? true : false

function getPixelIndex(M::Matrix{Bool}, pos::CartesianIndex{2}, step::Int) :: Int
  s = size(M)
  num :: Vector{Bool} = []
  (pos[1] != 1   ) && (pos[2] != 1   ) ? push!(num, M[pos + CartesianIndex(-1,-1)]) : push!(num,infPixel(step))
  (pos[1] != 1   ) ? push!(num, M[pos + CartesianIndex(-1,0)]) : push!(num,infPixel(step))
  (pos[1] != 1   ) && (pos[2] != s[2]) ? push!(num, M[pos + CartesianIndex(-1, 1)]) : push!(num,infPixel(step))

  (pos[2] != 1   ) ? push!(num, M[pos + CartesianIndex(0,-1)]) : push!(num,infPixel(step))
  push!(num, M[pos + CartesianIndex(0,0)])
  (pos[2] != s[2]) ? push!(num, M[pos + CartesianIndex(0,1)]) : push!(num,infPixel(step))

  (pos[1] != s[1]) && (pos[2] != 1   ) ? push!(num, M[pos + CartesianIndex( 1,-1)]) : push!(num,infPixel(step))
  (pos[1] != s[1]) ? push!(num, M[pos + CartesianIndex(1,0)]) : push!(num,infPixel(step))
  (pos[1] != s[1]) && (pos[2] != s[2]) ? push!(num, M[pos + CartesianIndex( 1, 1)]) : push!(num,infPixel(step))
  return toInt(num)
end

# Converts from binary to decimal (binary given as array of bools).
function toInt(v :: Array{Bool,1})
  num = 0
  for i in 1:length(v)
    num += (v[i] ? 1 : 0) * 2^(length(v)-i)
  end
  return num
end

getOutputPixel(algo::Vector{Bool}, M::Matrix{Bool}, pos::CartesianIndex{2}, step::Int) :: Bool = algo[getPixelIndex(M,pos,step)+1]

function getOutputImage(algo::Vector{Bool}, M::Matrix{Bool}, step::Int)
  return map(p -> getOutputPixel(algo,M,p,step), CartesianIndices(M))
end

function addPadding(M::Matrix{Bool}, n::Int) M::Matrix{Bool}
  s = size(M)
  res :: Matrix{Bool} = zeros(s[1] + 2*n, s[2] + 2*n)
  res[n+1:s[1]+n,n+1:s[2]+n] = M[:,:]
  return res
end

function parseLine(s::String) :: Union{Nothing,Vector{Bool}}
  (s == "") && (return nothing)
  return map(c -> (c == '#') ? true : false , collect(s))
end

function readInput()
  data = map(line -> parseLine(line), eachline("inputs/day20.txt"))

  algo  :: Vector{Bool} = data[1]
  image :: Matrix{Bool} = transpose(hcat(data[3:end]...))
  return (algo, image)
end

function problem1(algo::Vector{Bool}, image::Matrix{Bool})::Int
  ims = []
  im = addPadding(image,10)
  push!(ims, im)

  for i in 1:2
    push!(ims, getOutputImage(algo,ims[i],i))
  end
  return count(ims[3])
end

function problem2(algo::Vector{Bool}, image::Matrix{Bool})::Int
  ims = []
  im = addPadding(image,70)
  push!(ims, im)

  for i in 1:50
    push!(ims, getOutputImage(algo,ims[i],i))
  end
  return count(ims[51])
end

algo, image = readInput()

println("Problem 1:")
res = problem1(algo,image)
display(res)

println("Problem 2:")
res = problem2(algo,image)
display(res)


text:
using LinearAlgebra
using Statistics
using PyCall
using Plots
using NCFlex
using Optim
using JuLIP
using ASE
using FiniteDifferences

@pyimport ase.units as ase_units
@pyimport ase.lattice.cubic as ase_lattice_cubic
@pyimport ase.build as ase_build
@pyimport ase.constraints as ase_constraints
@pyimport quippy.potential as quippy_potential
@pyimport matscipy.fracture_mechanics.crack as crack
@pyimport matscipy.fracture_mechanics.clusters as clusters
@pyimport matscipy.elasticity as elasticity

# Parameters

r_III = 32.0 # radius of region III. Other regions are derived from this
crack_surface = [1, 1, 1] # y-diretion, the surface which is opened
crack_front = [1, -1, 0] # z-direction, crack front line
k = 1.0 # in units of K_G
relax_elasticity = false # if true, C_ij matrix computed with internal relaxation
relax_surface = false # if true, surface energy computed with internal relaxation

imsw = quippy_potential.Potential("IP SW", param_str="""
<SW_params n_types="1">
<per_type_data type="1" atomic_num="14" />
<per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
      p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />
<per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
      lambda="42.0" gamma="1.20" eps="2.1675" />
</SW_params>
""")

unitcell = ase_build.bulk("Si", cubic=true)
unitcell.calc = imsw
sf = ase_constraints.StrainFilter(unitcell)

# energy as a function of strain
f = x -> begin
            sf.set_positions(reshape(x, 2, 3))
            sf.get_potential_energy()
         end

# gradient of energy as a function of strain
g = x -> begin
            sf.set_positions(reshape(x, 2, 3))
            -reshape(sf.get_forces(), :)
         end
x0 = reshape(sf.get_positions(), :)
res = optimize(f, g, x0, inplace=false)

# relaxed lattice constant
alat = unitcell.cell[1, 1]
@show alat
unitcell = ase_build.bulk("Si", a=alat) # remove near-zero cell entries
unitcell.calc = imsw

# 6x6 elastic constant matrix
@pyimport ase.optimize as ase_optimize
C, C_err = elasticity.fit_elastic_constants(unitcell, symmetry="cubic",
                                            optimizer= relax_elasticity ? ase_optimize.LBFGS : nothing)
@show C[1, 1], C[1, 1] / ase_units.GPa
@show C[1, 2], C[1, 2] / ase_units.GPa
@show C[4, 4], C[4, 4] / ase_units.GPa

function surface_energy(; shift=2.0, relax=true)
   bulk = ase_lattice_cubic.Diamond(symbol="Si",
                                    latticeconstant=alat,
                                    directions=[[1,-1,0],[1,1,-2],[1,1,1]]) * (1, 1, 10)
   bulk.calc = imsw
   surface = bulk.copy()
   X = surface.positions
   X[:, 3] .+= shift
   surface.set_positions(X)
   surface.wrap()
   cell = surface.cell.array
   cell[3, :] += [0.0, 0.0, 10.0]
   surface.set_cell(cell)
   surface.calc = imsw
   area = norm(cross(bulk.cell.array[:, 1], bulk.cell.array[:, 2]))

   if relax
      surface.rattle()
      f = x -> begin
         surface.set_positions(reshape(x, length(surface), 3))
         surface.get_potential_energy()
      end

      # gradient of energy as a function of strain
      g = x -> begin
         surface.set_positions(reshape(x, length(surface), 3))
         -reshape(surface.get_forces(), :)
      end
      x0 = reshape(surface.get_positions(), :)
      print(optimize(f, g, x0, inplace=false))
   end
   γ = (surface.get_potential_energy() - bulk.get_potential_energy()) / (2 * area)
end

γ = surface_energy(relax=relax_surface)
@show γ, γ / (ase_units.J / ase_units.m^2)

# build the final crystal in correct orientation
r_cut = imsw._quip_atoms.cutoff # cutoff distance for the potential
r_I = r_III - 3 * r_cut

n = [2 * floor(Int, (r_III + r_cut)/ alat),
     2 * floor(Int, (r_III + r_cut)/ alat) - 1, 1]
@show n

cryst = clusters.diamond("Si", alat, n,
                         crack_surface=crack_surface,
                         crack_front=crack_front)

cluster = clusters.set_regions(cryst, r_I, r_cut, r_III)

X = cluster.positions
region = cluster.arrays["region"]

# object for computing the CLE displacments and gradients
crk = crack.CubicCrystalCrack(crack_surface, crack_front,
                              C11=C[1,1], C12=C[1,2], C44=C[4, 4])

k_G = crk.k1g(γ)

x0, y0, _ = diag(cluster.cell.array) ./ 2

function u_cle(α)
   u1, u2 = crk.displacements(X[:, 1], X[:, 2], x0 + α, y0, k * k_G)
   u_CLE = [u1 u2 zeros(length(u1))]
   return u_CLE
end

function ∇u_cle(α)
   Du = crk.deformation_gradient(X[:, 1], X[:, 2], x0 + α, y0, k * k_G)
   ∇u_CLE = hcat(-(Du[:, 1, 1] .- 1.0),
                 -Du[:, 1, 2],
                 zeros(length(u1)))
   return ∇u_CLE
end

# check gradient wrt finite differnces
u, ∇u = u_cle(0.0), ∇u_cle(0.0)
@assert maximum((central_fdm(2, 1))(u_cle, 0.0) - ∇u) < 1e-6

scatter(X[:, 1] + u[:, 1], X[:, 2] + u[:, 2],
        color=region, aspect_ratio=:equal, label=nothing)