text:
"""
map_rows(f::Function, a::Array)

Maps the function onto each row of the passed array
"""
function map_rows(f::Function, a::Array)
    n = size(a,1)
    vals = zeros(n)
    for i in 1:n
        vals[i] = f(a[i,:])
    end
    return vals
end

"""
generate_map_rows_function(f::Function, a::Array)

Returns a function that applies the function-generating function 'f' to
    the rows of 'a', which can then be called with a SCALAR argument

Example, where size(dat) == (3,1000), and ts==(time span)
    F = generate_map_rows_function(
                (x)->CubicSplineInterpolation(ts,x), dat)
    F(0) # size()=(3,1); interpolation of the data at this time
"""
function generate_map_rows_function(f::Function, a::Array)
    n = size(a,1)
    func_list = []
    for i in 1:n
        push!(func_list, f(a[i,:]))
    end
    function F(t)
        return t .|> func_list
    end
    return F
end

"""
drop_all_1dims(a)

Drops dimensions that are equal to one in multidimensional arrays
    from: https://stackoverflow.com/questions/46289554/drop-julia-array-dimensions-of-length-1
"""
function drop_all_1dims(a)
    return dropdims(a, dims=Tuple(findall(size(a).==1)))
end

"""
Get contiguous true values, returning all start and stop times as well as
    the list of block lengths.
"""
function calc_contiguous_blocks(dat::Vector;
                                minimum_length=1)
    dat = Bool.(dat)
    all_starts = []
    all_ends = []
    current_length = 0
    for i in 1:length(dat)
        if dat[i]
            if current_length == 0
                push!(all_starts, i)
            end
            current_length += 1
        else
            if current_length > 0
                push!(all_ends, i-1)
                current_length = 0
            end
        end
    end
    if current_length > 0
        # Ended on 'true'
        push!(all_ends, length(dat))
    end

    block_lengths = all_ends .- all_starts .+ 1
    # Remove short ones
    long_enough = (block_lengths .>= minimum_length)
    all_starts = all_starts[long_enough]
    all_ends = all_ends[long_enough]
    block_lengths = block_lengths[long_enough]

    # Get the "good" indices
    true_indices = []
    for (s, e) in zip(all_starts, all_ends)
        true_indices = vcat(true_indices, s:e)
    end

    return (true_indices=true_indices,
            block_lengths=block_lengths,
            all_starts=all_starts,
            all_ends=all_ends)
end

export map_rows, generate_map_rows_function, drop_all_1dims,
        calc_contiguous_blocks


text:
using Dash

test_cases = Dict(
    "not-boolean" => Dict(
        "fail"=> true,
        "name"=> "simple \"not a boolean\" check",
        "component"=> dcc_input,
        "props"=> (debounce = 0,),
    ),
    "missing-required-nested-prop"=> Dict(
        "fail"=> true,
        "name"=> "missing required \"value\" inside options",
        "component"=> dcc_checklist,
        "props"=> (options = [Dict("label" => "hello")], value = ["test"]),
    ),
    "invalid-nested-prop"=> Dict(
        "fail"=> true,
        "name"=> "invalid nested prop",
        "component"=> dcc_checklist,
        "props"=> (options = [Dict("label"=> "hello", "value"=> true)], value = ["test"]),
    ),
    "invalid-arrayOf"=> Dict(
        "fail"=> true,
        "name"=> "invalid arrayOf",
        "component"=> dcc_checklist,
        "props"=> (options = "test", value = []),
    ),
    "invalid-oneOf"=> Dict(
        "fail"=> true,
        "name"=> "invalid oneOf",
        "component"=> dcc_input,
        "props"=> (type = "test",),
    ),
    "invalid-oneOfType"=> Dict(
        "fail"=> true,
        "name"=> "invalid oneOfType",
        "component"=> dcc_input,
        "props"=> (max = true,),
    ),
    "invalid-shape-5"=> Dict(
        "fail"=> true,
        "name"=> "invalid not required key",
        "component"=> dcc_dropdown,
        "props"=> (options = [Dict("label"=> "new york", "value"=> "ny", "typo"=> "asdf")],),
    ),
    "string-not-list"=> Dict(
        "fail"=> true,
        "name"=> "string-not-a-list",
        "component"=> dcc_checklist,
        "props"=> (options = [Dict("label"=> "hello", "value"=> "test")], value = "test"),
    ),
    "no-properties"=> Dict(
        "fail"=> false,
        "name"=> "no properties",
        "component"=> dcc_input,
        "props"=> (),
    ),
    "nested-children"=> Dict(
        "fail"=> true,
        "name"=> "nested children",
        "component"=> html_div,
        "props"=> (children = [[1]],),
    ),
    "deeply-nested-children"=> Dict(
        "fail"=> true,
        "name"=> "deeply nested children",
        "component"=> html_div,
        "props"=> (children = html_div([html_div([3, html_div([[10]])])]),),
    ),
    "dict"=> Dict(
        "fail"=> true,
        "name"=> "returning a dictionary",
        "component"=> html_div,
        "props"=> (children = Dict("hello" => "world"),),
    ),
    "allow-null-2"=> Dict(
        "fail"=> false,
        "name"=> "allow null as value",
        "component"=> dcc_dropdown,
        "props"=> (value = nothing,),
    ),
    "allow-null-3"=> Dict(
        "fail"=> false,
        "name"=> "allow null in properties",
        "component"=> dcc_input,
        "props"=> (value = nothing,),
    ),
    "allow-null-4"=> Dict(
        "fail"=> false,
        "name"=> "allow null in oneOfType",
        "component"=> dcc_store,
        "props"=> (id = "store", data = nothing),
    ),
    "long-property-string"=> Dict(
        "fail"=> false,
        "name"=> "long property string with id",
        "component"=> html_div,
        "props"=> (id = "pink div", style = repeat("color=> hotpink; ", 1000)),
    ),
    "multiple-wrong-values"=> Dict(
        "fail"=> true,
        "name"=> "multiple wrong props",
        "component"=> dcc_dropdown,
        "props"=> (id = "dropdown", value = 10, options = "asdf"),
    ),
    "boolean-html-properties"=> Dict(
        "fail"=> true,
        "name"=> "dont allow booleans for dom props",
        "component"=> html_div,
        "props"=> (contentEditable = true,),
    ),
    "allow-exact-with-optional-and-required-1"=> Dict(
        "fail"=> false,
        "name"=> "allow exact with optional and required keys",
        "component"=> dcc_dropdown,
        "props"=> (options = [Dict("label"=> "new york", "value"=> "ny", "disabled"=> false)],),
    ),
    "allow-exact-with-optional-and-required-2"=> Dict(
        "fail"=> false,
        "name"=> "allow exact with optional and required keys 2",
        "component"=> dcc_dropdown,
        "props"=> (options = [Dict("label"=> "new york", "value"=> "ny")],),
    )
)
app = dash()

app.layout = html_div() do
    html_div(id="content"),
    dcc_location(id="location")
end

callback!(app, Output("content","children"), Input("location","pathname")) do pathname
    if isnothing(pathname) || pathname == "/"
        return "Initial state"
    end
    test_case = test_cases[strip(pathname, '/')]
    return html_div(
        id="new-component", children=test_case["component"](;test_case["props"]...)
    )
end

run_server(app, debug = true, dev_tools_hot_reload = false)

ext: jl


text:
module Equations

    export Equation, ODE, IODE, PODE, SODE, VODE, DAE, HDAE, IDAE, PDAE
    export SDE, PSDE, SPSDE

    include("equations/equations.jl")

    include("equations/ode.jl")
    include("equations/iode.jl")
    include("equations/pode.jl")
    include("equations/sode.jl")
    include("equations/vode.jl")
    include("equations/dae.jl")
    include("equations/hdae.jl")
    include("equations/idae.jl")
    include("equations/pdae.jl")
    include("equations/sde.jl")
    include("equations/psde.jl")
    include("equations/spsde.jl")

end


text:
# --------------  Keyword arguments

# test this +++++

"""
	$(SIGNATURES)

Object to hold default values for keyword arguments.
"""
mutable struct KwArgs
	d :: Dict{Symbol, Any}
end

"""
	$(SIGNATURES)

Retrieve a default value by name.
"""
function default_value(name :: Symbol, kwS :: KwArgs)
	@assert has_default(name, kwS)  "$name not found"
	return kwS.d[name]
end

function has_default(name :: Symbol, kwS :: KwArgs)
	return haskey(kwS.d, name)
end

"""
	$(SIGNATURES)

Retrieve either the default value or the override in `kwargs`.
"""
function kw_arg(name :: Symbol, kwS :: KwArgs; kwargs...)
	if has_default(name, kwS)
		defaultVal = default_value(name, kwS);
		found = true;
	else
		defaultVal = nothing
        found = false;
	end

	if haskey(kwargs, name)
		found = true;
		outVal = kwargs[name];
	else
		outVal = defaultVal;
	end
	@assert found  "Argument $name not found."
	return outVal
end

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


text:
using Utils

using Chain
using CompositeStructs
import IterTools

const file = find_input(@__FILE__)
const input = file |> slurp

function bitvector(input::AbstractString)
    @assert length(input) % 2 == 0 "input has length $(length(input))"
    bytes = @chain input IterTools.partition(_, 2) map(c -> parse(UInt8, join(c), base=16), _)
    re = BitVector((byte & (1 << i)) != 0 for byte=bytes for i=7:-1:0)
    append!(re, [false, false, false])
    re
end
@assert bitvector("D2")[1:8] == [true, true, false, true, false, false, true, false]
@assert bitvector("FE")[1:8] == [true, true, true, true, true, true, true, false]
@assert bitvector("FED2")[1:16] == [true, true, true, true, true, true, true, false, true, true, false, true, false, false, true, false]

decimal(bits::AbstractArray{Bool}) = foldl((l,r) -> (l << 1) | r, bits, init=0)
@assert decimal([true,true,false]) == 6

struct Packet
    version::UInt8
    type::UInt8
    content::Union{Int,Vector{Packet}}
end

function decode(bv::BitVector, pos=1; limit=-1)
    found = Packet[]
    while pos <= length(bv) && limit != 0
        version = decimal(bv[pos:pos+2])
        pos += 3
        type = decimal(bv[pos:pos+2])
        pos += 3
        if type == 4
            num = 0
            while true
                chunk = decimal(bv[pos:pos+4])
                pos += 5
                num = (num << 4) | (chunk & 0b1111)
                ((chunk & 0b10000) == 0) && break
            end
            push!(found, Packet(version, type, num))
        else
            if bv[pos]
                len = decimal(bv[pos+1:pos+11])
                pos += 1 + 11
                inc, children = decode(bv[pos:end]; limit=len)
            else
                len = decimal(bv[pos+1:pos+15])
                pos += 1 + 15
                inc, children = decode(bv[pos:pos+len-1])
            end
            return pos + inc - 1, [Packet(version, type, children)]
        end
        if limit > 0
            limit -= 1
        end
    end
    pos, found
end
decode(string::AbstractString) = string |> bitvector |> decode
decode("620080001611562C8802118E34")
decode("620080001611562C8802118E34" |> bitvector)

decode("EE00D40C823060")
decode("38006F45291200")
decode("C0015000016115A2E0802F182340")

function sumversions(paket::Packet)
    re = paket.version
    if !isa(paket.content, Number)
        re += map(sumversions, paket.content) |> sum
    end
    re
end

sumversions(pakets::Vector{Packet}) = map(sumversions, pakets)

function part1(input)
    @chain input begin
        decode
        _[2]
        sumversions
        sum
        convert(Int, _)
    end
end

@assert part1("8A004A801A8002F478") == 16
@assert part1("620080001611562C8802118E34") == 12
@assert part1("C0015000016115A2E0802F182340") == 23
@assert part1("A0016C880162017C3686B18A3D4780") == 31

function part2(input)
end

assertequal(part2(input))


text:
# Rational Example

struct RationalNumber
    n::Integer
    d::Integer
    function RationalNumber(n, d)
        (n, d) = _normalize(n, d)
        return new(n, d)
    end
end

function _normalize(n::Integer, d::Integer)
    g = gcd(n, d)
    m = d < 0 ? -1 : 1
    (m * n ÷ g, m * d ÷ g)
end

Base.:+(r1::RationalNumber, r2::RationalNumber) = RationalNumber(r1.n+r2.n, r1.d+r2.d)
Base.:-(r1::RationalNumber, r2::RationalNumber) = RationalNumber(r1.n-r2.n, r1.d-r2.d)
Base.:*(r1::RationalNumber, r2::RationalNumber) = RationalNumber(r1.n*r2.n, r1.d*r2.d)
Base.:/(r1::RationalNumber, r2::RationalNumber) = RationalNumber(r1.n*r2.d, r1.d*r2.n)

str(r::RationalNumber) = "$(r.n)/$(r.d)"

# Demonstration

x = half = RationalNumber(5, 10)
println("One half can be expressed as $(half.n)/$(half.d)")

y = two_thirds = RationalNumber(2,3)

ystr = str(y)
println("two thirds is $(ystr)")
println("two thirds is $(str(y))")

println("x+y = $(str(x+y))")
println("x-y = $(str(x-y))")
println("x*y = $(str(x*y))")
println("x/y = $(str(x/y))")

max_stars_repo_name: AdamVStephen/learning-julia

max_issues_repo_name: AdamVStephen/learning-julia

max_forks_repo_name: AdamVStephen/learning-julia


text:
module TokamakSmallCylindrical

    import ElectromagneticFields.AxisymmetricTokamakCylindrical

    export pauli_particle_3d_pode, hamiltonian, toroidal_momentum

    AxisymmetricTokamakCylindrical.@code() # inject magnetic field code

    const xᵢ = [1.05, 0.0, 0.0]
    const qᵢ = from_cartesian(0, xᵢ)
    const vᵢ = DF̄(0, xᵢ) * [2.1E-3, 4.3E-4, 0.0]

    include("pauli_particle_3d.jl")

end


text:
using OrdinaryDiffEq, Test
using MultiScaleArrays

struct CellGenotype{B<:Float64} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
    Iam::Symbol
end

struct Host{T<:AbstractMultiScaleArray,B<:Float64} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

host = construct(Host,[
        CellGenotype(# Genotype present at start for symbiont 1
        [Float64(1)],
        :aSymbiontGenotype, # I am
        ),
        CellGenotype(# Genotype present at start for symbiont 1
        [Float64(1)],
        :aSymbiontGenotype, # I am
        )]
    )

evolve = function (du,u::Host,p,t)
  for i in 1:length(u)
    du[i] = 2*u[i]
  end
end

similar(host)
similar(host.nodes[1])
similar(host.nodes[1],Float64)
similar(host,Float64)

## Solve
tspan = (0.0,50)
prob = ODEProblem(evolve, host, tspan)
sol = solve(prob,Tsit5())

struct Host2{T<:AbstractMultiScaleArray,B<:Float64} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
    Iam::Symbol
end

host2 = construct(Host2,[
        CellGenotype(# Genotype present at start for symbiont 1
        [Float64(1)],
        :aSymbiontGenotype, # I am
        ),
        CellGenotype(# Genotype present at start for symbiont 1
        [Float64(1)],
        :aSymbiontGenotype, # I am
        )],Float64[],
        :my_cell_type
    )

evolve = function (du,u::Host2,p,t)
  for i in 1:length(u)
    du[i] = 2*u[i]
  end
end

similar(host2)
similar(host2.nodes[1])
similar(host2.nodes[1],Float64)
similar(host2,Float64)

## Solve
tspan = (0.0,50)
prob = ODEProblem(evolve, host2, tspan)
sol = solve(prob,Tsit5())


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/DynamicBoundsBase.jl
# Main module code.
#############################################################################

module DynamicBoundsBase

using DocStringExtensions, ElasticArrays
using ForwardDiff: Chunk, Dual, Partials, construct_seeds, single_seed
using SciMLBase: remake, AbstractODEProblem, ODEProblem, AbstractContinuousCallback,
                 solve, ODEFunction, wrapfun_iip
using DiffEqSensitivity: extract_local_sensitivities, ODEForwardSensitivityProblem

import Base.get

export HasStateBounds, HasConstantStateBounds, HasVariableStateBounds,
       HasUserJacobian, ConstantStateBounds, VariableStateBounds,
       PolyhedralConstraint, AbstractDERelaxProblem, AbstractRelaxProblemAttribute,
       AbstractDERelaxIntegrator, AbstractIntegratorAttribute, AbstractDynamicIndex
include(joinpath(@__DIR__, "problem.jl"))

export Nominal, Lower, Upper, Undefined, IntegratorName, Value, Gradient,
       ConstantParameterValue, Subgradient, Bound, Relaxation,IsNumeric,
       IsSolutionSet, TerminationStatus, ParameterValue, ParameterBound,
       SupportSet, TerminationStatusCode, relax!,
       integrate!, make, AbstractDERelaxIntegator, AbstractODERelaxIntegrator,
       AbstractBoundLoc, TimeIndex, COMPLETED, EMPTY, RELAXATION_NOT_CALLED,
       NUMERICAL_ERROR, INVALID_OPTION, OTHER_ERROR, LIMIT_EXCEEDED,
       ParameterNumber, StateNumber, SupportNumber, LocalSensitivityOn,
       LocalIntegrator, AttachedProblem

include(joinpath(@__DIR__, "integrator.jl"))

export set!, setall!, get, getall, getall!, make, supports
include(joinpath(@__DIR__, "attribute_handlers.jl"))

export IntegratorStates
include(joinpath(@__DIR__, "utilities.jl"))

export ODERelaxProb, AbstractODERelaxProblem, ODELocalIntegrator
include(joinpath(@__DIR__, "problems", "ODERelaxProb.jl"))

end


text:
# Tests for the conditional mutual information function

entropy3 = get_entropy(arr3)
entropy13 = get_entropy(arr1, arr3)
entropy23 = get_entropy(arr2, arr3)
entropy123 = get_entropy(arr1, arr2, arr3)
cmi = entropy13 + entropy23 - entropy123 - entropy3

# Test conditional mutual information
@test get_conditional_mutual_information(arr1, arr2, arr3) ≈ cmi
println("Conditional mutual information passed.")

# Test conditional mutual information formula
@test apply_conditional_mutual_information_formula(entropy13, entropy23, entropy123, entropy3) ≈ cmi
println("Conditional mutual information formula passed.")


text:
# Model definition
n = eternity2_generator.n
m = eternity2_generator.m

activation = Flux.leakyrelu
approximator_GNN = GeometricFlux.GraphConv(64 => 64, activation)
target_approximator_GNN = GeometricFlux.GraphConv(64 => 64, activation)
gnnlayers = 1

approximator_model = SeaPearl.FullFeaturedCPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(numInFeatures=>12, activation),GeometricFlux.GraphConv(12=>12, activation)
    ),
    nodeChain = Flux.Chain(
    Flux.Dense(12, 12, activation),
    ),
    outputChain = Flux.Chain(
        Flux.Dense(24, 32, activation),
        Flux.Dense(32, 1),
    ),
) #|> gpu
target_approximator_model = SeaPearl.FullFeaturedCPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(numInFeatures=>12, activation),
        GeometricFlux.GraphConv(12=>12, activation)
    ),
    nodeChain = Flux.Chain(
    Flux.Dense(12, 12, activation),
    ),
    outputChain = Flux.Chain(
        Flux.Dense(24, 32, activation),
        Flux.Dense(32, 1),
    ),
) #|> gpu

#rng = MersenneTwister(33)

agent = RL.Agent(
    policy = RL.QBasedPolicy(
        learner = RL.DQNLearner(
            approximator = RL.NeuralNetworkApproximator(
                model = approximator_model,
                optimizer = ADAM()
            ),
            target_approximator = RL.NeuralNetworkApproximator(
                model = target_approximator_model,
                optimizer = ADAM()
            ),
            loss_func = Flux.Losses.huber_loss,
            γ = 0.9f0,
            batch_size = 8, #32,
            update_horizon = 7,
            min_replay_history = 8,
            update_freq = 8,
            target_update_freq = 100,
            #rng = rng,
        ),
        explorer = RL.EpsilonGreedyExplorer(
            ϵ_stable = 0.01,
            #kind = :exp,
            decay_steps = nbEpisodes,
            step = 1,
        )
    ),
    trajectory = RL.CircularArraySLARTTrajectory(
        capacity = 200,
        state = SeaPearl.DefaultTrajectoryState[] => (),
        legal_actions_mask = Vector{Bool} => (n*m, ),
    )
)

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
# these do most of the work in Schur decomposition, SVD, etc.
using LinearAlgebra: Givens

function _floatmin2(::Type{T}) where {T}
    twopar = 2one(T)
    twopar^trunc(Integer,log(floatmin(T)/eps(T))/log(twopar)/twopar)
end
const d32fm2 = _floatmin2(Double32)
const d64fm2 = _floatmin2(Double64)
LinearAlgebra.floatmin2(::Type{Double32}) = d32fm2
LinearAlgebra.floatmin2(::Type{Double64}) = d64fm2

@noinline function lmul!(G::Givens{DT}, A::AbstractVecOrMat{DT}
                       ) where {DT <: Complex{DoubleFloat{T}}} where T
    has_offset_axes(A) && throw(ArgumentError("not implemented "
                                              * "for offset axes"))
    m, n = size(A, 1), size(A, 2)
    if G.i2 > m
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    a1 = Vector{DT}(undef,n)
    a2 = Vector{DT}(undef,n)
    @inbounds begin
        for i=1:n
            a2[i] = A[G.i2,i]
        end
        for i=1:n
            a1[i] = A[G.i1,i]
        end
        rmul!(view(A,G.i1,:),G.c)
        axpy!(G.s,a2,view(A,G.i1,:))
        rmul!(view(A,G.i2,:),G.c)
        axpy!(-conj(G.s),a1,view(A,G.i2,:))
    end
    return A
end

@noinline function lmul!(G::Givens{DoubleFloat{T}},
                         A::StridedMatrix{DoubleFloat{T}}) where T
    has_offset_axes(A) && throw(ArgumentError("not implemented "
                                              * "for offset axes"))
    m, n = size(A, 1), size(A, 2)
    if G.i2 > m
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    _lmul!(G,view(A,G.i1,:),view(A,G.i2,:),Vec{Npref,T})
end
function _lmul!(G::Givens{DoubleFloat{T}}, xv::StridedVector{DoubleFloat{T}},
                yv::StridedVector{DoubleFloat{T}},::Type{Vec{N,T}}) where {N,T}
    n = size(xv,1)
    nd,nr = divrem(n, N)
    ahi, alo = HILO(G.c)
    chi = Vec{N,T}(ahi)
    clo = Vec{N,T}(alo)
    ahi, alo = HILO(G.s)
    shi = Vec{N,T}(ahi)
    slo = Vec{N,T}(alo)
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N

            xhi = vgethi(xv,i0,Vec{N,T})
            xlo = vgetlo(xv,i0,Vec{N,T})
            yhi = vgethi(yv,i0,Vec{N,T})
            ylo = vgetlo(yv,i0,Vec{N,T})

            uhi, ulo = dfvmul(xhi, xlo, chi, clo)
            vhi, vlo = dfvmul(yhi, ylo, shi, slo)
            zhi, zlo = dfvadd(uhi, ulo, vhi, vlo)
            vputhilo!(xv,i0,zhi,zlo)

            vhi, vlo = dfvmul(xhi, xlo, shi, slo)
            uhi, ulo = dfvmul(yhi, ylo, chi, clo)
            zhi, zlo = dfvsub(uhi, ulo, vhi, vlo)
            vputhilo!(yv,i0,zhi,zlo)
        end
    end
    (nr == 0) && return yv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            t1 = xv[i]
            t2 = yv[i]
            xv[i] = G.c*xv[i] + G.s*t2
            yv[i] = G.c*yv[i] - G.s*t1 # conj(G.s) == G.s for Real
        end
    end
    nothing
end

#@inline function rmul!(A::AbstractMatrix{DT}, G::Givens{DT}
@noinline function rmul!(A::AbstractMatrix{Complex{DoubleFloat{T}}},
                         G::Givens{Complex{DoubleFloat{T}}}
                       ) where T
    has_offset_axes(A) && throw(ArgumentError("not implemented "
                                              * "for offset axes"))
    m, n = size(A, 1), size(A, 2)
    if G.i2 > n
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    a1 = Vector{Complex{DoubleFloat{T}}}(undef,m)
    a2 = Vector{Complex{DoubleFloat{T}}}(undef,m)
    @inbounds begin
        for i=1:m
            a1[i] = A[i,G.i1]
            a2[i] = A[i,G.i2]
        end
        rmul!(view(A,:,G.i1),G.c)
        rmul!(view(A,:,G.i2),G.c)
        axpy!(-conj(G.s),a2,view(A,:,G.i1))
        axpy!(G.s,a1,view(A,:,G.i2))
    end
    return A
end

@noinline function rmul!(A::StridedMatrix{DoubleFloat{T}},
                         G::Givens{DoubleFloat{T}}
                         ) where T
    has_offset_axes(A) && throw(ArgumentError("not implemented "
                                              * "for offset axes"))
    m, n = size(A, 1), size(A, 2)
    if G.i2 > m
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    _rmul!(view(A,:,G.i1),view(A,:,G.i2),G,Vec{Npref,T})
end
function _rmul!(xv::StridedVector{DoubleFloat{T}},
                yv::StridedVector{DoubleFloat{T}},
                G::Givens{DoubleFloat{T}}, ::Type{Vec{N,T}}) where {N,T}
    n = size(xv,1)
    nd,nr = divrem(n, N)
    ahi, alo = HILO(G.c)
    chi = Vec{N,T}(ahi)
    clo = Vec{N,T}(alo)
    ahi, alo = HILO(G.s)
    shi = Vec{N,T}(ahi)
    slo = Vec{N,T}(alo)
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N

            xhi = vgethi(xv,i0,Vec{N,T})
            xlo = vgetlo(xv,i0,Vec{N,T})
            yhi = vgethi(yv,i0,Vec{N,T})
            ylo = vgetlo(yv,i0,Vec{N,T})

            uhi, ulo = dfvmul(xhi, xlo, chi, clo)
            vhi, vlo = dfvmul(yhi, ylo, shi, slo)
            zhi, zlo = dfvsub(uhi, ulo, vhi, vlo)
            vputhilo!(xv,i0,zhi,zlo)

            vhi, vlo = dfvmul(xhi, xlo, shi, slo)
            uhi, ulo = dfvmul(yhi, ylo, chi, clo)
            zhi, zlo = dfvadd(uhi, ulo, vhi, vlo)
            vputhilo!(yv,i0,zhi,zlo)
        end
    end
    (nr == 0) && return yv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            a1 = xv[i]
            a2 = yv[i]
            xv[i] = G.c*xv[i] - G.s*a2
            yv[i] = G.c*yv[i] + G.s*a1
        end
    end
    nothing
end


text:
module FinalExample

greet() = print("Hello World!")

include("extrafile_final.jl")

export my_fu

end # module


text:
struct Scenario
    mass_emission_rate::Union{Missing,Number}
    release_duration::Union{Missing,Number}
    jet_diameter::Union{Missing,Number}
    jet_velocity::Union{Missing,Number}
    jet_density::Union{Missing,Number}
    release_pressure::Union{Missing,Number}
    release_temperature::Union{Missing,Number}
    release_height::Union{Missing,Number}
    windspeed::Union{Missing,Number}
    ambient_density::Union{Missing,Number}
    ambient_pressure::Union{Missing,Number}
    ambient_temperature::Union{Missing,Number}
    pasquill_gifford::Union{Missing,String}
end

function Scenario(d::AbstractDict)
    return Scenario(
    get(Missing, d, :mass_emission_rate),
    get(Missing, d, :release_duration),
    get(Missing, d, :jet_diameter),
    get(Missing, d, :jet_velocity),
    get(Missing, d, :jet_density),
    get(Missing, d, :release_pressure),
    get(Missing, d, :release_temperature),
    get(Missing, d, :release_height),
    get(Missing, d, :windspeed),
    get(Missing, d, :ambient_density),
    get(Missing, d, :ambient_pressure),
    get(Missing, d, :ambient_temperature),
    get(Missing, d, :pasquill_gifford),
    )
end

function Scenario(base_scenario::Scenario; kwargs...)
    d = Dict([ key => getproperty(base_scenario, key) for key in fieldnames(Scenario)])
    for key in filter(k -> k ∈ fieldnames(Scenario), keys(kwargs))
        d[key] = kwargs[key]
    end
    return Scenario(d)
end

units = Dict{Symbol,String}([
    :mass_emission_rate => "kg/s",
    :release_duration => "s",
    :jet_diameter => "m",
    :jet_velocity => "m/s",
    :jet_density => "kg/m^3",
    :release_pressure => "Pa",
    :release_temperature => "K",
    :release_height => "m",
    :windspeed => "m/s",
    :ambient_density => "kg/m^3",
    :ambient_pressure => "Pa",
    :ambient_temperature => "K",
    :pasquill_gifford => ""
])

function Base.show(io::IO, ::MIME"text/plain", s::Scenario)
    print(io, "Release scenario:\n    ")
    for key in fieldnames(Scenario)
        val =  getproperty(s, key)
        unit = units[key]
        print(io, "$key: $val $unit \n    ")
    end
end


text:
module PlanktonIndividuals

if VERSION < v"1.6"
    error("This version of PlanktonIndividuals.jl requires Julia v1.6 or newer.")
end

export
    # Architectures
    Architecture, GPU, CPU,

    # Grids
    RegularRectilinearGrid, RegularLatLonGrid,
    VerticallyStretchedLatLonGrid, LoadVerticallyStretchedLatLonGrid,
    Periodic, Bounded,

    # Model
    PlanktonModel,
    CarbonMode, QuotaMode, MacroMolecularMode,

    # BoundaryConditions
    set_bc!,

    # Simulation
    PlanktonSimulation, update!, vel_copy!,

    # Parameters
    default_PARF, default_temperature,
    update_bgc_params, update_phyt_params,
    bgc_params_default, phyt_params_default,

    # Biogeochemistry
    generate_nutrients, default_nut_init,

    # Output
    PlanktonDiagnostics, PlanktonOutputWriter, interior,

    # Units
    second, minute, hour, meter, kilometer,
    seconds, minutes, hours, meters, kilometers,
    KiB, MiB, GiB, TiB

using CUDA
using Pkg.Artifacts

import Base: show

p=dirname(pathof(PlanktonIndividuals))
artifact_toml = joinpath(p, "../Artifacts.toml")
surface_mixing_vels_hash = artifact_hash("surface_mixing_vels", artifact_toml)
surface_mixing_vels = joinpath(artifact_path(surface_mixing_vels_hash)*"/velocities.jld2")
global_vels_hash = artifact_hash("OCCA_FlowFields", artifact_toml)
global_vels = joinpath(artifact_path(global_vels_hash)*"/OCCA_FlowFields.jld2")

"""
    AbstractMode
Abstract type for phytoplankton physiology modes supported by PlanktonIndividuals.
"""
abstract type AbstractMode end

"""
    CarbonMode <: AbstractMode
Type for the phytoplankton physiology mode which only resolves carbon quota.
"""
struct CarbonMode <: AbstractMode end

"""
    QuotaMode <: AbstractMode
Type for the phytoplankton physiology mode which resolves carbon, nitrogen, and phosphorus quotas.
"""
struct QuotaMode <: AbstractMode end

"""
    MacroMolecularMode <: AbstractMode
Type for the phytoplankton physiology mode which resolves marco-molecules.
"""
struct MacroMolecularMode <: AbstractMode end

include("Architectures.jl")
include("Units.jl")
include("Grids/Grids.jl")
include("Parameters/Parameters.jl")
include("Fields/Fields.jl")
include("Biogeochemistry/Biogeochemistry.jl")
include("Diagnostics/Diagnostics.jl")
include("Plankton/Plankton.jl")
include("Model/Model.jl")
include("Output/Output.jl")
include("Simulation/Simulation.jl")

using .Architectures
using .Grids
using .Parameters
using .Fields
using .Biogeochemistry
using .Diagnostics
using .Plankton
using .Model
using .Output
using .Simulation
using .Units

function __init__()
    if CUDA.has_cuda()
        CUDA.allowscalar(false)
    end
end

end # module


text:
""" Quantization: Interval clustering based on maximum entropy method.  """
function maxentropy(r::Int,data::Array{Float64,1})

    N = length(data)
    n = N ÷ r       # number of data points per chunk

    # throw away random data for equally sized chunks of data
    data = shuffle(data[:])[1:n*r]
    sort!(data)

    # reshape data into a matrix, each chunk one column
    dataview = reshape(view(data,:),(n,:))

    nodes = Array{Float64,1}(undef,r)

    for i in 1:r
        # arithmetic mean of min and max within chunk
        nodes[i] = (dataview[1,i] + dataview[end,i])/2
    end

    return nodes
end


text:
B_data = [
    -2.5   0.5     0.0   0.0   0.0   0.0   0.0   0.0   0.5      0.0    0.0    0.0     0.0     0.0       0.0   0.0
    1.0  -3.0     0.0   0.0   0.0   0.0   0.0   0.0   0.5      0.0    0.0    0.0     0.0     0.0       0.0   0.0
    0.0   0.0    -2.5   0.0   0.0   0.0   0.0   0.0   2.0      0.0    0.0    0.5     0.0     0.0       0.0   0.0
    0.0   0.0     0.0  -2.5   0.0   0.0   0.0   0.0   0.0      2.0    0.0    0.0     0.5     0.0       0.0   0.0
    0.0   0.0     0.0   0.0  -2.5   0.0   0.0   0.0   0.0      0.0    2.0    0.0     0.0     0.5       0.0   0.0
    0.0   0.0     0.0   0.0   0.0  -2.5   0.0   0.0   2.0      0.0    0.0    0.5     0.0     0.0       0.0   0.0
    0.0   0.0     0.0   0.0   0.0   0.0  -2.5   0.0   0.0      2.0    0.0    0.0     0.5     0.0       0.0   0.0
    0.0   0.0     0.0   0.0   0.0   0.0   0.0  -2.5   0.0      0.0    2.0    0.0     0.0     0.5       0.0   0.0
    0.0   0.0     1.0   0.0   0.0   0.0   0.0   0.0  -2.1875   0.25  -0.25   1.0     0.0     0.0       0.0   0.75
    0.0   0.0     0.0   1.0   0.0   0.0   0.0   0.0  -0.375   -2.0    1.0    0.0     1.0     0.0       0.0  -2.5
    0.0   0.0     0.0   0.0   1.0   0.0   0.0   0.0   0.0625  -0.25  -2.75   0.0     0.0     1.0       0.0   3.75
    0.0   5.625   0.0   0.0   0.0   1.0   0.0   0.0   2.0      0.0    0.0   -4.125  -0.375   0.09375   0.0   0.0
    0.0  -3.75    0.0   0.0   0.0   0.0   1.0   0.0   0.0      2.0    0.0    1.5    -3.0    -0.5625    0.0   0.0
    0.0   1.125   0.0   0.0   0.0   0.0   0.0   1.0   0.0      0.0    2.0   -0.375   0.375  -3.28125   0.0   0.0
    0.0   0.0     0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0    0.0    0.0     0.0     0.125    -2.5   2.0
    0.0   0.0     0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0    0.0    0.0     0.0     0.25      1.0  -2.0
]


text:
using StatsBase
using Distributions
using DelimitedFiles
using Combinatorics: with_replacement_combinations

import LinearAlgebra
import LinearAlgebra: diag
import Statistics: mean

mutable struct FactorGraph{}
    order_list::Array{Int64,1}
    variable_count::Int
    n_alphabets::Int
    terms::Dict{Tuple,Array} # TODO, would be nice to have a stronger tuple type
    #TODO In the sanity check make sure that an kth order hyper edge is mapped to  a k rank tensor
    #variable_names::Union{Vector{String}, Nothing}
    #FactorGraph(a,b,c,d,e) = check_model_data(a,b,c,d,e) ? new(a,b,c,d,e) : error("generic init problem")
end

permutations(items, order::Int; asymmetric::Bool = false) = sort(permutations([], items, order, asymmetric))

function permutations(partial_perm::Array{Any,1}, items, order::Int, asymmetric::Bool)
    """
    All possible permutations of a given size.
    If asymmetric is false, then it returns combinations of items of the given order
    If asymmetric is true it returns all possible tuples of the size given by order from items
    """
    if order == 0
        return [tuple(partial_perm...)]
    else
        perms = []
        for item in items
            if !asymmetric && length(partial_perm) > 0
                if partial_perm[end] >= item
                    continue
                end
            end
            perm = permutations(vcat(partial_perm, item), items, order-1, asymmetric)
            append!(perms, perm)
        end
        return perms
    end
end

function raw_sampler_potts(H::FactorGraph, n_samples::Int, centered::Bool)
    """
    Given the FactorGraph, return samples according to its Gibbs distribution

    """

    n = H.variable_count
    q = H.n_alphabets
    n_config = q^n
    configs = [ digits(i,base=q, pad=n) .+ 1 for i = 0:n_config-1]
    weights = [ exp(Energy_Potts(K, H, centered)) for K in configs ]
    #print(configs, weights/sum(weights))
    raw_samples =  wsample(configs, weights, n_samples)
    return raw_samples

   end

function Energy_Potts(state::Array{Int64,  1},H::FactorGraph, cent::Bool)
    """
    Given a state and a FactorGraph, return its energy
    """
    q =  H.n_alphabets
    b =  -1.0/q
    a = 1.0 -(1.0/q)
    if !(cent)
     E = 0.0

     for (e, theta) in H.terms
        edge=Any[]

        [push!(edge, state[j] ) for j in e]

        E += theta[edge...]
     end
    return E
    end

    if cent
        E = 0.0

     for (e, theta) in H.terms
        clrs=Any[]
        r = length(e) #order of interaction
        alphabet_keys =  permutations(Array(1:q), r, asymmetric=true) #No need to generate this everytime

        [push!(clrs, state[j] ) for j in e]
        ct =  Tuple(clrs)

        [E += a^(sum( ct.==c )) * b^(r -  sum(ct.==c)) *theta[c...] for c in alphabet_keys]
     end
    return E
    end

end

function TVD(truth::Dict{}, est::Dict{}, n_samples::Int)
    """
    Total variation distance between two distributions.
    """
    s = 0.0
    for (k ,v) in est
        if haskey(truth, k)
            s+= abs( v -  truth[k])
        else
            s+= v
        end
    end

    for (k,v) in truth
        if !haskey(est, k)
            s+=v
        end
    end

    return s/(2*n_samples)

end

function conditional_energy(u::Int, state::Array{Int64,  1},H::FactorGraph, cent::Bool)
"""
Given a state and a FactorGraph, return its energy
"""
q =  H.n_alphabets
b =  -1.0/q
a = 1.0 -(1.0/q)
if !(cent)
    E = 0.0

    for (e, theta) in H.terms
        if u in e
            edge=Any[]

            [push!(edge, state[j] ) for j in e]

            E += theta[edge...]
        end
    end
    return E
end

if cent
    E = 0.0

    for (e, theta) in H.terms
        if u in e
            clrs=Any[]
            r = length(e) #order of interaction
            alphabet_keys =  permutations(Array(1:q), r, asymmetric=true) #No need to generate this everytime

            [push!(clrs, state[j] ) for j in e]
            ct =  Tuple(clrs)

            [E += a^(sum( ct.==c )) * b^(r -  sum(ct.==c)) *theta[c...] for c in alphabet_keys]
        end
    end
    return E
    end

end

function pth_order_tensor(r)
    a =  Tuple([2 for i in 1:r])
    M =  zeros(a)
    for i in CartesianIndices(M)
        s = prod(2*[Tuple(i)...] .- 3)
        M[i] = s
    end
    return M
end


text:
using BlackBoxOptim

problem = BlackBoxOptim.example_problems["Ackley"]
numrepeats = 5
dim = 100
methods = [
  :generating_set_search,
  :probabilistic_descent,
  :adaptive_de_rand_1_bin_radiuslimited,
  :random_search,
  :separable_nes]
max_time = 10.0
ftol = 1e-5

BlackBoxOptim.repeated_bboptimize(numrepeats, problem, dim, methods, max_time, ftol)


text:
using Test
using SCIP
using SCIP_jll
using SCIP_PaPILO_jll

@show(@eval(SCIP, libscip) == SCIP_jll.libscip)
@show(SCIP_PaPILO_jll.is_available() && @eval(SCIP, libscip) == SCIP_PaPILO_jll.libscip)

@testset "MathOptInterface nonlinear expressions" begin
    include("MOI_nonlinear_exprs.jl")
end

@testset "direct library calls" begin
    include("direct_library_calls.jl")
end

@testset "managed memory" begin
    include("scip_data.jl")
end

# new type definitions in module (needs top level)
include("conshdlr_support.jl")

@testset "constraint handlers" begin
    include("conshdlr.jl")
end

include("sepa_support.jl")
@testset "separators" begin
    include("sepa.jl")
end
@testset "cut callbacks" begin
    include("cutcallback.jl")
end

const MOI_BASE_EXCLUDED = [
    "Semicontinuous",
    "ScalarAffineFunction_Semiinteger",
    "ScalarQuadraticFunction_Semiinteger",
    "VectorAffineFunction_GeometricMeanCone",
    "Indicator_GreaterThan",
    "Indicator_LessThan",
    "Indicator_ACTIVATE_ON_ZERO", # odd MOI bug?
    "test_constraint_get_ConstraintIndex", # accessing constraint from string name
    "BoundAlreadySet", # see TODO,
    "ScalarAffineFunction_ConstraintName", # get(::SCIP.Optimizer, ::Type{MathOptInterface.ConstraintIndex}, ::String)
    "duplicate_VariableName", # two identical variable names should error
    "test_conic_empty_matrix",
    "test_modification_delete_variables_in_a_batch",
    "test_modification_set_function_single_variable",
    "test_modification_set_scalaraffine_",
    "test_modification_set_singlevariable_",
    "test_modification_transform_",
    "test_nonlinear_", # None of tests provide expression graphs in the evaluator.
    "ObjectiveFunction_ScalarAffineFunction", # requires conversion of objective function
    "test_objective_set_via_modify", # ListOfModelAttributesSet
]

@testset "MathOptInterface tests (direct)" begin
    include("MOI_wrapper_direct.jl")
end

@testset "MathOptInterface additional tests" begin
    include("MOI_additional.jl")
end

@testset "MathOptInterface tests (bridged)" begin
    include("MOI_wrapper_bridged.jl")
end

@testset "MathOptInterface tests (bridged & cached)" begin
    include("MOI_wrapper_cached.jl")
end

include("MOI_conshdlr.jl")


text:
mock_run(cmd::Cmd) = cmd

Mocking.override(run, mock_run) do
    activate(TEST_CONFIG; dir=joinpath(TEST_TMP_DIR, "test-playground"))
    activate(TEST_CONFIG; name="myproject")
    activate(TEST_CONFIG; name="otherproject")
    activate(TEST_CONFIG)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -4.116064765345723 <= q <= -0.12537054293327043)

                     add_NL_constraint(m, :(sigmoid(-0.6620590692147021 + -0.9580027962881332*sigmoid(-0.38090481240177176 + -0.02939531016329644*sigmoid(0.2067810723024852 + 0.7469716803367321*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.9724733578803808*sigmoid(-0.4274287128450034 + -0.22403315015562253*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))) + -0.699289524220764*sigmoid(0.4498619401955626 + 0.5463276179102481*sigmoid(0.2067810723024852 + 0.7469716803367321*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.02415027804893688*sigmoid(-0.4274287128450034 + -0.22403315015562253*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + sigmoid(-0.8107501452143855 + 0.008470632510253484*sigmoid(-0.38090481240177176 + -0.02939531016329644*sigmoid(0.2067810723024852 + 0.7469716803367321*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.9724733578803808*sigmoid(-0.4274287128450034 + -0.22403315015562253*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))) + 0.39739497887500974*sigmoid(0.4498619401955626 + 0.5463276179102481*sigmoid(0.2067810723024852 + 0.7469716803367321*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.02415027804893688*sigmoid(-0.4274287128450034 + -0.22403315015562253*sigmoid(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*sigmoid(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
IP_back_legends_en = Dict("title" => "Information Plane (Backpropagation)", "cbartitle" => "Epochs")
IP_back_legends_pt = Dict("title" => "Plano de informação (Backpropagation)", "cbartitle" => "Épocas")

IP_legends_en = Dict("title" => "Information Plane", "cbartitle" => "Epochs")
IP_legends_pt = Dict("title" => "Plano de informação", "cbartitle" => "Épocas")

function plot_IP(model::MLPClassifier, input, output; imagename::String = "IP.png", nbins::Int = 10,  language::Symbol=:en)
    legends = language == :pt ? IP_legends_pt : IP_legends_en
    gr()
    positions = evaluate_information_plane(model, input, output, nbins)
    p = plot(positions, color=:purple, marker = (:circle, 6, 0.8, Plots.stroke(1, :gray)), legend=false,
             xlabel = L"I(X,Z)", ylabel = L"I(Z,Y)", title = legends["title"])
    savefig(p, imagename)
end

function get_net_sets(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, get_train_data, path)
    net_sets = map(1:nets_number) do i
        net = MLPClassifier(nfeatures, labels, hidden_layers)
        input_train, output_train = get_train_data(path, shuffle=true)
        net, input_train, output_train
    end
    net_sets
end

function plot_IP_average(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, get_train_data, path, epochs::Int,
        step_len::Int, imagename::String = "IP-average.png", nbins::Int=10; language::Symbol=:en)
    net_sets = get_net_sets(nfeatures, labels, hidden_layers, nets_number, get_train_data, path)
    generate_IP_average_plot(net_sets, nets_number, epochs, step_len, nbins, imagename, language=language)
end

function get_net_sets(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, data)
    net_sets = map(1:nets_number) do i
        net = MLPClassifier(nfeatures, labels, hidden_layers)
        train_data = shuffle(collect(data))
        (input_train, output_train) = getindex.(train_data, 1), getindex.(train_data, 2)
        net, input_train, output_train
    end
    net_sets
end

function plot_IP_average(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, data, epochs::Int,
        step_len::Int, imagename::String = "IP-average.png", nbins::Int=10; language::Symbol=:en)
    net_sets = get_net_sets(nfeatures, labels, hidden_layers, nets_number, data)
    generate_IP_average_plot(net_sets, nets_number, epochs, step_len, nbins, imagename, language=language)
end

function get_marker_vector(size::Int)
    if size == 0
        return []
    elseif size == 1
        return [:star6]
    elseif size >= 2
        markers = [:circle for i in 1:size]
        markers[end] = :star6
        return markers
    end
end

function generate_IP_average_plot(net_sets, nets_number::Int, epochs::Int, step_len::Int,
                                  nbins::Int=10, imagename::String = "IP-average.png";
                                   language::Symbol=:en)
    legends = language == :pt ? IP_back_legends_pt : IP_back_legends_en
    nets = getindex.(net_sets, 1)
    inputs = getindex.(net_sets, 2)
    outputs = getindex.(net_sets, 3)
    gr()
    #step_len = 10
    for step = 0:step_len:epochs
        map((n, in, out) -> fit!(n, in, out; max_epochs=step_len), nets, inputs, outputs)
        IP_points = map((n, in, out) -> evaluate_information_plane(n, in, out, nbins), nets, inputs, outputs)
        layer_points = collect.(hcat(IP_points...))'
        average_points = ((1/nets_number)*ones(Float64, (1,nets_number)))*layer_points |> adjoint
        x, y = separate_columns(average_points[:,1])
        z_colors = [cgrad(:inferno, [0.01, 0.99])[z] for z ∈ range(0.0, 1.0, length = Int((epochs/step_len)+1))]
        color_index = Int((step/step_len)+1)
        markers = get_marker_vector(length(x))
        if step == 0 #step_len
            plot(x, y, line_z = step, #=marker = (:circle, 6, 0.8, :inferno, Plots.stroke(1, :gray)),=# cbar = true, legend=false,
                 xlabel = L"I(X,Z)", ylabel = L"I(Z,Y)", title = legends["title"], cbartitle = legends["cbartitle"])
            Plots.scatter!(x, y, color = z_colors[color_index], m = markers, markersize = 7)
        else
            plot!(x, y, line_z = step)
            Plots.scatter!(x, y, color = z_colors[color_index], m = markers, markersize = 7)
            #Plots.scatter!(x[1:end-1], y[1:end-1], color = z_colors[color_index], markersize = 7)
            #Plots.scatter!([x[end]], [y[end]], color = z_colors[color_index], m = :hexagon, markersize = 7)
        end
    end
    savefig(imagename)
end

function calculate_IP_average_points(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, get_train_data, path, epochs::Int,
        step_len::Int, nbins::Int=10)
    net_sets = get_net_sets(nfeatures, labels, hidden_layers, nets_number, get_train_data, path)
    IP_average_points(net_sets, nets_number, epochs, step_len, nbins)
end

function calculate_IP_average_points(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, data, epochs::Int,
        step_len::Int, nbins::Int=10)
    net_sets = get_net_sets(nfeatures, labels, hidden_layers, nets_number, data)
    IP_average_points(net_sets, nets_number, epochs, step_len, nbins)
end

function IP_average_points(net_sets, nets_number::Int, epochs::Int, step_len::Int, nbins::Int=10)
    nets = getindex.(net_sets, 1)
    inputs = getindex.(net_sets, 2)
    outputs = getindex.(net_sets, 3)
    points = []
    #step_len = 10
    for step = 0:step_len:epochs
        println("Step: ", step)
        map((n, in, out) -> fit!(n, in, out; max_epochs=step_len), nets, inputs, outputs)
        IP_points = map((n, in, out) -> evaluate_information_plane(n, in, out, nbins), nets, inputs, outputs)
        layer_points = collect.(hcat(IP_points...))'
        average_points = ((1/nets_number)*ones(Float64, (1,nets_number)))*layer_points |> adjoint
        push!(points, average_points)
    end
    points
end

function generate_plot_IP_average(points, nets_number::Int, epochs::Int, step_len::Int,
                                  imagename::String = "IP-average.png"; language::Symbol=:en,
                                  xlim=nothing, ylim=nothing)
    legends = language == :pt ? IP_back_legends_pt : IP_back_legends_en
    gr()
    local p
    z_colors = [cgrad(:inferno, [0.01, 0.99])[z] for z ∈ range(0.0, 1.0, length = Int((epochs/step_len)+1))]
    for step = 0:step_len:epochs
        #step_index = Int(step/step_len)
        color_index = Int((step/step_len)+1)
        average_points = points[color_index]
        x, y = separate_columns(average_points[:,1])
        markers = get_marker_vector(length(x))
        if step == 0
            p = plot(x, y, line_z = step, #=marker = (:circle, 6, 0.8, :inferno, Plots.stroke(1, :gray)),=# cbar = true, legend=false,
                 xlabel = L"I(X,Z)", ylabel = L"I(Z,Y)", title = legends["title"], cbartitle = legends["cbartitle"])
            if xlim != nothing plot!(p, xlim=xlim) end
            if ylim != nothing plot!(p, ylim=ylim) end
            Plots.scatter!(p, x, y, color = z_colors[color_index],  m = markers, markersize = 7)

        else
            plot!(p, x, y, line_z = step)
            Plots.scatter!(p, x, y, color = z_colors[color_index], m = markers, markersize = 7)
        end
    end
    savefig(p, imagename)
end

function calculate_MI_output(model::MLPClassifier, input, Ŷ, epochs::Int, step::Int=1)
    epochs_array = 0:step:epochs
    num_epochs = length(epochs_array)
    MI_array = Vector{Float64}(undef, num_epochs)
    for epoch = 1:num_epochs
        if epoch > 0 fit!(model, input, Ŷ; max_epochs=step) end
        model_output = model.net.(input) |> flat_vector
        X_flatten = input |> flat_vector
        Ŷ_flatten =  convert_to_onehot(Ŷ, model.labels) |> flat_vector
        MI_array[epoch] = mutual_information(calculate_markov_chain_probability(Ŷ_flatten, X_flatten, model_output)...)
    end
    (epochs_array, MI_array)
end

function calculate_MI_hiddenlayer(model::MLPClassifier, input, output, hidden_pos::Int, epochs::Int, step::Int=1)
    epochs_array = 0:step:epochs
    num_epochs = length(epochs_array)
    MI_array = Vector{Float64}(undef, num_epochs)
    for epoch = 1:num_epochs
        if epoch > 0 fit!(model, input, output; max_epochs=step) end
        subnet = model.net[1:layer]
        subnet_output = subnet.(input) |> flat_vector
        X_flatten = input |> flat_vector
        MI_array[epoch] = mutual_information([X_flatten], [subnet_output])
    end
    (epochs_array, MI_array)
end

function calculate_MI_hiddenlayer(model::MLPClassifier, input, output, epochs::Int, step::Int=1)
    epochs_array = 0:step:epochs
    layers = number_hidden_layers(model)
    num_epochs = length(epochs_array)
    MI_array = zeros(Float64, (layers, num_epochs))
    for epoch = 1:num_epochs
        if epoch > 0 fit!(model, input, output; max_epochs=step) end
        for layer = 1:layers
            subnet = model.net[1:layer+1]
            subnet_output = subnet.(input) |> flat_vector
            X_flatten = input |> flat_vector
            MI_array[layer, epoch] = mutual_information([X_flatten], [subnet_output])
        end
    end
    (epochs_array, [MI_array[l,:] for l in 1:size(MI_array,1)]...,)
end

function plot_MI_hiddenlayers(model::MLPClassifier, input, output, epochs::Int=1, step::Int=1; imagename::String = "MI-hiddenlayers.png")
    gr()
    data = calculate_MI_hiddenlayer(model, input, output, epochs, step)
    p = plot(xlabel = "Epochs", ylabel = L"I(X,Z)", title = L"Epochs \times I(X,Z)",  legend=true)
    colors = [cgrad(:inferno, [0.01, 0.99])[z] for z ∈ range(0.0, 1.0, length = max(length(data)-1, 2))]
    for layer = 1:length(data)-1
        plot!(p, data[1], data[layer+1], color = colors[layer], marker = (6, 0.8, Plots.stroke(1, :gray)), label="I(X,Z$layer)")
    end
    savefig(imagename)
end

function plot_MI_output(model::MLPClassifier, input, output, epochs::Int=1, step::Int=1; imagename::String = "MI-output.png")
    gr()
    data = calculate_MI_output(model, input, output, epochs, step)
    p = plot(ylims = (0,1), xlabel = "Epochs", ylabel = L"I(X,Z)", title = "Epochs x I(Y,Ŷ)",  legend=true)
    p = plot!(p, data[1], data[2], color = :green, marker = (6, 0.8, Plots.stroke(1, :gray)),  label=L"I(Y,Ŷ)")
    savefig(p, imagename)
end

using Statistics

function plot_MLP_error_rate(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, nets_number::Int, train_data, test_data, epochs::Int,
        step_len::Int, imagename::String = "MLP-error_rate.png")

    net_sets = get_net_sets(nfeatures, labels, hidden_layers, nets_number, train_data)

    nets = getindex.(net_sets, 1)
    inputs = getindex.(net_sets, 2)
    outputs = getindex.(net_sets, 3)
    gr()

    steps = 0:step_len:epochs

    mean_train_error_rate = Vector{Float64}(undef, length(steps))
    mean_test_error_rate = Vector{Float64}(undef, length(steps))
    #step_len = 10

    i =  1
    for step = steps
        if step != 0 map((n, in, out) -> fit!(n, in, out; max_epochs=step_len), nets, inputs, outputs) end
        mean_train_error_rate[i] = map(n -> error_rate(train_data, n), nets) |> mean
        mean_test_error_rate[i] = map(n -> error_rate(test_data, n), nets) |> mean
        i += 1
    end

    println("Taxa erro final média (treinamento): ", mean_train_error_rate[end])
    println("Taxa erro final média (test): ", mean_test_error_rate[end])

    p = plot(title="Taxa de erro", xlabel="Época", ylabel="Taxa de erro")
    plot!(p, steps, mean_train_error_rate, label="Conjunto de treinamento", color=:green, line = 2)
    plot!(p, steps, mean_test_error_rate, label="Conjunto de teste", color=:blue, line = 2)
    savefig(p, imagename)
end

function plot_IP_network(nfeatures::Int, labels::Vector{<:Union{Int, Float64}},
        hidden_layers::Vector{Int}, data, epochs::Int, imagename::String = "IP-MLP-network.png",
        nbins::Int=10; language::Symbol=:en)

    net_sets = get_net_sets(nfeatures, labels, hidden_layers, 1, data)

    legends = language == :pt ? IP_legends_pt : IP_legends_en
    net = getindex.(net_sets, 1)
    input = getindex.(net_sets, 2)
    output = getindex.(net_sets, 3)
    gr()
    #step_len = 10

    if epochs > 0
        fit!(net[1], input[1], output[1]; max_epochs=epochs)
    end
    IP_points = evaluate_information_plane(net[1], input[1], output[1], nbins)
    layer_points = collect.(hcat(IP_points...))'
    #average_points = (ones(Float64, (1,1)))*layer_points |> adjoint
    x, y = getindex.(layer_points,1), getindex.(layer_points,2)

    color = cgrad(:inferno, [0.01, 0.99])[0.1]
    markers = get_marker_vector(length(x))

    plot(x, y, legend=false, color = color,
         xlabel = L"I(X,Z)", ylabel = L"I(Z,Y)", title = legends["title"])
    Plots.scatter!(x, y, color = color, m = markers, markersize = 7)

    savefig(imagename)
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
## Can we make a JavaScript-like object that natively supports a.b[1].c notation?

import JSON

export JSObject, stripobject

struct JSObject{T}
    data::T
end

const JSDict = Dict{String, Any}
const JSArray = Vector{Any}
const JSElement = Union{Integer, AbstractFloat, AbstractString}

JSObject(data::AbstractDict) = objectify(data)
JSObject(data::AbstractVector) = objectify(data)
JSObject(data::T) where T <: JSElement = JSObject{T}(data)

objectify(dict::AbstractDict) = JSObject{JSDict}(JSDict(String(key) => objectify(value) for (key, value) in dict))
objectify(array::AbstractVector) = JSObject{JSArray}(Any[objectify(item) for item in array])
objectify(x) = x

Base.length(object::JSObject) = length(object.data)
Base.size(object::JSObject) = size(object.data)
Base.getindex(object::JSObject, i...) = getindex(getfield(object, :data), i...)
Base.setindex!(object::JSObject, value, index...) = setindex!(getfield(object, :data), objectify(value), index...)
Base.show(io::IO, object::JSObject) = print(io, JSON.json(stripobject(object)))

Base.getproperty(object::JSObject{JSDict}, key::Symbol) = getindex(getfield(object, :data), String(key))
Base.setproperty!(object::JSObject{JSDict}, key::Symbol, value) = setindex!(getfield(object, :data), objectify(value), String(key))

## The reverse operations, going from JSObject to plain Julia dicts/vectors
Base.convert(::Type{Dict}, x::JSObject{JSDict}) = stripobject(x)
Base.convert(::Type{Array}, x::JSObject{JSArray}) = stripobject(x)

stripobject(x::JSObject{JSDict}) = JSDict(key => stripobject(value) for (key, value) in getfield(x, :data))
stripobject(x::JSObject{JSArray}) = Any[stripobject(element) for element in getfield(x, :data)]
stripobject(x::JSObject{T}) where T <: JSElement = getfield(x, :data)
stripobject(x) = x


text:
# The test is based on a modelling set-up designed for the
# 8th International Cloud Modelling Workshop
# (ICMW, Muhlbauer et al., 2013, case 1, doi:10.1175/BAMS-D-12-00188.1)
#
# See chapter 2 in Arabas et al 2015 for setup details:
#@Article{gmd-8-1677-2015,
#AUTHOR = {Arabas, S. and Jaruga, A. and Pawlowska, H. and Grabowski, W. W.},
#TITLE = {libcloudph++ 1.0: a single-moment bulk, double-moment bulk, and particle-based warm-rain microphysics library in C++},
#JOURNAL = {Geoscientific Model Development},
#VOLUME = {8},
#YEAR = {2015},
#NUMBER = {6},
#PAGES = {1677--1707},
#URL = {https://www.geosci-model-dev.net/8/1677/2015/},
#DOI = {10.5194/gmd-8-1677-2015}
#}

using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VTK

using LinearAlgebra
using StaticArrays
using Printf

using CLIMA.PlanetParameters
using CLIMA.MoistThermodynamics
using CLIMA.Microphysics

const _nstate = 5
const _ρ, _ρu, _ρw, _ρe_tot, _ρq_tot = 1:_nstate
const stateid = (ρid = _ρ, ρuid = _ρu, ρwid = _ρw,
                 ρe_tot_id = _ρe_tot, ρq_tot_id = _ρq_tot)
const statenames = ("ρ", "ρu", "ρw", "ρe_tot", "ρq_tot")

const _nauxcstate = 3
const _c_z, _c_x, _c_p = 1:_nauxcstate

# preflux computation for wavespeed function
@inline function preflux(Q)
  @inbounds begin
    # unpack all the state variables
    ρ, ρu, ρw, ρe_tot, ρq_tot = Q[_ρ], Q[_ρu], Q[_ρw], Q[_ρe_tot], Q[_ρq_tot]
    u, w, e_tot, q_tot = ρu / ρ, ρw / ρ, ρe_tot / ρ, ρq_tot / ρ

    return (u, w, ρ, q_tot, e_tot)
  end
end

# boundary condition
@inline function bcstate!(QP, VFP, auxP, nM, QM, VFM, auxM, bctype, t)
  @inbounds begin
    ρu_M, ρw_M, ρe_tot_M, ρq_tot_M = QM[_ρu], QM[_ρw], QM[_ρe_tot], QM[_ρq_tot]

    ρu_nM = nM[1] * ρu_M + nM[2] * ρw_M

    QP[_ρu] = ρu_M - 2 * nM[1] * ρu_nM
    QP[_ρw] = ρw_M - 2 * nM[2] * ρu_nM

    QP[_ρe_tot], QP[_ρq_tot] = ρe_tot_M, ρq_tot_M

    auxM .= auxP
  end
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t)
  u, w, ρ, q_tot, e_tot = preflux(Q)
  @inbounds abs(n[1] * u + n[2] * w)
end

@inline function constant_auxiliary_init!(aux, x, z, _...)
  @inbounds begin
    aux[_c_z] = z  # for gravity
    aux[_c_x] = x

    FT = eltype(aux)

    # initial condition
    θ_0::FT    = 289         # K
    p_0::FT    = 101500      # Pa
    p_1000::FT = 100000      # Pa
    qt_0::FT   = 7.5 * 1e-3  # kg/kg
    z_0::FT    = 0           # m

    R_m, cp_m, cv_m, γ = gas_constants(PhasePartition(qt_0))

    # Pressure profile assuming hydrostatic and constant θ and qt profiles.
    # It is done this way to be consistent with Arabas paper.
    # It's not neccesarily the best way to initialize with our model variables.
    p = p_1000 * ((p_0 / p_1000)^(R_d / cp_d) -
                R_d / cp_d * grav / θ_0 / R_m * (z - z_0)
               )^(cp_d / R_d)

    aux[_c_p] = p  # for prescribed pressure gradient (kinematic setup)
  end
end

# physical flux function
@inline function eulerflux!(F, Q, QV, aux, t)
  u, w, ρ, q_tot, e_tot = preflux(Q)
  @inbounds begin
    p = aux[_c_p]

    F .= 0
    # advect the moisture and energy
    F[1, _ρq_tot], F[2, _ρq_tot] = u *  ρ * q_tot,      w *  ρ * q_tot
    F[1, _ρe_tot], F[2, _ρe_tot] = u * (ρ * e_tot + p), w * (ρ * e_tot + p)
    # don't advect momentum (kinematic setup)
  end
end

# initial condition
const w_max = .6    # m/s
const Z_max = 1500. # m
const X_max = 1500. # m

@inline function single_eddy!(Q, t, x, z, _...)
  FT = eltype(Q)

  # initial condition
  θ_0::FT    = 289         # K
  p_0::FT    = 101500      # Pa
  p_1000::FT = 100000      # Pa
  qt_0::FT   = 7.5 * 1e-3  # kg/kg
  z_0::FT    = 0           # m

  R_m, cp_m, cv_m, γ = gas_constants(PhasePartition(qt_0))

  @inbounds begin
    # Pressure profile assuming hydrostatic and constant θ and qt profiles.
    # It is done this way to be consistent with Arabas paper.
    # It's not neccesarily the best way to initialize with our model variables.
    p = p_1000 * ((p_0 / p_1000)^(R_d / cp_d) -
                R_d / cp_d * grav / θ_0 / R_m * (z - z_0)
               )^(cp_d / R_d)
    T::FT = θ_0 * exner_given_pressure(p, PhasePartition(qt_0))
    ρ::FT = p / R_m / T

    # TODO should this be more "grid aware"?
    # the velocity is calculated as derivative of streamfunction
    ρu::FT = w_max * X_max/Z_max * cos(π * z/Z_max) * cos(2*π * x/X_max)
    ρw::FT = 2*w_max * sin(π * z/Z_max) * sin(2*π * x/X_max)
    u = ρu / ρ
    w = ρw / ρ

    ρq_tot::FT = ρ * qt_0

    e_int  = internal_energy(T, PhasePartition(qt_0))
    ρe_tot = ρ * (grav * z + (1//2)*(u^2 + w^2) + e_int)

    Q[_ρ], Q[_ρu], Q[_ρw], Q[_ρe_tot], Q[_ρq_tot] = ρ, ρu, ρw, ρe_tot, ρq_tot
  end
end

function main(mpicomm, FT, topl::AbstractTopology{dim}, N, timeend,
              ArrayType, dt) where {dim}

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )
  numflux!(x...) = NumericalFluxes.rusanov!(x...,
                                            eulerflux!,
                                            wavespeed
                                           )
  numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(x...,
                                                            eulerflux!,
                                                            bcstate!,
                                                            wavespeed
                                                           )

  # spacedisc = data needed for evaluating the right-hand side function
  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = _nstate,
                           flux! = eulerflux!,
                           numerical_flux! = numflux!,
                           numerical_boundary_flux! = numbcflux!,
                           auxiliary_state_length = _nauxcstate,
                           auxiliary_state_initialization! =
                             constant_auxiliary_init!)

  # This is a actual state/function that lives on the grid
  initialcondition(Q, x...) = single_eddy!(Q, FT(0), x...)
  Q = MPIStateArray(spacedisc, initialcondition)

  npoststates = 9
  v_q_liq, v_q_vap, v_q_tot, v_p, v_T, v_e_tot, v_e_int, v_e_kin, v_e_pot =
    1:npoststates
  postnames = ("q_liq", "q_vap","q_tot", "p", "T",
               "e_tot", "e_int", "e_kin", "e_pot")
  postprocessarray = MPIStateArray(spacedisc; nstate=npoststates)

  writevtk("initial_condition", Q, spacedisc, statenames)

  lsrk = LSRK54CarpenterKennedy(spacedisc, Q; dt = dt, t0 = 0)

  io = MPI.Comm_rank(mpicomm) == 0 ? stdout : devnull
  eng0 = norm(Q)
  @printf(io, "----\n")
  @printf(io, "||Q||₂ (initial) =  %.16e\n", eng0)

  # Set up the information callback
  timer = [time_ns()]
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(10, mpicomm) do (s=false)
    if s
      timer[1] = time_ns()
    else
      run_time = (time_ns() - timer[1]) * 1e-9
      (min, sec) = fldmod(run_time, 60)
      (hrs, min) = fldmod(min, 60)
      @printf(io, "----\n")
      @printf(io, "simtime =  %.16e\n", ODESolvers.gettime(lsrk))
      @printf(io, "runtime =  %03d:%02d:%05.2f (hour:min:sec)\n", hrs, min, sec)
      @printf(io, "||Q||₂  =  %.16e\n", norm(Q))
    end
    nothing
  end

  step = [0]
  mkpath("vtk")

  cbvtk = GenericCallbacks.EveryXSimulationSteps(60) do (init=false)
    DGBalanceLawDiscretizations.dof_iteration!(postprocessarray, spacedisc,
                                               Q) do R, Q, QV, aux
      @inbounds begin
        u, w, ρ, q_tot, e_tot = preflux(Q)
        z = aux[_c_z]
        p = aux[_c_p]

        e_int = e_tot - 1//2 * (u^2 + w^2) - grav * z
        ts = PhaseEquil(e_int, ρ, q_tot)  # saturation adjustment happens here
        pp = PhasePartition(ts)
        R[v_T] = ts.T
        R[v_p] = p

        R[v_q_tot] = q_tot
        R[v_q_vap] = q_tot - pp.liq
        R[v_q_liq] = pp.liq

        R[v_e_tot] = e_tot
        R[v_e_int] = e_int
        R[v_e_kin] = 1//2 * (u^2 + w^2)
        R[v_e_pot] = grav * z

      end
    end

    outprefix = @sprintf("vtk/ex_1_microphysics_sat_adj_%dD_mpirank%04d_step%04d",
                         dim, MPI.Comm_rank(mpicomm), step[1])
    @printf(io, "----\n")
    @printf(io, "doing VTK output =  %s\n", outprefix)
    writevtk(outprefix, Q, spacedisc, statenames, postprocessarray, postnames)
    step[1] += 1
    nothing
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  Qe = MPIStateArray(spacedisc,
                    (Q, x...) -> single_eddy!(Q, FT(timeend), x...))

  # Print some end of the simulation information
  engf = norm(Q)
  @printf(io, "----\n")
  @printf(io, "||Q||₂ ( final ) = %.16e\n", engf)
end

function run(dim, Ne, N, timeend, FT)

  CLIMA.init()
  ArrayType = CLIMA.array_type()

  mpicomm = MPI.COMM_WORLD

  brickrange = ntuple(j->range(FT(0); length=Ne[j]+1, stop=Z_max), 2)

  topl = BrickTopology(mpicomm, brickrange, periodicity=(true, false))
  dt = 1

  main(mpicomm, FT, topl, N, timeend, ArrayType, dt)

end

using Test
let
  timeend = 30 * 60
  numelem = (75, 75)
  lvls = 3
  dim = 2
  FT = Float64

  polynomialorder = 4

  run(dim, ntuple(j->numelem[j], dim), polynomialorder, timeend, FT)
end

nothing

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
function FDivGrad2Vec!(F,cCG,CG,Param)
nz=Param.Grid.nz;
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;

D1cCG = Param.CacheC1
D2cCG = Param.CacheC2
grad1CG = Param.CacheC3
grad2CG = Param.CacheC4
D1gradCG = Param.CacheC1
D2gradCG = Param.CacheC2
vC1 = Param.CacheC3
vC2 = Param.CacheC4

mul!(reshape(D1cCG,OP,OP*NF*nz),CG.DS,reshape(cCG,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(D2cCG,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(cCG,(2,1,3,4)),OP,OP*nz*NF))

grad1CG .= Param.dXdxIC11.*D1cCG .+ Param.dXdxIC21.*D2cCG
grad2CG .= Param.dXdxIC12.*D1cCG .+ Param.dXdxIC22.*D2cCG

D1gradCG .= Param.dXdxIC11.*grad1CG .+ Param.dXdxIC12.*grad2CG
D2gradCG .= Param.dXdxIC21.*grad1CG .+ Param.dXdxIC22.*grad2CG

mul!(reshape(vC1,OP,OP*NF*nz),CG.DW,reshape(D1gradCG,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(vC2,(2,1,3,4)),OP,OP*NF*nz),CG.DW,reshape(PermutedDimsArray(D2gradCG,(2,1,3,4)),OP,OP*nz*NF))
F .= F .- Param.HyperDDiv .* (vC1 .+ vC2) ./ Param.JC
end

function [div]=FDivGrad2Vec(cCG,CG,Param)
nz=Param.Grid.nz;
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
D1cCG=reshape(...
  CG.DS*reshape(cCG,OP,OP*NF*nz)...
  ,OP,OP,NF,nz);
D2cCG=permute(reshape(...
  CG.DS...
  *reshape(...
  permute(...
  reshape(cCG,OP,OP,NF,nz)...
  ,[2 1 3 4])...
  ,OP,OP*NF*nz)...
  ,OP,OP,NF,nz)...
  ,[2 1 3 4]);
gradCG1=Param.dXdxIC(:,:,:,:,1,1).*D1cCG...
  +Param.dXdxIC(:,:,:,:,2,1).*D2cCG;
gradCG2=Param.dXdxIC(:,:,:,:,1,2).*D1cCG...
  +Param.dXdxIC(:,:,:,:,2,2).*D2cCG;
div=(reshape(...
  CG.DW*(reshape(gradCG1.*Param.dXdxIC(:,:,:,:,1,1)...
  +gradCG2.*Param.dXdxIC(:,:,:,:,1,2)...
  ,OP,OP*NF*nz))...
  ,OP,OP,NF,nz)...
  +...
  permute(...
  reshape(...
  CG.DW*reshape(...
  permute(...
  gradCG1.*Param.dXdxIC(:,:,:,:,2,1)...
  +gradCG2.*Param.dXdxIC(:,:,:,:,2,2)...
  ,[2 1 3 4])...
  ,OP,OP*NF*nz)...
  ,OP,OP,NF,nz)...
  ,[2 1 3 4]))...
  ./Param.JC;
end


text:
using AssemblyCalculus
using Graphs
using Plots

function time_sim(assem_size::Int; mean_degree=100)
    n = assem_size^2
    p = mean_degree / n
    β = 0.1
    ba = Brain(num_areas=1, n=n, k=assem_size, p=p, β=β)
    stims = [rand_stim(ba[1])]
    assems = Assembly{Float64}[]
    timesteps = 50
    td = @timed as, sp = simulate!(stims, assems, timesteps)
    return td.time
end

# Precompile
time_sim(5)

assem_sizes = collect(100:50:350)
times = time_sim.(assem_sizes)
neurons = 2 .* assem_sizes .^ 2

plot(neurons, times, linewidth=3, color=:teal, legend=:bottomright)
xticks!(neurons, string.(neurons))
xlabel!("Number of neurons")
ylabel!("Seconds to run 50 iterations (PartialNeuralArea)")
title!("Time Complexity of Generative Synapses")


text:
#!/usr/bin/env julia
using GZip
using BioSequences
using ArgParse.ArgParseSettings
using ArgParse.@add_arg_table
using ArgParse.parse_args
using DataStructures
using DataFrames
using CSV

#=
  Calculates adapter-dimer and adapter ratio in total reads.
=#

function adapter_counter(fastq_file::String, csv_out::String, fw_adapt::String)
    adapt_array = []

    f = open(fw_adapt)
    lines = readlines(f)

    for i in lines
        push!(adapt_array, split(i, r"[\s\t]+")[1])
    end

    csv_data::Array{String} = []

    # Adds header.
    data = DataFrame(sample=String[], total_reads=Int64[], adapter=Float64[],
                     adapter_dimer=Int64[])

    # Reads every FASTQ file and calculates its occurence and number of total
    # reads.
    read_counter = 0
    adapter_counter = 0
    adapter_dimer_counter = 0
    stream = GZip.open(fastq_file, "r")

     while !eof(stream)
        name = rstrip(readline(stream),'\n')
        sequence = rstrip(readline(stream),'\n')
        strand = rstrip(readline(stream),'\n')
        quality = rstrip(readline(stream),'\n')
        read_seq = DNASequence(String(sequence))

        approx = 0
        for i in adapt_array
            approx2 = approxsearchindex(read_seq, DNASequence(String(i)), 1)

            if approx < approx2
                approx = approx2
            end
        end

        if approx  > 0
            adapter_counter += 1
        end

        if  approx in 1:4
            adapter_dimer_counter += 1
        end

        read_counter += 1
    end
    close(stream)

    push!(data, [splitext(basename(fastq_file))[1], read_counter,
                 adapter_counter, adapter_dimer_counter])
    CSV.write(csv_out, data)
end

function main(args)
    arg_parse_settings = ArgParseSettings(description="Calculates aligned bases.")
    @add_arg_table arg_parse_settings begin
        "--fastq", "-f"
            arg_type = String
            help = "fastq.gz file."
            required = true
            dest_name = "fastq_f"
        "--out", "-o"
            arg_type = String
            help = "directory and file name as one string."
            required = false
            dest_name = "outfile"
            default = "./output/adapter_metrics_all_samples.csv"
        "--fwd", "-w"
            arg_type = String
            help = "File for adapters."
            required = true
            dest_name = "fwd"
    end

    parsed_args = parse_args(arg_parse_settings)
    adapter_counter(parsed_args["fastq_f"], parsed_args["outfile"],
                    parsed_args["fwd"])

end

main(ARGS)

max_stars_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_issues_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_forks_repo_name: KarolisM/rnd-it-molbio-fidelity-umi


text:
using Plots

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

# Ch 7.1 Principles of Regression
# Julia code to fit data points using a straight line
N = 50
x = rand(N)
a = 2.5                         # true parameter
b = 1.3                         # true parameter
y = a*x .+ b + 0.2*rand(N)      # Synthesize training data

X = [x ones(N)]                 # construct the X matrix
θ = X\y                         # solve y = X*θ

t    = range(0,stop=1,length=200)
yhat = θ[1]*t .+ θ[2]                  # fitted values at t

p1 = scatter(x,y,makershape=:circle,label="data",legend=:topleft)
plot!(t,yhat,linecolor=:red,linewidth=4,label="best fit")
display(p1)
##------------------------------------------------------------------------------

# Julia code to fit data using a quadratic equation
N = 50
x = rand(N)
a = -2.5
b = 1.3
c = 1.2
X = x.^[0 1 2]                #same as [ones(N) x x.^2]
y = X*[c,b,a] + rand(N)

θ    = X\y
t    = range(0,stop=1,length=200)
yhat = (t.^[0 1 2])*θ               #same as (t.^collect(0:2)')*θ

p2 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft)
plot!(t,yhat,linecolor=:red,linewidth=4,label="fitted curve")
display(p2)
##------------------------------------------------------------------------------

# Julia code to fit data using Legendre polynomials

using LegendrePolynomials

N = 50
x = rand(N)*2 .- 1
a = [-0.001, 0.01, 0.55, 1.5, 1.2]
X = hcat([Pl.(x,p) for p=0:4]...)   # same as [Pl.(x,0) Pl.(x,1) Pl.(x,2) Pl.(x,3) Pl.(x,4)]
y = X*a + 0.5*randn(N)
θ = X\y

t    = range(-1,stop=1,length=200)
yhat = hcat([Pl.(t,p) for p=0:4]...)*θ

p3 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft)
plot!(t,yhat,linecolor=:red,linewidth=4,label="fitted curve")
display(p3)
##------------------------------------------------------------------------------

# Julia code for auto-regressive model
using ToeplitzMatrices

N = 500
y = cumsum(0.2*randn(N)) + 0.05*randn(N)        # generate data

L = 100                                         # use previous 100 samples
c = [0; y[1:400-1]]
r = zeros(L)
X = Matrix(Toeplitz(c,r))                       # Toeplitz matrix, converted
θ = X\y[1:400]                                  # solve y = X*θ
yhat = X*θ                                      # prediction

p4 = scatter(y[1:400],makershape=:circle,label="data",legend=:bottomleft)
plot!(yhat[1:400],linecolor=:red,linewidth=4,label="fitted curve")
display(p4)
##------------------------------------------------------------------------------

# Julia code to demonstrate robust regression

using LinearAlgebra, LegendrePolynomials, Convex, SCS
import MathOptInterface
const MOI = MathOptInterface

#----------------------------------------------------------
"""
    linprog_Convex(c,A,sense,b)

A wrapper for doing linear programming with Convex.jl/SCS.jl.
It solves, for instance, `c'x` st. `A*x<=b`
"""
function linprog_Convex(c,A,sense,b)
    n  = length(c)
    x  = Variable(n)
    #c1 = sense(A*x,b)             #restriction, for Ax <= b, use sense = (<=)
    c1 = A*x <= b
    prob = minimize(dot(c,x),c1)
    solve!(prob,()->SCS.Optimizer(verbose = false))
    x_i = prob.status == MOI.OPTIMAL ? vec(evaluate(x)) : NaN
    return x_i
end
#----------------------------------------------------------

N = 50
x = range(-1,stop=1,length=N)
a = [-0.001, 0.01, 0.55, 1.5, 1.2]
y = hcat([Pl.(x,p) for p=0:4]...)*a  + 0.05*randn(N)        # generate data
idx     = [10, 16, 23, 37, 45]
y[idx] .= 5

X = x.^[0 1 2 3 4]
A = [X -I; -X -I]
b = [y; -y]
c = [zeros(5);ones(N)]

θ = linprog_Convex(c,A,(<=),b)

t    = range(-1,stop=1,length=200)
yhat = (t.^[0 1 2 3 4])*θ[1:5]

p5 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft)
plot!(t,yhat,linecolor=:red,linewidth=4,label="fitted curve")
display(p5)
##------------------------------------------------------------------------------

# Ch 7.2 Overfitting
# Julia: An overfitting example

using LegendrePolynomials

N = 20
x = sort(rand(N)*2 .- 1)
a = [-0.001, 0.01, 0.55, 1.5, 1.2]
y = hcat([Pl.(x,p) for p=0:4]...)*a + 0.1*randn(N)

P = 20
X = hcat([Pl.(x,p) for p=0:P]...)

β = X\y

t    = range(-1,stop=1,length=50)
Xhat = hcat([Pl.(t,p) for p=0:P]...)
yhat = Xhat*β

p6 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft,ylims = (-3,3))
plot!(t,yhat,linecolor=:blue,linewidth=4,label="fitted curve")
display(p6)
##------------------------------------------------------------------------------

using Statistics

Nset = round.( Int,exp10.(range(1,stop=3,length=20)) )  #10,13,16,21,...

E_train = zeros(length(Nset))
E_test  = zeros(length(Nset))
a = [1.3, 2.5]
for j = 1:length(Nset)
    local N,x,E_train_temp,E_test_temp,X
    N = Nset[j]
    x = range(-1,stop=1,length=N)
    E_train_temp = zeros(1000)
    E_test_temp  = zeros(1000)
    X = [ones(N) x]
    for i = 1:1000
        local y,y1,θ,yhat
        y  = a[1] .+ a[2]*x + randn(N)
        y1 = a[1] .+ a[2]*x + randn(N)
        θ = X\y
        yhat = θ[1] .+ θ[2]*x
        E_train_temp[i] = mean(abs2,yhat-y)
        E_test_temp[i]  = mean(abs2,yhat-y1)
    end
    E_train[j] = mean(E_train_temp)
    E_test[j]  = mean(E_test_temp)
end

p7 = scatter(Nset,E_train,xscale=:log10,markershape=:x,markercolor=:black,label="Training Error")
scatter!(Nset,E_test,xscale=:log10,markershape=:circle,markercolor=:red,label="Testing Error")
plot!(Nset,1.0 .- 2.0./Nset,xscale=:log10,linestyle=:solid,color=:black,label="")
plot!(Nset,1.0 .+ 2.0./Nset,xscale=:log10,linestyle=:solid,color=:red,label="")
display(p7)
##------------------------------------------------------------------------------

# Ch 7.3 Bias and Variance
# Julia code to visualize the average predictor

N = 20
a = [5.7, 3.7, -3.6, -2.3, 0.05]
x = range(-1,stop=1,length=N)
X = x.^[0 1 2 3 4]

t    = range(-1,stop=1,length=50)
tMat = t.^[0 1 2 3 4]

yhat = zeros(50,100)                      #50x100 instead of 100x50
for i = 1:100
    local y,θ
    y = X*a + 0.5*randn(N)
    θ = X\y
    yhat[:,i] = tMat*θ
end

p8 = plot(t,yhat,linecolor=:gray,label="")
plot!(t,mean(yhat,dims=2),linecolor=:red,linewidth=4,label="")
display(p8)
##------------------------------------------------------------------------------

# Ch 7.4 Regularization
#  Julia code to demonstrate a ridge regression example

#----------------------------------------------------------
"""
    RidgeRegression(Y,X,λ,β₀=0)

Calculate ridge regression estimate with targetr vector β₀.
"""
function RidgeRegression(Y,X,λ,β₀=0)
    K = size(X,2)
    isa(β₀,Number) && (β₀=fill(β₀,K))
    b = (X'X+λ*I)\(X'Y+λ*β₀)      #same as inv(X'X+λ*I)*(X'Y+λ*β₀)
    return b
end
#----------------------------------------------------------

N = 20                                 # Generate data
x = range(-1,stop=1,length=N)
a = [0.5, -2, -3, 4, 6]
y = x.^[0 1 2 3 4]*a + 0.25*randn(N)

λ = 0.1                                # Ridge regression
d = 20
X = x.^collect(0:d-1)'
A = [X; sqrt(λ)I]
b = [y; zeros(d)]
θ = A\b
θ_alt = RidgeRegression(y,X,λ)         # alternative ridge regression
display([θ θ_alt])

t    = range(-1,stop=1,length=500)
Xhat = t.^collect(0:d-1)'
yhat = Xhat*θ

p9 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft)
plot!(t,yhat,linecolor=:blue,linewidth=4,label="fitted curve")
display(p9)
##------------------------------------------------------------------------------

#----------------------------------------------------------
"""
    LassoEN(Y,X,λ,δ,β₀)

Do Lasso (set λ>0,δ=0), ridge (set λ=0,δ>0) or elastic net regression (set λ>0,δ>0).
Setting β₀ allows us to specify the target level.

"""
function LassoEN(Y,X,λ,δ=0.0,β₀=0)

    K = size(X,2)
    isa(β₀,Number) && (β₀=fill(β₀,K))

    #b_ls = X\Y                    #LS estimate of weights, no restrictions

    Q  = X'X
    c  = X'Y                      #c'b = Y'X*b

    b  = Variable(K)              #define variables to optimize over
    L1 = quadform(b,Q)            #b'Q*b
    L2 = dot(c,b)                 #c'b
    L3 = norm(b-β₀,1)             #sum(|b-β₀|)
    L4 = sumsquares(b-β₀)         #sum((b-β₀)^2)

    Sol = minimize(L1-2*L2+λ*L3+δ*L4)      #u'u + λ*sum(|b|) + δ*sum(b^2), where u = Y-Xb
    solve!(Sol,()->SCS.Optimizer(verbose = false))
    b_i = vec(evaluate(b))

    return b_i

end
#----------------------------------------------------------

using LinearAlgebra, DelimitedFiles, Convex, SCS
import MathOptInterface
const MOI = MathOptInterface

data  = readdlm("dataset/ch7_data_crime.txt")
y     = data[:,1]        # the observed crime rate
X     = data[:,3:end]    # feature vectors
(N,d) = size(X)

(y,X) = (y.-mean(y),X.-mean(X,dims=1))           # de-meaning all variables

λset = exp10.(range(-2,stop=4,length=50))        # adjusted range of λ values since X/100 is used
(θ_store,θ_store2,θ_store3) = [zeros(50,d) for _=1:3]   # 50 x d instead of d x 50
for i = 1:length(λset)
    θ_store[i,:]  = LassoEN(y/100,X/100,λset[i])    # LASSO, /100 improves the numerical stability
    θ_store2[i,:] = LassoEN(y/100,X/100,0,λset[i])  # ridge regression
    #θ_store3[i,:] = RidgeRegression(y/100,X/100,λset[i])   # alternative ridge regression
end
#display(θ_store2-θ_store3)                                 # compare solutions

p10 = plot(λset,θ_store,xscale=:log10,title="LASSO",xlabel="λ",ylabel="feature attribute",
           label = ["funding" "% high" "% no high" "% college" "% graduate"],
            linecolor = [:blue :red :yellow3 :magenta :green],ylims=(-5,15))
display(p10)

p10b = plot(λset,θ_store2,xscale=:log10,title="ridge regression",xlabel="λ",ylabel="feature attribute",
           label = ["funding" "% high" "% no high" "% college" "% graduate"],
            linecolor = [:blue :red :yellow3 :magenta :green],ylims=(-5,15))
display(p10b)
##------------------------------------------------------------------------------

# Julia code to demonstrate overfitting and LASSO

using LinearAlgebra, LegendrePolynomials, Convex, SCS
import MathOptInterface
const MOI = MathOptInterface

N = 20                                          # Generate data
x = range(-1,stop=1,length=N)
a = [1, 0.5, 0.5, 1.5, 1]
y = hcat([Pl.(x,p) for p=0:4]...)*a + 0.25*randn(N)

d = 20                                          # Solve LASSO using Convex.jl
X = hcat([Pl.(x,p) for p=0:d-1]...)
λ = 2
θ = LassoEN(y,X,λ)

t    = range(-1,stop=1,length=200)
Xhat = hcat([Pl.(t,p) for p=0:d-1]...)
yhat = Xhat*θ

p11 = scatter(x,y,makershape=:circle,label="data",legend=:bottomleft)
plot!(t,yhat,linecolor=:blue,linewidth=4,label="fitted curve")
display(p11)
#------------------------------------------------------------------------------


text:
## Functions to pull GeoFRED data
function fred_api_buildstr(series, apikey, date)
    url = "https://api.stlouisfed.org/geofred/series/data?"
    url = url * "series_id=$(series)&"; # Add series name
    url = url * "api_key=$(apikey)&"; # Add your API key
    url = url * "file_type=json&date=$(date)" # add file type and date
    return url
end

function fred_json_df(json_str, date_str)
    # Clean Fred specific JSON file to DataFrame of state level
    json_mat = JSON.parse(json_str);                          # Parse out json format.
    json_mat = json_mat["meta"]["data"][date_str]             # Pull data
    res_mat = DataFrame();
    for i = 1:length(json_mat)
        if i == 1
            res_mat = DataFrame(json_mat[i])
        end
        res_mat = vcat(res_mat, DataFrame(json_mat[i]));
    end
    res_mat["code"] = parse.(Float64, res_mat["code"])
    res_mat["value"] = parse.(Float64, res_mat["value"])
    return res_mat
end

# Unemployment, leveraging cal as the index, ***************************************************
unemp_map = "CAUR"

# Unemployment, baseline (2020 January) --------------------------------------------------------
# API request address build.
idx_url = fred_api_buildstr(unemp_map, my_api_key, date_baseline)
# API request
unemp_mapdata = HTTP.request("GET", idx_url);
unemp_mapdata = String(unemp_mapdata.body);
# Clean up data
map_data_unemp_start = fred_json_df(unemp_mapdata, date_baseline_txt)
# Unemployment, latest -------------------------------------------------------------------------
# API request address build.
idx_url = fred_api_buildstr(unemp_map, my_api_key, date_last)
# API request
unemp_mapdata = HTTP.request("GET", idx_url);
unemp_mapdata = String(unemp_mapdata.body);
# Clean up data
map_data_unemp_end = fred_json_df(unemp_mapdata, date_last_txt)

# Calculate the difference ---------------------------------------------------------------------
map_data_unemp = leftjoin(map_data_unemp_start, map_data_unemp_end,
                            on = :region, makeunique = true
                        );
map_data_unemp["Recovery"] = map_data_unemp["value_1"] - map_data_unemp["value"];

# Map the difference
map_unemp = @vlplot(
    :geoshape,
    width=500, height=300,
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    transform=[{
        lookup=:id,
        from={
            data=map_data_unemp,
            key=:code,
            fields=["Recovery"]
        }
    }],
    projection={
        type=:albersUsa
    },
    color= {
        "Recovery:q",
        legend = {
            title = "$(date_last_txt)"
        }

    },
    title = "Percentage point change in unemployment rate since $(date_baseline_txt)",
    config = {
        view = {
            stroke = "transparent"
        }
    }
)
VegaLite.save("map_unempdel.svg", map_unemp)

# Labor force participation rate, leveraging cal as the index, *********************************
lfptc_map = "LBSSA06"

# Labor force participation rate, baseline (2020 January) --------------------------------------
# API request address build.
idx_url = fred_api_buildstr(lfptc_map, my_api_key, date_baseline)
# API request
lfptc_mapdata = HTTP.request("GET", idx_url);
lfptc_mapdata = String(lfptc_mapdata.body);
# Clean up data
map_data_lfptc_start = fred_json_df(lfptc_mapdata, date_baseline_txt)
# Labor force participation rate, latest -------------------------------------------------------
# API request address build.
idx_url = fred_api_buildstr(lfptc_map, my_api_key, date_last)
# API request
lfptc_mapdata = HTTP.request("GET", idx_url);
lfptc_mapdata = String(lfptc_mapdata.body);
# Clean up data
map_data_lfptc_end = fred_json_df(lfptc_mapdata, date_last_txt)

# Calculate the difference ---------------------------------------------------------------------
map_data_lfptc = leftjoin(map_data_lfptc_start, map_data_lfptc_end,
                            on = :region, makeunique = true
                        );
map_data_lfptc["Recovery"] = map_data_lfptc["value_1"] - map_data_lfptc["value"];

# Map the difference
map_lfptc = @vlplot(
    :geoshape,
    width=500, height=300,
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    transform=[{
        lookup=:id,
        from={
            data=map_data_lfptc,
            key=:code,
            fields=["Recovery"]
        }
    }],
    projection={
        type=:albersUsa
    },
    color= {
        bin = true,
        field = "Recovery",
        legend = {
            title = "$(date_last_txt)"
        }
    },
    title = "Percentage point change in labor force participation rate since $(date_baseline_txt)",
    config = {
        view = {
            stroke = "transparent"
        }
    }
)
VegaLite.save("map_lfptcdel.svg", map_lfptc)


text:
using LightGraphs
using Test
using Pkg

Pkg.activate(".")
using VertexEliminationOrder

# The file to read the graph from.
graph_file = "test/example_graphs/sycamore_53_20.gr"

ENV["JULIA_DEBUG"]=VertexEliminationOrder
g = smallgraph("house")
iterative_dissection(g)

g = graph_from_gr(graph_file)
times = Vector{Float64}()
for i in 1:10
    tic = time()
    @info "timing main algorithm iterative_dissection"
    iterative_dissection(g)
    push!(times, time() - tic)
end
@info "average time = $(round(sum(times) / 10; digits=2))"

@info "new wrapper single threaded"
order_tw_by_dissections_simple(g, 30)


text:
#include("./src/SeisBeamforming.jl")
include("/Users/kurama/Documents/kurama/research/SCECInternship2019/Beamforming/SeisBeamforming/src/SeisBeamforming.jl")
using .SeisBeamforming
using SeisIO, Noise, JLD2, FileIO, Dates, Printf, HDF5
#==================================================#
# Input Parameters

finame      = "./dataset/BP_beamform.jld2"
fodir   = "./dataset"
foname  = "BP_polarplot" # data is saved at ./dataset/$foname.jld2

#==================================================#

mkpath(fodir)
fopath=joinpath(fodir, foname*".h5")

t = jldopen(finame)

DLtimestamplist = t["info/DLtimestamplist"]

h5open(fopath, "w") do file
    write(file, "info/DLtimestamplist", string.(DLtimestamplist))  # alternatively, say "@write file A"
end

for tstamp = DLtimestamplist

    Beam = t[joinpath("BeamData", tstamp)]

    # obtain examined frequency vector
    fvec = Float64[]
    fvecstamplist = String[]
    for freqid = 1:length(Beam)
        push!(fvec, Beam[freqid].f[1])
        push!(fvecstamplist, @sprintf("f:%04.2fHz", Beam[freqid].f[1]))
    end

    group1 = joinpath("BeamPower", tstamp)
    h5open(fopath, "r+") do file
        write(file, group1*"/fvec",  fvec)
        write(file, group1*"/fvecstamplist",  fvecstamplist)
    end

    for freqid = 1:length(Beam)
        B = Beam[freqid]
        cvec = B.c
        θvec = round.(rad2deg.(B.θ), digits=3)
        freq = B.f
        time = string.(u2d.(B.t))

        group1 = joinpath("BeamPower", tstamp, fvecstamplist[freqid])
        h5open(fopath, "r+") do file
            write(file, group1*"/cvec",  cvec)
            write(file, group1*"/thetavec",  θvec)
            write(file, group1*"/f",  fvec[freqid])
            write(file, group1*"/time",  time)
        end

        for tid = 1:length(time)
            bp = B.bp[:, : , tid]

            group1 = joinpath("BeamPower", tstamp, fvecstamplist[freqid], time[tid])
            h5open(fopath, "r+") do file
                write(file, group1*"/bp",  abs.(bp))
            end
        end
    end

end

max_stars_repo_name: kura-okubo/SeisBeamforming

max_issues_repo_name: kura-okubo/SeisBeamforming

max_forks_repo_name: kura-okubo/SeisBeamforming


text:
function convert_xyz_latlon_to_enu(lat_lons::AbstractArray{<:AbstractFloat, 2},
                                   xyz_array::AbstractArray{<:AbstractFloat, 2})
    num_points = size(lat_lons, 2)
    out = Array{eltype(lat_lons), 2}(undef, (3, num_points))
    for i = 1:num_points
        out[:, i] = rotate_xyz_to_enu(xyz_array[:, i], lat_lons[:, i]...)
    end
    return out
end

"""Rotates a vector in XYZ coords to ENU

Args:
    xyz: 3xk array of xyz vectors
    lat: latitude (deg) of point to rotate into
    lon: longitude (deg) of point to rotate into

Reference: https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates

"""
function rotate_xyz_to_enu(xyz, lat, lon)
    # Rotate about axis 3 with longitude, then axis 1 with latitude
    R3 = rot(90 + lon, 3, in_degrees=true)
    R1 = rot(90 - lat, 1, in_degrees=true)
    return R1 * (R3 * xyz)
end

"""
Find a 3x3 euler rotation matrix given an angle and axis.

Rotation matrix used for rotating a vector about a single axis.

Args:
    `angle` in degrees to rotate
    axis: 1, 2 or 3 for R1, R2, or R3
    in_degrees (bool): specify the angle in degrees. if false, using radians
"""
function rot(angle, axis; in_degrees=true)
    R = Matrix{Float64}(I, 3, 3)
    if in_degrees
        angle = deg2rad(angle)
    end
    cang = cos(angle)
    sang = sin(angle)
    if axis == 1
        R[2, 2] = cang
        R[3, 3] = cang
        R[2, 3] = sang
        R[3, 2] = -sang
    elseif axis == 2
        R[1, 1] = cang
        R[3, 3] = cang
        R[1, 3] = -sang
        R[3, 1] = sang
    elseif axis == 3
        R[1, 1] = cang
        R[2, 2] = cang
        R[2, 1] = -sang
        R[1, 2] = sang
    else
        throw("axis must be 1, 2 or 2")
    end
    return R
end

"""Find magnitude of an ENU vector in the LOS direction

Rotates the line of sight vector to ENU coordinates at
(lat, lon), then dots with the enu data vector

Args:
    enu (list[float], ndarray[float]): E,N,U coordinates, either
        as list of 3, or a (3, k) array of k ENU vectors
    los_vec (ndarray[float]) length 3 line of sight, in XYZ frame
    lat (float): degrees latitude of los point
    lon (float): degrees longitude of los point
    enu_coeffs (ndarray) size 3 array of the E,N,U coefficients
    of a line of sight vector. Comes from `find_enu_coeffs`.
        If this arg is used, others are not needed

Returns:
    ndarray: magnitudes same length as enu input, (k, 1)

Examples:
>>> print('%.2f' % project_enu_to_los([1,2,3],[1, 0, 0], 0, 0))
-2.00
>>> print('%.2f' % project_enu_to_los([1,2,3],[0, 1, 0], 0, 0))
-3.00
>>> print('%.2f' % project_enu_to_los([1,2,3],[0, 0, 1], 0, 0))
1.00
"""
function project_enu_to_los(enu, los_vec=nothing, lat=nothing, lon=nothing, enu_coeffs=nothing)
    if isnothing(enu_coeffs)
        los_hat = los_vec / norm(los_vec)
        enu_coeffs = rotate_xyz_to_enu(los_hat, lat, lon)
    end
    return enu_coeffs * enu
end


text:
module ChangeRespawn

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Change Respawn" => Ahorn.EntityPlacement(
        Maple.ChangeRespawnTrigger,
        "rectangle"
    )
)

function Ahorn.nodeLimits(trigger::Maple.ChangeRespawnTrigger)
    return 0, 1
end

end


text:
module Relaxation

	export longitudinal, transverse
	export fit_inversion_recovery, inversion_recovery_model
	export fit_transverse

	using LinearAlgebra
	using LsqFit

	@inline function longitudinal(t::Real, R1::Real, M::Real)::Real
		1.0 - (1.0 - M) * exp(-R1 * t)
	end
	@inline function transverse(t::Real, R2::Real)::Real
		exp(-R2 * t)
	end

	@inline function inversion_recovery_model(
		signal::AbstractVector{<: Real},
		t::AbstractVector{<: Real},
		(R1, Minv, M0)::AbstractVector{<: Real}
	)::Vector{Real}
		@. signal = M0 * abs(longitudinal(t, R1, Minv))
	end

	"""
		TODO: Two pool models?

	"""
	function fit_inversion_recovery(
		Tinv::AbstractVector{<: Real},
		signal::AbstractVector{<: Real},
		Δsignal::AbstractVector{<: Real},
		T1::Real,
		Minv::Real,
		M0::Real
	)::NTuple{6, Real}

		NaNs = ntuple(i -> NaN,  Val(6))

		result = curve_fit(
			inversion_recovery_model,
			Tinv,
			signal,
			1.0 ./ Δsignal.^2,
			[1.0 / T1, Minv, M0]; # Initial guess
			autodiff=:forwarddiff,
			inplace=true
		)
		!result.converged && return NaNs

		# Get parameters
		R1, Minv, M0 = coef(result)

		# Compute std-errors
		ΔR1, ΔMinv, ΔM0 = let
			J = result.jacobian
			Q, R = qr(J)
			det(R) == 0 && return NaNs
			Rinv = inv(R)
			covar = Rinv * Rinv'
			var = diag(covar)
			any(var .< 0.0) && return NaNs
			sqrt.(var)
		end

		# Get T1 from R1
		T1 = 1.0 / R1
		ΔT1 = T1^2 * ΔR1

		return T1, ΔT1, Minv, ΔMinv, M0, ΔM0
	end

	"""
		fit_transverse(TE, signal)

	Find the transverse relaxation time.

	Fits an exponential decay to the absolute of the transverse magnetisation `signal`
	measured at echo times `TE`.

	The exponential decay

		s(t) = M0 * exp(-R2 * t)

	with R2 = 1/T2, can be rewritten as

		log(s(t)) = log(M0) - R2 * t

	making it linear in log(M0) and R2.

	Errors for the fit are computed via Gaussian error propagation

		Δ(log(s(t))) = |∂log(s)/∂s| * Δs =  1/s * Δs

	where Δs is assumed constant. Note that s >=0 (absolute!).

	Errors of the parameter estimates are computed using the same technique

		ΔT2 = |∂T2/∂R2| * ΔR2 = T2² * ΔR2

		ΔM0 = |∂M0/∂log(M0)| * ΔM0 = M0 * Δlog(M0)

	signal must be > 0
	"""
	function fit_transverse(
		TE::AbstractVector{<: Real},
		signal::AbstractVector{<: Real},
		Δsignal::AbstractVector{<: Real}
	)::NTuple{4, Real}

		c, Δc = linreg(
			[ TE ones(length(TE)) ], # Use TE and constant term, the latter fitting for M0
			log.(signal),
			signal ./ Δsignal # 1 / Propagated error of log.(signal)
		)
		isnan(c[1]) && return NaN, NaN, NaN, NaN

		R2 = -c[1]
		T2 = 1.0 / R2
		ΔT2 = T2^2 * Δc[1]

		M0 = exp(c[2])
		ΔM0 = M0 * Δc[2]

		return T2, ΔT2, M0, ΔM0
	end
	function linreg( # TODO: outsource
		X::Matrix{T},
		y::Vector{T},
		w::Vector{T}
	)::Tuple{Vector{T}, Vector{T}} where T <: Union{Float64, ComplexF64}
		# Make note about upper case variables
		W = diagm(w.^2)
		A = X' * W
		H = A * X
		if det(H) == 0
			θ = Vector{T}(undef, size(X, 1))
			fill!(θ, NaN)
			Δθ = θ
		else
			θ = H \ (A * y)
			Δθ = sqrt.(diag(inv(H)))
		end
		return θ, Δθ
	end

end