text:
## Uniform ##

struct TuringUniform{T} <: ContinuousUnivariateDistribution
    a::T
    b::T
end
TuringUniform() = TuringUniform(0.0, 1.0)
function TuringUniform(a::Int, b::Int)
    return TuringUniform{Float64}(Float64(a), Float64(b))
end
function TuringUniform(a::Real, b::Real)
    T = promote_type(typeof(a), typeof(b))
    return TuringUniform{T}(T(a), T(b))
end
Distributions.logpdf(d::TuringUniform, x::Real) = uniformlogpdf(d.a, d.b, x)

Base.minimum(d::TuringUniform) = d.a
Base.maximum(d::TuringUniform) = d.b

function uniformlogpdf(a, b, x)
    diff = b - a
    if a <= x <= b
        return -log(diff)
    else
        return log(zero(diff))
    end
end

## PoissonBinomial ##

struct TuringPoissonBinomial{T<:Real, TV1<:AbstractVector{T}, TV2<:AbstractVector} <: DiscreteUnivariateDistribution
    p::TV1
    pmf::TV2
end

# use the faster `poissonbinomial_pdf`
function TuringPoissonBinomial(p::AbstractArray{<:Real}; check_args = true)
    pb = Distributions.poissonbinomial_pdf(p)
    ϵ = eps(eltype(pb))
    check_args && @assert all(x -> x >= -ϵ, pb) && isapprox(sum(pb), 1; atol=ϵ)
    return TuringPoissonBinomial(p, pb)
end

function logpdf(d::TuringPoissonBinomial{T}, k::Int) where T<:Real
    insupport(d, k) ? log(d.pmf[k + 1]) : -T(Inf)
end
quantile(d::TuringPoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1
Base.minimum(d::TuringPoissonBinomial) = 0
Base.maximum(d::TuringPoissonBinomial) = length(d.p)


text:
abstract type AbstractMetric end

_string_M⁻¹(mat::AbstractMatrix, n_chars::Int=32) = _string_M⁻¹(diag(mat), n_chars)
function _string_M⁻¹(vec::AbstractVector, n_chars::Int=32)
    s_vec = string(vec)
    l = length(s_vec)
    s_dots = " ...]"
    n_diag_chars = n_chars - length(s_dots)
    return s_vec[1:min(n_diag_chars,end)] * (l > n_diag_chars ? s_dots : "")
end

struct UnitEuclideanMetric{T,A<:Union{Tuple{Int},Tuple{Int,Int}}} <: AbstractMetric
    M⁻¹::UniformScaling{T}
    size::A
end

UnitEuclideanMetric(::Type{T}, sz) where {T} = UnitEuclideanMetric(UniformScaling{T}(one(T)), sz)
UnitEuclideanMetric(sz) = UnitEuclideanMetric(Float64, sz)
UnitEuclideanMetric(::Type{T}, dim::Int) where {T} = UnitEuclideanMetric(UniformScaling{T}(one(T)), (dim,))
UnitEuclideanMetric(dim::Int) = UnitEuclideanMetric(Float64, (dim,))

renew(ue::UnitEuclideanMetric, M⁻¹) = UnitEuclideanMetric(M⁻¹, ue.size)

Base.size(e::UnitEuclideanMetric) = e.size
Base.size(e::UnitEuclideanMetric, dim::Int) = e.size[dim]
Base.show(io::IO, uem::UnitEuclideanMetric) = print(io, "UnitEuclideanMetric($(_string_M⁻¹(ones(uem.size))))")

struct DiagEuclideanMetric{T,A<:AbstractVecOrMat{T}} <: AbstractMetric
    # Diagnal of the inverse of the mass matrix
    M⁻¹     ::  A
    # Sqare root of the inverse of the mass matrix
    sqrtM⁻¹ ::  A
    # Pre-allocation for intermediate variables
    _temp   ::  A
end

function DiagEuclideanMetric(M⁻¹::AbstractVecOrMat{T}) where {T<:AbstractFloat}
    return DiagEuclideanMetric(M⁻¹, sqrt.(M⁻¹), similar(M⁻¹))
end
DiagEuclideanMetric(::Type{T}, sz) where {T} = DiagEuclideanMetric(ones(T, sz...))
DiagEuclideanMetric(sz) = DiagEuclideanMetric(Float64, sz)
DiagEuclideanMetric(::Type{T}, dim::Int) where {T} = DiagEuclideanMetric(ones(T, dim))
DiagEuclideanMetric(dim::Int) = DiagEuclideanMetric(Float64, dim)

renew(ue::DiagEuclideanMetric, M⁻¹) = DiagEuclideanMetric(M⁻¹)

Base.size(e::DiagEuclideanMetric, dim...) = size(e.M⁻¹, dim...)
Base.show(io::IO, dem::DiagEuclideanMetric) = print(io, "DiagEuclideanMetric($(_string_M⁻¹(dem.M⁻¹)))")

struct DenseEuclideanMetric{
    T,
    AV<:AbstractVecOrMat{T},
    AM<:Union{AbstractMatrix{T},AbstractArray{T,3}},
    TcholM⁻¹<:UpperTriangular{T},
} <: AbstractMetric
    # Inverse of the mass matrix
    M⁻¹::AM
    # U of the Cholesky decomposition of the mass matrix
    cholM⁻¹::TcholM⁻¹
    # Pre-allocation for intermediate variables
    _temp::AV
end

# TODO: make dense mass matrix support matrix-mode parallel
function DenseEuclideanMetric(M⁻¹::Union{AbstractMatrix{T},AbstractArray{T,3}}) where {T<:AbstractFloat}
    _temp = Vector{T}(undef, Base.front(size(M⁻¹)))
    return DenseEuclideanMetric(M⁻¹, cholesky(Symmetric(M⁻¹)).U, _temp)
end
DenseEuclideanMetric(::Type{T}, D::Int) where {T} = DenseEuclideanMetric(Matrix{T}(I, D, D))
DenseEuclideanMetric(D::Int) = DenseEuclideanMetric(Float64, D)
DenseEuclideanMetric(::Type{T}, sz::Tuple{Int}) where {T} = DenseEuclideanMetric(Matrix{T}(I, first(sz), first(sz)))
DenseEuclideanMetric(sz::Tuple{Int}) = DenseEuclideanMetric(Float64, sz)

renew(ue::DenseEuclideanMetric, M⁻¹) = DenseEuclideanMetric(M⁻¹)

Base.size(e::DenseEuclideanMetric, dim...) = size(e._temp, dim...)
Base.show(io::IO, dem::DenseEuclideanMetric) = print(io, "DenseEuclideanMetric(diag=$(_string_M⁻¹(dem.M⁻¹)))")

# `rand` functions for `metric` types.

function _rand(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    metric::UnitEuclideanMetric{T}
) where {T}
    r = randn(rng, T, size(metric)...)
    return r
end

function _rand(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    metric::DiagEuclideanMetric{T}
) where {T}
    r = randn(rng, T, size(metric)...)
    r ./= metric.sqrtM⁻¹
    return r
end

function _rand(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    metric::DenseEuclideanMetric{T}
) where {T}
    r = randn(rng, T, size(metric)...)
    ldiv!(metric.cholM⁻¹, r)
    return r
end

Base.rand(rng::AbstractRNG, metric::AbstractMetric) = _rand(rng, metric)    # this disambiguity is required by Random.rand
Base.rand(rng::AbstractVector{<:AbstractRNG}, metric::AbstractMetric) = _rand(rng, metric)
Base.rand(metric::AbstractMetric) = rand(GLOBAL_RNG, metric)

ext: jl


text:
using PerlaTonettiWaugh
using Test, LinearAlgebra, Statistics, Compat
using Sundials, BenchmarkTools, Interpolations, Parameters, SimpleDifferentialOperators
using DataFrames, DataFramesMeta

@elapsed begin
  @time @testset "Simple Stationary" begin include("simple/stationarytest.jl") end
  @time @testset "Full Stationary" begin include("full/stationarytest.jl") end
  @time @testset "Simple Dynamic" begin include("simple/dynamictest.jl") end
  @time @testset "Full Dynamic" begin include("full/dynamictest.jl") end
  @time @testset "Type Stability and Speed" begin include("benchmarking/inference.jl") end
end


text:
using OptimizationParameters
using Test

@testset "OptimizationParameters.jl" begin

    # dictionary read and write
    dict = read_parameters("example.csv", dict=true)
    write_parameters("tmp.csv", dict)
    dict2 = read_parameters("tmp.csv", dict=true)
    for key in keys(dict)
        for field in fieldnames(OptimizationParameter)
            @test getproperty(dict[key], field) == getproperty(dict2[key], field)
        end
    end
    write_parameters("example.csv", "tmp.csv", dict)
    dict3 = read_parameters("tmp.csv", dict=true)
    for key in keys(dict)
        for field in fieldnames(OptimizationParameter)
            @test getproperty(dict[key], field) == getproperty(dict3[key], field)
        end
    end

    # named tuple read and write
    nt = read_parameters("example.csv", dict=false)
    write_parameters("tmp.csv", nt)
    nt2 = read_parameters("tmp.csv", dict=false)
    for key in keys(nt)
        for field in fieldnames(OptimizationParameter)
            @test getproperty(nt[key], field) == getproperty(nt2[key], field)
        end
    end
    write_parameters("example.csv", "tmp.csv", nt)
    nt3 = read_parameters("tmp.csv", dict=false)
    for key in keys(nt)
        for field in fieldnames(OptimizationParameter)
            @test getproperty(nt[key], field) == getproperty(nt3[key], field)
        end
    end

    # remove temporary file
    rm("tmp.csv")

    # get and set methods for OptimizationParameter
    optparam = OptimizationParameter(1)
    optparam = set_x0(optparam, 5)
    optparam = set_lb(optparam, 0)
    optparam = set_ub(optparam, 10)
    optparam = set_scaling(optparam, 0.1)
    optparam = set_dv(optparam, true)
    optparam = set_description(optparam, "test variable")

    @test get_x0(optparam) == 5
    @test get_lb(optparam) == 0
    @test get_ub(optparam) == 10
    @test get_scaling(optparam) == 0.1
    @test get_dv(optparam) == true
    @test get_description(optparam) == "test variable"

    # get and set methods for dictionaries
    tmp1 = copy(dict)
    set_x0!(tmp1, :scalar1, 5)
    set_lb!(tmp1, :scalar1, 0)
    set_ub!(tmp1, :scalar1, 10)
    set_scaling!(tmp1, :scalar1, 0.1)
    set_dv!(tmp1, :scalar1, true)
    set_description!(tmp1, :scalar1, "test variable")

    @test get_x0(tmp1, :scalar1) == 5
    @test get_lb(tmp1, :scalar1) == 0
    @test get_ub(tmp1, :scalar1) == 10
    @test get_scaling(tmp1, :scalar1) == 0.1
    @test get_dv(tmp1, :scalar1) == true
    @test get_description(tmp1, :scalar1) == "test variable"

    tmp2 = copy(dict)
    tmp2 = set_x0(tmp2, :scalar1, 5)
    tmp2 = set_lb(tmp2, :scalar1, 0)
    tmp2 = set_ub(tmp2, :scalar1, 10)
    tmp2 = set_scaling(tmp2, :scalar1, 0.1)
    tmp2 = set_dv(tmp2, :scalar1, true)
    tmp2 = set_description(tmp2, :scalar1, "test variable")

    @test get_x0(tmp2, :scalar1) == 5
    @test get_lb(tmp2, :scalar1) == 0
    @test get_ub(tmp2, :scalar1) == 10
    @test get_scaling(tmp2, :scalar1) == 0.1
    @test get_dv(tmp2, :scalar1) == true
    @test get_description(tmp2, :scalar1) == "test variable"

    # get and set methods for named tuples
    tmp3 = nt
    tmp3 = set_x0(tmp3, :scalar1, 5)
    tmp3 = set_lb(tmp3, :scalar1, 0)
    tmp3 = set_ub(tmp3, :scalar1, 10)
    tmp3 = set_scaling(tmp3, :scalar1, 0.1)
    tmp3 = set_dv(tmp3, :scalar1, true)
    tmp3 = set_description(tmp3, :scalar1, "test variable")

    @test get_x0(tmp3, :scalar1) == 5
    @test get_lb(tmp3, :scalar1) == 0
    @test get_ub(tmp3, :scalar1) == 10
    @test get_scaling(tmp3, :scalar1) == 0.1
    @test get_dv(tmp3, :scalar1) == true
    @test get_description(tmp3, :scalar1) == "test variable"

    # update parameters (still need proper tests here)
    p_dict = get_values(dict)
    p_nt = get_values(nt)
    tmp = update_parameters(nt, p_dict)
    tmp = update_parameters(dict, p_nt)
    tmp = update_parameters(dict, p_dict)
    tmp = update_parameters(nt, p_nt)
    tmp = update_parameters!(dict, p_dict)

    # update_design variables (still need proper tests here
    tmp = update_design_variables(nt, dict)
    tmp = update_design_variables(dict, nt)
    tmp = update_design_variables(dict, dict2)
    tmp = update_design_variables(nt, nt2)
    tmp = update_design_variables!(dict, dict2)

    # assemble input (still need proper tests here)
    x0, lb, ub = assemble_input(nt)
    x0, lb, ub = assemble_input(dict)

    # extract parameter values (still need proper test here)
    parameters = get_values(nt, x0)
    parameters = get_values(dict, x0)

    # test that offset works properly
    lb = 0
    ub = 10
    optparam = OptimizationParameter(0; lb=lb, ub=ub, dv=true)
    @test OptimizationParameters.get_value(optparam, 0) == (lb + ub)/2

end


text:
function ga_v2(;
        model::Model,
        nth_param_set::Int64,
        max_generation::Int64,
        n_population::Int64,
        n_gene::Int64,
        initial_threshold::Float64,
        allowable_error::Float64,
        local_search_method::String,
        n_children::Int64,
        maxiter::Int64)::Nothing
    """
        1. Initialization
            As an initial population, create np individuals randomly.
            ga_v2 also represents individuals as n-dimensional real number
            vectors, where n is the dimension of the search space. Set
            Generation to 0, and set the iteration number of converging
            operations Niter to 1.
        2. Selection for reproduction
            As parents for the recombination operator, ENDX, select m
            individuals, p1, p2, . . . ,pm, without replacement from the
            population.
        3. Generation of offsprings
            Generate Nc children by applying ENDX to the selected parents.
            This algorithm assigns the worst objective value to the children.
        4. Local Search
            Apply the local search method to the best individual in a family
            consisting of the two parents, i.e., p1 and p2, and their children.
            Note here that the children are assumed to have the worst objective
            value. Thus, whenever the objective values of the two parents have
            been actually computed in previous generations, the algorithm
            applies the local search to either of the parents. When all of the
            individuals in the family have the same objective value, on the
            other hand, the local search is applied to a randomly selected
            individual from the family.
        5. Selection for survival
            Select two individuals from the family. The first selected
            individual should be the individual with the best objective value,
            and the second should be selected randomly. Then, replace the two
            parents (p1 and p2) with the selected individuals. Note that the
            individual to which the local search has been applied in the
            previous step is always selected as the best.
        6. Application of ENDX/MGG
            To achieve a good search performance, ga_v2 optimizes a function,
            gradually narrowing the search space. For this purpose, the
            converging phase slightly converges the population by repeating the
            following procedure Niter times.
            (i) Select m individuals without replacement from the population.
                The selected individuals, expressed here as p1, p2, . . . , pm,
                are used as the parents for an extended normal distribution
                crossover (ENDX) applied in the next step.
            (ii) Generate Nc children by applying ENDX to the parents selected
                in the previous step. To reduce the computational cost, ga_v2
                forgoes any computation of the objective values of the Nc
                individuals generated here. Instead, the algorithm assigns the
                newly generated children a single objective value, one which is
                inferior to the objective values of any of the possible
                candidate solutions.
            (iii) Select two individuals from a family containing the two
                parents, i.e., p1 and p2, and their children. The first
                selected individual should be the one with the best objective
                value, and the second should be selected randomly. Then,
                replace the two parents with the selected individuals.
        7. Adaptation of Niter
            If the best individual has not improved during the last np
            generations, Niter <- 2 * Niter. Otherwise, set Niter to 1.
        8. Termination
            Stop if the halting criteria are satisfied.
            Otherwise, Generation <- Generation + 1, and return to the step 2.
    """
    if n_population < n_gene + 2
        error("n_population must be larger than $(n_gene + 2)")
    end

    n_iter::Int64 = 1
    N0::Vector{Float64} = zeros(3 * n_population)

    population::Matrix{Float64} = get_initial_population(
        model, nth_param_set, n_population, n_gene, initial_threshold
    )
    N0[1] = population[1, end]
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ), "a"
    ) do f
        write(f,
            @sprintf(
                "Generation%d: Best Fitness = %.6e\n",
                1, population[1, end]
            )
        )
    end
    best_indiv::Vector{Float64} = model.gene2val(population[1, 1:n_gene])
    best_fitness::Float64 = population[1, end]

    open(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "fit_param1.dat"
        ), "w"
    ) do f
        for val in best_indiv
            write(f, @sprintf("%.6e\n", val))
        end
    end
    open(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "generation.dat"
        ), "w"
    ) do f
        write(f, @sprintf("%d", 1))
    end
    open(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "best_fitness.dat"
        ), "w"
    ) do f
        write(f, @sprintf("%.6e", best_fitness))
    end

    if population[1, end] <= allowable_error
        return
    end

    generation::Int64 = 2
    while generation <= max_generation
        ip::Vector{Int} = sample(
            collect(1:n_population), n_gene + 2, replace=false
        )
        population .= converging!(
            model.obj_func, ip, population, n_population, n_gene
        )
        population .= local_search!(
            model.obj_func, ip, population, n_population, n_gene,
            method=local_search_method, n_children=n_children, maxiter=maxiter
        )
        if n_iter > 1
            for _ in 1:n_iter
                ip .= sample(
                    collect(1:n_population), n_gene + 2, replace=false
                )
                population .= converging!(
                    model.obj_func, ip, population, n_population, n_gene
                )
            end
        end

        if generation % length(N0) == 0
            N0[end] = population[1, end]
            if N0[1] == N0[end]
                n_iter *= 2
            else
                n_iter = 1
            end
        else
            N0[generation % length(N0)] = population[1, end]
        end
        open(
            joinpath(
                model.path,
                "logs",
                "$nth_param_set.log"
            ), "a"
        ) do f
            write(f,
                @sprintf(
                    "Generation%d: Best Fitness = %.6e\n",
                    generation, population[1, end]
                )
            )
        end
        best_indiv .= model.gene2val(population[1, 1:n_gene])
        if population[1, end] < best_fitness
            open(
                joinpath(
                    model.path,
                    "fitparam",
                    "$nth_param_set",
                    "fit_param$generation.dat"
                ), "w"
            ) do f
                for val in best_indiv
                    write(f, @sprintf("%.6e\n", val))
                end
            end
            open(
                joinpath(
                    model.path,
                    "fitparam",
                    "$nth_param_set",
                    "generation.dat"
                ), "w"
            ) do f
                write(f, @sprintf("%d", generation))
            end
        end
        best_fitness = population[1, end]

        open(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "best_fitness.dat"
            ), "w"
        ) do f
            write(f, @sprintf("%.6e", best_fitness))
        end

        open(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "count_num.dat"
            ), "w"
        ) do f
            write(f, @sprintf("%d", generation))
        end

        if population[1, end] <= allowable_error
            break
        end

        generation += 1
    end

    return
end

function ga_v2_continue(;
        model::Model,
        nth_param_set::Int64,
        max_generation::Int64,
        n_population::Int64,
        n_gene::Int64,
        initial_threshold::Float64=1e12,
        allowable_error::Float64,
        local_search_method::String,
        n_children::Int64,
        maxiter::Int64,
        p0_bounds::Vector{Float64})::Nothing
    if n_population < n_gene + 2
        error("n_population must be larger than $(n_gene + 2)")
    end

    n_iter::Int64 = 1
    N0::Vector{Float64} = zeros(3 * n_population)

    count::Int64 = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "count_num.dat"
        )
    )[1, 1]
    best_generation::Int64 = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "generation.dat"
        )
    )[1, 1]
    best_indiv::Vector{Float64} = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "fit_param$best_generation.dat"
        )
    )[:, 1]
    best_indiv_gene::Vector{Float64} = model.val2gene(best_indiv)
    best_fitness::Float64 = model.obj_func(best_indiv_gene)

    population::Matrix{Float64} = get_initial_population_continue(
        model, nth_param_set, n_population, n_gene, initial_threshold, p0_bounds
    )
    if best_fitness < population[1, end]
        for i in 1:n_gene
            @inbounds population[1, i] = best_indiv_gene[i]
        end
        population[1, end] = best_fitness
    else
        best_indiv = model.gene2val(population[1, 1:n_gene])
        best_fitness = population[1, end]
        open(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "fit_param$count.dat"
            ), "w"
        ) do f
            for i = 1:n_gene
                write(f, @sprintf("%.6e", best_indiv[i]))
            end
        end
    end

    N0[1] = population[1, end]
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ), "a"
    ) do f
        write(f,
            @sprintf(
                "Generation%d: Best Fitness = %.6e\n",
                count + 1, population[1, end]
            )
        )
    end

    if population[1, end] <= allowable_error
        return
    end

    generation::Int64 = 2 + count
    while generation <= max_generation
        ip::Vector{Int} = sample(
            collect(1:n_population), n_gene + 2, replace=false
        )
        population .= converging!(
            model.obj_func, ip, population, n_population, n_gene
        )
        population .= local_search!(
            model.obj_func, ip, population, n_population, n_gene,
            method=local_search_method, n_children=n_children, maxiter=maxiter
        )
        if n_iter > 1
            for _ in 1:n_iter
                ip .= sample(
                    collect(1:n_population), n_gene + 2, replace=false
                )
                population .= converging!(
                    model.obj_func, ip, population, n_population, n_gene
                )
            end
        end

        if generation % length(N0) == 0
            N0[end] = population[1, end]
            if N0[1] == N0[end]
                n_iter *= 2
            else
                n_iter = 1
            end
        else
            N0[generation % length(N0)] = population[1, end]
        end
        open(
            joinpath(
                model.path,
                "logs",
                "$nth_param_set.log"
            ), "a"
        ) do f
            write(f,
                @sprintf(
                    "Generation%d: Best Fitness = %.6e\n",
                    generation, population[1, end]
                )
            )
        end
        best_indiv .= model.gene2val(population[1, 1:n_gene])
        if population[1, end] < best_fitness
            open(
                joinpath(
                    model.path,
                    "fitparam",
                    "$nth_param_set",
                    "fit_param$generation.dat"
                ), "w"
            ) do f
                for val in best_indiv
                    write(f, @sprintf("%.6e\n", val))
                end
            end
            open(
                joinpath(
                    model.path,
                    "fitparam",
                    "$nth_param_set",
                    "generation.dat"
                ), "w"
            ) do f
                write(f, @sprintf("%d", generation))
            end
        end
        best_fitness = population[1, end]

        open(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "best_fitness.dat"
            ), "w"
        ) do f
            write(f, @sprintf("%.6e", best_fitness))
        end

        open(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "count_num.dat"
            ), "w"
        ) do f
            write(f, @sprintf("%d", generation))
        end

        if population[1, end] <= allowable_error
            break
        end

        generation += 1
    end

    return
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    partition!(mesh::ParallelTreeMesh, allow_coarsening=true)

Partition `mesh` using a static domain decomposition algorithm
based on leaf cell count and tree structure.
If `allow_coarsening` is `true`, the algorithm will keep leaf cells together
on one rank when needed for local coarsening (i.e. when all children of a cell are leaves).
"""
function partition!(mesh::ParallelTreeMesh; allow_coarsening=true)
  # Determine number of leaf cells per rank
  leaves = leaf_cells(mesh.tree)
  @assert length(leaves) > mpi_nranks() "Too many ranks to properly partition the mesh!"
  n_leaves_per_rank = OffsetArray(fill(div(length(leaves), mpi_nranks()), mpi_nranks()),
                                  0:(mpi_nranks() - 1))
  for d in 0:(rem(length(leaves), mpi_nranks()) - 1)
    n_leaves_per_rank[d] += 1
  end
  @assert sum(n_leaves_per_rank) == length(leaves)

  # Assign MPI ranks to all cells such that all ancestors of each cell - if not yet assigned to a
  # rank - belong to the same rank
  mesh.first_cell_by_rank = similar(n_leaves_per_rank)
  mesh.n_cells_by_rank = similar(n_leaves_per_rank)

  leaf_count = 0
  mesh.first_cell_by_rank[0] = 1
  # Iterate over all ranks
  for d in 0:(mpi_nranks() - 1)
    leaf_count += n_leaves_per_rank[d]
    last_id = leaves[leaf_count]
    parent_id = mesh.tree.parent_ids[last_id]

    # Check if all children of the last parent are leaves
    if allow_coarsening &&
        all(id -> is_leaf(mesh.tree, id), @view mesh.tree.child_ids[:, parent_id]) &&
        d < length(n_leaves_per_rank) - 1

      # To keep children of parent together if they are all leaves,
      # all children are added to this rank
      additional_cells = (last_id+1):mesh.tree.child_ids[end, parent_id]
      if length(additional_cells) > 0
        last_id = additional_cells[end]

        additional_leaves = count(id -> is_leaf(mesh.tree, id), additional_cells)
        leaf_count += additional_leaves
        # Add leaves to this rank, remove from next rank
        n_leaves_per_rank[d] += additional_leaves
        n_leaves_per_rank[d+1] -= additional_leaves
      end
    end

    @assert all(n -> n > 0, n_leaves_per_rank) "Too many ranks to properly partition the mesh!"

    mesh.n_cells_by_rank[d] = last_id - mesh.first_cell_by_rank[d] + 1
    mesh.tree.mpi_ranks[mesh.first_cell_by_rank[d]:last_id] .= d

    # Set first cell of next rank
    if d < length(n_leaves_per_rank) - 1
      mesh.first_cell_by_rank[d+1] = mesh.first_cell_by_rank[d] + mesh.n_cells_by_rank[d]
    end
  end

  @assert all(x->x >= 0, mesh.tree.mpi_ranks[1:length(mesh.tree)])
  @assert sum(mesh.n_cells_by_rank) == length(mesh.tree)

  return nothing
end

function get_restart_mesh_filename(restart_filename, mpi_parallel::Val{true})
  # Get directory name
  dirname, _ = splitdir(restart_filename)

  if mpi_isroot()
    # Read mesh filename from restart file
    mesh_file = ""
    h5open(restart_filename, "r") do file
      mesh_file = read(attributes(file)["mesh_file"])
    end

    buffer = Vector{UInt8}(mesh_file)
    MPI.Bcast!(Ref(length(buffer)), mpi_root(), mpi_comm())
    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
  else # non-root ranks
    count = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())
    buffer = Vector{UInt8}(undef, count[])
    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
    mesh_file = String(buffer)
  end

  # Construct and return filename
  return joinpath(dirname, mesh_file)
end

end # @muladd


text:
include("bool_cube.jl")
include("bool_fun.jl")
include("ANF.jl")
include("bool_gen.jl")
include("boolean_utils.jl")

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
"""
Apply different weighting functions to `TD`. Use this method to create a
data preconditioning matrix

# Arguments Modified

* `dw::TD` :

# Keyword Arguments

* `offsetlim::Vector{Float64}=[-Inf,Inf]` : [xoffsetlim, zoffsetlim], where the records with offsets > `offsetlim` are given zero weight
* `tlim::Vector{Float64}=[dw.tgrid[1], dw.tgrid[end]]` : [tminimum, tmaximum], time mute window
* `offsetpow::Vector{Float64}=[0.0,0.0]` :
* `tpow::Float64=0.0` :
* `ttaperperc::Float64=0.` : taper window percentage for time window

* NOTE: if more than one simultaneous source are present, their mean position is considered to calculate offset.
"""
function TD_weight!(
		    dw::TD;
		    offsetlim::Vector{Float64}=[-Inf,Inf],
		    tlim::Vector{Float64}=[dw.tgrid[1], dw.tgrid[end]],
		    offsetpow::Vector{Float64}=[0.0,0.0],
		    tpow::Float64=0.0,
		    ttaperperc::Float64=0.,
		  )

	tvecexp = dw.tgrid - dw.tgrid[1]
	tmaxI = maximum(inv.(abs.(tvecexp)))
	nt = length(dw.tgrid)
	fields=dw.fields
	nss=dw.ageom.nss
	rx=dw.ageom.r[:x]; rz=dw.ageom.r[:z]
	sx=dw.ageom.s[:x]; sz=dw.ageom.s[:z]
	nr=dw.ageom.nr; ns=dw.ageom.ns

	itlim = sort(broadcast(indmin,[abs.(dw.tgrid-tlim[i]) for i in 1:2]))
	twin=zeros(nt)
	twin[itlim[1] : itlim[2]] = Utils.taper(ones(itlim[2]-itlim[1]+1),ttaperperc)

	for ifield = 1:length(fields), iss = 1:nss
		zo = sqrt.((rz[iss][:]-mean(sz[iss])).^2) # offsets computed using mean of source position
		xo = sqrt.((rx[iss][:]-mean(sx[iss])).^2)
		zomaxI = maximum(zo)^(-1)
		xomaxI = maximum(xo)^(-1)
		for ir = 1:nr[iss]
			inoffsetlim = ( (abs.(zo[ir]) < abs.(offsetlim[2])) & (abs.(xo[ir]) < abs.(offsetlim[1])) )
			#if(inoffsetlim)
			#	# apply xoffset weighting
			#	if(!isinf(xomaxI))
			#		dw.d[iss,ifield][:,ir] .*= exp(xo[ir]*offsetpow[1]*xomaxI)
			#	end
			#	# apply zoffset weighting
			#	if(!isinf(zomaxI))
			#		dw.d[iss,ifield][:,ir] .*= exp(zo[ir]*offsetpow[2]*zomaxI)
			#	end
			#else
			#	dw.d[iss,ifield][:,ir] = 0.
			#end
			for it=1:nt
				dw.d[iss,ifield][it,ir] *= twin[it] #* exp(tvecexp[it]*tpow*tmaxI)
			end
		end
	end

end


text:
import Base: |>
export button,
       slider,
       checkbox,
       togglebutton,
       textinput,
       dropdown,
       labelfor

# A widget can signal some state
abstract Widget <: Tile

## Button

immutable Button <: Widget
    label::Tile
    raised::Bool
    noink::Bool
end

button(label; name=:_button, raised=false, noink=false) =
    clickable(Button(label, raised, noink), name=name)

render(b::Button) =
    Elem("paper-button", render(b.label),
        raised=boolattr(b.raised, "raised"), noink=boolattr(b.raised, "raised"))

## Slider

immutable Slider{T <: Real} <: Widget
    value::T
    range::Range{T}
    editable::Bool
    pin::Bool
    disabled::Bool
    secondaryprogress::T
end

slider{T}(range::Range{T};
          name=:_slider,
          value=first(range),
          editable=true,
          pin=false,
          disabled=false,
          secondaryprogress=zero(T)) =
    hasstate(
        Slider(convert(T, value), range, editable, pin, disabled, secondaryprogress),
        name=name)

render(s::Slider) =
    Elem("paper-slider",
        min=first(s.range),
        max=last(s.range),
        step=step(s.range),
        value=s.value,
        editable=s.editable,
        disabled=s.disabled,
        secondaryProgress=s.secondaryprogress)

## Boolean widgets: Checkbox and Toggle Button

immutable BoolWidget{typ} <: Widget
    value::Bool
    label::String
    disabled::Bool
end

checkbox(;
         name=:_checkbox,
         value=false,
         label="",
         disabled=false) =
    hasstate(
        BoolWidget{:checkbox}(value, label, disabled),
        name=name, attr="checked", trigger="change")

render(c::BoolWidget{:checkbox}) =
    Elem("paper-checkbox",
        checked=c.value,
        disabled=boolattr(c.disabled, "disabled"))

togglebutton(;
             name=:_togglebutton,
             value=false,
             label="",
             disabled=false) =
    hasstate(
        BoolWidget{:toggle}(value, label, disabled),
        name=name, attr="checked", trigger="change")

render(t::BoolWidget{:toggle}) =
    Elem("paper-toggle-button",
        checked=t.value,
        disabled=boolattr(t.disabled, "disabled"))

## Text input

immutable TextInput <: Widget
    value::String
    label::String
    floatinglabel::Bool
    disabled::Bool
end

textinput(value::String="";
          label="",
          name=:_textinput,
          floatinglabel=false,
          disabled=false) =
    hasstate(
        TextInput(value, label, floatinglabel, disabled),
        name=name, attr="value", trigger="keyup")

render(t::TextInput) =
    Elem("paper-input",
        label=t.label,
        value=t.value,
        floatingLabel=boolattr(t.floatinglabel, "floatingLabel"),
        disabled=boolattr(t.disabled, "disabled"))

## Dropdown

immutable SelectionItem{T} <: Tile
    value::T
    item::Tile
end

render(t::SelectionItem) =
    Elem("paper-item", render(t.tile), value=t.value)

immutable Dropdown <: Widget
    value::String
    label::String
    items::Vector{SelectionItem}
    disabled::Bool
end

render(d::Dropdown) =
    Elem("paper-dropdown-menu",
        value=d.value,
        label=d.label,
        floatingLabel=boolattr(d.floatinglabel, "floatingLabel"),
        disabled=boolattr(d.disabled, "disabled")) |>
    (wrap -> reduce(<<, wrap, map(render, d.items)))

## TODO: Look into this stuff again!
makeitems(xs) =
    [SelectionItem(k, v) for (k, v) in xs]

dropdown(items::AbstractArray;
         name=:_dropdown,
         value=first(items)) =
    hasstate(Dropdown(makeitems(items), value), name=name)

# label a widget
immutable Label <: Widget
    target::Symbol
    label::Tile
end


text:
using Base.Test
# symbols.jl
import WaveletScattering: appendsymbols

@test_throws AssertionError appendsymbols(Symbol[], 1)
@test_throws AssertionError appendsymbols([:time], 0)
@test appendsymbols([:time], 1) == [:time]
@test appendsymbols([:time], 2) == [:time, :var1]
@test appendsymbols([:time], 3) == [:time, :var1, :var2]
@test appendsymbols([:time, :space], 2) == [:time, :space]
@test appendsymbols([:time, :space], 3) == [:time, :space, :var1]


text:
import Random

"""
  RandomSampler(dataset)

Random Sampler that enables to get index without any modifications of your dataset.

`dataset` should satisfy `IndexStyle(dataset) == IndexLinear()`. `getindex(dataset, index::Int)` and `size(dataset)` is required too.
"""
struct RandomSampler
  dataset
  transform::Array{Int64, 1}

  function RandomSampler(::IndexLinear, dataset)
    instance = new(dataset, collect(1:length(dataset)))
    Random.shuffle!(instance)
    instance
  end
end

function RandomSampler(dataset)
  RandomSampler(IndexStyle(dataset), dataset)
end

function Base.getindex(sampler::RandomSampler, index::Int)
  sampler.dataset[sampler.transform[index]]
end

Base.size(sampler::RandomSampler) = (length(sampler),)

Base.length(sampler::RandomSampler) = length(sampler.dataset)

Base.IndexStyle(::RandomSampler) = IndexLinear()

function Random.shuffle!(sampler::RandomSampler)
  Random.shuffle!(sampler.transform)
  nothing
end

function Random.shuffle!(rng::Random.AbstractRNG, sampler::RandomSampler)
  Random.shuffle!(rng, sampler.transform)
  nothing
end


text:
using DataDeps

register(DataDep(
  "heart-disease-hungarian",
  "http://archive.ics.uci.edu/ml/datasets/Heart+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.hungarian.data",
  "d1ad108f785768cd3d7e82dc522e6f5a61eea93cccfb3a46ee8076f73fc3d796",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=14,
        feature_indices=1:13,
        categoric_indices=[3; 7; 11; 13],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# This file is part of Kpax3. License is MIT.

function reorderunits(R::Vector{Int},
                      P::Matrix{Float64},
                      clusterorder::Vector{Int})
  n = length(R)
  k = maximum(R)

  M = zeros(Float64, n)
  v = zeros(Float64, n)
  for j in 1:(n - 1), i in (j + 1):n
    if R[i] == R[j]
      M[i] += P[i, j]
      v[i] += 1

      M[j] += P[i, j]
      v[j] += 1
    end
  end

  M ./= v

  neworder = zeros(Int, n)
  midpoint = zeros(Float64, k)
  seppoint = zeros(Float64, k + 1)

  h = 1
  u = 1
  for g in 1:k
    idx = findall(R .== clusterorder[g])
    u = length(idx)
    ord = sortperm(M[idx], rev=true)

    copyto!(neworder, h, idx[ord], 1, u)
    midpoint[g] = (2 * h + u - 1) / 2
    seppoint[g] = h - 0.5

    h += u
  end

  seppoint[k + 1] = n + 0.5

  (neworder, midpoint, seppoint)
end

function expandsquarediag(j::Int,
                          n::Int,
                          val::String,
                          z::Vector{String})
  # probability matrix is symmetric, so we only need to scan either the columns
  # or the rows. We scan the column because elements are sequential in
  # variable z
  h = j
  expand = true

  while expand && h < n
    idx = LinearIndices((n, n))[j, h + 1]

    i = j
    while expand && i <= h + 1
      if z[idx] != val
        expand = false
      end

      idx += 1
      i += 1
    end

    if expand
      h += 1
    end
  end

  h
end

function expandsquare(i::Int,
                      j::Int,
                      ni::Int,
                      nj::Int,
                      val::String,
                      z::Vector{String},
                      processed::BitArray{2})
  imax = i
  jmax = j
  expand = true

  while expand && imax < ni && jmax < nj
    # check next column from imin to imax + 1
    imin = i
    idx = LinearIndices((ni, nj))[imin, jmax + 1]

    while expand && imin <= imax + 1
      if processed[idx] || z[idx] != val
        expand = false
      end

      idx += 1
      imin += 1
    end

    # check next row from j to jmax + 1 only if column was ok
    if expand
      jmin = j
      while expand && jmin <= jmax + 1
        idx = LinearIndices((ni, nj))[imax + 1, jmin]

        if processed[idx] || z[idx] != val
          expand = false
        end

        jmin += 1
      end
    end

    if expand
      imax += 1
      jmax += 1
    end
  end

  (imax, jmax)
end

function expandrecthoriz(i::Int,
                         j::Int,
                         imax::Int,
                         jmax::Int,
                         ni::Int,
                         nj::Int,
                         val::String,
                         z::Vector{String},
                         processed::BitArray{2})
  expand = true

  while expand && jmax < nj
    idx = LinearIndices((ni, nj))[i, jmax + 1]

    imin = i
    while expand && imin <= imax
      if processed[idx] || z[idx] != val
        expand = false
      end

      idx += 1
      imin += 1
    end

    if expand
      jmax += 1
    end
  end

  jmax
end

function expandrectverti(i::Int,
                         j::Int,
                         imax::Int,
                         jmax::Int,
                         ni::Int,
                         nj::Int,
                         val::String,
                         z::Vector{String},
                         processed::BitArray{2})
  expand = true

  while expand && imax < ni
    jmin = j
    while expand && jmin <= jmax
      idx = LinearIndices((ni, nj))[imax + 1, jmin]

      if processed[idx] || z[idx] != val
        expand = false
      end

      jmin += 1
    end

    if expand
      imax += 1
    end
  end

  imax
end

function expandrect(i::Int,
                    j::Int,
                    ni::Int,
                    nj::Int,
                    val::String,
                    z::Vector{String},
                    processed::BitArray{2})
  (imax, jmax) = expandsquare(i, j, ni, nj, val, z, processed)

  h = expandrecthoriz(i, j, imax, jmax, ni, nj, val, z, processed)
  v = expandrectverti(i, j, imax, jmax, ni, nj, val, z, processed)

  if v - imax > h - jmax
    h = jmax
  elseif h - jmax > v -imax
    v = imax
  end

  (v, h)
end


text:
using JuMP, LightGraphs, GLPK, GLPKMathProgInterface

time_factor = 1.0

solver_instance = with_optimizer(GLPK.Optimizer, tm_lim=timeout*1000); time_param="tm_lim"; time_factor=1000.0
# using Gurobi; const GRB_ENV = Gurobi.Env(); solver_instance = with_optimizer(Gurobi.Optimizer, GRB_ENV, OutputFlag=0, LazyConstraints=1, TimeLimit=timeout); time_param = "TimeLimit"
# using CPLEX, MathProgBase; solver_instance = with_optimizer(CPLEX.Optimizer, CPX_PARAM_TILIM=timeout)

include("basic.jl")
include("intermediate.jl")
include("parallel.jl")
include("transport.jl")

"""
chain_cycle_match
Returns a maximum-weight matching in the form of three elements:
- An array of matched vertices
- An array of matched edges
- The corresponding match value.
"""
function chain_cycle_match(input_data::Input, objective_fn::String = "basic")
    if input_data.max_cycle_length == 0
        input_data.max_cycle_length = 1
    end
    while true
        t0 = time()
        @time begin
            if input_data.max_cycle_length <= 0 || input_data.max_chain_length < 0
                println("Error in parameters:")
                println("\tMaximum cycle length: $(input_data.max_cycle_length)")
                println("\tMaximum chain length: $(input_data.max_chain_length)")
                exit(1)
            end
            input_data.solver_instance = solver_instance
            input_data.time_param = time_param
            input_data.time_factor = time_factor
            input_data.start_time = t0

            println("Solving for parameters:")
            println("\tMaximum cycle length: $(input_data.max_cycle_length)")
            println("\tMaximum chain length: $(input_data.max_chain_length)")

            if objective_fn == "basic"
                problem = Basic(input_data)
            elseif objective_fn == "pctsp"
                problem = Intermediate(input_data)
            else
                problem = Parallel(input_data)
            end

            output_data = solve(problem)
        end

        if !isnothing(output_data) && output_data.optimal
            match_vertices = output_data.match_vertices
            match_edges = output_data.match_edges
            match_cycles = output_data.match_cycles
            graph_cycles = output_data.graph_cycles
            value = round(output_data.value, digits=5)

            println("\nObjective value: $value")
            return (match_vertices, match_edges, match_cycles, graph_cycles, value)
        end
        println()
        println("Solution is not optimal, decreasing maximum lengths")
        input_data.max_cycle_length = input_data.max_cycle_length-1
    end
end

"""
generate_graph
Creates a LightGraph directed graph from an adjacency edge-weight matrix.
"""
function generate_graph(weights::Array{Float64,2})
    @assert size(weights)[1] == size(weights)[2]
    n = size(weights)[1]
    graph = DiGraph(n)
    for i in 1:n
        for j in 1:n
            if weights[i,j] > 0
                add_edge!(graph, i, j)
            end
        end
    end
    return graph
end

function get_match_list(match_edges, match_cycles, graph_cycles, graph::DiGraph)
    n = nv(graph)
    cycle_chains_list = []
    flag = zeros(Int64, n)
    for i in 1:n
        if flag[i] == 0
            current_vertex = i
            cycle_chain_length = 0
            cycle_chain = []
            for k in 1:n  # replaces while true in case the match_edges is not a valid match.
                next_vertex = find_recipient(current_vertex, match_edges, graph)
                if next_vertex == 0 # no recipient, end of a chain
                    if cycle_chain_length > 0
                        push!(cycle_chain, current_vertex)
                        push!(cycle_chains_list, ("chain", cycle_chain))
                    end
                    break
                end

                if next_vertex == i  # closed the cycle
                    push!(cycle_chain, current_vertex)
                    push!(cycle_chains_list, ("cycle", cycle_chain))
                    break
                end
                push!(cycle_chain, current_vertex)
                flag[next_vertex] = 1
                cycle_chain_length += 1
                current_vertex = next_vertex
            end
        end
    end
    if !isnothing(match_cycles)
        num_cycles = size(graph_cycles)[1]
        for i in 1:num_cycles
            if match_cycles[i] > 0.5
                push!(cycle_chains_list, ("cycle", graph_cycles[i]))
            end
        end
    end
    return cycle_chains_list
end


text:
@doc """
Heavy tailed distribution used for proposals
""" ->
function randh()
	t = randn()/sqrt(-log(rand()))
	return 10.0^(1.5 - 3*abs(t))*randn()
end

@doc """
log(sum(exp(x)))
""" ->
function logsumexp(x::Array{Float64, 1})
	top = maximum(x)
	y = exp(x - top)
	return log(sum(y)) + top
end

@doc """
log(exp(a) - exp(b))
""" ->
function logdiffexp(a::Float64, b::Float64)
	@assert a > b
	return b + log(exp(a - b) - 1.0)
end


text:
module ImplicitDifferentiation

using ChainRulesCore
using Krylov
using LinearOperators
using SparseArrays

include("implicit_function.jl")
include("simplex.jl")

export ImplicitFunction
export simplex_projection

end


text:
module mholestr1
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, 1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  GaussRule(2, 2)),  material)

    K = stiffness(femm,  geom,  u)
    U=  K\(Fm)
    scattersysvec!(u, U[:])

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mholestr1-s1.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 2.210557410276065)/2.210557410276065 <= 1.0e-4
    # File =  "mholestr1-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr1
mholestr1.test()

module mholestr2
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, 1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStrain
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  GaussRule(2, 2)),  material)

    K = stiffness(femm,  geom,  u)
    U=  K\(Fm)
    scattersysvec!(u, U[:])

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mholestr2-s1.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 2.2102738887214257)/2.2102738887214257 <= 1.0e-4
    # File =  "mholestr2-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr2
mholestr2.test()

module mholestr3
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, -1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  GaussRule(2, 2)),  material)

    K = stiffness(femm,  geom,  u)
    U=  K\(Fm)
    scattersysvec!(u, U[:])

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    File =  "mholestr3-s1.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 5.921999943843146)/5.921999943843146 <= 1.0e-4
    # File =  "mholestr3-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr3
mholestr3.test()


text:
# get values for 10 best parameter sets of each run
# get number of saved values
using DelimitedFiles
using Plots
using StatsPlots
using CSV
using Plots.PlotMeasures
using DocStringExtensions
# pyplot()
# Plots.PyPlotBackend()

"""
This function retrieves the best parameter sets and saves them in a seperate csv file

$(SIGNATURES)

"""
function getbest_parametersets(path_to_file, number_best)
    calibration = readdlm(path_to_file, ',')
    # sort the calibration according to the euclidean distance
    calibration_sorted = sortslices(calibration, dims=1)#5:end,:]
    #number_best = 10
    calibration_best = calibration_sorted[1:number_best,:]
    print(size(calibration_best))
    writedlm(path_to_file[1:end-4]*"_best_"*string(number_best)*".csv", calibration_best, ',')
end

"""
Combines Calibration Results of several files into one file.

$(SIGNATURES)

The functions takes a path with data and saves the combined data to the path_to_save
"""
function combine_calibrations(path, path_to_save)
    all_calibrations = Array{Float64,2}[]
    total_saved = 0
    all_calibrations = zeros((1,30))
    #files = readdir(path)
    files = filter(name -> endswith(name, ".csv"), readdir(path))
    for i in 1: length(files)
        calibration = readdlm(path*files[i], ',')
        #adds the new array to former array
        all_calibrations = vcat(all_calibrations, calibration)
        total_saved+= size(calibration)[1]
    end

    #return all_calibrations[2:end, :], total_saved
    writedlm(path_to_save, all_calibrations, ',')
end

function calibration_statistics(path_to_file, number_best, lower_boundary_y_axis)
    max_Obj = Float64[]
    max_NSE = Float64[]
    max_NSElog = Float64[]
    max_VE = Float64[]
    max_NSE_FDC = Float64[]
    max_Reative_Error_AC_1day = Float64[]
    max_NSE_AC_90day = Float64[]
    max_Relative_Error_Runoff = Float64[]
    max_Snow_Cover = Float64[]
    end_file = findlast(isequal('t'), path_to_file)
    names_obj = ["NSE", "NSElog", "VE", "NSE_FDC", "Reative_Error_AC_1day", "NSE_AC_90day", "NSE_Runoff", "Snow_Cover"]
    Parameters = ["beta_Bare", "beta_Forest", "beta_Grass", "beta_Rip", "Ce", "Interceptioncapacity_Forest", "Interceptioncapacity_Grass", "Interceptioncapacity_Rip", "Kf_Rip", "Kf", "Ks", "Meltfactor", "Mm", "Ratio_Pref", "Ratio_Riparian", "Soilstoaragecapacity_Bare", "Soilstoaragecapacity_Forest", "Soilstoaragecapacity_Grass", "Soilstoaragecapacity_Rip", "Temp_Thresh"]
    #Objective_Functions = [max_NSE, max_NSElog, max_VE, max_NSE_FDC, max_Reative_Error_AC_1day, max_NSE_AC_90day, max_Relative_Error_Runoff, max_Snow_Cover]
    # get array with all calibtation data
    calibration = readdlm(path_to_file, ',')
    # sort the calibration according to the euclidean distance
    calibration_sorted = sortslices(calibration, dims=1)
    #number_best = 10
    calibration_best = calibration_sorted[1:number_best,:]
    ED_best = calibration_best[:,1]
    plots_obj = []
    for i in 1:8
        scatter(ED_best, calibration_best[:,i+1], xlabel = "Euclidean Distance", ylabel= names_obj[i])
        push!(plots_obj, scatter(ED_best, calibration_best[:,i+1], xlabel = "Euclidean Distance", ylabel= names_obj[i]))
    end
    plot(plots_obj[1], plots_obj[2], plots_obj[3], plots_obj[4], plots_obj[5], plots_obj[6], plots_obj[7], plots_obj[8], layout= (2,4), legend = false, size=(1400,800))
    #ylims!(lower_boundary_y_axis,1)
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_objective_functions.png")
    #plot the parameter distribution
    plots_par = []
    for i in 1:20
        scatter(ED_best, calibration_best[:,i+9], xlabel = "Euclidean Distance", ylabel= string(Parameters[i]))
        push!(plots_par, scatter(ED_best, calibration_best[:,i+9], xlabel = "Euclidean Distance", ylabel= Parameters[i]))
    end
    print(size(plots_par), typeof(plots_par))
    plot(plots_par[1], plots_par[2], plots_par[3], plots_par[4], plots_par[5], plots_par[6], plots_par[7], plots_par[8], plots_par[9], plots_par[10], plots_par[11], plots_par[12], layout= (3,4), legend=false, size=(1400,1000))
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_parameters1.png")
    plot(plots_par[13], plots_par[14], plots_par[15], plots_par[16], plots_par[17], plots_par[18], plots_par[19], plots_par[20], layout= (2,4), legend=false, size=(1400,1000))
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_parameters2.png")
    # scatter(ED_best, calibration_best[:,30], xlabel = "Euclidean Distance", ylabel= "loss parameter")
    # savefig(path_to_file[1:end_file+1]*string(number_best)*"_loss_parameters.png")
end

function calibration_statistics_pitztal(path_to_file, number_best, lower_boundary_y_axis)
    max_Obj = Float64[]
    max_NSE = Float64[]
    max_NSElog = Float64[]
    max_VE = Float64[]
    max_NSE_FDC = Float64[]
    max_Reative_Error_AC_1day = Float64[]
    max_NSE_AC_90day = Float64[]
    max_Relative_Error_Runoff = Float64[]
    max_Snow_Cover = Float64[]
    end_file = findlast(isequal('t'), path_to_file)
    names_obj = ["NSE", "NSElog", "VE", "NSE_FDC", "Reative_Error_AC_1day", "NSE_AC_90day", "NSE_Runoff", "Snow_Cover"]
    Parameters = ["beta_Bare", "beta_Forest", "beta_Grass", "beta_Rip", "Ce", "Interceptioncapacity_Forest", "Interceptioncapacity_Grass", "Interceptioncapacity_Rip", "Kf_Rip", "Kf", "Ks", "Meltfactor", "Mm", "Ratio_Pref", "Ratio_Riparian", "Soilstoaragecapacity_Bare", "Soilstoaragecapacity_Forest", "Soilstoaragecapacity_Grass", "Soilstoaragecapacity_Rip", "Temp_Thresh", "loss_term"]
    #Objective_Functions = [max_NSE, max_NSElog, max_VE, max_NSE_FDC, max_Reative_Error_AC_1day, max_NSE_AC_90day, max_Relative_Error_Runoff, max_Snow_Cover]
    # get array with all calibtation data
    calibration = readdlm(path_to_file, ',')
    # sort the calibration according to the euclidean distance
    calibration_sorted = sortslices(calibration, dims=1)
    #number_best = 10
    calibration_best = calibration_sorted[1:number_best,:]
    ED_best = calibration_best[:,1]
    plots_obj = []
    for i in 1:8
        scatter(ED_best, calibration_best[:,i+1], xlabel = "Euclidean Distance", ylabel= names_obj[i])
        push!(plots_obj, scatter(ED_best, calibration_best[:,i+1], xlabel = "Euclidean Distance", ylabel= names_obj[i]))
    end
    plot(plots_obj[1], plots_obj[2], plots_obj[3], plots_obj[4], plots_obj[5], plots_obj[6], plots_obj[7], plots_obj[8], layout= (2,4), legend = false, size=(1400,800))
    #ylims!(lower_boundary_y_axis,1)
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_objective_functions.png")
    #plot the parameter distribution
    plots_par = []
    for i in 1:21
        scatter(ED_best, calibration_best[:,i+9], xlabel = "Euclidean Distance", ylabel= string(Parameters[i]))
        push!(plots_par, scatter(ED_best, calibration_best[:,i+9], xlabel = "Euclidean Distance", ylabel= Parameters[i]))
    end
    print(size(plots_par), typeof(plots_par))
    plot(plots_par[1], plots_par[2], plots_par[3], plots_par[4], plots_par[5], plots_par[6], plots_par[7], plots_par[8], plots_par[9], plots_par[10], plots_par[11], plots_par[12], layout= (3,4), legend=false, size=(1400,1000))
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_parameters1.png")
    plot(plots_par[13], plots_par[14], plots_par[15], plots_par[16], plots_par[17], plots_par[18], plots_par[19], plots_par[20], layout= (2,4), legend=false, size=(1400,1000))
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_parameters2.png")
    plot(plots_par[21])
    savefig(path_to_file[1:end_file+1]*string(number_best)*"_loss_parameter.png")
    # scatter(ED_best, calibration_best[:,30], xlabel = "Euclidean Distance", ylabel= "loss parameter")
    # savefig(path_to_file[1:end_file+1]*string(number_best)*"_loss_parameters.png")
end

function projections_statistics(path_to_file)
    names_obj = ["NSE", "NSElog", "VE","NSE_FDC", "Reative_Error_AC_1day", "NSE_AC_90day", "NSE_Runoff", "Snow_Cover"]
    #Parameters = ["beta_Bare", "beta_Forest", "beta_Grass", "beta_Rip", "Ce", "Interceptioncapacity_Forest", "Interceptioncapacity_Grass", "Interceptioncapacity_Rip", "Kf_Rip", "Kf", "Ks", "Meltfactor", "Mm", "Ratio_Pref", "Ratio_Riparian", "Soilstoaragecapacity_Bare", "Soilstoaragecapacity_Forest", "Soilstoaragecapacity_Grass", "Soilstoaragecapacity_Rip", "Temp_Thresh"]
    #Objective_Functions = [max_NSE, max_NSElog, max_VE, max_NSE_FDC, max_Reative_Error_AC_1day, max_NSE_AC_90day, max_Relative_Error_Runoff, max_Snow_Cover]
    # get array with all calibtation data
    projections = readdlm(path_to_file, ',')
    number_best = size(projections)[1]
    runs = collect(1:number_best)
    print(size(projections))
    plots_obj = []
    for i in 1:8
        #scatter(ED_best, calibration_best[:,i+1], xlabel = "Euclidean Distance", ylabel= names_obj[i])
        # xlabel!("Euclidean Distance")
        # ylabel!(names_obj[i])
        #savefig(names_obj[i]*".png")
        push!(plots_obj, scatter(runs, projections[:,i+1], xlabel = "Runs", ylabel= names_obj[i]))
    end
    plot(plots_obj[1], plots_obj[2], plots_obj[3], plots_obj[4], plots_obj[5], plots_obj[6], plots_obj[7], plots_obj[8], layout= (2,4), legend = false, size=(1400,800))
    savefig("Gailtal/Projections/objbestfit_"*string(number_best)*".png")
end

function boxplot_projection(path_to_calibration, path, number_best, Catchment_Name)
    #path = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/"
    # 14 different projections
    Name_Projections = readdir(path)
    plots_obj = []
    # run the model for all projections using the best 100 parameter sets
    All_Obj_Functions = Array{Float64,2}[]
    for (i, name) in enumerate(Name_Projections)
            getData = readdlm(path*name*"/"*Catchment_Name*"/300_model_results_snow_redistr_1986_2005.csv",',')
            push!(All_Obj_Functions, getData[1:number_best,1:9])
    end

    Calibration = readdlm(path_to_calibration, ',')[:,1:9]

    names_obj = ["Euclidean Distance", "NSE", "NSElog", "VE","NSE_FDC", "Reative_Error_AC_1day", "NSE_AC_90day", "NSE_Runoff", "Snow_Cover"]
    #number = collect(1:size(All_Obj_Functions)[1])
    for obj in 1:size(names_obj)[1]
        plot()
        for i in 1:size(All_Obj_Functions)[1]
            box = boxplot!(["Proj " *string(i)],All_Obj_Functions[i][1:number_best,obj],leg = false)

        end
        box = boxplot!(["Calibration"],Calibration[:,obj],leg = false, color="darkgrey")
        ylabel!(names_obj[obj])
        push!(plots_obj, box)
        #savefig("/home/sarah/Master/Thesis/Results/Projektionen/10best/"*names_obj[obj]*"rcp45_comparison_Calibration_"*string(number_best)*"best.png")
    end

    plot(plots_obj[2], plots_obj[3], plots_obj[4], plots_obj[5], plots_obj[6], plots_obj[7], plots_obj[8], plots_obj[9], layout= (2,4), legend = false, size=(2000,1000), left_margin = [5mm 0mm], bottom_margin = 20px, xrotation = 60)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/Comparison_Real_Proj/rcp45_comparison_Calibration_"*string(number_best)*"best.png")

    plot(plots_obj[1], left_margin = [5mm 0mm], bottom_margin = 20px, xrotation = 60)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/Comparison_Real_Proj/rcp45_comparison_Calibration_ED_"*string(number_best)*"best.png")
end

function plot_FDC(path)
    Name_Projections = readdir(path)
    All_Discharges = Array{Float64,2}[]
    for (i, name) in enumerate(Name_Projections)
            getDischarge = readdlm(path*name*"/Gailtal/100_model_results_85_05_discharge.csv",',')
            push!(All_Discharges, getDischarge)
    end
    #observed discharge

    Discharge = CSV.read("/home/sarah/HBVModel/Gailtal/Q-Tagesmittel-212670.csv", header= false, skipto=23, decimal=',', delim = ';', types=[String, Float64])
    Discharge = convert(Matrix, Discharge)
    startindex = findfirst(isequal("01.10.1985 00:00:00"), Discharge)
    endindex = findfirst(isequal("30.09.2005 00:00:00"), Discharge)
    Observed_Discharge = Array{Float64,1}[]
    push!(Observed_Discharge, Discharge[startindex[1]:endindex[1],2])
    #Observed_Discharge = Observed_Discharge[1]
    observed_FDC = flowdurationcurve(log.(Observed_Discharge[1]))
    plot()

    #print(size(All_Discharges[1][1,:]))
    for proj in 1:length(Name_Projections)
        NSE_FDC_observations = Float64[]
        for i in 1:size(All_Discharges[1])[1]
            modeled_FDC = flowdurationcurve(log.(All_Discharges[proj][i,:]))
            #plot!(modeled_FDC[2], modeled_FDC[1], color="black", legend=false, size=(1400,800))
            NSE_FDC = nse(observed_FDC[1], modeled_FDC[1])
            append!(NSE_FDC_observations, NSE_FDC)
        end
        # plot!(observed_FDC[2], observed_FDC[1], color="red", size=(1400,800))
        # title!(Name_Projections[proj])
        # xlabel!("Exceedance Probability")
        # ylabel!("Discharge [m3/s]")
        # savefig("/home/sarah/Master/Thesis/Results/Projektionen/FDC_"*Name_Projections[proj]*".png")
        boxplot!(NSE_FDC_observations, leg=false)
    end
    Modelled_Discharge_Observations = readdlm("Gailtal/Calibration_8.05/Discharges_best100.csv", '\t')
    #print(size(Modelled_Discharge_Observations)[1])
    NSE_FDC_observations = Float64[]
    for i in 1:size(Modelled_Discharge_Observations)[2]
        modeled_FDC = flowdurationcurve(log.(Modelled_Discharge_Observations[:,i]))
        #plot!(modeled_FDC[2], modeled_FDC[1], color="black", legend=false, size=(1400,800))
        NSE_FDC = nse(observed_FDC[1], modeled_FDC[1])
        append!(NSE_FDC_observations, NSE_FDC)
    end
    boxplot!(["obs"],NSE_FDC_observations, leg=false)
    # plot!(observed_FDC[2], observed_FDC[1], color="red")
    # title!("Observed Data")
    # xlabel!("Exceedance Probability")
    # ylabel!("Discharge [m3/s]")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/FDC_Boxplot.png")
    plot()
    observed_FDC = flowdurationcurve((Observed_Discharge[1]))
    for proj in 1:length(Name_Projections)
        NSE_FDC_observations = Float64[]
        for i in 1:size(All_Discharges[1])[1]
            modeled_FDC = flowdurationcurve(All_Discharges[proj][i,:])
            #plot!(modeled_FDC[2], modeled_FDC[1], color="black", legend=false, size=(1400,800))
            NSE_FDC = lognse(observed_FDC[1], modeled_FDC[1])
            append!(NSE_FDC_observations, NSE_FDC)
        end
        # plot!(observed_FDC[2], observed_FDC[1], color="red", size=(1400,800))
        # title!(Name_Projections[proj])
        # xlabel!("Exceedance Probability")
        # ylabel!("Discharge [m3/s]")
        # savefig("/home/sarah/Master/Thesis/Results/Projektionen/FDC_"*Name_Projections[proj]*".png")
        boxplot!(NSE_FDC_observations, leg=false)
    end
    Modelled_Discharge_Observations = readdlm("Gailtal/Calibration_8.05/Discharges_best100.csv", '\t')
    #print(size(Modelled_Discharge_Observations)[1])
    NSE_FDC_observations = Float64[]
    for i in 1:size(Modelled_Discharge_Observations)[2]
        modeled_FDC = flowdurationcurve(Modelled_Discharge_Observations[:,i])
        #plot!(modeled_FDC[2], modeled_FDC[1], color="black", legend=false, size=(1400,800))
        NSE_FDC = lognse(observed_FDC[1], modeled_FDC[1])
        append!(NSE_FDC_observations, NSE_FDC)
    end
    boxplot!(["obs"],NSE_FDC_observations, leg=false)
    # plot!(observed_FDC[2], observed_FDC[1], color="red")
    # title!("Observed Data")
    # xlabel!("Exceedance Probability")
    # ylabel!("Discharge [m3/s]")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/FDC_old_Boxplot.png")

end

boxplot_projection("/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/Pitztal_Parameterfit_All_runs_snow_redistr_best_300.csv", "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/", 300, "Pitztal")
#Modelled_Discharge_Observations = plot_FDC("/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/")

# #----------------- COMBINE RESULTS OF ONE DEVICE-------------
#combine_calibrations("/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/", "/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/Pitztal_Parameterfit_All_runs_snow_redistr.csv")
# --------------- STORE BEST PARAMETER SETS ---------------------
#getbest_parametersets("/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/Pitztal_Parameterfit_All_runs_snow_redistr.csv", 300)
# -------------- GET STATISTICS -------

#calibration_statistics_pitztal("/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/Pitztal_Parameterfit_All_runs_snow_redistr_best_10000.csv", 10000, 0.7)

# Defreggen_parameters = readdlm("/home/sarah/Master/Thesis/Calibrations/Pitztal_Snow_Redistribution/Pitztal_Parameterfit_All_runs_snow_redistr_best_1000.csv", ',')
# #
# Params_unique = unique(DataFrame(Defreggen_parameters))
# test2 = convert(Vector, Params_unique[:,1])
# sort!(test2)
# # writedlm("/home/sarah/Master/Thesis/Calibrations/Paltental_less_dates/Paltental_Parameterfit_All_less_dates_unique2.csv", test2, ',')

function EC_calibration(path_to_file, lower_threshold, upper_threshold, nr_runs)
    calibration = readdlm(path_to_file, ',')
    # sort the calibration according to the euclidean distance
    calibration_sorted = sortslices(calibration, dims=1)
    #number_best = 10
    #calibration_best = calibration_sorted[1:number_best,:]
    number_values = Float64[]
    threshold_values = collect(lower_threshold:0.001:upper_threshold)
    print(size(threshold_values))
    for threshold in threshold_values
        all_values_below_threshold = length(findall(x -> x < threshold, calibration_sorted[:,1]))
        append!(number_values, all_values_below_threshold)
    end
    print(number_values/nr_runs * 100)
    scatter(threshold_values, number_values/nr_runs * 100 , size=(1400,800))
    xlabel!("Euclidean Distance")
    ylabel!("Percent of Runs below the Euclidean Distance")
    savefig(path_to_file[1:end-4]* "_compare_ED.png")
end

#EC_calibration("/home/sarah/Master/Thesis/Calibrations/Pitztal_loss_less_dates/Pitztal_Parameterfit_All_runs.csv", 0.09, 0.3, 3000002)

#projections_statistics("Gailtal/Projections/Gailtal_Parameterfit_best100_projection1.csv")

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
module TestLinqMacro

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta

srand(1)
n = 100
df = DataFrame(a = rand(1:3, n),
               b = ["a","b","c","d"][rand(1:4, n)],
               x = rand(n))

x = @where(df, :a .> 2)
x = @transform(x, y = 10 * :x)
x = @by(x, :b, meanX = mean(:x), meanY = mean(:y))
x = @orderby(x, :b, -:meanX)
x = @select(x, var = :b, :meanX, :meanY)

x1 = @linq transform(where(df, :a .> 2), y = 10 * :x)
x1 = @linq by(x1, :b, meanX = mean(:x), meanY = mean(:y))
x1 = @linq select(orderby(x1, :b, -:meanX), var = :b, :meanX, :meanY)

## chaining
xlinq = @linq df  |>
    where(:a .> 2)  |>
    transform(y = 10 * :x)  |>
    by(:b, meanX = mean(:x), meanY = mean(:y))  |>
    orderby(:b, -:meanX)  |>
    select(var = :b, :meanX, :meanY)

@test x == x1
@test x == xlinq

end


text:
# ROUTINES FOR COMPUTING STOCHASTIC GRADIENT OF THE OBJECTIVE FUNCTION

# get the partial derivatives of recourse variables wrt the first-stage variables for a given scenario
function getSensitivities(alpha::Array{Float64},pt::Array{Float64},omega_d::Array{Float64},tau_sat::Float64)

# store values of the derivative of the saturation function with respect to its argument
	dsat_dpt = zeros(Float64,numGen)
	for i = 1:numGen
		dsat_dpt[i] = smoothSaturationDeriv(pt[i],pmin[i],pmax[i],tau_sat)
	end
	den_s::Float64 = sum(alpha.*dsat_dpt)

# first, get partial derivatives wrt the p0s
	deriv_s_p0 = -dsat_dpt/den_s
	dp_dp0 = (dsat_dpt.*alpha)*(deriv_s_p0') + diagm(dsat_dpt)

	return dp_dp0
end

# compute a stochastic gradient of the objective function wrt the first-stage variables
function getStochasticGradient(p0::Array{Float64},rp::Array{Float64},rm::Array{Float64},alpha::Array{Float64},
								numSamples::Int64,tau_sat::Float64,gamma_res::Float64,tau_pos::Float64,
								gamma_line::Float64,delta_line::Float64)

# generate scenarios for computing the stochastic gradient
	omega_d, omega_w = generateScenarios(numSamples)

	stochgrad_p0 = zeros(Float64,numGen)
	stochgrad_rp = zeros(Float64,numGen)
	stochgrad_rm = zeros(Float64,numGen)
	stochgrad_alpha = zeros(Float64,numGen)

	for iter = 1:numSamples
	# first, solve the recourse system. this routine also updates the wind generator levels at the current scenario
		pt, p, slack, theta = solveRecourseProblem(p0,alpha,omega_d[:,iter],omega_w[:,iter],true,tau_sat)

	# then, compute partial derivatives of recourse variables wrt first-stage variables at above solution
	# ignore the ones we know to be zero
		dp_dp0 = getSensitivities(alpha,pt,omega_d,tau_sat)

	# compute stochastic gradient contribution from this scenario
	# add on the derivatives from the penalties for positive, negative reserves for generator i
		derivPosViol = zeros(Float64,numGen)
		derivNegViol = zeros(Float64,numGen)
		for i = 1:numGen
			posResViol::Float64 = p[i]-p0[i]-rp[i]
			derivPosViol[i] = reserveSmoothPenDeriv(posResViol,posResCost[i],gamma_res,tau_pos)

			negResViol::Float64 = p0[i]-p[i]-rm[i]
			derivNegViol[i] = reserveSmoothPenDeriv(negResViol,negResCost[i],gamma_res,tau_pos)
		end

		sumDemandDeviations::Float64 = sum(omega_d[j,iter] for j = 1:numLoads) - sum(meanDemand[j] for j = 1:numLoads)
		stochgrad_p0_iter1 = (dp_dp0')*(derivPosViol-derivNegViol)
		stochgrad_p0 += stochgrad_p0_iter1 + (derivNegViol - derivPosViol)
		stochgrad_alpha += stochgrad_p0_iter1*(sumDemandDeviations + slack)
		stochgrad_rp -= derivPosViol
		stochgrad_rm -= derivNegViol

	# add on the derivatives from the line limit penalties for line k
		derivLineViol = zeros(Float64,numLines)
		for k = 1:numLines
			lineFlow::Float64 = beta[k]*(theta[lines[k][1]] - theta[lines[k][2]])
			derivLineViol[k] = linePenDeriv(lineFlow,flow_max[k],gamma_line,delta_line)
		end

		stochgrad_p0_iter2 = (dp_dp0')*(Line_mat_times_B_over_G_trans*(derivLineViol.*beta))
		stochgrad_p0 += stochgrad_p0_iter2
		stochgrad_alpha += stochgrad_p0_iter2*(sumDemandDeviations + slack)
	end

	stochgrad_p0 /= numSamples
	stochgrad_rp /= numSamples
	stochgrad_rm /= numSamples
	stochgrad_alpha /= numSamples

# add the derivatives from the first-stage generation costs
	stochgrad_p0 += genCost
	stochgrad_rp += posResCost
	stochgrad_rm += negResCost

	return stochgrad_p0, stochgrad_rp, stochgrad_rm, stochgrad_alpha
end

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c51f94fc-4b81-11eb-1bb4-e988fea88167
md"## _If Statements_
"

# ╔═╡ 429023cc-4b82-11eb-06f8-cd8b628a57fb

# ╔═╡ 89797bda-4b82-11eb-192b-a179248f700e
num1 =  `$(@bind num1 html"<input type=text >")`

# ╔═╡ 6d9f5fb0-4b82-11eb-3949-cb21832e9734
md"Enter the first number: $num1 "

# ╔═╡ b7598e14-4b82-11eb-1bfc-9f9c12c10464

# ╔═╡ bf81ea96-4b82-11eb-20f7-c9c2a151b433

# ╔═╡ c5c10e3c-4b82-11eb-27cb-7168044bbdc8
begin
	if num1 <= num2
		lowest = num1
	else
		lowest = num2
	end
	md"The lowest number is $lowest"
end

# ╔═╡ Cell order:
# ╠═c51f94fc-4b81-11eb-1bb4-e988fea88167
# ╟─429023cc-4b82-11eb-06f8-cd8b628a57fb
# ╠═6d9f5fb0-4b82-11eb-3949-cb21832e9734
# ╠═89797bda-4b82-11eb-192b-a179248f700e
# ╟─b7598e14-4b82-11eb-1bfc-9f9c12c10464
# ╟─bf81ea96-4b82-11eb-20f7-c9c2a151b433
# ╠═c5c10e3c-4b82-11eb-27cb-7168044bbdc8


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Abstract base type for time integration schemes of storage class `3S*`
abstract type SimpleAlgorithm3Sstar end

"""
    HypDiffN3Erk3Sstar52()

Five stage, second-order acurate explicit Runge-Kutta scheme with stability region optimized for
the hyperbolic diffusion equation with LLF flux and polynomials of degree polydeg=3.
"""
struct HypDiffN3Erk3Sstar52 <: SimpleAlgorithm3Sstar
  gamma1::SVector{5, Float64}
  gamma2::SVector{5, Float64}
  gamma3::SVector{5, Float64}
  beta::SVector{5, Float64}
  delta::SVector{5, Float64}
  c::SVector{5, Float64}

  function HypDiffN3Erk3Sstar52()
    gamma1 = SVector(0.0000000000000000E+00, 5.2656474556752575E-01,  1.0385212774098265E+00, 3.6859755007388034E-01, -6.3350615190506088E-01)
    gamma2 = SVector(1.0000000000000000E+00, 4.1892580153419307E-01, -2.7595818152587825E-02, 9.1271323651988631E-02,  6.8495995159465062E-01)
    gamma3 = SVector(0.0000000000000000E+00, 0.0000000000000000E+00,  0.0000000000000000E+00, 4.1301005663300466E-01, -5.4537881202277507E-03)
    beta   = SVector(4.5158640252832094E-01, 7.5974836561844006E-01,  3.7561630338850771E-01, 2.9356700007428856E-02,  2.5205285143494666E-01)
    delta  = SVector(1.0000000000000000E+00, 1.3011720142005145E-01,  2.6579275844515687E-01, 9.9687218193685878E-01,  0.0000000000000000E+00)
    c      = SVector(0.0000000000000000E+00, 4.5158640252832094E-01,  1.0221535725056414E+00, 1.4280257701954349E+00,  7.1581334196229851E-01)

    new(gamma1, gamma2, gamma3, beta, delta, c)
  end
end

"""
    ParsaniKetchesonDeconinck3Sstar94()

Parsani, Ketcheson, Deconinck (2013)
  Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems
[DOI: 10.1137/120885899](https://doi.org/10.1137/120885899)
"""
struct ParsaniKetchesonDeconinck3Sstar94 <: SimpleAlgorithm3Sstar
  gamma1::SVector{9, Float64}
  gamma2::SVector{9, Float64}
  gamma3::SVector{9, Float64}
  beta::SVector{9, Float64}
  delta::SVector{9, Float64}
  c::SVector{9, Float64}

  function ParsaniKetchesonDeconinck3Sstar94()
    gamma1 = SVector(0.0000000000000000E+00, -4.6556413837561301E+00, -7.7202649689034453E-01, -4.0244202720632174E+00, -2.1296873883702272E-02, -2.4350219407769953E+00, 1.9856336960249132E-02, -2.8107894116913812E-01, 1.6894354373677900E-01)
    gamma2 = SVector(1.0000000000000000E+00, 2.4992627683300688E+00, 5.8668202764174726E-01, 1.2051419816240785E+00, 3.4747937498564541E-01, 1.3213458736302766E+00, 3.1196363453264964E-01, 4.3514189245414447E-01, 2.3596980658341213E-01)
    gamma3 = SVector(0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 7.6209857891449362E-01, -1.9811817832965520E-01, -6.2289587091629484E-01, -3.7522475499063573E-01, -3.3554373281046146E-01, -4.5609629702116454E-02)
    beta   = SVector(2.8363432481011769E-01, 9.7364980747486463E-01, 3.3823592364196498E-01, -3.5849518935750763E-01, -4.1139587569859462E-03, 1.4279689871485013E+00, 1.8084680519536503E-02, 1.6057708856060501E-01, 2.9522267863254809E-01)
    delta  = SVector(1.0000000000000000E+00, 1.2629238731608268E+00, 7.5749675232391733E-01, 5.1635907196195419E-01, -2.7463346616574083E-02, -4.3826743572318672E-01, 1.2735870231839268E+00, -6.2947382217730230E-01, 0.0000000000000000E+00)
    c      = SVector(0.0000000000000000E+00, 2.8363432481011769E-01, 5.4840742446661772E-01, 3.6872298094969475E-01, -6.8061183026103156E-01, 3.5185265855105619E-01, 1.6659419385562171E+00, 9.7152778807463247E-01, 9.0515694340066954E-01)

    new(gamma1, gamma2, gamma3, beta, delta, c)
  end
end

"""
    ParsaniKetchesonDeconinck3Sstar32()

Parsani, Ketcheson, Deconinck (2013)
  Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems
[DOI: 10.1137/120885899](https://doi.org/10.1137/120885899)
"""
struct ParsaniKetchesonDeconinck3Sstar32 <: SimpleAlgorithm3Sstar
  gamma1::SVector{3, Float64}
  gamma2::SVector{3, Float64}
  gamma3::SVector{3, Float64}
  beta::SVector{3, Float64}
  delta::SVector{3, Float64}
  c::SVector{3, Float64}

  function ParsaniKetchesonDeconinck3Sstar32()
    gamma1 = SVector(0.0000000000000000E+00, -1.2664395576322218E-01,  1.1426980685848858E+00)
    gamma2 = SVector(1.0000000000000000E+00,  6.5427782599406470E-01, -8.2869287683723744E-02)
    gamma3 = SVector(0.0000000000000000E+00,  0.0000000000000000E+00,  0.0000000000000000E+00)
    beta   = SVector(7.2366074728360086E-01,  3.4217876502651023E-01,  3.6640216242653251E-01)
    delta  = SVector(1.0000000000000000E+00,  7.2196567116037724E-01,  0.0000000000000000E+00)
    c      = SVector(0.0000000000000000E+00,  7.2366074728360086E-01,  5.9236433182015646E-01)

    new(gamma1, gamma2, gamma3, beta, delta, c)
  end
end

mutable struct SimpleIntegrator3SstarOptions{Callback}
  callback::Callback # callbacks; used in Trixi
  adaptive::Bool # whether the algorithm is adaptive; ignored
  dtmax::Float64 # ignored
  maxiters::Int # maximal numer of time steps
  tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored
end

function SimpleIntegrator3SstarOptions(callback, tspan; maxiters=typemax(Int), kwargs...)
  SimpleIntegrator3SstarOptions{typeof(callback)}(
    callback, false, Inf, maxiters, [last(tspan)])
end

mutable struct SimpleIntegrator3Sstar{RealT<:Real, uType, Params, Sol, Alg, SimpleIntegrator3SstarOptions}
  u::uType #
  du::uType
  u_tmp1::uType
  u_tmp2::uType
  t::RealT
  dt::RealT # current time step
  dtcache::RealT # ignored
  iter::Int # current number of time step (iteration)
  p::Params # will be the semidiscretization from Trixi
  sol::Sol # faked
  alg::Alg
  opts::SimpleIntegrator3SstarOptions
  finalstep::Bool # added for convenience
end

# Forward integrator.destats.naccept to integrator.iter (see GitHub PR#771)
function Base.getproperty(integrator::SimpleIntegrator3Sstar, field::Symbol)
  if field === :destats
    return (naccept = getfield(integrator, :iter),)
  end
  # general fallback
  return getfield(integrator, field)
end

# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
function solve(ode::ODEProblem, alg::T;
               dt, callback=nothing, kwargs...) where {T<:SimpleAlgorithm3Sstar}
  u = copy(ode.u0)
  du = similar(u)
  u_tmp1 = similar(u)
  u_tmp2 = similar(u)
  t = first(ode.tspan)
  iter = 0
  integrator = SimpleIntegrator3Sstar(u, du, u_tmp1, u_tmp2, t, dt, zero(dt), iter, ode.p,
                  (prob=ode,), alg,
                  SimpleIntegrator3SstarOptions(callback, ode.tspan; kwargs...), false)

  # initialize callbacks
  if callback isa CallbackSet
    for cb in callback.continuous_callbacks
      error("unsupported")
    end
    for cb in callback.discrete_callbacks
      cb.initialize(cb, integrator.u, integrator.t, integrator)
    end
  elseif !isnothing(callback)
    error("unsupported")
  end

  solve!(integrator)
end

function solve!(integrator::SimpleIntegrator3Sstar)
  @unpack prob = integrator.sol
  @unpack alg = integrator
  t_end = last(prob.tspan)
  callbacks = integrator.opts.callback

  integrator.finalstep = false
  @trixi_timeit timer() "main loop" while !integrator.finalstep
    if isnan(integrator.dt)
      error("time step size `dt` is NaN")
    end

    # if the next iteration would push the simulation beyond the end time, set dt accordingly
    if integrator.t + integrator.dt > t_end || isapprox(integrator.t + integrator.dt, t_end)
      integrator.dt = t_end - integrator.t
      terminate!(integrator)
    end

    # one time step
    integrator.u_tmp1 .= zero(eltype(integrator.u_tmp1))
    integrator.u_tmp2 .= integrator.u
    for stage in eachindex(alg.c)
      t_stage = integrator.t + integrator.dt * alg.c[stage]
      prob.f(integrator.du, integrator.u, prob.p, t_stage)

      delta_stage   = alg.delta[stage]
      gamma1_stage  = alg.gamma1[stage]
      gamma2_stage  = alg.gamma2[stage]
      gamma3_stage  = alg.gamma3[stage]
      beta_stage_dt = alg.beta[stage] * integrator.dt
      @trixi_timeit timer() "Runge-Kutta step" begin
        @threaded for i in eachindex(integrator.u)
          integrator.u_tmp1[i] += delta_stage * integrator.u[i]
          integrator.u[i]       = (gamma1_stage * integrator.u[i] +
                                   gamma2_stage * integrator.u_tmp1[i] +
                                   gamma3_stage * integrator.u_tmp2[i] +
                                   beta_stage_dt * integrator.du[i])
        end
      end
    end
    integrator.iter += 1
    integrator.t += integrator.dt

    # handle callbacks
    if callbacks isa CallbackSet
      for cb in callbacks.discrete_callbacks
        if cb.condition(integrator.u, integrator.t, integrator)
          cb.affect!(integrator)
        end
      end
    end

    # respect maximum number of iterations
    if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep
      @warn "Interrupted. Larger maxiters is needed."
      terminate!(integrator)
    end
  end

  return TimeIntegratorSolution((first(prob.tspan), integrator.t),
                                (prob.u0, integrator.u),
                                integrator.sol.prob)
end

# get a cache where the RHS can be stored
get_du(integrator::SimpleIntegrator3Sstar) = integrator.du
get_tmp_cache(integrator::SimpleIntegrator3Sstar) = (integrator.u_tmp1, integrator.u_tmp2)

# some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
u_modified!(integrator::SimpleIntegrator3Sstar, ::Bool) = false

# used by adaptive timestepping algorithms in DiffEq
function set_proposed_dt!(integrator::SimpleIntegrator3Sstar, dt)
  integrator.dt = dt
end

# stop the time integration
function terminate!(integrator::SimpleIntegrator3Sstar)
  integrator.finalstep = true
  empty!(integrator.opts.tstops)
end

# used for AMR
function Base.resize!(integrator::SimpleIntegrator3Sstar, new_size)
  resize!(integrator.u, new_size)
  resize!(integrator.du, new_size)
  resize!(integrator.u_tmp1, new_size)
  resize!(integrator.u_tmp2, new_size)
end

end # @muladd


text:
function value(obj::RealFieldUnivariateFunction, arg0::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement,), arg0)
end


text:
using DataStructures

function ancestors(str1, str2) :: Vector{Pair{String, String}}
	f = open("7.txt")
	res = []
	for i=1:594
		l = readline(f)
		words = split(l, " ")
		for j=3:length(words)-1
			if words[j] == str1 && words[j+1] == str2
				push!(res, Pair(words[1], words[2]))
				#println(words[1], words[2])
				break
			end
		end
	end
	res
end

#returns solution of part 1
function x()
	q = Queue{Pair{String,String}}()
	s = Set{Pair{String, String}}()
	v = ancestors("shiny", "gold")
	for e in v
		enqueue!(q,e)
	end
	n = 0
	while !isempty(q)
		#println(q)
		y = dequeue!(q)
		if(in(y, s))
			continue
		end
		push!(s, y)
		n += 1
		v = ancestors(y...)
		for e in v
			enqueue!(q,e)
		end
	end
	n
end

function descendants(str1, str2) :: Vector{Pair{Pair{String, String}, Int32}}
	f = open("7.txt")
	res = []
	for i=1:594
		l = readline(f)
		words = split(l, " ")
		if words[1] == str1 && words[2] == str2
			for j = 3:length(words)-2
				n = 0
				try
					n = parse(Int, words[j])
				catch
					continue
				end
				push!(res, Pair(Pair(words[j+1], words[j+2]), n))
			end
			break
		end
	end
	res
end

#returns solution for part 2
function y()
	q = Queue{Pair{Pair{String, String}, Int32}}()
	v = descendants("shiny", "gold")
	for e in v
		enqueue!(q,e)
	end
	n = 0
	while !isempty(q)
		#println(q)
		y = dequeue!(q)
		n += y.second
		v = descendants(y.first...)
		for e in v
			enqueue!(q, Pair(e.first, e.second * y.second))
		end
	end
	n
end


text:
# f(x)=hcat(x...)
Fhcat(x...)=(hcat(x...),nothing)
Fhcat_inplace(value,auxvalue,x...)=copy!(value,hcat(x...))

function Dhcat(derivativeIDX,f_c,faux_c,grad_c,grad_n,x...)
    startind=1
    for i=1:length(x)
        endind=startind+length(x[i])-1
        if pointer(x[i])==pointer(x[derivativeIDX]) # use pointers so we can deal with hcat([A B A])
            axpy!(1.0,grad_c[startind:endind],grad_n)
        end
        startind=endind+1
    end
end

if PROC=="GPU"

    function Fhcat_inplace(value::CudaArray,auxvalue,x::CudaArray...) # inplace
        totallength=1
        for i in 1:length(x)
            copyinto!(value,x[i],totallength)
            totallength+=length(x[i])
        end
    end

    function Dhcat(derivativeIDX,f_c,faux_c,grad_c,grad_n::CudaArray,x...)
        startind=1
        for i=1:length(x)
            endind=startind+length(x[i])-1
            if pointer(x[i])==pointer(x[derivativeIDX]) # use pointers so we can deal with vcat([A B A])
                copyfrom_update!(grad_n,grad_c,startind)
            end
            startind=endind+1
        end
    end

end

Derivative[Fhcat]=Dhcat # Define dictionary lookup
Inplace[Fhcat]=Fhcat_inplace

import Base.hcat
hcat(n::ADnode)=ADnode(Fhcat,n)

hcat(n::Array{ADnode,1})=n
hcat(n::ArrayADnode)=ADnode(Fhcat,n)

export hcat

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
Analytic axisymmetric medium-size tokamak equilibrium in cartesian coordinates.
"""
module TokamakMediumCartesian

    import ElectromagneticFields.AxisymmetricTokamakCartesian

    export initial_conditions_barely_passing, initial_conditions_barely_trapped,
           initial_conditions_deeply_passing, initial_conditions_deeply_trapped

    export hamiltonian, toroidal_momentum

    AxisymmetricTokamakCartesian.@code(2., 5., 2.) # inject magnetic field code

    initial_conditions_barely_passing() = ([2.5, 0., 0., 3.425E-1], (μ = 1E-2,)) # Δt=2.5, nt=50
    initial_conditions_barely_trapped() = ([2.5, 0., 0., 3.375E-1], (μ = 1E-2,)) # Δt=3.0, nt=100
    initial_conditions_deeply_passing() = ([2.5, 0., 0., 5E-1], (μ = 1E-2,))     # Δt=2.5, nt=25
    initial_conditions_deeply_trapped() = ([2.5, 0., 0., 1E-1], (μ = 1E-2,))     # Δt=5.0, nt=50

    μ_loop() = 1E-3
    μ_surface() = 1E-3

    function f_loop(s)
        X0 = 1.75
        Y0 = 0.0
        Z0 = 0.0
        u0 = 0.5
        rx = 0.1
        ry = 0.1

        Xs = X0 + rx*cos(2π*s)
        Zs = Z0 + ry*sin(2π*s)

        qs = [Xs, Y0, Zs, u0]

        return qs
    end

    function f_surface(s,t)
        X0 = 1.75
        Y0 = 0.0
        Z0 = 0.0
        u0 = 0.5
        r0 = 0.1

        Xt = X0 + r0*(s-0.5)
        Zt = Z0 + r0*(t-0.5)

        qt = [Xt, Y0, Zt, u0]

        return qt
    end

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")
    include("guiding_center_4d_loop.jl")
    include("guiding_center_4d_surface.jl")
    include("guiding_center_4d_diagnostics.jl")

end


text:
struct Lesson
    prompt::String
    request::String
    parser::Evaluator
    hint::String
end

"""
    Lesson(prompt::String, request::String, parser::Evaluator)

Will create a Lesson instance without a hint.
Note: it is strongly preferred that Lesson's contain `hints` so please construct via
Lesson( prompt::String, request::String, parser::Evaluator, hint::String ) except when neccesary.
"""
Lesson(prompt::String, request::String, parser::Evaluator) = Lesson(prompt, request, parser, "")

mutable struct Tutor
    lesson_path::String
    current_lesson::Int16
    lesson_plan::Vector{Lesson}
end

"""
    display_prompt_and_request(t::Tutor)

Displays a lesson prompt and requests an action from the user.
"""
function display_prompt_and_request(t::Tutor)
    active_lesson = t.lesson_plan[t.current_lesson]
    inform(active_lesson)
    request(active_lesson)
end

"""
    display_hint(t::Tutor)

Displays a hint pertaining to the users current lesson, if it exists. See `hint(lesson::Lesson)`
"""
function display_hint(t::Tutor)
    active_lesson = t.lesson_plan[t.current_lesson]
    hint(active_lesson)
end

"""
    keywords(t, user_input)

Assesses if the user has asked for help, hint, or some other invention specific to JuliaTutor.
"""
function keywords(t, user_input)
    used_keyword = false
    #ToDo:: strip whitespace, add more keywords?
    if (user_input == "help") || (user_input == "help()")
        help()
        used_keyword = true
    elseif (user_input == "hint") || (user_input == "hint()")
        display_hint( t )
        used_keyword = true
    end
    return used_keyword
end

"""
    ( t::Tutor )( user_input::String )

This is the parser, and lesson plan orchestration mechanism of Julia tutor. This functor progresses the user
through lessons and parsers their inputs.
"""
function ( t::Tutor )( user_input::String )
    active_lesson = t.lesson_plan[ t.current_lesson ]
    if active_lesson.parser( user_input )
        if t.current_lesson < length(t.lesson_plan)
            t.current_lesson += 1
            display_prompt_and_request( t )
        else
            println( Crayon(foreground = :green, bold = true),
                    "\n Congratulations! " *
                    "\n You've completed this lesson. Care to try another? \n")
            print(Crayon( foreground = :white, bold = false), "" )
            menu()
        end
    end
end