text:
"""
    flux_balance_analysis_vec(args...)::Maybe{Vector{Float64}}

A variant of FBA that returns a vector of fluxes in the same order as reactions
of the model, if the solution is found.

Arguments are passed to [`flux_balance_analysis`](@ref).

This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
flux_balance_analysis_vec(args...; kwargs...)::Maybe{Vector{Float64}} =
    flux_vector(flux_balance_analysis(args...; kwargs...))

"""
    flux_balance_analysis_dict(model::MetabolicModel, args...)::Maybe{Dict{String, Float64}}

A variant of FBA that returns a dictionary assigning fluxes to reactions, if
the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref).

This function is kept for backwards compatibility, use [`flux_dict`](@ref)
instead.
"""
flux_balance_analysis_dict(
    model::MetabolicModel,
    args...;
    kwargs...,
)::Maybe{Dict{String,Float64}} =
    flux_dict(model, flux_balance_analysis(model, args...; kwargs...))

"""
    flux_balance_analysis(
        model::M,
        optimizer;
        modifications = [],
    ) where {M<:MetabolicModel}

Run flux balance analysis (FBA) on the `model` optionally specifying
`modifications` to the problem.  Basically, FBA solves this optimization problem:
```
max cᵀx
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat
Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more
information.

The `optimizer` must be set to a `JuMP`-compatible optimizer, such as
`GLPK.Optimizer` or `Tulip.Optimizer`

Optionally, you may specify one or more modifications to be applied to the
model before the analysis, such as [`change_optimizer_attribute`](@ref),
[`change_objective`](@ref), and [`change_sense`](@ref).

Returns an optimized `JuMP` model.

# Example
```
model = load_model("e_coli_core.json")
solution = flux_balance_analysis(model, GLPK.optimizer)
value.(solution[:x])  # extract flux steady state from the optimizer

biomass_reaction_id = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")

modified_solution = flux_balance_analysis(model, GLPK.optimizer;
    modifications=[change_objective(biomass_reaction_id)])
```

"""
function flux_balance_analysis(
    model::M,
    optimizer;
    modifications = [],
) where {M<:MetabolicModel}

    opt_model = make_optimization_model(model, optimizer)

    for mod in modifications
        mod(model, opt_model)
    end

    optimize!(opt_model)
    return opt_model
end


text:
using Test
using ClassicAlgorithmsCollections

@testset "BIT Manipulation" begin
    # Tests for BIT Manipulation
    arr = [8, 1, 2, 12, 45, 100, 43]
    @test ClassicAlgorithmsCollections.maximum_xor_subarray(arr) == 110
    @test ClassicAlgorithmsCollections.sum_of_bit_difference(arr) == 104
    @test ClassicAlgorithmsCollections.magic_number(5) == 130
    @test ClassicAlgorithmsCollections.swapping_even_odd_bits(23) == 43
    arr_three_one = [12, 1, 12, 3, 12, 1, 1, 2, 3, 3]
    @test ClassicAlgorithmsCollections.get_single_integer(arr_three_one) == 2
    @test ClassicAlgorithmsCollections.int2binary(131) == ([1, 0, 0, 0, 0, 0, 1, 1], 1)
    @test ClassicAlgorithmsCollections.count_total_bits(17) == 35
end

@testset "Dynamic Programming" begin
    # Tests for BIT Manipulation
    arr = [8, 1, 2, 12, 45, 100, 43]
    @test ClassicAlgorithmsCollections.lcs_length("AGGTAB", "GXTXAYB") == 4
    arr = [10, 22, 9, 33, 21, 50, 41, 60]
    ClassicAlgorithmsCollections.lis_length(arr) == 5
    str1 = "sunday"
    str2 = "saturday"
    ClassicAlgorithmsCollections.minumum_operations(str1, str2) == 3
    ClassicAlgorithmsCollections.number_of_steps(8, 3) == 81
    set = [3, 34, 4, 12, 5, 2]
    sum_true = 9
    sum_false = 19
    ClassicAlgorithmsCollections.subset_sum_test(set, sum_true) == true
    ClassicAlgorithmsCollections.subset_sum_test(set, sum_false) == false
    val = [60, 100, 120, 230, 220]
    wt = [10, 20, 25, 30, 25]
    W_1 = 50
    W_2 = 10
    ClassicAlgorithmsCollections.knapsack_solver(val, wt, W_1) == 340
    ClassicAlgorithmsCollections.knapsack_solver(val, wt, W_2) == 60
end

@testset "Graph Theory" begin
    # Tests for Graph Theory
    graph = Dict(1 => [2, 3], 2 => [3], 3 => [1, 4], 4 => [4])
    @test ClassicAlgorithmsCollections.breadth_first_search(graph, 3) == [3, 1, 4, 2]
    @test ClassicAlgorithmsCollections.depth_first_search(graph, 3) == [3, 1, 2, 4]

    graph_with_weights_1 =
        Dict(1 => [(3, -2)], 2 => [(1, 4), (3, 3)], 3 => [(4, 2)], 4 => [(2, -1)])
    @test ClassicAlgorithmsCollections.shortest_path_tree(graph_with_weights_1, 2, 4) ==
          ([0 -1 -2 0; 4 0 2 4; 5 1 0 2; 3 -1 1 0], [2, 1, 3, 4])

    graph_with_weights_2 =
        Dict(1 => [(3, -2)], 2 => [(1, 4), (3, 3)], 3 => [(4, 2)], 4 => [(2, -1), (3, -1)])

    @test ClassicAlgorithmsCollections.shortest_path_tree(graph_with_weights_2) ==
          [0 -1 -2 0; 4 0 2 4; 5 1 0 2; 3 -1 -1 0]

    graph_cycle_true = Dict(1 => [2], 2 => [3], 3 => [1, 4])
    @test ClassicAlgorithmsCollections.graph_cycle_check(graph_cycle_true) == true

    graph_cycle_false = Dict(1 => [2], 2 => [5], 3 => [1, 4])
    @test ClassicAlgorithmsCollections.graph_cycle_check(graph_cycle_false) == false

    graph_with_spanning_tree =
        Dict(1 => [(2, 10), (3, 6), (4, 5)], 2 => [(4, 15)], 3 => [(4, 4)])
    @test ClassicAlgorithmsCollections.minimum_spanning_tree(graph_with_spanning_tree) ==
          [3 4 4; 1 4 5; 1 2 10]

    graph_bridge = Dict(1 => [2, 3, 4], 2 => [1, 3], 3 => [1, 2], 4 => [1, 5], 5 => [4])
    @test ClassicAlgorithmsCollections.graph_bridge_check(graph_bridge) == [4 5; 1 4]

    word_list = ["GEEKS", "FOR", "QUIZ", "GO"]
    graph_boogle = Dict(
        1 => ["G", "I", "Z"],
        2 => ["U", "E", "K"],
        3 => ["Q", "S", "E"],
        4 => ["D", "O", "P"],
        5 => ["F", "O", "R"],
    )

    @test ClassicAlgorithmsCollections.boogle_word_check(graph_boogle, word_list) ==
          ["GEEKS", "QUIZ", "FOR"]
end

@testset "Number Theory" begin
    # Tests for Number Theory
    @test ClassicAlgorithmsCollections.modular_exponentiation(2, 3, 15) == 8
    @test ClassicAlgorithmsCollections.modular_inverse(3, 11) == 4
    @test ClassicAlgorithmsCollections.modular_inverse(3, 1) == 1
    @test ClassicAlgorithmsCollections.euler_totient(20) ==
          [1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8]
    @test ClassicAlgorithmsCollections.sieve_of_eratosthenes(50) ==
          [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    num = [3, 4, 5]
    rem = [2, 3, 1]
    @test ClassicAlgorithmsCollections.chinese_remainder_theorem(num, rem) == 11
end

@testset "Sorting & Searching" begin
    # Test for sorting and searching
    array_sorting = [64, 34, 25, 12, 22, 11, 90]
    ref_sorting = [11, 12, 22, 25, 34, 64, 90]
    @test ClassicAlgorithmsCollections.bubble_sorting(array_sorting) == ref_sorting
    @test ClassicAlgorithmsCollections.insertion_sorting(array_sorting) == ref_sorting
    @test ClassicAlgorithmsCollections.merge_sorting(array_sorting) == ref_sorting
    @test ClassicAlgorithmsCollections.heap_sorting(array_sorting) == ref_sorting
    @test ClassicAlgorithmsCollections.quick_sorting(array_sorting) == ref_sorting

    array_searching = [10, 11, 12, 14, 16, 18, 19, 20, 21, 22, 23, 24, 33, 35, 42, 47]
    target_searching = 12
    ref_searching = 3
    @test ClassicAlgorithmsCollections.binary_search(array_searching, target_searching) ==
          ref_searching
    @test ClassicAlgorithmsCollections.interpolation_search(
        array_searching,
        target_searching,
    ) == ref_searching

    array_searching_pair = [10, 22, 28, 29, 30, 40]
    target_searching_pair = 56
    ref_searching_pair = (28, 29)
    @test ClassicAlgorithmsCollections.closest_pair_search(
        array_searching_pair,
        target_searching_pair,
    ) == ref_searching_pair

    array_unsorted_searching = [2, 5, 4, 7, 2, 8, 9, 3, 10, 2]
    target_unsorted_searching = 3
    ref_unsorted_searching = 8
    @test ClassicAlgorithmsCollections.binary_pivot_search(
        array_unsorted_searching,
        target_unsorted_searching,
    ) == ref_unsorted_searching
end

@testset "String & Array" begin
    old_string = "a!!!b.c79.d,e'f,ghi3###"
    reversed_string = "i!!!h.g79.f,e'd,cba3###"
    @test ClassicAlgorithmsCollections.reverse_pure_string(old_string) == reversed_string

    array_before_zigzag = [4, 3, 7, 8, 6, 2, 1, 10, 13, 3]
    array_after_zigzag = [3, 7, 4, 8, 2, 6, 1, 13, 3, 10]
    @test ClassicAlgorithmsCollections.zigzag_ordering(array_before_zigzag) ==
          array_after_zigzag

    arr = [5, 1, 3, 4, 7]
    sum = 12
    @test ClassicAlgorithmsCollections.count_triplet_elements(arr, sum) == 2

    array_pythagorean = [5, 1, 3, 4, 17, 8, 15, 2, 2, 13, 12]
    @test ClassicAlgorithmsCollections.count_pythagorean_elements(array_pythagorean) == 3

    arr_1 = [10, 15, 25]
    arr_2 = [5, 20, 30]
    sorted_combinations = [
        [10, 20],
        [10, 20, 25, 30],
        [10, 30],
        [15, 20],
        [15, 20, 25, 30],
        [15, 30],
        [25, 30],
    ]
    @test ClassicAlgorithmsCollections.combinations_of_2arrays(arr_1, arr_2) ==
          sorted_combinations
    subarray_in_array = [10, 12, 11, 9, 13, 14, 17, 18, 15]
    @test ClassicAlgorithmsCollections.find_maxlength_subarray(subarray_in_array) == 5

    @test ClassicAlgorithmsCollections.smallest_subset4sum(subarray_in_array, 52) == 4
    @test ClassicAlgorithmsCollections.smallest_subset4sum(subarray_in_array, 2) == 1
    @test ClassicAlgorithmsCollections.smallest_subset4sum(subarray_in_array, 5555) ==
          nothing

    sort_array = [1, 2, 5, 10, 20, 40]
    @test ClassicAlgorithmsCollections.find_smallest_nonelement(sort_array) == 4

    gradients_in_array = [100, 180, 260, 310, 40, 535, 695]
    @test ClassicAlgorithmsCollections.sum_of_postive_gradients(gradients_in_array) == 865
end

ext: jl

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
abstract type Transform end

struct pipeline
    transforms
    inplace::Bool
end

"""
    Pipeline(Transforms)

Constructs a transformation pipeline from vector/tuple of `Transforms`. The Transforms vector are effectively a vector of functions which transform data.
"""
Pipeline(Transforms) = pipeline(Transforms, false)

"""
    PipelineInPlace( X, FnStack...)

Construct a pipeline object from vector/tuple of `Transforms`. The Transforms vector are effectively a vector of functions which transform data.
This function makes "inplace" changes to the Array `X` as though it has been sent through the pipeline.
This is more efficient if memory is a concern, but can irreversibly transform data in memory depending on the transforms in the pipeline.
"""
function PipelineInPlace( X, FnStack...)
    pipe = Array{Any,1}(undef, length(FnStack))
    for (i, fn) in enumerate( FnStack )
        pipe[i] = isa(fn, Function) ? fn : fn(X)
        X .= pipe[i]( X )
    end
    return pipeline(Tuple(pipe), true)
end

"""
    Pipeline( X, FnStack... )

Construct a pipeline object from vector/tuple of `Transforms`. The Transforms vector are effectively a vector of functions which transform data.
"""
function Pipeline( X, FnStack... )
    pipe = Array{Any,1}( undef, length( FnStack ))
    for (i, fn) in enumerate( FnStack )
        pipe[ i ] = isa( fn, Function ) ? fn : fn( X )
        X = pipe[ i ]( X )
    end
    return pipeline(Tuple( pipe ), false)
end

"""
    (P::pipeline)(X; inverse = false)

Applies the stored transformations in a pipeline object `P` to data in X.
The inverse flag can allow for the transformations to be reversed provided they are invertible functions.
"""
function (P::pipeline)(X; inverse = false)
    if inverse
        if P.inplace
            for fn in reverse( P.transforms ); X .= fn( X; inverse = true ) ; end
        else
            foldr( ( p, X ) -> p(X; inverse = true), P.transforms, init = X)
        end
    else
        if P.inplace
            for fn in enumerate( X.transforms ); X .= fn( X ) ; end
        else
            foldl( ( X, p ) -> p(X), P.transforms, init = X)
        end
    end
end

struct QuantileTrim <: Transform
    Quantiles::Array
    invertible::Bool
end

"""
    QuantileTrim(Z; quantiles::Tuple{Float64,Float64} = (0.05, 0.95) )

Trims values above or below the specified columnwise quantiles to the quantile values themselves.
"""
function QuantileTrim(Z; quantiles::Tuple{Float64,Float64} = (0.05, 0.95) )
    @assert length(quantiles) == 2
    return QuantileTrim( EmpiricalQuantiles(Z, quantiles), false )
end

"""
    (T::QuantileTrim)(X, inverse = false)

Trims data in array `X` columns wise according to learned quantiles in QuantileTrim object `T`
This function does NOT have an inverse.
"""
function (T::QuantileTrim)(X, inverse = false)
    if inverse == false
        for c in 1:size(X)[2]
            if T.Quantiles[ 1, c ] != T.Quantiles[ 2, c ]
                lt = X[ : , c ] .< T.Quantiles[ 1, c ]
                gt = X[ : , c ] .> T.Quantiles[ 2, c ]
                X[ lt, c ] .= T.Quantiles[ 1, c ]
                X[ gt, c ] .= T.Quantiles[ 2, c ]
            end
        end
    else
        println("QuantileScale does not provide an inverse, skipping operation.")
    end
    return X
end

struct Center{B} <: Transform
    Mean::B
    invertible::Bool
end

"""
    Center(Z)

Acquires the mean of each column in `Z` provided and returns a transform that will subtract those column means from any future data.
"""
Center(Z::Array) = Center( StatsBase.mean(Z, dims = 1), true )

"""
    (T::Center)(Z; inverse = false)

Centers data in array `Z` column-wise according to learned mean centers in Center object `T`.
"""
(T::Center)(Z; inverse = false) = (inverse) ? (Z .+ T.Mean) : (Z .- T.Mean)

struct Scale{B} <: Transform
    StdDev::B
    invertible::Bool
end

"""
    Scale(Z)

Acquires the standard deviation of each column in `Z` provided and returns a transform that will divide those column-wise standard deviation from any future data.
"""
Scale(Z::Array) = Scale( StatsBase.std(Z, dims = 1),  true )

"""
    (T::Scale)(Z; inverse = false)

Scales data in array `Z` column-wise according to learned standard deviations in Scale object `T`.
"""
(T::Scale)(Z; inverse = false) = (inverse) ? (Z .* T.StdDev) : (Z ./ T.StdDev)

struct CenterScale{B,C} <: Transform
    Mean::B
    StdDev::C
    invertible::Bool
end

"""
    CenterScale(Z)

This is a composition of Center and Scale (in that order).
"""
function CenterScale(Z::Array)
    mu = StatsBase.mean(Z, dims = 1)
    stdev = StatsBase.std(Z, dims = 1)
    CenterScale( mu, stdev, true)
end

"""
    (T::CenterScale)(Z; inverse = false)

Centers and Scales data in array `Z` column-wise according to learned measures of central tendancy in Scale object `T`.
"""
(T::CenterScale)(Z; inverse = false) = (inverse) ? ((Z .* T.StdDev) .+ T.Mean) : ((Z .- T.Mean) ./ T.StdDev)

struct RangeNorm{B,C} <: Transform
    Mins::B
    Maxes::C
end

"""
    RangeNorm( Z )

Acquires the minimum and maximum of each column in `Z` provided and returns a transform that performs the following operation (Z - min(X))/(max(X) - min(X)) on any future data. This has the important effect of scaling all values observed in the range of `Z` to be between 0 and 1 with respect to each column.
"""
function RangeNorm( Z::Array )
    mins = minimum(Z, dims = 1)
    maxes = maximum(Z, dims = 1)
    return RangeNorm( mins, maxes)
end

"""
    (T::RangeNorm)(Z; inverse = false)

Scales and shifts data in array `Z` column-wise according to learned min-maxes in RangeNorm object `T`.
"""
(T::RangeNorm)(Z; inverse = false) = (inverse) ? ( (Z .* ( T.Maxes .- T.Mins ) ) .+ T.Mins) : (Z .- T.Mins) ./ ( T.Maxes .- T.Mins )

struct BoxCox <: Transform
    innercall
end

"""
    BoxCox(lambda)

Returns a BoxCox transform operator/function. To be used in a pipeline.
"""
BoxCox(lambda::Float64) = return BoxCox(X; inverse = false) = begin
    Z = zeros(size(X))
    if inverse
        if lambda != 0.0
            Z =  ((X .* lambda).+ 1.0) .^ (1/lambda)
        else
            Z = exp.(X)
        end
    else
        if lambda != 0.0
            Z = ((X .^ lambda) .- 1.0) ./ lambda
        else
            Z = log.(X)
        end
    end
    return Z
end

"""
    Logit(Z; inverse = false)

Logit transforms (```ln( X / (1 - X) ))```) every element in `Z`. The inverse may also be applied.
*Warning: This can return Infs and NaNs if elements of Z are not suited to the transform*
"""
function Logit(Z; inverse = false)
    if inverse
        return (exp.(Z) ./ (1.0 .+ exp.(Z)))
    else
        return log.( (Z ./ (1.0 .- Z)))
    end
end


text:
module plotCircuits

using PyCall
using PyPlot
using Juqst

export scheduleIt,plotIt

pushfirst!(PyVector(pyimport("sys")."path"), "./additionalPythonCode/PlotQCircuit/")

pqc = pyimport("plot_quantum_circuit")

function scheduleIt(t::Tableau)
    (li,la)=drawCircuit(t)
    pqc.plot_quantum_schedule(pqc.make_schedule(li),labels=la)
end

function plotIt(t::Tableau)
    (li,la)=drawCircuit(t)
    pqc.plot_quantum_circuit(li,labels=la)
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LibYAML_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LibYAML")
JLLWrappers.@generate_main_file("LibYAML", UUID("f531010a-71ee-5b11-8e71-2f13f29d8a82"))
end  # module LibYAML_jll


text:
@testset "Levels & States" begin
    @testset "Construction" begin
        csfs_3s_3p = csfs(c"3s 3p")
        csf_1 = first(csfs_3s_3p)
        @test_throws ArgumentError Level(csf_1, hi"2")
        @test_throws ArgumentError State(Level(csf_1, hi"1"), hi"2")

        @test string(Level(csf_1, hi"1")) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩"
        @test string(State(Level(csf_1, hi"1"), hi"-1")) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = -1⟩"

        @test sort([State(Level(csf_1, hi"1"), M_J) for M_J ∈ reverse(hi"-1":hi"1")]) ==
            [State(Level(csf_1, hi"1"), M_J) for M_J ∈ hi"-1":hi"1"]

        @test states.(csfs_3s_3p) == [[[State(Level(csf_1, hi"1"), M_J) for M_J ∈ hi"-1":hi"1"]],
                                      [[State(Level(csfs_3s_3p[2], J), M_J) for M_J ∈ -J:J] for J ∈ hi"0":hi"2" ]]

        rcsfs_3s_3p = csfs(rc"3s" ⊗ rcs"3p")

        @test states.(rcsfs_3s_3p)== [[[State(Level(rcsfs_3s_3p[1], hi"0"), hi"0")]],
                                      [[State(Level(rcsfs_3s_3p[2], hi"1"), M_J) for M_J ∈ hi"-1":hi"1"]],
                                      [[State(Level(rcsfs_3s_3p[3], hi"1"), M_J) for M_J ∈ hi"-1":hi"1"]],
                                      [[State(Level(rcsfs_3s_3p[4], hi"2"), M_J) for M_J ∈ hi"-2":hi"2"]]]
    end
end


text:
using SignalTemps
using Test

@testset "SignalTemps.jl" begin
    # Write your tests here.
end


text:
abstract type Setter <: KaleidoLens end

Setfield.get(::Any, setter::Setter) = error("Setters do not support `get`.")

struct NullSetter <: Setter end
Setfield.set(obj, ::NullSetter, ::Any) = obj

"""
    nullsetter :: Setter

A setter that does nothing; i.e., `set(x, nullsetter, y) === x` for any
`x` and `y`.

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> set(1, nullsetter, 2)
1
```
"""
const nullsetter = NullSetter()

Base.show(io::IO, ::NullSetter) = print(io, "nullsetter")
# TODO: handle prefix

"""
    ToField(f) :: Setter

Apply `f` when setting.  Use `x -> get(x, f)` if `f` is a `Lens`.

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> setter = (@lens _.x) ∘ ToField(@lens _.a)
(@lens _.x) ∘ (←(@lens _.a)|❌→)

julia> set((x = 1, y = 2), setter, (a = 10, b = 20))
(x = 10, y = 2)
```
"""
struct ToField{F} <: Setter
    f::F
end

_eval(f, x) = f(x)
_eval(f::Lens, x) = get(x, f)

Setfield.set(::Any, setter::ToField, x) = _eval(setter.f, x)


text:
# Aliev_Panfilov_2D_compute_v4.jl
# ================================   u.p. 5.8.21  / 5.9.21

# Aliev-Panfilov model
# with ODE-solver
# compute a save solution

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using PyPlot, OrdinaryDiffEq
using JLD
using Random
using Statistics
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

function Aliev_Panfilov_ODE(du,u,p,t)
    d, k, a, epsilon, mu1, mu2, hsq6  = p

    uu = u[1,:,:]
    vv = u[2,:,:]
    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]

    # 9-point stencil
    diff_term = d .* (
                4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+
                      uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+
                    uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+
                    uuu[1:end-2,1:end-2] .- 20 .*uu  )  ./ hsq6

    du[1,:,:] = k .* uu .*(1 .-uu).*(uu.-a) .- uu.*vv .+ diff_term
    du[2,:,:] = (epsilon .+ (mu1 .* vv) ./ (mu2 .+ uu))  .* (.- vv .- k .* uu .* (uu .- a .- 1 ) )
end

@time begin

# new simulation starting from initial conditions
# -----------------------------------------------
L = 100  #  size of domain
N = 400  #  no. of grid points
h = L/(N-1) # spatial steps size
hsq6 = 6*h*h

# parameters
k = 8
a = 0.05
mu1 = 0.2
mu2 = 0.3
d = 0.2   #  diffusion
epsilon = 0.002

# initial values
# ---------------
uu = zeros(N,N)  #  u steady state
vv = zeros(N,N)  #  v steady state
uu[1:end,1:20] .= 0.9   # N = 400
# uu[1:end,1:10] .= 0.9   # N = 200

u0 = zeros(2,N,N)
u0[1,:,:] = uu
u0[2,:,:] = vv

p = d, k, a, epsilon, mu1, mu2, hsq6
saveat = (0.0, 120.0) # [0. , 20. , 120.]
tspan = (0., saveat[end])

prob = ODEProblem(Aliev_Panfilov_ODE, u0, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat = saveat,
progress = true, progress_steps = 100,
)
uout = sol.u

# perturbation
# -------------
ures = uout[end]
uu = ures[1,:,:]
vv = ures[2,:,:]

# uu[1:30,:] .= 1.  # N = 200
uu[1:60,:] .= 1.  # N = 400

u0[1,:,:] = uu
u0[2,:,:] = vv

# evolve forward
# --------------
# saveat = [700. , 710., 720., 730., 740., 750., 760., 770., 780.]
# tspan = (0., saveat[end])
saveat = [600]
tspan = (0., saveat[end])

prob = ODEProblem(Aliev_Panfilov_ODE, u0, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat = saveat,
progress = true, progress_steps = 100,
)

# compute snapshots
# -----------------
uout = sol.u
ures = uout[end]

saveat = 1.0
tspan = (0., 200)

prob = ODEProblem(Aliev_Panfilov_ODE, ures, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat = saveat,
progress = true, progress_steps = 100,
)

end # @time

tvec = sol.t
uout = sol.u
nts = length(tvec)

# plot
# ----
fig = plt.figure("AP_2D",figsize=(figx,2*figy))
iplt = 0

for n = 1:9
    iplt = iplt + 1
    ures = uout[n]
    subplot(3,3,iplt)
    im = plt.pcolormesh(ures[1,:,:],cmap = "inferno")  # "gnuplot" ) #"inferno")

    ax = gca()
    ax.set_xticks([])
    ax.set_yticks([])

    plt.title(L"t = "*string(floor(Int,tvec[n])), fontsize = 24)
end

fig.tight_layout(pad=0.3)
path = datadir("Aliev_Panfilov_snapshots.jld")
# save(path, "d", d, "k", k, "a", a, "epsilon", epsilon, "mu1", mu1, "mu2", mu2, "L", L, "N", N, "tvec", tvec, "uout", uout)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
include("../matplotlib_settings.jl")

using GeometricExamples

include("guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

equ = problem.guiding_center_4d_iode()

include("guiding_center_4d.jl")

run_tableau_list(tableau_list, equ, Δt, ntime)

generate_html(tableau_list, figures, prefix, fig_suff, fig_dir, out_dir, html_file;
              title=coll_title, params=get_params(Δt, ntime, problem))


text:
"""
    sphericalhankel2(nu, x)

Spherical Hankel function of the second kind, order `nu`
and argument `x`.
"""
function sphericalhankel2(nu, x)
    jn = SpecialFunctions.sphericalbesselj(nu, x)
    yn = SpecialFunctions.sphericalbessely(nu, x)
    return jn - im*yn
end

"""
    sphericalhankel1(nu, x)

Spherical Hankel function of the first kind, order `nu`
and argument `x`.
"""
function sphericalhankel1(nu, x)
    return conj(sphericalhankel2(nu, x))
end

"""
    sphericalhankel2_and_derivatives(nu::UnitRange, x)
    sphericalhankel2_and_derivatives(nu::Number, x)

Evaluates the spherical Hankel function of the second kind, its first and
second derivative, for orders `nu::UnitRange` and argument `x`.
"""
function sphericalhankel2_and_derivatives(nu::UnitRange, x)
    start = nu.start
    stop = nu.stop
    len = stop-start+1   # number of terms to compute
    sph = sphericalhankel2.(start:stop+2, x)  # compute two more terms
    # compute first derivative using recurrence
    # h_n'(x) = -h_{n+1}(x) + n/x*h_n(z)
    sph_d = similar(sph, len+1)   # compute one more term
    n = stop+1   # order
    for index in len+1:-1:1
        sph_d[index] = -sph[index+1] + n/x*sph[index]
        n -= 1
    end
    # compute second derivative using recurrence
    # h_n''(x) = n / x^2 * (x * h_n'(x) - h_n(x)) - h_{n+1}'(x)
    sph_dd = similar(sph, len)
    n = stop   # order
    for index in len:-1:1
        sph_dd[index] = n / x^2 * (x*sph_d[index] - sph[index]) - sph_d[index+1]
        n -= 1
    end
    return sph[1:len], sph_d[1:len], sph_dd
end
function sphericalhankel2_and_derivatives(nu::Number, x)
    sph, sph_d, sph_dd = sphericalhankel2_and_derivatives(nu:nu, x)
    return first(sph), first(sph_d), first(sph_dd)
end

"""
    sphericalhankel1_and_derivatives(nu::UnitRange, x)
    sphericalhankel1_and_derivatives(nu::Number, x)

Evaluates the spherical Hankel function of the first kind, its first and
second derivative, for orders `nu` and argument `x`.
"""
function sphericalhankel1_and_derivatives(nu, x)
    return conj.(sphericalhankel2_and_derivatives(nu, x))
end

"""
    riccatihankel2_and_derivatives(nu::UnitRange, x)

Evaluates the spherical Riccati-Bessel Hankel function of the second kind, its first and
second derivative, for orders `nu::UnitRange` and argument `x`.
"""
function riccatihankel2_and_derivatives(nu::UnitRange, x)
    sph, sph_d, sph_dd = sphericalhankel2_and_derivatives(nu, x)
    rh = x.*sph
    rh_d = sph + x.*sph_d
    rh_dd = 2*sph_d + x.*sph_dd
    return rh, rh_d, rh_dd
end

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
using BenchmarkTools
using DataFrames
using Transducers
using VegaLite

resultpath = joinpath(@__DIR__, "result.json")
result, = BenchmarkTools.load(resultpath)

df_raw =
    BenchmarkTools.leaves(result) |>
    Map() do ((basesize, needleloc, ex), trial)
        (
            basesize = parse(Int, basesize),
            needleloc = parse(Int, needleloc),
            executor = Symbol(ex),
            trial = trial,
        )
    end |>
    DataFrame
#-

begin
    df_tmp = select(df_raw, Not(:trial))
    df_tmp[!, :minimum] = map(trial -> minimum(trial).time, df_raw.trial)
    df_tmp[!, :median] = map(trial -> median(trial).time, df_raw.trial)
    df_tmp[!, :memory] = map(trial -> trial.memory, df_raw.trial)
    df_stats = stack(
        df_tmp,
        [:minimum, :median],
        variable_name = :time_stat,
        value_name = :time_ns,
    )
end
#-

begin
    df = transform(groupby(df_stats, [:basesize, :needleloc, :time_stat])) do group
        baseline = only(eachrow(filter(:executor => ==(:ThreadedEx), group)))
        (speedup = baseline.time_ns ./ group.time_ns,)
    end
    filter!(:executor => !(==(:ThreadedEx)), df)
    df
end
#-

plt1 = @vlplot(
    facet = {row = {field = :executor, type = :nominal}, column = {field = :time_stat, type = :nominal}},
    spec = {
        layer = [
            {
                mark = {type = :line, point = true},
                encoding = {
                    x = {field = :needleloc, scale = {type = :log, base = 2}},
                    y = {field = :speedup, axis = {title = "Speedup (T_default / T_\$executor)"}},
                    color = {field = :basesize, type = :nominal},
                },
            },
            {
                mark = :rule,
                encoding = {y = {datum = 1}},
            },
        ],
    },
    data = df,
)
#-

plt2 = @vlplot(
    mark = {type = :line, point = true},
    x = {field = :needleloc, scale = {type = :log, base = 2}},
    y = {field = :time_ns, axis = {title = "Time [ns]"}},
    color = {field = :basesize, type = :nominal},
    row = :executor,
    column = :time_stat,
    resolve = {scale = {y = "independent"}},
    data = df_stats,
)
#-

plt3 = @vlplot(
    mark = {type = :line, point = true, clip = true},
    x = {field = :needleloc, scale = {type = :log, base = 2}},
    y = {field = :time_ns, axis = {title = "Time [ns]"}, scale = {domain = [0, 800_000]}},
    color = {field = :basesize, type = :nominal},
    row = :executor,
    column = :time_stat,
    data = df_stats,
)
#-


text:
# error codes
const NO_ERROR = 0
const ERROR_INTERNAL = 1	        # unexpected error
# parse errors
const ERROR_BAD_ESCAPE = 2	        # bad escape sequence
const ERROR_BAD_CHAR_CLASS = 3	    # bad character class
const ERROR_BAD_CHAR_RANGE = 4	    # bad character class range
const ERROR_MISSING_BRACKET = 5	    # missing closing ]
const ERROR_MISSING_PAREN = 6	    # missing closing )
const ERROR_TRAILING_BACKSLASH = 7  # trailing \ at end of regexp
const ERROR_REPEAT_ARGUMENT = 8	    # repeat argument missing, e.g. "*"
const ERROR_REPEAT_SIZE = 9	        # bad repetition argument
const ERROR_REPEAT_OP = 10		    # bad repetition operator
const ERROR_BAD_PERL_OP = 11	    # bad perl operator
const ERROR_BAD_UTF8 = 12		    # invalid UTF-8 in regexp
const ERROR_BAD_NAMED_CAPTURE = 13	# bad named capture group
const ERROR_PATTERN_TOO_LARGE = 14	# pattern too large (compile failed)

# compile options
const UTF8 = 1
const Latin1 = 2

# bitmasks for compile options
const LOG_ERRORS =      UInt32(0x00000001)
const LATIN1 =          UInt32(0x00000002)
const LONGEST_MATCH =   UInt32(0x00000004)
const LITERAL =         UInt32(0x00000008)
const NEVER_NL =        UInt32(0x00000010)
const DOTALL =          UInt32(0x00000020)
const NEVER_CAPTURE =   UInt32(0x00000040)
const CASELESS =        UInt32(0x00000080)
const POSIX_SYNTAX =    UInt32(0x00000100)
const PERL_CLASSES =    UInt32(0x00000200)
const WORD_BOUNDARY =   UInt32(0x00000400)
const MULTILINE =       UInt32(0x00000800)
const COMPILE_MASK =    UInt32(0x00000fff)

# matching options
const UNANCHORED   = 1
const ANCHOR_START = 2
const ANCHOR_BOTH  = 3

# bitmasks for match options
const ANCHORED =        UInt32(0x00000001)
const ENDANCHORED =     UInt32(0x00000002)
const EXECUTE_MASK =    UInt32(0x00000003)

const DEFAULT_MAX_MEM = UInt32(0x00800000)  # 8 MB


text: USE_HISTORIC = false


text:
module JuliaForHEP

const __revise_mode__ = :sigs

using Flux, Zygote, CUDA, WeightedOnlineStats, ProgressMeter, Statistics
using Pandas: read_hdf
import Random

export accuracy, extract_cols, scale_weights, split_indices, split_datasets, weighted_mean,
    step!, prefix_labels

"""
    accuracy(ŷ::AbstractMatrix, y::AbstractMatrix{Bool})

Calculate accuracy of predictions `ŷ` against ground truth `y`.
"""
function accuracy(ŷ::AbstractMatrix{T}, y::AbstractMatrix{Bool}) where {T}
    return mean(Flux.onecold(ŷ) .== Flux.onecold(y))
end

"""
    extract_cols(h5files, column_names)

Extract features given by `column_names` from all files in `h5files` and concatenate them.
"""
function extract_cols(h5files, column_names, idx=:)
    return mapreduce(hcat, h5files) do f
        df = read_hdf(f)
        PermutedDimsArray(Array(df[column_names][idx]), (2, 1))
    end
end

"""
    scale_weights(weights, outputs)

Scale `weights`, so that the sum over all weights for each output feature is always one.
"""
function scale_weights(weights, outputs)
    sums_weights = Dict{Symbol,Float32}()
    for (i, class) in pairs(outputs)
        sums_weights[class] = get(sums_weights, class, 0f0) + weights[i]
    end
    return map(CartesianIndices(weights)) do i
        weights[i] / sums_weights[outputs[i]]
    end
end

function split_indices(indices, proportions; shuffle=true, rng=Random.GLOBAL_RNG)
    if shuffle
        indices = Random.shuffle(rng, indices)
    end
    split = [0; round.(Int, cumsum(proportions) .* length(indices))]

    if split[end] > length(indices)
        throw(ArgumentError("Invalid proportions $proportions, sum is > 1"))
    elseif split[end] < length(indices)
        @warn "$(length(indices) - slit[end]) indices not assigned to any dataset"
    end

    return ntuple(length(split) - 1) do i
        indices[split[i]+1:split[i+1]]
    end
end

function split_datasets(datasets, proportions; kwargs...)
    return split_datasets(identity, datasets, proportions; kwargs...)
end
function split_datasets(f, datasets, proportions; kwargs...)
    if !all(ds -> axes(ds)[end] == axes(datasets[1])[end], datasets[2:end])
        throw(ArgumentError("Last axes of datasets don't match"))
    end
    return map(split_indices(axes(datasets[1])[end], proportions; kwargs...)) do idx
        map(datasets) do ds
            f(ds[ntuple(_ -> (:), ndims(ds)-1)..., idx])
        end
    end
end

"""
    weighted_mean(x, weights)

AD-friendly weighted mean. `weights` is always normalized.
"""
weighted_mean(x, weights) = sum(x .* weights ./ sum(weights))

function step!(model, ds_train, ds_tests...; loss_function, optimizer, measures=(;))
    ps = Flux.params(model)

    trainmode!(model)

    measures_train = map(_ -> WeightedMean(), measures)
    p = Progress(div(ds_train.nobs, ds_train.batchsize, RoundUp))
    for (x, y, weights) in ds_train
        local ŷ
        loss, pb = Zygote.pullback(ps) do
            ŷ = model(x)
            return loss_function(ŷ, y; weights)
        end
        gs = pb(one(loss))
        Flux.update!(optimizer, ps, gs)

        state = (; x, y, weights, loss, ŷ, mode=:train)
        for (k, f) in pairs(measures)
            fit!(measures_train[k], f(state), size(y, 2))
        end
        next!(p)
    end

    testmode!(model)

    measures_tests = map(enumerate(ds_tests)) do (i, ds_test)
        measures_test = map(_ -> WeightedMean(), measures)
        for (x, y, weights) in ds_test
            ŷ = model(x)
            loss = loss_function(ŷ, y; weights)

            state = (; x, y, weights, loss, ŷ, mode=Symbol(:test_, i))
            for (k, f) in pairs(measures)
                fit!(measures_test[k], f(state), size(y, 2))
            end
        end
        return measures_test
    end

    return map.(value, (measures_train, measures_tests...))
end

function prefix_labels(nt::NamedTuple{keys}, prefix) where {keys}
    return NamedTuple{Symbol.(prefix, keys)}(Tuple(nt))
end

export CoLa, Linear, LoLa, LorentzSidechain

using Tullio, CUDA, KernelAbstractions, Adapt, Compat, LinearAlgebra
using StaticArrays
using ChainRulesCore

# 🏴☠
Flux.params!(p::Zygote.Params, k::AbstractArray{<:SArray}, seen=IdSet()) = push!(p, k)
Flux.fmap(f, a::AbstractArray{<:SArray}) = f(a)
using ChainRulesCore: backing
Zygote.accum(x::Union{NamedTuple,Composite}, y::Union{NamedTuple,Composite}) = Zygote.accum(backing(x), backing(y))

function ChainRulesCore.rrule(::typeof(reinterpret), ::typeof(reshape), ::Type{T}, x::AbstractArray) where {T}
    function reinterpret_reshape_pullback(Δ)
        return NO_FIELDS, NO_FIELDS, NO_FIELDS, reinterpret(reshape, eltype(x), Δ)
    end
    return reinterpret(reshape, T, x), reinterpret_reshape_pullback
end

#include("reinterpret_reshape.jl")

#macro tullio(x...)
#    return Tullio.var"@tullio"(__source__, __module__, :(verbose=2), x...)
#end

include("cola.jl")
include("lola.jl")
include("lorentz_sidechain.jl")

end


text:
using Statistics
using Test
using GeoPhyInv

# testing attenuation
#=
np2 = nextpow(2, 2 * length(tgrid));
fnpow2grid = FFTW.fftfreq(np2, inv(step(tgrid)));
fill!(medium.m[:Q], 100)
medium.fc = NamedArray([20.0, 100 * 2 * pi], ([:freqmin, :freqmax],))
medium.ic = NamedArray([5], ([:nsls],))

vp0 = medium.ref[:vp]
rho0 = medium.ref[:rho]
Ku = rho0 * vp0 * vp0
tau_epsilon = medium[:tau_epsilon]
tau_sigma = medium[:tau_sigma]

Kc = GeoPhyInv.complexK(Ku, 2.0 * pi .* fnpow2grid, tau_sigma, tau_epsilon)
=#

function check()
    # least-squares misfit
    paerr = GeoPhyInv.VNamedD_misfit(rec1, pa.c.data[1])
    err = GeoPhyInv.func_grad!(paerr)

    # normalize error
    error = err[1] / paerr.ynorm

    # desired accuracy?
    @test error < 1e-2
end

medium = Medium(:elastic_homo3D);
ageom=AGeom(medium.mgrid, SSrcs(2), [Srcs(3), Srcs(2)], [Recs(10), Recs(20)])
wav, tgrid=ricker(medium, 4, 1, 0.4)
srcwav = SrcWav(tgrid, ageom, [:p])
update!(srcwav, [:p], wav)

pa = SeisForwExpt(
    FdtdElastic(),
    npw = 1,
    medium = medium,
    ageom = [ageom],
    srcwav = [srcwav],
    sflags = [2],
    rflags = [1],
    tgrid = tgrid,
    verbose = true,
);

@time update!(pa);

vp0 = mean(medium[:vp])
rho0 = mean(medium[:rho])
rec1 = GeoPhyInv.Born.mod(
    medium,
    medium_pert = medium,
    ageom = ageom,
    srcwav = srcwav,
    tgridmod = tgrid,
    src_flag = 2,
)

check()


text:
"""
    ComputeLimberArray!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities)

This function compute the Limber grid. In the Limber approximation, ``k`` and
``z`` are related by the following relation:
```math
k_\\ell(z)=\\frac{\\ell+1/2}{r(z)}.
```
"""
function ComputeLimberArray!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities)
    for (idx_l, myl) in enumerate(CosmologicalGrid.ℓBinCenters)
        CosmologicalGrid.KLimberArray[idx_l, :] =  (myl + 0.5) ./
        BackgroundQuantities.rZArray
    end
end

"""
    InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
    ::InterpolationMethod)

This function interpolates the Power Spectrum on the ``k-z`` grid and evaluates
it on the Limber grid.
"""
function InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
    ::RectBivSplineDierckx)
    InterpPmm = Dierckx.Spline2D(log10.(CosmologicalGrid.KArray),
    CosmologicalGrid.ZArray, log10.(PowerSpectrum.PowerSpectrumNonlinArray);
    kx=5, ky=5, s=0.0)
    for (idx_l, myl) in enumerate(CosmologicalGrid.ℓBinCenters)
        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
        CosmologicalGrid.ZArray))
    end
end

function InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
    ::BSplineCubic)
    x = LinRange(log10(first(CosmologicalGrid.KArray)),
    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
    length(CosmologicalGrid.ZArray))
    InterpPmm = Interpolations.interpolate(
    log10.(PowerSpectrum.PowerSpectrumNonlinArray),
    BSpline(Cubic(Line(OnGrid()))))
    InterpPmm = scale(InterpPmm, x, y)
    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
    for (idx_l, myl) in enumerate(CosmologicalGrid.ℓBinCenters)
        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
        CosmologicalGrid.ZArray))
    end
end

function InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
    ConvolvedDensity::ConvolvedDensity)
    x = LinRange(log10(first(CosmologicalGrid.KArray)),
    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
    length(CosmologicalGrid.ZArray))
    InterpPmm = Interpolations.interpolate(
    log10.(PowerSpectrum.PowerSpectrumLinArray),
    BSpline(Cubic(Line(OnGrid()))))
    InterpPmm = scale(InterpPmm, x, y)
    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
    for iidx in 1:length(ConvolvedDensity.ZBinArray)-1
        for lidx in 1:length(CosmologicalGrid.ℓBinCenters)
            PowerSpectrum.InterpolatedPowerSpectrumBeyondLimber[lidx, :] =
            10 .^(InterpPmm.(log10.(
            CosmologicalGrid.KBeyondLimberArray[lidx, :]),
            CosmologicalGrid.ZArray[1]))
        end
    end
end

function ExtractGrowthFactor!(PowerSpectrum::PowerSpectrum)
    PowerSpectrum.GrowthFactor = zeros(
        length(PowerSpectrum.PowerSpectrumLinArray[1,:]))
    for zidx in 1:length(PowerSpectrum.GrowthFactor)
        PowerSpectrum.GrowthFactor[zidx] =
        sqrt(PowerSpectrum.PowerSpectrumLinArray[1,zidx]./
        PowerSpectrum.PowerSpectrumLinArray[1,1])
    end
end


text:
using Photometry.Aperture: bounds, center,
                           oblique_coefficients

@testset "Apertures" begin
    ap_ellipse = EllipticalAperture(0, 0, 20, 10, 0)
    @test center(ap_ellipse) == (0, 0)
    @test bounds(ap_ellipse) == (-20, 20, -10, 10)
    @test size(ap_ellipse) == (41, 21)
    @test size(ap_ellipse, 1) == 41
    @test all(axes(ap_ellipse) .== (-20:20, -10:10))
    @test eachindex(ap_ellipse) == CartesianIndex(-20, -10):CartesianIndex(20, 10)

    @test EllipticalAperture([0, 0], 20, 10, 0) == ap_ellipse

    ap_ellipse = EllipticalAperture(0, 0, 2, 1, 45)
    @test bounds(ap_ellipse) == (-1, 1, -1, 1)
end

@testset "Elliptical Aperture" begin
    e = EllipticalAperture(0, 0, 20, 10, 0)
    @test sprint(show, e) == "EllipticalAperture(0, 0, a=20, b=10, θ=0°)"
end

@testset "oblique_coefficients" begin
    @test all(oblique_coefficients(2, 2, 0) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 2, 90) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 1, 30) .≈ (7 / 16, 13 / 16, -6sqrt(3) / 16))
end

@testset "Elliptical Annulus" begin
    e0 = EllipticalAnnulus(0, 0, 8, 16, 4, 45)
    @test sprint(show, e0) == "EllipticalAnnulus(0.0, 0.0, a_in=8.0, a_out=16.0, b_in=2.0, b_out=4.0, θ=45.0°)"
    @test EllipticalAnnulus([0, 0], 8, 16, 4, 45) == EllipticalAnnulus(0, 0, 8, 16, 4, 45)
end

@testset "Elliptical Annulus bounding box" begin
    e = EllipticalAnnulus(0, 0, 8, 16, 4, 0)
    @test center(e) == (0, 0)
    @test bounds(e) == (-16, 16, -4, 4)
    @test size(e) == (33, 9)
    @test all(axes(e) .== (-16:16, -4:4))
    @test eachindex(e) == CartesianIndex(-16, -4):CartesianIndex(16, 4)
end

@testset "regression circular ellipse" begin
    # some weird bug where centered on-grid past a certain size (and rotated) would fail
    e = EllipticalAperture(5.5, 5.5, 4, 4, 20)
    @test sum(e) ≈ photometry(e, ones(9, 9)).aperture_sum # just a test that no errors occur
end


text:
TypeObject = Union{DataType, Union, UnionAll}

"""representing the in-matching object in pattern compile time
"""
struct Target{IsComplex}
    repr::Any
    type::Ref{TypeObject}
end

@nospecialize

function target_method end

function Base.getproperty(target::Target, sym::Symbol)
    target_method(target, Val(sym))
end

function target_method(target::Target, ::Val{:type_narrow!})
    function (ty::TypeObject)
        tyref = getfield(target, :type)
        if tyref[] <: ty
        else
            tyref[] = ty
        end
    end
end

function target_method(target::Target, ::Val{:type})
    getfield(target, :type)[]
end

function target_method(target::Target, ::Val{:repr})
    getfield(target, :repr)
end

function target_method(target::Target{IsC}, ::Val{:with_repr}) where {IsC}
    function ap(repr::Any)
        Target{IsC}(repr, getfield(target, :type))
    end
    function ap(repr::Any, ::Val{IsC′}) where {IsC′}
        Target{IsC′}(repr, getfield(target, :type))
    end
    ap
end

function target_method(target::Target{IsC}, ::Val{:with_type}) where {IsC}
    function ap(ty::TypeObject)
        Target{IsC}(target.repr, Ref{TypeObject}(ty))
    end
    function ap(ty::TypeObject, ::Val{IsC′}) where {IsC′}
        Target{IsC′}(target.repr, Ref{TypeObject}(ty))
    end
    ap
end

function target_method(target::Target{IsC}, ::Val{:clone}) where {IsC}
    Target{IsC}(target.repr, Ref{TypeObject}(target.type))
end

@specialize


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

using Markdown
using InteractiveUtils

# ╔═╡ 4c508010-c5a2-11eb-1fc1-19d8e9de32c4
md"""
### Test Code for Module NB
"""

# ╔═╡ 942227ca-947b-497e-a20b-c7a251e72f74
module script include("./script.jl") end

# ╔═╡ da372bf3-dd21-40fe-84af-eee8722a401a
NB = script.NB

# ╔═╡ bf0086c8-58b1-442b-ac75-f61a30349c63
module tools include("../tools.jl") end

# ╔═╡ 1d91c553-ba63-4c45-97c8-c14effbe5712
JuTools = tools.JuTools

# ╔═╡ 13e3bf10-217d-40d4-a23d-1f4f6b450b1e
import Plots

# ╔═╡ 22fe721d-390b-42a0-9abe-44ab67afdd38
import Random

# ╔═╡ d90f92cf-1ab3-4a7f-a1ba-c32a2ca76a04
begin
X_data1, Y_data1 = JuTools.data_generate_cluster_2d(pos1=(20.0, 20.0), pos2=(50.0, 50.0),
    radius1=5.0, radius2=5.0, random_scale=8.0, data_size=1000)
size(X_data1), size(Y_data1)
end

# ╔═╡ 23f676f9-430e-4264-b059-2494c00bb3ba
let
trained = NB.train(X_data1, Y_data1)
plot_X1 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_Y1 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_X2 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
plot_Y2 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
X_test_1 = 0:70
X_test_2 = 0:70
X_test_contour = hcat(vec(repeat(reshape(X_test_1, :, 1), 1, 71)'), vec(repeat(X_test_2, 71, 1)))
X_test_contour = reshape(NB.predict(X_test_contour, trained), 71, :)
Plots.gr()
Plots.contour(X_test_1, X_test_2, X_test_contour, fill=true, c=Plots.cgrad(:matter, rev=false, categorical=false),
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.xlims!((0.0, 70.0))
Plots.ylims!((0.0, 70.0))
Plots.title!("Naive Bayes (Gaussian)")
end

# ╔═╡ Cell order:
# ╟─4c508010-c5a2-11eb-1fc1-19d8e9de32c4
# ╠═942227ca-947b-497e-a20b-c7a251e72f74
# ╠═da372bf3-dd21-40fe-84af-eee8722a401a
# ╠═bf0086c8-58b1-442b-ac75-f61a30349c63
# ╠═1d91c553-ba63-4c45-97c8-c14effbe5712
# ╠═13e3bf10-217d-40d4-a23d-1f4f6b450b1e
# ╠═22fe721d-390b-42a0-9abe-44ab67afdd38
# ╠═d90f92cf-1ab3-4a7f-a1ba-c32a2ca76a04
# ╠═23f676f9-430e-4264-b059-2494c00bb3ba

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using LsqFit

model(t, p) = p[1] * exp.(-p[2] * t)

tdata = range(0,10, length=20)
ydata = model(tdata, [1.0, 2.0]) + 0.01*randn(length(tdata))

p0 = [0.5, 0.5]

fit = curve_fit(model, tdata, ydata, p0)

using Plots
plotly()

pt = plot(title="Nonlinear Regression", size=(600,400))
plot!(pt, tdata, ydata, seriestype=:scatter)

a = fit.param[1]
b = fit.param[2]
fitmodel(t) = a * exp.(-b*t)
#plot!(pt, fitmodel, extrema(tdata)..., seriestype=:line)
x = range(0,10,length=100)
y = @. fitmodel(x)
plot!(pt, x, y, seriestype=:line)
gui(pt)
nothing


text:
include("FEM1D.jl")
include("FEM2D.jl")

function Laplace(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    F = [x -> -1, x -> 0, x -> 0, x -> x^2]
    BCs_ = [[1, BCs[1]], [2, BCs[2]]]
    A, _, _ = BVP_1D(L, F, BCs_)
    return A
end

function Laplace_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    F = [x -> -1, x -> 0, x -> -1, x -> 0, x -> 0, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Helmholtz(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    κ = parameters[1]
    F = [x -> -1, x -> 0, x -> -κ^2*1, x -> x^2]
    BCs_ = [[1, BCs[1]], [2, BCs[2]]]
    A, _, _ = BVP_1D(L, F, BCs_)
    return A
end

function Helmholtz_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    κ = parameters[1]
    F = [x -> -1, x -> 0, x -> -1, x -> 0, x -> 0, x -> -κ^2*1, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Convection_Diffusion(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    v = parameters[1]
    F = [x -> -1, x -> -v*1, x -> 0, x -> x^2]
    BCs_ = [[1, BCs[1]], [2, BCs[2]]]
    A, _, _ = BVP_1D(L, F, BCs_)
    return A
end

function Convection_Diffusion_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    v1 = parameters[1]
    v2 = parameters[2]
    F = [x -> -1, x -> 0, x -> -1, x -> v1*1, x -> v2*1, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Disc_Diffusion(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    σ = parameters[1]
    F = [x -> -1-σ*(x>0.5), x -> 0, x -> 0, x -> x^2]
    BCs_ = [[1, BCs[1]], [2, BCs[2]]]
    A, _, _ = BVP_1D(L, F, BCs_)
    return A
end

function Disc_Diffusion_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    σ = parameters[1]
    g(x) = -1*(x[1]<=0.5)*(x[2]>=0.5)-1*(x[1]>=0.5)*(x[2]<=0.5)-σ*(x[1]>0.5)*(x[2]>0.5)-1/σ*(x[1]<0.5)*(x[2]<0.5)
    F = [g, x -> 0, g, x -> 0, x -> 0, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Implicit_Diffusion(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    μ = parameters[1]
    F = [x -> μ*1, x -> 0, x -> 0, x -> x^2]
    BCs_ = [[1, BCs[1]], [2, BCs[2]]]
    A, _, _ = BVP_1D(L, F, BCs_)
    A = I - A
    return A
end

function Implicit_Diffusion_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    μ = parameters[1]
    F = [x -> -1, x -> 0, x -> -1, x -> 0, x -> 0, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    A = I - μ*A
    return A
end

function Anisotropic_Laplace_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    ϵ = parameters[1]
    F = [x -> -1, x -> 0, x -> -ϵ*1, x -> 0, x -> 0, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Mixed_Derivative_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    ϵ = parameters[1]
    F = [x -> -1, x -> -ϵ*1, x -> -1, x -> 0, x -> 0, x -> 0, x -> 0]
    A, _ = BVP_2D(L, F, BCs)
    return A
end

function Biharmonic_equation_2D(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    rows, columns, values = [], [], Array{Float64,1}(undef, 0)
    lex(i, j) = i + N*(j-1)
    for (i, j) = collect(Iterators.product(1:N, 1:N))
        for (p, q) = collect(Iterators.product([+1, -1], [+1, -1]))
            if (1<=i+p<=N) & (1<=j+q<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i+p, j+q))
                push!(values, +2)
            end
        end
        for p = [+1, -1]
            if (1<=i+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i+p, j))
                push!(values, -8)
            end
            if (1<=j+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i, j+p))
                push!(values, -8)
            end
        end
        for p = [+2, -2]
            if (1<=i+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i+p, j))
                push!(values, +1)
            end
            if (1<=j+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i, j+p))
                push!(values, +1)
            end
        end
        if (2<=i<=(N-1)) & (2<=j<=(N-1))
            push!(rows, lex(i, j))
            push!(columns, lex(i, j))
            push!(values, +20)
        else
            push!(rows, lex(i, j))
            push!(columns, lex(i, j))
            if ((i == 1) & (j == 1)) || ((i == N) & (j == N)) || ((i == N) & (j == 1)) || ((i == 1) & (j == N))
                push!(values, +22)
            else
                push!(values, +21)
            end
        end
    end
    A = sparse(rows, columns, values)
end

function Mehrstellen(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    rows, columns, values = [], [], Array{Float64,1}(undef, 0)
    lex(i, j) = i + N*(j-1)
    for (i, j) = collect(Iterators.product(1:N, 1:N))
        for (p, q) = collect(Iterators.product([+1, -1], [+1, -1]))
            if (1<=i+p<=N) & (1<=j+q<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i+p, j+q))
                push!(values, -1)
            end
        end
        for p = [+1, -1]
            if (1<=i+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i+p, j))
                push!(values, -4)
            end
            if (1<=j+p<=N)
                push!(rows, lex(i, j))
                push!(columns, lex(i, j+p))
                push!(values, -4)
            end
        end
        push!(rows, lex(i, j))
        push!(columns, lex(i, j))
        push!(values, +20)
    end
    A = sparse(rows, columns, values)
end

function fourth_Laplace_2D_1(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    rows, columns, values = [], [], Array{Float64,1}(undef, 0)
    lex(i, j) = i + N*(j-1)
    for (i, j) = collect(Iterators.product(1:N, 1:N))
        for p = [+2, -2]
            if (i ≠ 1) & (i ≠ N)
                if (1<=i+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i+p, j))
                    push!(values, 1)
                end
            end
            if (j ≠ 1) & (j ≠ N)
                if (1<=j+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i, j+p))
                    push!(values, 1)
                end
            end
        end
        for p = [+1, -1]
            if (i ≠ 1) & (i ≠ N)
                if (1<=i+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i+p, j))
                    push!(values, -16)
                end
            else
                if (1<=i+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i+p, j))
                    push!(values, -12)
                end
            end
            if (j ≠ 1) & (j ≠ N)
                if (1<=j+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i, j+p))
                    push!(values, -16)
                end
            else
                if (1<=j+p<=N)
                    push!(rows, lex(i, j))
                    push!(columns, lex(i, j+p))
                    push!(values, -12)
                end
            end
        end
        if (2<=i<=(N-1)) & (2<=j<=(N-1))
            push!(rows, lex(i, j))
            push!(columns, lex(i, j))
            push!(values, +60)
        else
            push!(rows, lex(i, j))
            push!(columns, lex(i, j))
            if ((i == 1) & (j == 1)) || ((i == N) & (j == N)) || ((i == N) & (j == 1)) || ((i == 1) & (j == N))
                push!(values, +48)
            else
                push!(values, +54)
            end
        end
    end
    A = sparse(rows, columns, values)
end

function fourth_Laplace_2D_2(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    A = spdiagm(0 => 30*ones(N), 1 => -16*ones(N-1), -1 => -16*ones(N-1), 2 => ones(N-2), -2 => ones(N-2))
    A[1, 1] = 15; A[1, 2] = 4; A[1, 3] = -14; A[1, 4] = 6; A[1, 5] = -1
    A[end, end] = 15; A[end, end-1] = 4; A[end, end-2] = -14; A[end, end-3] = 6; ; A[end, end-4] = -1
    D = spdiagm(0 => ones(N))
    A = kron(A, D) + kron(D, A)
end

function naive_Biharmonic(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    # A = spdiagm(0 => 6*ones(N), -1 => -4*ones(N-1), 1 => -4*ones(N-1), 2 => ones(N-2), -2 => ones(N-2))
    A = spdiagm(0 => (28/3)*ones(N), -1 => (-13/2)*ones(N-1), 1 => (-13/2)*ones(N-1), 2 => 2*ones(N-2), -2 => 2*ones(N-2), 3 => (-1/6)*ones(N-3), -3 => (-1/6)*ones(N-3))
    return A
end

function Biharmonic_equation(L::Integer, BCs::Array{String, 1}, parameters::Array{<:Real, 1})
    N = 2^L-1
    rows, columns, values = [], [], Array{Float64,1}(undef, 0)
    for i = 1:N
        for p = [+1, -1]
            if (1<=i+p<=N)
                push!(rows, i)
                push!(columns, i+p)
                push!(values, -8)
            end
        end
        for p = [+2, -2]
            if (1<=i+p<=N)
                push!(rows, i)
                push!(columns, i+p)
                push!(values, +1)
            end
        end
        if (2<=i<=(N-1))
            push!(rows, i)
            push!(columns, i)
            push!(values, +10)
        else
            push!(rows, i)
            push!(columns, i)
            push!(values, +11)
        end
    end
    A = sparse(rows, columns, values)
end

function test_Biharmonic()
    errors, H = [], []
    for L = 2:6
        N = 2^L-1
        h = 2.0^(-L)
        A = Biharmonic_equation_2D(L, ["D", "D"], [10.0, ])/h^4
        rhs = zeros(N^2)
        x = LinRange(h, 1-h, N)
        lex(i, j) = i + N*(j-1)
        u_exact = Array{Float64, 1}(undef, N^2)
        for (j, i) = collect(Iterators.product(1:N, 1:N))
            u_exact[lex(i, j)] = (x[i]*x[j]*(1-x[i])*(1-x[j]))^2
            rhs[lex(i, j)] = 8*(6*x[i]^2-6*x[i]+1)*(6*x[j]^2-6*x[j]+1)+24*((x[i]*(1-x[i]))^2+(x[j]*(1-x[j]))^2)
        end
        push!(H, h)
        push!(errors, norm(inv(Array(A))*rhs - u_exact, Inf))
    end
    errors = log10.(errors)
    H = log10.(H)
    return (errors[1] - errors[end])/(H[1] - H[end])
end

function test_Laplace()
    errors, H = [], []
    for L = 3:6
        N = 2^L-1
        h = 2.0^(-L)
        A = fourth_Laplace_2D_2(L, ["D", "D"], [10.0, ])/(12*h^2)
        rhs = zeros(N^2)
        x = LinRange(h, 1-h, N)
        lex(i, j) = i + N*(j-1)
        u_exact = Array{Float64, 1}(undef, N^2)
        for (j, i) = collect(Iterators.product(1:N, 1:N))
            u_exact[lex(i, j)] = sin(π*x[i])*sin(π*x[j])
            rhs[lex(i, j)] = 2*π^2*sin(π*x[i])*sin(π*x[j])
        end
        push!(H, h)
        push!(errors, norm(inv(Array(A))*rhs - u_exact, Inf))
    end
    errors = log10.(errors)
    H = log10.(H)
    return (errors[1] - errors[end])/(H[1] - H[end])
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
using Eirene

# Translate kwarg names.
function _eirene(data; dim_max=1, threshold=nothing, model)
    if isnothing(threshold)
        return eirene(data; maxdim=dim_max, model=model)
    else
        return eirene(data; maxdim=dim_max, maxrad=threshold, model=model)
    end
end

function eirene_benchmark(b::Benchmark)
    if b.data isa Vector
        data = Ripserer.to_matrix(b.data)
        model = "pc"
    else
        data = b.data
        model = "vr"
    end
    return @benchmarkable _eirene(
        $data; $(b.kwargs)..., threshold=$(b.threshold), model=$model
    ) samples=5 seconds=1000 gcsample=true
end


text:
#=== elastic ===#
#=== write Recording Components ===#
function wCT(fid::IOStream,OptCpnt::Array{String},nwf::nelwf2d)

    SVx = nwf.BDnvx[1]*nwf.BDnvx[2];
    SVz = nwf.BDnvz[1]*nwf.BDnvz[2];
    STxx = nwf.BDntpp[1]*nwf.BDntpp[2];
    STzz = nwf.BDntpp[1]*nwf.BDntpp[2];
    STxz = nwf.BDntxz[1]*nwf.BDntxz[2];
    #SWy = nwf.N_BDWy[1]*nwf.N_BDWy[2];

    SizeWFLib = Dict(
    "V"=>SVx+SVz,
    "Vx"=>SVx,
    "Vz"=>SVz,
    "Txx"=>STxx,
    "Tzz"=>STzz,
    "Tii"=>STxx+STzz,
    "Txz"=>STxz,
    #"Wy"=>SWy,
    "full"=>SVx+SVz+STxx+STzz+STxz)#+SWy)

    if "full" in OptCpnt
        Opt = Dict("full"=>"[Vx;Vz;Txx;Tzz;Txz]\n")
        write(fid,Opt["full"])
        NumCpnt = 5;
        SizeWFOpt = SizeWFLib["full"];
    else
        NumCpnt = 0
        SizeWFOpt = 0
        Opt = Dict(
        "V"=>["Vx","Vz"],
        "Tii"=>["Txx","Tzz"])

        for Cpnt in ["V","Tii"]
            if Cpnt in OptCpnt
                filter!(x->!(x in Opt[Cpnt]),OptCpnt)
            end
        end

        wOpt = Dict(
        "V"=>"Vx;Vz",
        "Tii"=>"Txx;Tzz")

        write(fid,"[")
        NCpnt = length(OptCpnt)-1==0?1:length(OptCpnt)-1
        for nCpnt in 1:NCpnt
            if OptCpnt[nCpnt] in ["V","Tii"]
                write(fid, wOpt[OptCpnt[nCpnt]])
                NumCpnt = NumCpnt + 2
            else write(fid, OptCpnt[nCpnt])
                NumCpnt = NumCpnt + 1
            end
            if length(OptCpnt) > 1
                write(fid,";")
            end
            SizeWFOpt = SizeWFOpt + SizeWFLib[OptCpnt[nCpnt]]
        end
        if length(OptCpnt) > 1
            if OptCpnt[length(OptCpnt)] in ["V","Tii"]
                write(fid, wOpt[OptCpnt[length(OptCpnt)]])
                NumCpnt = NumCpnt + 2
            else write(fid, OptCpnt[length(OptCpnt)])
                NumCpnt = NumCpnt + 1
            end
            SizeWFOpt = SizeWFOpt + SizeWFLib[OptCpnt[length(OptCpnt)]]
        end
        write(fid,"]\n")
    end
    return NumCpnt, SizeWFOpt
end

#=== acoustic ===#
#=== write Recording Components ===#
function wCT(fid::IOStream,OptCpnt::Array{String},nwf::nacwf2d)

    SVx = nwf.BDnvx[1]*nwf.BDnvx[2];
    SVz = nwf.BDnvz[1]*nwf.BDnvz[2];
    STxx = nwf.BDntpp[1]*nwf.BDntpp[2];
    STzz = nwf.BDntpp[1]*nwf.BDntpp[2];
    #SWy = nwf.N_BDWy[1]*nwf.N_BDWy[2];

    SizeWFLib = Dict(
    "V"=>SVx+SVz,
    "Vx"=>SVx,
    "Vz"=>SVz,
    "Txx"=>STxx,
    "Tzz"=>STzz,
    "Tii"=>STxx+STzz,
    #"Wy"=>SWy,
    "full"=>SVx+SVz+STxx+STzz)#+SWy)

    if "full" in OptCpnt
        Opt = Dict("full"=>"[Vx;Vz;Txx;Tzz]\n")
        write(fid,Opt["full"])
        NumCpnt = 4;
        SizeWFOpt = SizeWFLib["full"];
    else
        NumCpnt = 0
        SizeWFOpt = 0
        Opt = Dict(
        "V"=>["Vx","Vz"],
        "Tii"=>["Txx","Tzz"])

        for Cpnt in ["V","Tii"]
            if Cpnt in OptCpnt
                filter!(x->!(x in Opt[Cpnt]),OptCpnt)
            end
        end

        wOpt = Dict(
        "V"=>"Vx;Vz",
        "Tii"=>"Txx;Tzz")

        write(fid,"[")
        NCpnt = length(OptCpnt)-1==0?1:length(OptCpnt)-1
        for nCpnt in 1:NCpnt
            if OptCpnt[nCpnt] in ["V","Tii"]
                write(fid, wOpt[OptCpnt[nCpnt]])
                NumCpnt = NumCpnt + 2
            else write(fid, OptCpnt[nCpnt])
                NumCpnt = NumCpnt + 1
            end
            if length(OptCpnt) > 1
                write(fid,";")
            end
            SizeWFOpt = SizeWFOpt + SizeWFLib[OptCpnt[nCpnt]]
        end
        if length(OptCpnt) > 1
            if OptCpnt[length(OptCpnt)] in ["V","Tii"]
                write(fid, wOpt[OptCpnt[length(OptCpnt)]])
                NumCpnt = NumCpnt + 2
            else write(fid, OptCpnt[length(OptCpnt)])
                NumCpnt = NumCpnt + 1
            end
            SizeWFOpt = SizeWFOpt + SizeWFLib[OptCpnt[length(OptCpnt)]]
        end
        write(fid,"]\n")
    end
    return NumCpnt, SizeWFOpt
end

#=== write Recording Parameters ===#
function wRPT(fid::IOStream,nwf::Union{nelwf2d,nacwf2d},rec::receiver,
    OptCpnt::Array{String},medium::Union{elastic2d,acoustic2d})
    write(fid,"recinfo.bin: ") #   Recordings information file
    write(fid,"Simulated recordings information file, structured as\n")
    write(fid,"[Total Receiver Number; Total Time Sampling; Receiver Location Coordinates [Nz Nx] ]\n\n")

    write(fid,"rec.bin: ") #   Recordings data file
    write(fid,"recordings data file, structured as ")

    NumCpnt,SizeWFOpt = wCT(fid,OptCpnt,nwf)
    write(fid,"File Size $(round(rec.nr*NumCpnt*medium.nT*8/1024/1024,2))MB\n\n")
end

#=== write WaveField Parameters ===#
function wWFPT(fid::IOStream,nwf::Union{nelwf2d,nacwf2d},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}},
    OptCpnt::Array{String})
    write(fid,"wfinfo.bin: ") # Wavefields.bin
    write(fid,"Simulated wavefields information file, structured as\n")
    write(fid,"[Tpp Model Size(PML) [NZ NX]; Total Time Sampling Number; Corresponding Sampling Time ]\n\n")

    write(fid,"wf.bin: ")
    write(fid,"Simulated wavefields data file, structured as ")
    NumCpnt,SizeWFOpt = wCT(fid,OptCpnt,nwf)
    datasize = round(SizeWFOpt*length(Slices)*8/1024/1024/1024,2)
    write(fid,"File Size $(datasize)GB\n\n")
end

#=== write medium parameters ===#
function wOPT(fid::IOStream)
    wSPT(fid)
    write(fid,"No Recordings\n") # Recordings
    write(fid,"No Wavefields\n\n") # Wavefields
end

#=== write medium and recording parameters ===#
function wOPT(
    fid::IOStream,
    model::Union{elmod2d,acmod2d},
    rec::receiver,
    OptCpnt::Array{String})

    wSPT(fid)
    wRPT(fid,model.nwf,rec,OptCpnt,model.medium)
    write(fid,"No Wavefields\n\n") # Wavefields
end

#=== write medium and wavefield parameters ===#
function wOPT(
    fid::IOStream,
    model::Union{elmod2d,acmod2d},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}},
    OptCpnt::Array{String})

    wSPT(fid)
    write(fid,"No Recordings\n") # Recordings
    wWFPT(fid,model.nwf,Slices,OptCpnt)
end

#=== write medium, recording and wavefield parameters ===#
function wOPT(
    fid::IOStream,
    model::Union{elmod2d,acmod2d},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}},
    rec::receiver,
    OptCpnt::Array{String})

    wSPT(fid)
    wRPT(fid,model.nwf,rec.nr,OptCpnt,medium)
    wWFPT(fid,model.nwf,Slices,OptCpnt)
end


text:
# Pkg.build("QuantEcon")
import Dolo

path = Dolo.pkg_path
using AxisArrays

filename = joinpath(path,"examples","models","rbc_mc.yaml")
model = Dolo.yaml_import(filename)
@time dr = Dolo.time_iteration(model, verbose=true, maxit=10000)
@time dr = Dolo.improved_time_iteration(model, verbose=true, maxit=10000)

N=10
T= 50
Index_mc = Dolo.simulate(model.exogenous, N, T, 1)

s0=model.calibration[:states]
s0=repeat(s0,N)
typeof(s0)<:AbstractVector
Index_mc.data

dp_process=model.exogenous

sim = Dolo.simulate(model, dr.dr, Index_mc, model.exogenous)
sim = Dolo.simulate(model, dr.dr; i0=2)
sim = Dolo.simulate(model, dr.dr)

Tvalues = linspace(1, sim[Axis{:T}][end], sim[Axis{:T}][end])

a_sort = sort(collect(sim[Axis{:V}(:k)]),1)
a_median = a_sort[Int(size(sim,1)*0.5),:]

import PyPlot
plt = PyPlot;
fig = plt.figure("Markov Chaine")
for i in 1:sim[Axis{:N}][end]
    plt.plot(Tvalues, sim[Axis{:N}(i), Axis{:V}(:k)], color="red", alpha = .05)
end
plt.title("Simulations");
plt.xlabel("Horizon");
plt.ylabel("Capital");

plt.plot(Tvalues, a_median, marker="o")

Shocks_mc= Dolo.simulate(model.exogenous, 10,50, 1; return_indexes=false)

s0=model.calibration[:states]

# Check it won't work with values
# Dolo.simulate(model, dr, s0, Shocks_mc)
# Check

# Check it works without providing States
sim_mc2 = Dolo.simulate(model, dr, Index_mc, model.exogenous)
sim_mc2==sim

# Check it works without providing a driving process
sim_mc3 = Dolo.simulate(model, dr, model.exogenous; N=10, T=50, m0=2)

# Check it works without providing a driving process
sim_mc4 = Dolo.simulate(model, dr, s0, model.exogenous)
Tvalues = linspace(1, sim_mc4[Axis{:T}][end], sim_mc4[Axis{:T}][end])

import PyPlot
plt = PyPlot;
fig = plt.figure("Markov Chaine check")
for i in 1:sim_mc4[Axis{:N}][end]
    plt.plot(Tvalues, sim_mc4[Axis{:N}(i), Axis{:V}(:k)], color="red", alpha = .05)
end
plt.title("Simulations");
plt.xlabel("Horizon");
plt.ylabel("Capital");

################################################################################
# Tabulation
s0=model.calibration[:states]
m0=1

Dolo.tabulate(model, dr, :k, s0, m0)
bounds = [4.0, 10.0]

Dolo.tabulate(model, dr, :k, bounds, model.calibration[:states], m0)

# Check Simulation work for AR1 process

####################################################################################

filename = joinpath(path,"examples","models","rbc.yaml")
model2 = Dolo.yaml_import(filename)
@time dr2 = Dolo.time_iteration(model2, verbose=true, maxit=10000)

driving_process = Dolo.simulate(model2.exogenous, N, T)
dprocess = model2.exogenous
sim_ar = Dolo.simulate(model2, dr2.dr, driving_process)
sim_ar = Dolo.simulate(model2, dr2.dr, dprocess; m0=[0.1], N=20)
sim_ar = Dolo.simulate(model2, dr2.dr)

Tvalues = linspace(1, sim_ar[Axis{:T}][end], sim_ar[Axis{:T}][end])

a_sort = sort(collect(sim_ar[Axis{:V}(:k)]),1)
a_median = a_sort[Int(size(sim_ar,1)*0.5),:]

import PyPlot
plt = PyPlot;
fig = plt.figure("Simulate a model with AR1")
for i in 1:sim_ar[Axis{:N}][end]
    plt.plot(Tvalues, sim_ar[Axis{:N}(i), Axis{:V}(:k)], color="red", alpha = .05)
end
plt.title("Simulations");
plt.xlabel("Horizon");
plt.ylabel("Capital");

plt.plot(Tvalues, a_median, marker="o")

sim_ar2 = Dolo.simulate(model2, dr2, driving_process)
sim_ar2==sim_ar

# Check it works without providing a driving process
s0=model2.calibration[:states]
sim_ar3 = Dolo.simulate(model2, dr2, s0)

# Check it works without providing a driving process
m0 = model2.calibration[:exogenous]
sim_ar4 = Dolo.simulate(model2, dr2, s0, m0; N=100)
Tvalues = linspace(1, sim_ar4[Axis{:T}][end], sim_ar4[Axis{:T}][end])

import PyPlot
plt = PyPlot;
fig = plt.figure("AR1 check")
for i in 1:sim_ar4[Axis{:N}][end]
    plt.plot(Tvalues, sim_ar4[Axis{:N}(i), Axis{:V}(:k)], color="red", alpha = .05)
end
plt.title("Simulations");
plt.xlabel("Horizon");
plt.ylabel("Capital");

################################################################################
# Tabulation
s0=model2.calibration[:states]
m0 = model2.calibration[:exogenous]

Tab= Dolo.tabulate(model2, dr2, :k, s0, m0)
Dolo.tabulate(model2, dr2, :k)

Tab[:n]

# Check Simulation work for AR1 process when discretized with MC

####################################################################################

filename = joinpath(path,"examples","models","rbc.yaml")
model3 = Dolo.yaml_import(filename)

n_states=5
mc_ar = Dolo.discretize_mc(model3.exogenous;N=n_states)

@time dr_armc = Dolo.time_iteration(model3, mc_ar, verbose=true, maxit=10000)

# If you solve a model with a continuous exogenous process using MC disretization, then to simulate you need first proved a driving process originated form this new MC. model.exogenous
# At this point there is no option in the code. If you do not specify any driving porcess(series of shocks), then Dolo takes model.exogenous to do it itself, but then it can't combine dr(i,s) with values of shocks. Shall it be added?
Index_mc = Dolo.simulate(mc_ar, N, T, 1)

s0=model3.calibration[:states]
sim_armc = Dolo.simulate(model3, dr_armc, s0, Index_mc, mc_ar)
sim_armc_2 = Dolo.simulate(model3, dr_armc, s0, mc_ar)
sim_armc_3 = Dolo.simulate(model3, dr_armc, mc_ar)


text:
function test_gradient(model, parameters; rtol = 1e-10, atol = 0)
    true_grad = FiniteDiff.finite_difference_gradient(x -> objective!(model, x)[1], parameters)

    # F and G
    gradient(model) .= 0
    model(parameters, true, true, false)
    correct1 = isapprox(gradient(model), true_grad; rtol = rtol, atol = atol)

    # only G
    gradient(model) .= 0
    model(parameters, false, true, false)
    correct2 = isapprox(gradient(model), true_grad; rtol = rtol, atol = atol)

    return correct1 & correct2
end

function test_hessian(model, parameters; rtol = 1e-4, atol = 0)
    true_hessian = FiniteDiff.finite_difference_hessian(x -> objective!(model, x)[1], parameters)

    # F and H
    hessian(model) .= 0
    model(parameters, true, false, true)
    correct1 = isapprox(hessian(model), true_hessian; rtol = rtol, atol = atol)

    # G and H
    hessian(model) .= 0
    model(parameters, false, true, true)
    correct2 = isapprox(hessian(model), true_hessian; rtol = rtol, atol = atol)

    # only H
    hessian(model) .= 0
    model(parameters, false, false, true)
    correct3 = isapprox(hessian(model), true_hessian; rtol = rtol, atol = atol)

    return correct1 & correct2 & correct3
end

fitmeasure_names_ml = Dict(
    :AIC => "aic",
    :BIC => "bic",
    :df => "df",
    :χ² => "chisq",
    :p_value => "pvalue",
    :n_par => "npar",
    :RMSEA => "rmsea",
)

fitmeasure_names_ls = Dict(
    :df => "df",
    :χ² => "chisq",
    :p_value => "pvalue",
    :n_par => "npar",
    :RMSEA => "rmsea",
)

function test_fitmeasures(measures, measures_lav; rtol = 1e-4, atol = 0, fitmeasure_names = fitmeasure_names_ml)
    correct = []
    for key in keys(fitmeasure_names)
        measure = measures[key]
        measure_lav = measures_lav.x[measures_lav[:, 1] .==  fitmeasure_names[key]][1]
        push!(correct, isapprox(measure, measure_lav; rtol = rtol, atol = atol))
    end
    return correct
end

function compare_estimates(ens_partable::EnsembleParameterTable, partable_lav;
        rtol = 1e-10, atol = 0, col = :estimate, lav_col = :est,
        lav_groups)

    correct = []

    for key in keys(ens_partable.tables)

        group = lav_groups[key]
        partable = ens_partable.tables[key]

        for i in findall(partable.columns[:free])

            from = partable.columns[:from][i]
            to = partable.columns[:to][i]
            type = partable.columns[:parameter_type][i]
            estimate = partable.columns[col][i]

            if from == Symbol("1")

                lav_ind = findall(
                    (partable_lav.lhs .== String(to)) .&
                    (partable_lav.op .== "~1") .&
                    (partable_lav.group .== group))

                if length(lav_ind) == 0
                    throw(ErrorException("At least one parameter could not be found in the lavaan solution"))
                elseif length(lav_ind) > 1
                    throw(ErrorException("At least one parameter was found twice in the lavaan solution"))
                else
                    is_correct = isapprox(estimate, partable_lav[:, lav_col][lav_ind[1]]; rtol = rtol, atol = atol)
                    push!(correct, is_correct)
                end

            else

                if type == :↔
                    type = "~~"
                elseif type == :→
                    if (from ∈ partable.variables[:latent_vars]) & (to ∈ partable.variables[:observed_vars])
                        type = "=~"
                    else
                        type = "~"
                        from, to = to, from
                    end
                end

                lav_ind = findall(
                    (partable_lav.lhs .== String(from)) .&
                    (partable_lav.rhs .== String(to)) .&
                    (partable_lav.op .== type).&
                    (partable_lav.group .== group))

                if length(lav_ind) == 0
                    throw(ErrorException("At least one parameter could not be found in the lavaan solution"))
                elseif length(lav_ind) > 1
                    throw(ErrorException("At least one parameter was found twice in the lavaan solution"))
                else
                    is_correct = isapprox(estimate, partable_lav[:, lav_col][lav_ind[1]]; rtol = rtol, atol = atol)
                    push!(correct, is_correct)
                end

            end

        end

    end

    return all(correct)
end


text:
# This code by Tara Abrishami

export transitive_orientation

"""
`transitive_orientation(G)` finds a transitive orientation of the
simple graph `G`. The result is a `SimpleDigraph`. An error is raised
if `G` does not have a transitive orientation.
"""
function transitive_orientation(G::SimpleGraph)
  err_msg = "This graph does not have a transitive orientation"
  vertices = deepcopy(vlist(G))
  edges = deepcopy(elist(G))
  V = vertex_type(G)
  diredges = Tuple{V,V}[]
  D = SimpleDigraph{V}()
  while length(diredges) != 0 || length(edges) != 0
    if length(diredges) == 0
      e = popfirst!(edges)
      add!(D, e[1], e[2])
    else
      e = popfirst!(diredges)
    end
    v1 = e[1]
    v2 = e[2]
    vs = V[]
    append!(vs,G[v1])
    append!(vs,G[v2])
    for v in vs
      if has(G, v1, v) && !has(D, v1, v) && !has(G, v, v2)
        if has(D, v, v1)
          error(err_msg)
        end
        add!(D, v1, v)
        edg = (v1, v)
        swapedg = (v, v1)
        pushfirst!(diredges, edg)
        edges = filter(x -> x != edg && x != swapedg,edges)
      elseif has(G, v2, v) && !has(D, v, v2) && !has(G, v, v1)
        if has(D, v2, v)
          error(err_msg)
        end
        add!(D, v, v2)
        edg = (v, v2)
        swapedg = (v2, v)
        pushfirst!(diredges, edg)
        edges = filter(x -> x != edg && x != swapedg,edges)
      end
    end
  end
  if length(elist(D)) == length(elist(G))
    return D
  else
    error(err_msg)
  end
end

export num_trans_orientations

"""
`num_trans_orientations(G)` returns the number of transitive
orientations of the graph `G`.
"""
function num_trans_orientations(G2::SimpleGraph)
  if cache_check(G2,:num_trans_orientations)
    return cache_recall(G2,:num_trans_orientations)
  end
  G = deepcopy(G2)
  V = vertex_type(G)
  col = Dict{Tuple{V,V}, Int}()
  try
    col = makeColorClass(G)
  catch
    return 0
  end
  multiplexes = BigInt[]
  while length(elist(G)) != 0
    makeSimplex!(G, multiplexes, col)
  end
  ans = 1
  for m in multiplexes
    ans = ans*factorial(m)
  end
  cache_save(G2,:num_trans_orientations,ans)
  return ans
end

function makeSimplex!(G::SimpleGraph, multiplexes::Array, col::Dict)
  edge = elist(G)[1]
  vert1 = edge[1]
  vert2 = edge[2]
  colorSet = Set()
  T = vertex_type(G)
  simp = Set{T}()
  push!(simp, vert1)
  push!(simp, vert2)
  push!(colorSet, col[edge])
  for v in neighbors(G, vert1)
    colorSetTest = deepcopy(colorSet)
    isEdge = true
    for vo in simp
      if has(G, vo, v) && (((haskey(col, (v, vo)) && !(col[(v, vo)] in colorSetTest)) || (haskey(col, (vo,v)) &&!(col[(vo, v)] in colorSetTest))))
        if haskey(col, (v, vo))
          push!(colorSetTest, col[(v, vo)])
        else
          push!(colorSetTest, col[(vo, v)])
        end
      else
        isEdge = false
        break
      end
    end
    if isEdge == true
      push!(simp, v)
      colorSet = deepcopy(colorSetTest)
    end
  end
  pushfirst!(multiplexes, length(simp))
  for e in elist(G)
    if col[e] in colorSet
      delete!(G, e[1], e[2])
    end
  end
end

function makeColorClass(G1::SimpleGraph)
  err_msg = "error"
  G = deepcopy(G1)
  vertices = deepcopy(vlist(G))
  edge = deepcopy(elist(G))
  V = vertex_type(G)
  diredge = Tuple{V,V}[]
  D = SimpleDigraph{V}()
  E = SimpleGraph{V}()
  classes = SimpleGraph{V}[]
  while length(diredge) != 0 || length(edge) != 0
    if length(diredge) == 0
      if length(elist(E)) != 0
        pushfirst!(classes,E)
        E1 = SimpleGraph{V}()
        E = E1
      end
      e = popfirst!(edge)
      add!(D, e[1], e[2])
      add!(E, e[1], e[2])
    else
      e = popfirst!(diredge)
    end
    v1 = e[1]
    v2 = e[2]
    vs = V[]
    append!(vs,G[v1])
    append!(vs,G[v2])
    for v in vs
      if has(G, v1, v) && !has(D, v1, v) && !has(G, v, v2)
        if has(D, v, v1)
          error(err_msg)
        end
        add!(D, v1, v)
        add!(E, v1, v)
        edg = (v1, v)
        swapedg = (v, v1)
        pushfirst!(diredge, edg)
        edge = filter(x -> x != edg && x != swapedg,edge)
      elseif has(G, v2, v) && !has(D, v, v2) && !has(G, v, v1)
        if has(D, v2, v)
          error(err_msg)
        end
        add!(D, v, v2)
        add!(E, v, v2)
        edg = (v, v2)
        swapedg = (v2, v)
        pushfirst!(diredge, edg)
        edge = filter(x -> x != edg && x != swapedg,edge)
      end
    end
  end
  if length(elist(E)) != 0
    pushfirst!(classes,E)
  end
  col = Dict{Tuple{V,V}, Int}()
  i = 1
  for q in classes
    for x in elist(q)
      col[x] = i
    end
    i = i + 1
  end
  return col
end


text:
# compute the dynamic change in mean over time
function dμ_dt(network::DiffEqBase.AbstractReactionNetwork, momentstruct::MomentStruct; algebra=sympy())
    species_dict = speciesmap(network)
    dμ = repeat([zero(algebra.Sym)], momentstruct.num_species)
    for reaction in network.reactions
        for reactant in union(getfield.(reaction.products, :reactant), getfield.(reaction.substrates, :reactant))
            propensity = mass_rate_MA(reaction; algebra=algebra)
            dμ_i = algebra.expand(taylor_expand(propensity, momentstruct; algebra=algebra))
            dμ_i *= DiffEqBiological.get_stoch_diff(reaction, reactant)
            dμ[species_dict[reactant]] += dμ_i
        end
    end
    dμ
end


text:
using CooperationPolarization
using BenchmarkTools
using Profile
using PyPlot: plt, savefig
using Dates
using Random

const b = 1.0 # benefit to cooperating
const c = 0.2 # cost to cooperating

const N = 10 # population size
const M = 3 # number of sets
const K = 2 # number of issues to care about
const P = 2 # number of parties

const β = 0.000    # selection strength
const u = 0.001    # strategy mutation rate
const v = 0.001    # set mutation rate
const p = 0.75     # party bias (0 = low, 1 = high)
const ϵ = 1.       # damping factor (0 = no bias in mutation, 1 = max bias (= p) in mutation)
# note: α = 0 automatically

# simulation parameters
const verbose       = false
const verbose_track = false
const generations   = 10000000
const testtype      = "test_neutral"

# plotting parameters
const gen_skip  = 1 # max( div(generations, 10000), 1)
const gen_start = max( 1, generations-1000000 )
const gen_end   = generations

const start_time = Dates.format(now(),"yymmdd-HHMMSS.sss")

println("parameters defined")

sets      = random_sets(N, M, K, P)
game      = Game(b, c, β, u, v, p, ϵ)
pop       = Population(sets, game, verbose)
n_tracker = NeutralTracker(pop, generations, verbose_track)

# test
for i in 1:N pop.sets.h[i,:] = shuffle( vcat( zeros(Int64,M-K), rand([-1,1],K) ) ) end

println("Payoff matrix: \t $(game.A)")
println("Payoffs:\t $(pop.payoffs)")
println("Verbose?:\t $(pop.verbose)")
println("Initial:\t strategies_changed = $(pop.strategies_changed), sets_changed = $(pop.sets_changed)")

@time for gen::Int64 in 1:generations
    if generations > 10 && (gen-1) % (generations ÷ 10) == 0  # set up to print 100 generations
        println("\tinitiating generation $(n_tracker.pop.generation)")
    end
    # update population
    evolve_neutral!(n_tracker.pop, 1) # evolve without selection
    track_neutral!(n_tracker, gen)

    if gen == generations
        println("replicate complete, ran $(generations) generations:\n")
        println("\tnum_imitating:\t\t\t$(n_tracker.pop.num_imitating)")
        println("\tnum_not_imitating:\t\t$(n_tracker.pop.num_not_imitating)")
        println("\tnum_strategies_changed:  \t$(n_tracker.pop.num_strategies_changed)")
        println("\t\tnum_random_strategies: \t$(n_tracker.pop.num_random_strategies)")
        println("\t\t\tsa_flipped: \t$(n_tracker.pop.num_random_sa_flipped)")
        println("\t\t\tsd_flipped: \t$(n_tracker.pop.num_random_sd_flipped)")
        println("\tnum_sets_changed: \t\t$(n_tracker.pop.num_sets_changed)")
        println("\t\tnum_random_opinions: \t$(n_tracker.pop.num_random_opinions)")
        println("\t\t\tnum_flipped: \t$(n_tracker.pop.num_random_opinions_flipped)")
        println("\t\tnum_biased_opinions: \t$(n_tracker.pop.num_biased_opinions)")
        println("\t\t\tnum_flipped: \t$(n_tracker.pop.num_biased_opinions_flipped)")
        println("\tnum_both_changed: \t\t$(n_tracker.pop.num_both_changed)\n")
    end
end

# plotting
# fig = plotter( n_tracker, gen_start, gen_skip, gen_end )

fig, axs = plt.subplots(4, 4, figsize = (15,8))

plot_total_opn( n_tracker.pop, n_tracker.ens_y_running, axs[1], "y, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_z_running, axs[2], "z, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_g_running, axs[3], "g, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_h_running, axs[4], "h, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_sia_sjd_running, axs[9],  "<s_ia s_jd>, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_sia_sid_running, axs[10], "<s_ia s_id>, runnig avg", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_y, axs[5], "y", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_z, axs[6], "z", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_g, axs[7], "g", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_h, axs[8], "h", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_sia_sjd, axs[13], "<s_ia s_jd>", gen_start, gen_skip, gen_end)
plot_total_opn( n_tracker.pop, n_tracker.ens_sia_sid, axs[14], "<s_ia s_id>", gen_start, gen_skip, gen_end)

plot_strat_dist( n_tracker.pop, n_tracker.ens_strat_dist, axs[15], "abundance", gen_start, gen_skip, gen_end)
plot_strat_dist( n_tracker.pop, n_tracker.ens_strat_running, axs[11], "abundance", gen_start, gen_skip, gen_end)

plt.suptitle("Parameters: b=$b, c=$c, N=$N, M=$M, K=$K, β=$β, u=$u, v=$v, p=$p, ϵ=$ϵ", fontsize=12)
fig.tight_layout(rect=[0, 0.03, 1, 0.97])

display(fig)

# save file
# savefig("plots/$(testtype)_$(filename).pdf")


text:
module EmailScraper

using Cascadia
using Gumbo
using HTTP

include("TLDs.jl")
include("domain.jl")
include("email.jl")

"""
    scrape_url(url::String)

Scrapes a `url` and returns all email adresses found in links and in text.

# Arguments
- `url::String`: the desired url to scrape.

# Examples
```julia-repl
julia> scrape_url("julialang.org/about/help/")
2-element Vector{String}:
 "contact@julialang.org"
 "logan@julialang.org"
```
"""
function scrape_url(url::String)
    body = request_body(url)
    # links = get_links(body) # in the future this will scrap links
    emails = find_emails(body)
    return emails
end

export scrape_url

end # module


text:
# Test methods for reduced infinite variables
@testset "Reduced Variable" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 1 <= par1 <= 2)
    @infinite_parameter(m, 1 <= par2 <= 2)
    @infinite_variable(m, 0 <= inf1(par1, par2) <= 1, Int)
    @infinite_variable(m, inf2(par1, par2) == 1, Bin, start = 0)
    index = m.next_var_index + 1
    rvref1 = ReducedInfiniteVariableRef(m, index)
    m.reduced_info[index] = ReducedInfiniteInfo(inf1, Dict(1 => 1))
    rvref2 = ReducedInfiniteVariableRef(m, index + 1)
    m.reduced_info[index + 1] = ReducedInfiniteInfo(inf2, Dict(2 => 1))
    # test _reduced_info
    @testset "_reduced_info" begin
        @test InfiniteOpt._reduced_info(rvref1).infinite_variable_ref == inf1
        @test InfiniteOpt._reduced_info(rvref2).infinite_variable_ref == inf2
        @test InfiniteOpt._reduced_info(rvref1).eval_supports == Dict(1 => 1)
        @test InfiniteOpt._reduced_info(rvref2).eval_supports == Dict(2 => 1)
    end
    # test infinite_variable_ref
    @testset "infinite_variable_ref" begin
        @test infinite_variable_ref(rvref1) == inf1
        @test infinite_variable_ref(rvref2) == inf2
    end
    # test eval_supports
    @testset "eval_supports" begin
        @test eval_supports(rvref1) == Dict(1 => 1)
        @test eval_supports(rvref2) == Dict(2 => 1)
    end
    # test name
    @testset "JuMP.name" begin
        @test name(rvref1) == "inf1(1, par2)"
        @test name(rvref2) == "inf2(par1, 1)"
    end
    # parameter_refs (reduced infinite variable)
    @testset "parameter_refs (reduced infinite)" begin
        # test the references
        @test parameter_refs(rvref1) == (par2, )
        @test parameter_refs(rvref2) == (par1, )
    end
    # has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(rvref1)
        @test !has_lower_bound(rvref2)
    end
    # lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(rvref1) == 0
        @test_throws ErrorException lower_bound(rvref2)
    end
    # _lower_bound_index
    @testset "JuMP._lower_bound_index" begin
        @test JuMP._lower_bound_index(rvref1) == JuMP._lower_bound_index(inf1)
        @test_throws ErrorException JuMP._lower_bound_index(rvref2)
    end
    # LowerBoundRef
    @testset "JuMP.LowerBoundRef" begin
        @test LowerBoundRef(rvref1) == LowerBoundRef(inf1)
        @test_throws ErrorException LowerBoundRef(rvref2)
    end
    # has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(rvref1)
        @test !has_upper_bound(rvref2)
    end
    # upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(rvref1) == 1
        @test_throws ErrorException upper_bound(rvref2)
    end
    # _upper_bound_index
    @testset "JuMP._upper_bound_index" begin
        @test JuMP._upper_bound_index(rvref1) == JuMP._upper_bound_index(inf1)
        @test_throws ErrorException JuMP._upper_bound_index(rvref2)
    end
    # UpperBoundRef
    @testset "JuMP.UpperBoundRef" begin
        @test UpperBoundRef(rvref1) == UpperBoundRef(inf1)
        @test_throws ErrorException UpperBoundRef(rvref2)
    end
    # is_fixed
    @testset "JuMP.is_fixed" begin
        @test is_fixed(rvref2)
        @test !is_fixed(rvref1)
    end
    # fix_value
    @testset "JuMP.fix_value" begin
        @test fix_value(rvref2) == 1
        @test_throws ErrorException fix_value(rvref1)
    end
    # _fix_index
    @testset "JuMP._fix_index" begin
        @test JuMP._fix_index(rvref2) == JuMP._fix_index(inf2)
        @test_throws ErrorException JuMP._fix_index(rvref1)
    end
    # FixRef
    @testset "JuMP.FixRef" begin
        @test FixRef(rvref2) == FixRef(inf2)
        @test_throws ErrorException FixRef(rvref1)
    end
    # start_value
    @testset "JuMP.start_value" begin
        @test start_value(rvref2) == 0
    end
    # is_binary
    @testset "JuMP.is_binary" begin
        @test is_binary(rvref2)
        @test !is_binary(rvref1)
    end
    # _binary_index
    @testset "JuMP._binary_index" begin
        @test JuMP._binary_index(rvref2) == JuMP._binary_index(inf2)
        @test_throws ErrorException JuMP._binary_index(rvref1)
    end
    # BinaryRef
    @testset "JuMP.BinaryRef" begin
        @test BinaryRef(rvref2) == BinaryRef(inf2)
        @test_throws ErrorException BinaryRef(rvref1)
    end
    # is_integer
    @testset "JuMP.is_integer" begin
        @test is_integer(rvref1)
        @test !is_integer(rvref2)
    end
    # _integer_index
    @testset "JuMP._integer_index" begin
        @test JuMP._integer_index(rvref1) == JuMP._integer_index(inf1)
        @test_throws ErrorException JuMP._integer_index(rvref2)
    end
    # IntegerRef
    @testset "JuMP.IntegerRef" begin
        @test IntegerRef(rvref1) == IntegerRef(inf1)
        @test_throws ErrorException IntegerRef(rvref2)
    end
    # used_by_constraint
    @testset "used_by_constraint" begin
        # test not used
        @test !used_by_constraint(rvref1)
        # prepare use case
        m.reduced_to_constrs[JuMP.index(rvref1)] = [1]
        # test used
        @test used_by_constraint(rvref1)
    end
    # used_by_measure
    @testset "used_by_measure" begin
        # test not used
        @test !used_by_measure(rvref1)
        # prepare use case
        m.reduced_to_meas[JuMP.index(rvref1)] = [1]
        # test used
        @test used_by_measure(rvref1)
    end
    # test is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, rvref1)
        @test is_valid(m, rvref2)
        @test !is_valid(m, ReducedInfiniteVariableRef(m, 100))
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module HopfieldNetworks

using Lazy
using Quaternions
using Compat

include("base.jl")
include("learn_hebb.jl")
include("learn_projection.jl")
include("learn_local.jl")
include("associate.jl")
include("associate_cue.jl")
include("activatefuncs.jl")
include("overlap.jl")

end # module


text:
module MathTools
  # this module contains algebra and math functions
  # Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
  #
  #  MIT license.
  #
  # Permission is hereby granted, free of charge, to any person
  # obtaining a copy of this software and associated documentation files (the
  # "Software"), to deal in the Software without restriction, including without
  # limitation the rights to use, copy, modify, merge, publish, distribute,
  # sublicense, and/or sell copies of the Software, and to permit persons to whom
  # the Software is furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included in all
  # copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  # SOFTWARE.

  export numel
  export max

  # wrapper for the numel function
  function numel(matrix)
    numberOfElements = length(matrix)
    return numberOfElements
  end

  # wrapper for the max function
  function max(vector)
    maxOfVector = maximum(vector)
    return maxOfVector
  end

end #End of StringTools


text:
using BinaryBuilder

name = "LibUV"
version = v"2"

# Collection of sources required to build libuv
sources = [
    GitSource("https://github.com/JuliaLang/libuv.git",
              "67d1dc2ca2289f25a23b016a0bc427d45301ec05"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libuv/

# Touch some files so that the build system doesn't attempt to re-run `autoconf`:
touch -c aclocal.m4
touch -c Makefile.in
touch -c configure

# `--with-pic` isn't enough; we really really need -fPIC and -DPIC everywhere...
# everywhere, especially on FreeBSD
./configure --prefix=$prefix --host=$target --with-pic CFLAGS="${CFLAGS} -DPIC -fPIC" CXXFLAGS="${CXXFLAGS} -DPIC -fPIC"
make -j${nproc} V=1
make install
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libuv", :libuv),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
version = v"2.0.1"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
# ---
# title: 1626. Best Team With No Conflicts
# id: problem1626
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/best-team-with-no-conflicts/description/>
# hidden: true
# ---
#
# You are the manager of a basketball team. For the upcoming tournament, you
# want to choose the team with the highest overall score. The score of the team
# is the **sum** of scores of all the players in the team.
#
# However, the basketball team is not allowed to have **conflicts**. A
# **conflict** exists if a younger player has a **strictly higher** score than
# an older player. A conflict does **not** occur between players of the same
# age.
#
# Given two lists, `scores` and `ages`, where each `scores[i]` and `ages[i]`
# represents the score and age of the `ith` player, respectively, return _the
# highest overall score of all possible basketball teams_.
#
#
#
# **Example 1:**
#
#
#
#     Input: scores = [1,3,5,10,15], ages = [1,2,3,4,5]
#     Output: 34
#     Explanation:  You can choose all the players.
#
#
# **Example 2:**
#
#
#
#     Input: scores = [4,5,6,5], ages = [2,1,2,1]
#     Output: 16
#     Explanation:  It is best to choose the last 3 players. Notice that you are allowed to choose multiple people of the same age.
#
#
# **Example 3:**
#
#
#
#     Input: scores = [1,2,3,5], ages = [8,9,10,1]
#     Output: 6
#     Explanation:  It is best to choose the first 3 players.
#
#
#
#
# **Constraints:**
#
#   * `1 <= scores.length, ages.length <= 1000`
#   * `scores.length == ages.length`
#   * `1 <= scores[i] <= 106`
#   * `1 <= ages[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function mandelbrot_set(xmin,xmax,ymin,ymax,width,height,maxiter)
    r1 = float_range(xmin, xmax, width)
    r2 = float_range(ymin, ymax, height)

    cartesian =  Array{Float64, 3}(undef, 3, height, width)

    for y = 1:height, x = 1:width
	point = mandel(ComplexF64(r1[x], r2[y]), maxiter)
	if point == maxiter
		point = 0
	else
		point /= convert(Float64, maxiter)

	end

 	red_color, blue_color, green_color = 1.0, 1.0, 1.0
    	if point < 0.01
            red_color, blue_color, green_color = 1.0, 1.0, 1.0
    	else
            red_color = point * 0.80
            blue_color = point * 0.80
            green_color = point * 0.15
    	end

	cartesian[1,y,x] = red_color
	cartesian[2,y,x] = blue_color
	cartesian[3,y,x] = green_color
    end

    return cartesian
end

# the mandelbrot iteration function
function mandel(c, maxiter::Int64)
    z=0+0im
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

function float_range(xmin, xmax, width)
	r = Array{Float64, 1}(undef, width)

	delta = abs(xmax - xmin)/(width - 1)

	for i = 1:width
		r[i] = xmin + (i-1)*delta
	end

	return r
end

max_stars_repo_name: rabarar/juliaFractals

max_issues_repo_name: rabarar/juliaFractals

max_forks_repo_name: rabarar/juliaFractals


text:
using Documenter
using TcpInstruments

makedocs(
    sitename = "TcpInstruments",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
             "Home" => "index.md",
             "General" => "general.md",
             "Supported Instruments" => "instruments.md",
            ],
    modules = [TcpInstruments]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/Orchard-Ultrasound-Innovation/TcpInstruments.jl.git"
)