text:
"""
    MeanVariance <: AbstractMeanVariance

Mean Variance Formulation (Markowitz).

PS.: Equivalent to a single point uncertainty set.

Atributes:
- `predicted_mean::Array{Float64,1}`: Predicted mean of returns.
- `predicted_covariance::Array{Float64,2}`: Predicted covariance of returns.
"""
struct MeanVariance <: AbstractMeanVariance
    predicted_mean::Array{Float64,1}
    predicted_covariance::Array{Float64,2}
    number_of_assets::Int
end

function MeanVariance(;
    predicted_mean::Array{Float64,1},
    predicted_covariance::Array{Float64,2},
)
    number_of_assets = size(predicted_mean, 1)
    if number_of_assets != size(predicted_covariance, 1)
        error("size of predicted mean ($number_of_assets) different to size of predicted covariance ($(size(predicted_covariance, 1)))")
    end

    @assert issymmetric(predicted_covariance)
    return MeanVariance(
        predicted_mean, predicted_covariance, number_of_assets
    )
end

"""
    predicted_portfolio_return!(model::JuMP.Model, w, formulation::AbstractMeanVariance)

Return predicted portfolio return used in the formulation.

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: Portfolio optimization investment variable ("weights").
"""
function predicted_portfolio_return!(::JuMP.Model, w, formulation::AbstractMeanVariance)
    return sum(formulation.predicted_mean'w)
end

"""
    portfolio_return!(model::JuMP.Model, w, formulation::AbstractMeanVariance; kwargs...)

Return worst case portfolio return of the Uncertainty set defined by the formulation.

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: Portfolio optimization investment variable ("weights").
"""
function portfolio_return!(model::JuMP.Model, w, formulation::AbstractMeanVariance; kwargs...)
    return predicted_portfolio_return!(model, w, formulation; kwargs...)
end

"""
    predicted_portfolio_variance!(model::JuMP.Model, w, formulation::AbstractMeanVariance)

Return predicted portfolio variance used in the formulation:

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: Portfolio optimization investment variable ("weights").
"""
function predicted_portfolio_variance!(::JuMP.Model, w, formulation::AbstractMeanVariance)
    return sum(w'formulation.predicted_covariance * w)
end

"""
    portfolio_variance!(model::JuMP.Model, w, formulation::AbstractMeanVariance; kwargs...)

Return worst case portfolio variance of the Uncertainty set defined by the formulation.

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: Portfolio optimization investment variable ("weights").
"""
function portfolio_variance!(model::JuMP.Model, w, formulation::AbstractMeanVariance; kwargs...)
    return predicted_portfolio_variance!(model, w, formulation; kwargs...)
end

function _po_min_variance_limit_return_latex()
    return """
        ```math
        \\begin{aligned}
        \\min_{w} \\quad & V \\\\
        s.t. \\quad & R = (\\min r'w \\; | \\; r \\in \\Omega) \\\\
        & V = (\\max w ' \\Sigma w  \\; | \\; \\Sigma \\in \\Omega) \\\\
        & R \\geq R * W_0 \\\\
        & w \\in \\mathcal{X} \\\\
        \\end{aligned}
        ```
        """
end

"""
    po_min_variance_limit_return(formulation::AbstractPortfolioFormulation, R)

Mean-Variance Portfolio Allocation (with a risk free asset). Posed as a quadratic convex problem.
Minimizes the worst case portfolio variance (``V``) and limits the worst case portfolio return (``R``) in the uncertainty set (``\\Omega``)
to a minimal return parameter (``R_0``) normalized by current wealth (``W_0``).

$(_po_min_variance_limit_return_latex())

Where ``\\mathcal{X}`` represents the additional constraints defined in the model by the user
(e.g. a limit on maximum invested money).

Arguments:
 - `formulation::AbstractPortfolioFormulation`: Struct containing attributes of the formulation.
 - `minimal_return::Real`: Minimal normalized return accepted.

Optional Keywork Arguments:
 - `rf::Real = 0.0`: Risk-free return of money not invested.
 - `current_wealth::Real = 1.0`: Current available wealth to be invested.
 - `model::JuMP.Model = base_model(...)`: JuMP upper level portfolio optimization model. Defaults to a generic market specification.
 - `w`: Portfolio optimization investment variable reference ("weights").
 - `portfolio_return::Function = portfolio_return!`: Function that calculates ``R``.
 - `portfolio_variance::Function = portfolio_variance!`: Function that calculates ``V``.
"""
function po_min_variance_limit_return(formulation::AbstractPortfolioFormulation, minimal_return::Real;
    rf::Real = 0.0, current_wealth::Real = 1.0,
    model::JuMP.Model = base_model(formulation.number_of_assets; current_wealth=current_wealth),
    w=model[:w],
    portfolio_return::Function = portfolio_return!,
    portfolio_variance::Function = portfolio_variance!
)
    # auxilary variables
    @variable(model, R)
    sum_invested = create_sum_invested_variable(model, w)

    # model
    @constraint(model, R == portfolio_return(model, w, formulation) + rf * (current_wealth - sum_invested))
    @constraint(model, R >= minimal_return * current_wealth)
    @objective(model, Min, portfolio_variance(model, w, formulation))
    return model
end

function _po_max_return_limit_variance_latex()
    return """
        ```math
        \\begin{aligned}
        \\max_{w} \\quad & R \\\\
        s.t. \\quad & R = (\\min r'w \\; | \\; r \\in \\Omega) \\\\
        & V \\leq V_0 * W_0\\\\
        & V = (\\max w ' \\Sigma w  \\; | \\; \\Sigma \\in \\Omega) \\\\
        & w \\in \\mathcal{X} \\\\
        \\end{aligned}
        ```
        """
end

"""
    po_max_return_limit_variance(formulation::AbstractPortfolioFormulation, V_0)

Mean-Variance Portfolio Allocation (with a risk free asset). Posed as a quadratic convex problem.
Maximizes the worst case portfolio return (``R``) and limits the worst case portfolio variance (``V``) in the uncertainty set (``\\Omega``)
to a maximal risk parameter (``V_0``) normalized by current wealth (``W_0``).

$(_po_max_return_limit_variance_latex())

Where ``\\mathcal{X}`` represents the additional constraints defined in the model by the user
(e.g. a limit on maximum invested money).

Arguments:
 - `formulation::AbstractPortfolioFormulation`: Struct containing attributes of formulation.
 - `max_risk::Real`: Maximal normalized variance accepted.

Optional Keywork Arguments:
 - `rf::Real = 0.0`: Risk-free return of money not invested.
 - `current_wealth::Real = 1.0`: Current available wealth to be invested.
 - `model::JuMP.Model = base_model(...)`: JuMP upper level portfolio optimization model. Defaults to a generic market specification.
 - `w`: Portfolio optimization investment variable reference ("weights").
 - `portfolio_return::Function = portfolio_return!`: Function that calculates ``R``.
 - `portfolio_variance::Function = portfolio_variance!`: Function that calculates ``V``.
"""
function po_max_return_limit_variance(formulation::AbstractPortfolioFormulation, max_risk::Real;
    rf::Real = 0.0, current_wealth::Real = 1.0,
    model::JuMP.Model = base_model(formulation.number_of_assets; current_wealth=current_wealth),
    w=model[:w],
    portfolio_return::Function = portfolio_return!,
    portfolio_variance::Function = portfolio_variance!
)
    # auxilary variables
    @variable(model, R)
    sum_invested = create_sum_invested_variable(model, w)

    # model
    @constraint(model, R == portfolio_return(model, w, formulation) + rf * (current_wealth - sum_invested))
    @constraint(model, portfolio_variance(model, w, formulation) <= max_risk * current_wealth)
    @objective(model, Max, R)
    return model
end

###################### Not updated ######################
"""
Mean-Variance Portfolio Allocation With no risk free asset. Analytical solution.
"""
function mean_variance_noRf_analytical(formulation::MeanVariance, R_0)
    # parameters
    r̄ = formulation.predicted_mean
    Σ = formulation.predicted_covariance
    numA = formulation.number_of_assets
    # model
    invΣ = pinv(Σ, 1E-25)
    x = zeros(numA)
    A = sum(invΣ'r̄)
    B = sum(r̄'r̄)
    C = sum(invΣ)
    one = ones(size(r̄, 1))
    D = sum(r̄' * (invΣ'r̄))
    mu = (R_0 * C - A) / (C * D - A)
    x = (1 / C) * (invΣ * one) + mu * invΣ * (r̄ - one * (A / C))
    return x
end


text:
@testset "Manifold" begin
    f_0 = 1575420e3
    c₀ = 299_792_458
    λ = c₀ / f_0
    ant_pos = SVector(λ / 4 * SVector(1, 1, 0), λ / 4 * SVector(-1, 1, 0), λ / 4 * SVector(1, -1, 0), λ / 4 * SVector(-1, -1 , 0))
    manifold = @inferred IdealManifold(f_0, ant_pos)

    @test @inferred(get_num_ants(manifold)) == 4

    num_ants = length(EXAMPLE_LUT)
    manifold_real = RealManifold(EXAMPLE_LUT)
    @test @inferred(get_num_ants(manifold_real)) == num_ants
end

ext: jl


text:
# Tests the load corpus function
include(Pkg.dir("KUparser/conll17/rnnparser.jl"))
include("set_testenv.jl")
include("../src/types.jl")
include("../src/preprocess.jl")
using JLD

# add-hoc solution for vocab type of the parser
function eqtypes{T}(a::T, b::T)
    for item in fieldnames(a)
        if item == :vocab
            x = eqtypes(getfield(a, item), getfield(b, item))
            (x != true) && error("field $item are not equal")
        else
            @assert(getfield(a, item) == getfield(b, item), "$item test failed in $(typeof(a))")
        end
    end
    return true
end

function test_loadcorpus()
    d = load(chmodel)
    vocab = create_vocab(d)
    corpus = loadcorpus(real_tdata, vocab)

    orig_corpora, orig_vocab = main("--load $chmodel --datafiles $real_tdata $real_ddata --otrain 0")
    orig_corpus = orig_corpora[1]

    for item in fieldnames(orig_vocab)
        @assert(getfield(orig_vocab, item) == getfield(vocab, item), "$item are not equal in those vocabularies")
    end
    info("Vocabulary assertion passed")

    for i=1:length(corpus)
        eqtypes(corpus[i], orig_corpus[i])
    end
    info("Datafiles read as the with original rnnparser")
    return vocab, corpus, orig_corpus, orig_vocab
end
!isinteractive() && test_loadcorpus()

# TODO locate en_chmodel.jld
# test with rnnparser.jl implementation and your implementation


text:
@testset "Core definitions" begin
    # Can't do much without a transformation, and the only one defined in
    # core.jl is IdentityTransformation... will have to test
    # ComposeTransformation in test/commontransformations.jl
    identity_trans = IdentityTransformation()

    @test inv(identity_trans) == identity_trans

    @test identity_trans ∘ identity_trans == identity_trans

    x = SVector(1.0, 2.0, 3.0)
    @test identity_trans(x) == x
end


text:
import SparseArrays.SparseVector
import NearestNeighbors: KDTree, nn, inrange
import LinearAlgebra.Vector

# Abstract type for super nodes
abstract type AbstractSuperNode{PT} end
# A struct that contains a supernode, the corresponding basis vectors, and their identity
struct SuperNodeBasis{PT, RT} <: AbstractSuperNode{PT}
    center::PT
    basis_functions::Vector{BasisFunction{PT, RT}}
end

# A struct that contains a supernode for the DOFs
struct SuperNodeDomain{PT} <: AbstractSuperNode{PT}
    center::PT
    domains::Vector{Domain{PT}}
end

function SuperNodeBasis(center::PT, basis_functions::Vector{BasisFunction{PT,RT}}) where {RT<:Real,PT<:AbstractArray{RT}}
    return SuperNodeBasis{PT, RT}(center, basis_functions)
end

function SuperNodeDomain(center::PT, domains::Vector{Domain{PT}}) where {PT<:AbstractArray{<:Real}}
    return SuperNodeDomain{PT}(center, domains)
end

function domains(node::SuperNodeDomain)
    return node.domains
end

function basis_functions(node::SuperNodeBasis)
    return node.basis_functions
end

function center(in::AbstractSuperNode)
    return in.center
end

# Construct supernodes about aggregation centers provided by the user
# centers is the vector of centers of basis functions, and
function construct_supernodes(aggregation_centers, basis_functions::AbstractVector{BasisFunction{PT,RT}}, tree_function) where {PT<:AbstractVector{<:Real},RT<:Real}
    @assert !isempty(aggregation_centers)
    # allocating output array
    # Constructing themembership lists of the different supernodes, by assigning them to the closest aggregation center
    member_lists = construct_member_lists(nn(tree_function(aggregation_centers), center.(basis_functions))[1])

    # allocating output array
    out = Vector{SuperNodeBasis{PT,RT}}(undef, length(member_lists))
    # list[1] contrains the id of the the aggregation center, list[2] contains the actual list of domains associated to it
    for (k, list) in enumerate(member_lists)
        out[k] = SuperNodeBasis(aggregation_centers[list[1]], basis_functions[list[2]])
    end
    return out
end

# Construct supernodes about aggregation centers provided by the user
# centers is the vector of centers of basis functions, and
function construct_supernodes(aggregation_centers, domains::AbstractVector{<:Domain}, tree_function)
    @assert !isempty(aggregation_centers)
    # Constructing themembership lists of the different supernodes, by assigning them to the closest aggregation center
    member_lists = construct_member_lists(nn(tree_function(aggregation_centers), center.(domains))[1])
    # Creating the new SuperNode

    # allocating output array
    out = Vector{SuperNodeDomain{eltype(aggregation_centers)}}(undef, length(member_lists))
    # list[1] contrains the id of the the aggregation center, list[2] contains the actual list of domains associated to it
    for (k, list) in enumerate(member_lists)
        out[k] = SuperNodeDomain(aggregation_centers[list[1]], domains[list[2]])
    end
    return out
end

function construct_multicolor_ordering(input_supernodes::AbstractArray{<:SuperNodeBasis}, ρh, tree_function)
    # Vector (colors) of Vectors of supernodes
    out = typeof(input_supernodes)[]
    assigned = falses(length(input_supernodes))
    # While not all nodes are assigned to a color
    tree = tree_function(center.(input_supernodes))
    while sum(length.(out)) < length(input_supernodes)
        ruled_out = falses(size(assigned))
        # add new color
        push!(out, typeof(input_supernodes)(undef, 0))
        for (k, node) in enumerate(input_supernodes)
            #Check that node was not ruled out or assigned yet
            if !(ruled_out[k] || assigned[k])
                # Add new node to present color
                push!(out[end], node)
                # Make note that node was assigned to color
                assigned[k] = true
                # Make note that its neighbors are not allowed to be assigned to the same color
                ruled_out[inrange(tree, center(node), 2 * ρh)] .= true
            end
        end
    end
    return out
end

# function that can directly take an array of arrays (corresponding to different scales) as supernodes
function construct_multicolor_ordering(input_supernodes::AbstractArray{<:AbstractArray{<:SuperNodeBasis}}, ρh::AbstractArray{<:Real}, tree_function)
    # Both lengths should be equal to total number of scales
    @assert length(ρh) == length(input_supernodes)
    q = length(ρh)
    out = Vector{typeof(input_supernodes)}(undef, q)
    for k = 1 : q
        out[k] = construct_multicolor_ordering(input_supernodes[k], ρh[k], tree_function)
    end
    return vcat(out...)
end


text:
abstract type PairKernel{K1<:Kernel,K2<:Kernel} <: CompositeKernel end

leftkern(k::PairKernel) = throw(MethodError(leftkern, (k,)))
rightkern(k::PairKernel) = throw(MethodError(rightkern, (k,)))
components(k::PairKernel) = [leftkern(k), rightkern(k)]

function Base.show(io::IO, pairkern::PairKernel, depth::Int = 0)
    pad = repeat(" ", 2 * depth)
    println(io, "$(pad)Type: $(typeof(pairkern))")
    show(io, leftkern(pairkern), depth+1)
    show(io, rightkern(pairkern), depth+1)
end

num_params(pairkern::PairKernel) = num_params(leftkern(pairkern))+num_params(rightkern(pairkern))
get_params(pairkern::PairKernel) = vcat(get_params(leftkern(pairkern)), get_params(rightkern(pairkern)))
get_param_names(pairkern::PairKernel) = composite_param_names([leftkern(pairkern), rightkern(pairkern)], :sk)

function set_params!(pairkern::PairKernel, hyp::Vector{Float64})
    npl = num_params(leftkern(pairkern))
    hyp_left = hyp[1:npl]
    hyp_right = hyp[npl+1:end]
    set_params!(leftkern(pairkern), hyp_left)
    set_params!(rightkern(pairkern), hyp_right)
end

##########
# Priors #
##########

function set_priors!(pairkern::PairKernel, priors::Array)
    npl = num_params(leftkern(pairkern))
    priors_left = priors[1:npl]
    priors_right = priors[npl+1:end]
    set_priors!(leftkern(pairkern), priors_left)
    set_priors!(rightkern(pairkern), priors_right)
end

get_priors(pairkern::PairKernel) = vcat(get_priors(leftkern(pairkern)), get_priors(rightkern(pairkern)))

#################
# PairData      #
#################

struct PairData{KD1 <: KernelData, KD2 <: KernelData} <: KernelData
    data1::KD1
    data2::KD2
end
function KernelData(pairkern::PairKernel, X::MatF64)
    kl = leftkern(pairkern)
    kr = rightkern(pairkern)
    # this is a bit broken:
    if kernel_data_key(kl, X) == kernel_data_key(kr, X)
        kdata = KernelData(kl, X)
        return PairData(kdata, kdata)
    else
        return PairData(
                KernelData(kl, X),
                KernelData(kr, X)
               )
    end
end

function kernel_data_key(pairkern::PairKernel, X::MatF64)
    kl = leftkern(pairkern)
    kr = rightkern(pairkern)
    @sprintf("PairData:%s+%s", kernel_data_key(kl, X), kernel_data_key(kr, X))
end


text:
module TestNoOpTask

using Base.Test
using SimpleTasksTests.Utils.TestTasks
using SimpleTasksTests.Utils.MockServices
using SimpleTasks.Examples.NoOpTask

import SimpleTasks.Tasks.DaemonTask
import SimpleTasks.Tasks.BasicTask
import SimpleTasks.Services.Datasource

function test_create()
    task = NoOpTaskDetails(make_valid_basic_info(), "TEST")
    task.basic_info.name = NoOpTask.NAME
    @test task != nothing
end

function test_run()
    task = NoOpTaskDetails(make_valid_basic_info(), "TEST")
    task.basic_info.name = NoOpTask.NAME
    datasource = MockDatasourceService()

    #prime the cache
    for test_input in TEST_INPUTS
        buffer = IOBuffer(test_input)
        seekstart(buffer)
        Datasource.put!(datasource, DaemonTask.make_key(task, test_input), buffer)
        Datasource.put!(datasource, DaemonTask.make_key(task, test_input), buffer)
    end

    DaemonTask.run(task, datasource)
end

function __init__()
    test_create()
    test_run()
end

end # module TestNoOpTask

max_forks_repo_name: seung-lab/SimpleTasks


text:
@testset "Arithmetic" begin
    @testset "Sequence" begin
        a = Sequence(ParameterSpace() × (Taylor(1) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1, collect(1.0:13.0))
        b = Sequence(ParameterSpace() × (Taylor(2) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1, collect(1.0:19.0))
        c = Sequence(ParameterSpace() × (Taylor(0) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1, collect(1.0:7.0))
        d = Sequence(ParameterSpace() × (Taylor(2) ⊗ Fourier(0, 1.0) ⊗ Chebyshev(1))^1, collect(1.0:7.0))

        @test +(4.0\(8.0*a*3.0))/3.0 == a + a == a - (-a) == a +̄ a == a -̄ (-a)

        @test a + b == a - (-b) == ladd!(a, copy(b)) == lsub!(a, -b)
        @test a +̄ b == a -̄ (-b) == radd!(copy(a), b) == rsub!(copy(a), -b)

        @test a + c == a - (-c) == radd!(copy(a), c) == rsub!(copy(a), -c)
        @test a +̄ c == a -̄ (-c) == ladd!(a, copy(c)) == lsub!(a, -c)

        @test a + d == a - (-d)
        @test a +̄ d == a -̄ (-d)
    end

    @testset "LinearOperator" begin
        𝒮₁ = ParameterSpace() × (Taylor(1) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1
        𝒮₂ = ParameterSpace() × (Taylor(2) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1
        𝒮₃ = ParameterSpace() × (Taylor(0) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1
        𝒮₄ = ParameterSpace() × (Taylor(2) ⊗ Fourier(0, 1.0) ⊗ Chebyshev(1))^1

        A = LinearOperator(𝒮₁, 𝒮₁, [i+j for i ∈ indices(𝒮₁), j ∈ indices(𝒮₁)])
        B = LinearOperator(𝒮₂, 𝒮₁, [i+j for i ∈ indices(𝒮₁), j ∈ indices(𝒮₂)])
        C = LinearOperator(𝒮₃, 𝒮₁, [i+j for i ∈ indices(𝒮₁), j ∈ indices(𝒮₃)])
        D = LinearOperator(𝒮₄, 𝒮₁, [i+j for i ∈ indices(𝒮₁), j ∈ indices(𝒮₄)])

        @test +(4.0\(8.0*A*3.0))/3.0 == A + A == A - (-A) == A +̄ A == A -̄ (-A)
        @test A == (A - I) + I == -(I - (I + A)) == radd!(rsub!(copy(A), I), I) == -lsub!(I, ladd!(I, copy(A)))

        @test A + B == A - (-B) == ladd!(A, copy(B)) == lsub!(A, -B)
        @test A +̄ B == A -̄ (-B) == radd!(copy(A), B) == rsub!(copy(A), -B)

        @test A + C == A - (-C) == radd!(copy(A), C) == rsub!(copy(A), -C)
        @test A +̄ C == A -̄ (-C) == ladd!(A, copy(C)) == lsub!(A, -C)

        @test A + D == A - (-D)
        @test A +̄ D == A -̄ (-D)

        @test A ∘ A == A*A == A^2
        @test A ∘ B == mul!(similar(B), A, B, true, false)
        @test B ∘ A == mul!(similar(A), B, A, true, false)
        @test A ∘ C == mul!(similar(C), A, C, true, false)
        @test C ∘ A == mul!(similar(A), C, A, true, false)
        @test A ∘ D == mul!(similar(D), A, D, true, false)
        @test D ∘ A == mul!(similar(A), D, A, true, false)
    end

    @testset "Convolution" begin
        function conv(a, b)
            n = fft_size(space(a), space(b))
            space_c = image(*, space(a), space(b))
            c = Sequence(space_c, Vector{Float64}(undef, dimension(space_c)))
            rifft!(c, fft(a, n) .* fft(b, n))
            return c
        end

        a = Sequence(Taylor(1), [1.0, 2.0])
        b = Sequence(Taylor(2), [1.0, 2.0, 3.0])
        @test fft_size(space(a), space(a)) == fft_size(space(a), 2)
        @test conv(a, a) ≈ a * a == a ^ 2
        @test conv(a, b) ≈ a * b
        @test a *̄ a == a ^̄ 2

        a = Sequence(Fourier(1, 1.0), [1.0, 2.0, 3.0])
        b = Sequence(Fourier(2, 1.0), [1.0, 2.0, 3.0, 4.0, 5.0])
        @test fft_size(space(a), space(a)) == fft_size(space(a), 2)
        @test conv(a, a) ≈ a * a == a ^ 2
        @test conv(a, b) ≈ a * b
        @test a *̄ a == a ^̄ 2

        a = Sequence(Chebyshev(1), [1.0, 2.0])
        b = Sequence(Chebyshev(2), [1.0, 2.0, 3.0])
        @test fft_size(space(a), space(a)) == fft_size(space(a), 2)
        @test conv(a, a) ≈ a * a == a ^ 2
        @test conv(a, b) ≈ a * b
        @test a *̄ a == a ^̄ 2

        a = Sequence(Taylor(1) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1), collect(1.0:12.0))
        b = Sequence(Taylor(2) ⊗ Fourier(0, 1.0) ⊗ Chebyshev(1), collect(1.0:6.0))
        @test fft_size(space(a), space(a)) == fft_size(space(a), 2)
        @test conv(a, a) ≈ a * a == a ^ 2
        @test conv(a, b) ≈ a * b
        @test a *̄ a == a ^̄ 2
    end
end


text:
# test/runtests.jl

using MIRT

using Plots: default

if !isinteractive() # for usual noninteractive test, do not show plots
	default(show=false, reuse=true) # trick from Plots/runtests.jl
	ENV["GKSwstype"] = "100" # also from Plots/runtests.jl
end

mirt_setup_plot() # see src/plot/setup.jl

test_all_mirt()


text:
function ei_global(g::SimpleGraph, groups::AbstractVector)
    L = ne(g)
    E = 0

    @inbounds for e in edges(g)
        groups[src(e)] == groups[dst(e)] && (E += 1)
    end

    I = L - E
    return (E - I) / (E + I)
end

function ei_group(g::SimpleGraph, groups::AbstractVector)

    ng = maximum(groups)
    E = zeros(ng)
    I = zeros(ng)

    @inbounds for v in vertices(g)
        @inbounds for nb in neighbors(g, v)
            if groups[v] == groups[nb]
                I[groups[v]] += 1
            else
                E[groups[v]] += 1
            end
        end
    end

    ei =  (E .- I ) ./ (E .+ I)

    return Dict(collect(1:ng) .=> ei)
end

function ei_indiv(g::SimpleGraph, groups::AbstractVector)

    I = zeros(nv(g))
    @inbounds for v in vertices(g)
        @inbounds for nb in neighbors(g, v)
            if groups[v] == groups[nb]
                I[v] += 1
            end
        end
    end

    deg = degree(g)
    E = deg .- I

    return  (E .- I) ./ deg
end

max_stars_repo_name: jfb-h/EIIndex

max_issues_repo_name: jfb-h/EIIndex

max_forks_repo_name: jfb-h/EIIndex


text:
"""
    Forcing(func; parameters=nothing, field_dependencies=(), discrete_form=false)

Returns a forcing function added to the tendency of an Oceananigans model field.

If `discrete_form=false` (the default), and neither `parameters` nor `field_dependencies`
are provided, then `func` must be callable with the signature

```
func(x, y, z, t)
```

where `x, y, z` are the east-west, north-south, and vertical spatial coordinates, and `t` is time.
Note that this form is also default in the constructor for `NonhydrostaticModel`, so that `Forcing` is
not needed.

If `discrete_form=false` (the default), and `field_dependencies` are provided,
the signature of `func` must include them. For example, if `field_dependencies=(:u, :S)`
(and `parameters` are _not_ provided), then `func` must be callable with the signature

```
func(x, y, z, t, u, S)`
```

where `u` is assumed to be the `u`-velocity component, and `S` is a tracer. Note that any field
which does not have the name `u`, `v`, or `w` is assumed to be a tracer and must be present
in `model.tracers`.

If `discrete_form=false` (the default) and `parameters` are provided, then the _last_ argument
to `func` must be `parameters`. For example, if `func` has no `field_dependencies` but does
depend on `parameters`, then it must be callable with the signature

```
func(x, y, z, t, parameters)
```

The object `parameters` is arbitrary in principle, however GPU compilation can place
constraints on `typeof(parameters)`.

With `field_dependencies=(:u, :v, :w, :c)` and `parameters`, then `func` must be
callable with the signature

```
func(x, y, z, t, u, v, w, c, parameters)
```

If `discrete_form=true` then `func` must be callable with the "discrete form"

```
func(i, j, k, grid, clock, model_fields)
```

where `i, j, k` is the grid point at which the forcing is applied, `grid` is `model.grid`,
`clock.time` is the current simulation time and `clock.iteration` is the current model iteration,
and `model_fields` is a `NamedTuple` with `u, v, w`, the fields in `model.tracers`,
and the fields in `model.diffusivity_fields`, each of which is an `OffsetArray`s (or `NamedTuple`s
of `OffsetArray`s depending on the turbulence closure) of field data.

When `discrete_form=true` and `parameters` _is_ specified, `func` must be callable with the signature

```
func(i, j, k, grid, clock, model_fields, parameters)
```

Examples
========

```jldoctest forcing
using Oceananigans

# Parameterized forcing
parameterized_func(x, y, z, t, p) = p.μ * exp(z / p.λ) * cos(p.ω * t)

v_forcing = Forcing(parameterized_func, parameters = (μ=42, λ=0.1, ω=π))

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}}
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()
```

Note that because forcing locations are regularized within the
`NonhydrostaticModel` constructor:

```jldoctest forcing
grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, forcing=(v=v_forcing,))

model.forcing.v

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}} at (Center, Face, Center)
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()
```

After passing through the constructor for `NonhydrostaticModel`, the `v`-forcing location
information is available and set to `Center, Face, Center`.

```jldoctest forcing
# Field-dependent forcing
growth_in_sunlight(x, y, z, t, P) = exp(z) * P

plankton_forcing = Forcing(growth_in_sunlight, field_dependencies=:P)

# output
ContinuousForcing{Nothing}
├── func: growth_in_sunlight (generic function with 1 method)
├── parameters: nothing
└── field dependencies: (:P,)
```

```jldoctest forcing
# Parameterized, field-dependent forcing
tracer_relaxation(x, y, z, t, c, p) = p.μ * exp((z + p.H) / p.λ) * (p.dCdz * z - c)

c_forcing = Forcing(tracer_relaxation,
                    field_dependencies = :c,
                            parameters = (μ=1/60, λ=10, H=1000, dCdz=1))

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :H, :dCdz), Tuple{Float64, Int64, Int64, Int64}}}
├── func: tracer_relaxation (generic function with 1 method)
├── parameters: (μ = 0.016666666666666666, λ = 10, H = 1000, dCdz = 1)
└── field dependencies: (:c,)
```

```jldoctest forcing
# Unparameterized discrete-form forcing function
filtered_relaxation(i, j, k, grid, clock, model_fields) =
    @inbounds - (model_fields.c[i-1, j, k] + model_fields.c[i, j, k] + model_fields.c[i+1, j, k]) / 3

filtered_forcing = Forcing(filtered_relaxation, discrete_form=true)

# output
DiscreteForcing{Nothing}
├── func: filtered_relaxation (generic function with 1 method)
└── parameters: nothing
```

```jldoctest forcing
# Discrete-form forcing function with parameters
masked_damping(i, j, k, grid, clock, model_fields, parameters) =
    @inbounds - parameters.μ * exp(grid.zᵃᵃᶜ[k] / parameters.λ) * model_fields.u[i, j, k]

masked_damping_forcing = Forcing(masked_damping, parameters=(μ=42, λ=π), discrete_form=true)

# output
DiscreteForcing{NamedTuple{(:μ, :λ), Tuple{Int64, Irrational{:π}}}}
├── func: masked_damping (generic function with 1 method)
└── parameters: (μ = 42, λ = π)
```
"""
function Forcing(func; parameters=nothing, field_dependencies=(), discrete_form=false)
    if discrete_form
        return DiscreteForcing(func; parameters=parameters)
    else
        return ContinuousForcing(func; parameters=parameters, field_dependencies=field_dependencies)
    end
end


text:
__precompile__()

"""
We implement the *COCOA* algorithm from the following reference with the
*smoothed hinge loss* function.

> Jaggi, M., Smith, V., Takác, M., Terhorst, J., Krishnan, S., Hofmann, T. and Jordan, M.I., 2014.
> *Communication-efficient distributed dual coordinate ascent*.
> In Advances in Neural Information Processing Systems (pp. 3068-3076).

"""
module COCOA

using ArgParse
using DistLearn

function parse_cmdline(cmdline :: Vector{String}, s :: ArgParseSettings)
  @add_arg_table s begin
    "--beta"
      arg_type = Float64
      default = 1.0
      help = "scaling parameter beta"
    "--smoothness"
      arg_type = Float64
      default = 1.0
      help = "the smoothness of the loss"
  end

  return parse_args(cmdline, s)
end

################################################################################
# Parameters
################################################################################
struct HyperParameters
  β  :: Float64
  λ  :: Float64  # regularizer coefficient
  γ  :: Float64  # smoothness
  R2 :: Float64  # max l2 norm^2 of train samples
  n  :: Int      # (avg) num train per worker
  m  :: Int      # num workers
end
function Base.show(io :: IO, hp :: HyperParameters)
  println(io, "COCOA Hyper Parameters")
  println(io, "----------------------")
  println(io, "    β: ", hp.β)
  println(io, "    λ: ", hp.λ)
  println(io, "    γ: ", hp.γ)
  println(io, "   R2: ", hp.R2)
  println(io, "    n: ", hp.n)
  println(io, "    m: ", hp.m)
end

mutable struct LocalParameters
  α :: Vector{Float64}  # local dual variables
end

mutable struct Parameters
  w :: Vector{Float64}  # primal weights

  Parameters(dim :: Int) = new(zeros(dim))
end

################################################################################
# Worker
################################################################################
"To be run on remote workers"
function make_worker(w_id :: Int, n_workers :: Int, dset_name :: String)
  dset_tr, dset_tt = load_datasets(dset_name, w_id, n_workers)
  worker = Worker{HyperParameters,LocalParameters}()
  worker.dset_tr = dset_tr
  worker.dset_tt = dset_tt

  n_dim, n_smp = get_num_features(dset_tr), get_num_samples(dset_tr)
  worker.lp = LocalParameters(zeros(n_smp))

  return worker
end

"To be run on remote workers"
function gather_local_statistics(worker_ref :: Future)
  worker = try_fetch(worker_ref)

  stats = Dict(:n_train => get_num_samples(worker.dset_tr),
               :n_test => get_num_samples(worker.dset_tt),
               :fea_dim => get_num_features(worker.dset_tr),
               :max_l2_norm => maximum(norm(worker.dset_tr.data[:, i])
                                       for i = 1:size(worker.dset_tr.data, 2)))
  return stats
end

###############################################################################
# COCOA Building Blocks
################################################################################
"""
    smoothed_hinge_gradient(pred, y)

Compute gradient w.r.t. pred for the smoothed hinge loss defined as:

```
loss(pred, y) =
  let margin = pred * y
    if margin >= 1
      0
    elseif margin <= 0
      0.5 - margin
    else
      0.5 * (1 - margin)^2
```

The smoothed hinge loss is a smoothed version of the hinge loss, with smoothness
parameter 1.0. If needed, one can define general `Loss` type to dispatch gradient
computation automatically for different loss.

Note because the COCOA here compute one example at a time, so we define this
function for scalar prediction and label.
"""
function smoothed_hinge_gradient(pred :: Float64, y :: Int)
  margin = pred * y
  if margin >= 1
    return 0
  elseif margin <= 0
    return -y
  else
    return - (1 - margin) * y
  end
end

"To be run on remote worker"
function sdca(worker_ref :: Future, w :: Vector{Float64})
  worker = try_fetch(worker_ref)

  hp = worker.hp
  X = worker.dset_tr.data
  y = worker.dset_tr.labels
  Δα = zeros(size(worker.lp.α))
  α = copy(worker.lp.α)
  idx_all = shuffle(1:length(Δα))

  N = hp.n * hp.m
  λNγ = hp.λ * N * hp.γ

  for idx in idx_all
    pred = dot(X[:, idx], w)
    grad = smoothed_hinge_gradient(pred, y[idx])

    Δ = λNγ / (hp.R2 + λNγ) * (-grad - α[idx])

    α[idx] += Δ
    Δα[idx] += Δ
    w += Δ * X[:, idx] / (hp.λ * N)
  end

  worker.lp.α += (hp.β / hp.m) * Δα
  Δw = X * Δα / (hp.λ * N)
  Δw *= (hp.β / hp.m)
  return Δw
end

function cocoa_epoch(workers, params :: Parameters)
  Δw_all = invoke_on_workers(sdca, workers, params.w)
  params.w .+= sum(Δw_all)
end

################################################################################
# COCOA Algorithm
################################################################################
function cocoa(args)
  println("Spawning workers and loading data...")
  workers = spawn_workers(make_worker, args["data"])

  println("Gathering statistics from workers...")
  stats_all = invoke_on_workers(gather_local_statistics, workers)
  for stats in stats_all
    @assert stats[:fea_dim] == stats_all[1][:fea_dim]
  end
  n_tr_tot = sum(stats[:n_train] for stats in stats_all)

  println("Initializing parameters")
  params = Parameters(stats_all[1][:fea_dim])

  println("Initializing hyper-parameters")
  R2 = maximum(stats[:max_l2_norm] for stats in stats_all)^2
  n = round(Int, mean(stats[:n_train] for stats in stats_all))
  h_params = HyperParameters(args["beta"], args["regu-coef"],
                             args["smoothness"], R2, n, length(workers))
  invoke_on_workers(workers, h_params) do worker_ref, hp
    try_fetch(worker_ref).hp = hp
  end
  println(h_params)

  println("----------------------")
  println("Entering training loop...")
  for epoch = 1:args["n-epoch"]
    met_tr, met_tt = accumulate_metric(workers, Accuracy, params.w)

    # show progress
    print(@sprintf("E%03d WT(%s): ", epoch, Dates.format(now(), "dd.HH:MM:SS.sss")))
    print("TRAIN-", met_tr)
    println(", TEST-", met_tt)

    cocoa_epoch(workers, params)
  end
end

end  # module COCOA


text:
using Printf
using Statistics
using Random
using Oceananigans
using Oceananigans.Units
using GLMakie

function baroclinic_adjustment(kskew, ksymmetric, tapering, scale; dt = 20minutes, gradient = "y")

# Architecture
architecture = CPU()

# Domain
Lx = 1000kilometers
Ly = 1000kilometers  # north-south extent [m]
Lz = 1kilometers     # depth [m]
Nx = 20
Ny = 20
Nz = 20
save_fields_interval = 0.5day
stop_time = 30days
Δt = dt

zfaces = -Lz .+ (-0.5 * (cos.(range(0.0, π, length=Nz + 1))) .+ 0.5) .* Lz # range(-Lz, 0.0, length = Nz+1)
zfaces = range(-Lz, 0.0, length=Nz + 1)

grid = RectilinearGrid(architecture;
                       topology = (Bounded, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       x = (-Lx/2, Lx/2),
                       y = (-Ly/2, Ly/2),
                       z = zfaces,
                       halo = (3, 3, 3))

coriolis = FPlane(latitude = -45)

println("The diffusive timescale is ", (zfaces[1] - zfaces[2])^2 / 1e3)

Δy = Ly/Ny
@show κh = νh = Δy^4 / 10days
vertical_closure = VerticalScalarDiffusivity(ν=1e-2, κ=1e-4)
horizontal_closure = HorizontalScalarBiharmonicDiffusivity(ν=νh, κ=κh)

gerdes_koberle_willebrand_tapering = FluxTapering(tapering)
gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew=kskew,
    κ_symmetric=ksymmetric,
    slope_limiter=gerdes_koberle_willebrand_tapering)

closures = (vertical_closure, horizontal_closure, gent_mcwilliams_diffusivity)

@info "Building a model..."

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    coriolis = coriolis,
                                    buoyancy = BuoyancyTracer(),
                                    closure = closures,
                                    tracers = (:b, :c),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    free_surface = ImplicitFreeSurface())

@info "Built $model."

"""
Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:

y < y₀           => ramp = 0
y₀ < y < y₀ + Δy => ramp = y / Δy
y > y₀ + Δy      => ramp = 1
"""
function ramp(x, y, Δ)
    gradient == "x" && return min(max(0, x / Δ + 1/2), 1)
    gradient == "y" && return min(max(0, y / Δ + 1/2), 1)
    gradient == "xy" && return 0.5 *( min(max(0, x / Δ + 1/2), 1) + min(max(0, y / Δ + 1/2), 1) )
end

# Parameters
N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

Δy = 100kilometers * 1.0
Δz = 100

Δc = 2Δy
Δb = Δy * M²
ϵb = 0e-2 * Δb # noise amplitude

bᵢ(x, y, z) = N² * z + scale * Δb * ramp(x, y, Δy)
cᵢ(x, y, z) = 00.25 # exp(-x^2 / 2Δc^2) * exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2)

set!(model, b=bᵢ, c=cᵢ)

#####
##### Simulation building
#####

simulation = Simulation(model; Δt, stop_time)
wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(50))

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)
println("done with gradient in ", gradient)

loss = extrema(interior(simulation.model.tracers.c)[:, :, :, end])
return loss, simulation
end

klist = [0.0, 1e2, 1e3, 2e3]
taperinglist = [0.0, 1e-2, 2e-2, 0.5e-2, 10.0]
scalelist = [0.5, 1.0, 2.0]
namelist = (; κskew = 1e3, κsymmetric = 1e3, tapering = 1e-2, scale = 2.0)

clist = []
plist = []
for k in klist, tapering in taperinglist, scale in scalelist
    namelist = (; κskew=k, κsymmetric=k, tapering=tapering, scale=scale)
    println("current on ", namelist)
    cextrema, simulation = baroclinic_adjustment(namelist...)
    push!(clist, cextrema)
    push!(plist, namelist)
end

for i in eachindex(clist)
    if (abs(clist[i]) > 0.3) | (abs(clist[i]) < 0.2)
        println("-----------------------")
        println("For parameter values ", plist[i])
        println("the simulation was not good")
        println("----------------------")
    end
end

# (; κskew = 2000.0, κsymmetric = 2000.0, tapering = 0.01, scale = 2.0) failed


text:
# plot kernek density estimation

module MarginalKde700
    using DataFrames, CSV, StatsPlots
    pyplot()

    function main(is_test=false)
        # input
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_700.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        # extract between 12 and 16 o'clock
        df_ext = df_org[(df_org.time.>=120000).&(df_org.time.<160000),:]

        # plot marginal kde
        marginalkde(df_ext.ir, df_ext.lidar, c=:ice,
                    xlabel="ir", ylabel="liDAR")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/multi_dim_gauss_dist/marginal_kde/marginal_kde_700.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module MyIOModule

using DelimitedFiles

export readtodict
export niceprint
export nicewrite
export replacelineswithstring

function readtodict(filename::String, allowed)
    println("")
    f = open(filename)
    vars = Dict{String,Any}()
    println("Reading config from $filename")
    for line in eachline(f)

        # skip comments
        if length(line) == 0 || line[1] == '#'
            #println("skip comment line")
            continue
        end
        print("\t")
        println(line)

        var, val = split(line, "=")
        stripped = strip(var)
        if stripped in allowed
            vars[stripped] = parse(Float64, val)
        end
    end
    println("")
    return vars
end

function niceprint(A)
    println("")
    #Base.showarray(STDOUT,A,false)
    println(A)
    println("")
end

function nicewrite(filename::String, A)
    writedlm(filename, A, ", ")
end

function replacelineswithstring(
    filename::String,
    newfilename::String,
    keywords::Array{String},
    replaceby::Array{String},
)
    # replaces every line,  that contains the "keyword" with a line that consists of "replaceby"
    # result is written into a new file called "newfilename"

    # file from which to read
    f = open(filename)
    lines = readlines(f)

    # file to which we write
    open(newfilename, "w") do tmpf
        for l in lines
            containssomething = false
            matchid = 0
            for k = 1:size(keywords, 1)
                if occursin(keywords[k], l)
                    matchid = k
                    containssomething = true
                end
            end
            if !(containssomething)
                write(tmpf, l * "\n")
            else
                write(tmpf, replaceby[matchid] * "\n")
            end
        end
    end

end

end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
"""
    BlockIndex{N}

A `BlockIndex` is an index which stores a global index in two parts: the block
and the offset index into the block.

It can be used to index into `BlockArrays` in the following manner:

```jldoctest; setup = quote using BlockArrays end
julia> arr = Array(reshape(1:25, (5,5)));

julia> a = PseudoBlockArray(arr, [3,2], [1,4])
2×2-blocked 5×5 PseudoBlockArray{Int64,2}:
 1  │   6  11  16  21
 2  │   7  12  17  22
 3  │   8  13  18  23
 ───┼────────────────
 4  │   9  14  19  24
 5  │  10  15  20  25

julia> a[BlockIndex((1,2), (1,2))]
11

julia> a[BlockIndex((2,2), (2,3))]
20
```
"""
struct BlockIndex{N}
    I::NTuple{N, Int}
    α::NTuple{N, Int}
end

@inline BlockIndex(a::Int, b::Int) = BlockIndex((a,), (b,))
@inline BlockIndex(a::Tuple, b::Int) = BlockIndex(a, (b,))
@inline BlockIndex(a::Int, b::Tuple) = BlockIndex((a,), b)

@inline BlockIndex(a::Block, b::Tuple) = BlockIndex(a.n, b)
@inline BlockIndex(a::Block, b::Int) = BlockIndex(a, (b,))

@generated function BlockIndex(I::NTuple{N, Int}, α::NTuple{M, Int}) where {M,N}
    @assert M < N
    α_ex = Expr(:tuple, [k <= M ? :(α[$k]) : :(1) for k = 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds α2 = $α_ex
        BlockIndex(I, α2)
    end
end

"""
    global2blockindex{N}(block_sizes::BlockSizes{N}, inds...) -> BlockIndex{N}

Converts from global indices `inds` to a `BlockIndex`.
"""
@generated function global2blockindex(block_sizes::AbstractBlockSizes{N}, i::NTuple{N, Int}) where {N}
    block_index_ex = Expr(:tuple, [:(_find_block(block_sizes, $k, i[$k])) for k = 1:N]...)
    I_ex = Expr(:tuple, [:(block_index[$k][1]) for k = 1:N]...)
    α_ex = Expr(:tuple, [:(block_index[$k][2]) for k = 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds block_index = $block_index_ex
        @inbounds I = $I_ex
        @inbounds α = $α_ex
        return BlockIndex(I, α)
    end
end

"""
    blockindex2global{N}(block_sizes::BlockSizes{N}, block_index::BlockIndex{N}) -> inds

Converts from a block index to a tuple containing the global indices
"""
@generated function blockindex2global(block_sizes::AbstractBlockSizes{N}, block_index::BlockIndex{N}) where {N}
    ex = Expr(:tuple, [:(cumulsizes(block_sizes, $k, block_index.I[$k]) + block_index.α[$k] - 1) for k = 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds v = $ex
        return $ex
    end
end

# I hate having these function definitions but the generated function above sometimes(!) generates bad code and starts to allocate
@inline function blockindex2global(block_sizes::AbstractBlockSizes{1}, block_index::BlockIndex{1})
    @inbounds v =(cumulsizes(block_sizes, 1, block_index.I[1]) + block_index.α[1] - 1,)
    return v
end

@inline function blockindex2global(block_sizes::AbstractBlockSizes{2}, block_index::BlockIndex{2})
    @inbounds v =(cumulsizes(block_sizes, 1, block_index.I[1]) + block_index.α[1] - 1,
                  cumulsizes(block_sizes, 2, block_index.I[2]) + block_index.α[2] - 1)
    return v
end

@inline function blockindex2global(block_sizes::AbstractBlockSizes{3}, block_index::BlockIndex{3})
    @inbounds v =(cumulsizes(block_sizes, 1, block_index.I[1]) + block_index.α[1] - 1,
                  cumulsizes(block_sizes, 2, block_index.I[2]) + block_index.α[2] - 1,
                  cumulsizes(block_sizes, 3, block_index.I[3]) + block_index.α[3] - 1)
    return v
end

##
# checkindex
##

@inline checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::Block{N}) where N = blockcheckbounds(Bool, A, I.n...)
@inline function checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::BlockIndex{N}) where N
    checkbounds(Bool, A, Block(I.I)) || return false
    @inbounds block = getblock(A, I.I...)
    checkbounds(Bool, block, I.α...)
end

checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::AbstractVector{BlockIndex{N}}) where N =
    all(checkbounds.(Bool, Ref(A), I))


text:
@doc """
### NavierStokesMod.calcElemFurfaceArea
This function calculates the wet area of each element. A weight of 2 is given to
faces with Dirichlet boundary conditions.
Arguments:
mesh: AbstractMesh
sbp: SBP operator
eqn: an implementation of NSData. Does not have to be fully initialized.
"""->
# used by NSData Constructor
function calcElemSurfaceArea(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator,
                             eqn::NSData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}
  nfaces = length(mesh.interfaces)
  nrm = zeros(Tmsh, Tdim, mesh.numNodesPerFace)
  area = zeros(Tmsh, mesh.numNodesPerFace)
  face_area = zero(Tmsh)
  sbpface = mesh.sbpface

  #
  # Compute the wet area of each element
  #
  for f = 1:nfaces
    face = mesh.interfaces[f]
    eL = face.elementL
    eR = face.elementR
    fL = face.faceL
    fR = face.faceR
    #
    # Compute the size of face
    face_area = 0.0

    for n = 1 : mesh.numNodesPerFace
      nrm_xy = ro_sview(mesh.nrm_face, :, n, f)
      area[n] = norm(nrm_xy)
      face_area += sbpface.wface[n]*area[n]
    end

    eqn.area_sum[eL] += face_area
    eqn.area_sum[eR] += face_area
  end

  for bc = 1 : mesh.numBC
    indx0 = mesh.bndry_offsets[bc]
    indx1 = mesh.bndry_offsets[bc+1] - 1

    for f = indx0 : indx1
      face = mesh.bndryfaces[f].face
      elem = mesh.bndryfaces[f].element

      # Compute the size of face
      face_area = 0.0
      for n=1:mesh.numNodesPerFace
        nrm_xy = ro_sview(mesh.nrm_bndry, :, n, f)
        area[n] = norm(nrm_xy)
        face_area += sbpface.wface[n]*area[n]
      end
      eqn.area_sum[elem] += 2.0*face_area
    end
  end
  return nothing
end

@doc """

Compute the constant coefficent in inverse trace ineqality, i.e.,
the largest eigenvalue of
B^{1/2} R H^{-1} R^{T} B^{1/2}

Input:
  sbp
Output:
  cont_tii
"""->

function calcTraceInverseInequalityConst(sbp::AbstractOperator{Tsbp},
                                         sbpface::AbstractFace{Tsbp}) where Tsbp
  R = sview(sbpface.interp, :,:)
  BsqrtRHinvRtBsqrt = Array{Tsbp}(sbpface.numnodes, sbpface.numnodes)
  perm = zeros(Tsbp, sbp.numnodes, sbpface.stencilsize)
  Hinv = zeros(Tsbp, sbp.numnodes, sbp.numnodes)
  Bsqrt = zeros(Tsbp, sbpface.numnodes, sbpface.numnodes)
  for s = 1:sbpface.stencilsize
    perm[sbpface.perm[s, 1], s] = 1.0
  end
  for i = 1:sbp.numnodes
    Hinv[i,i] = 1.0/sbp.w[i]
  end
  for i = 1:sbpface.numnodes
    Bsqrt[i,i] = sqrt(sbpface.wface[i])
  end

  BsqrtRHinvRtBsqrt = Bsqrt*R.'*perm.'*Hinv*perm*R*Bsqrt
  const_tii = eigmax(BsqrtRHinvRtBsqrt)

  return const_tii

end

@doc """

Compute derivative operators

Input:
  mesh
  sbp
  dxidx : derivatives of mapping, i.e., Jacobian matrix
  jac   : determinant of Jacobian

Input/Output:
  Dx    : derivative operators in physical domain, incling Dx, Dy
"""->
function calcDx(sbp::AbstractOperator,
                dxidx::AbstractArray{Tmsh, 3},
                jac::AbstractArray{Tmsh, 1},
                Dx::AbstractArray{Tmsh, 3}) where Tmsh
  @assert(size(Dx, 1) == sbp.numnodes)
  @assert(size(Dx, 1) == size(dxidx, 3))
  @assert(size(Dx, 2) == size(Dx, 1))
  @assert(size(Dx, 3) == size(dxidx, 1))

  Dx[:,:,:] = 0.0
  dim = size(Dx, 3)
  numNodes = sbp.numnodes

  for d=1:dim            # loop over direction in which derivative is computing
    for dd=1:dim
      for n1 = 1:numNodes
        for n2 = 1:numNodes
          # Since dxidx is scaled by 1/|J|, we need to get it back,
          # that's why jac is here
          Dx[n1, n2, d] += dxidx[dd, d, n1]*jac[n1]*sbp.Q[n1, n2, dd]
        end
      end
    end

    # Until here Dx stores Qx, we need to left multiply H^(-1)
    for n2=1:numNodes
      for n1 = 1:numNodes
        Dx[n1, n2, d] /= sbp.w[n1]
      end
    end
  end
  return nothing
end

@doc """

Compute derivative operators

Input:
  mesh
  sbp
  elem    : index of element of which we are computing the derivatives

Input/Output:
  Dx    : derivative operators in physical domain, incling Dx, Dy
"""->
function calcDx(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                elem::Integer,
                Dx::AbstractArray{Tmsh, 3}) where Tmsh
  @assert(size(Dx, 1) == mesh.numNodesPerElement)
  @assert(size(Dx, 2) == mesh.numNodesPerElement)
  @assert(size(Dx, 3) == size(mesh.dxidx, 1))
  dxidx = sview(mesh.dxidx, :,:,:,elem) # (dim, dim, numNodesPerElement)
  jac = sview(mesh.jac, :, elem)
  dim = size(Dx, 3)

  for i = 1 : length(Dx)
    Dx[i] = 0.0
  end

  for d = 1 : dim            # loop over direction in which derivative is computing
    for dd = 1 : dim
      for n1 = 1 : mesh.numNodesPerElement
        for n2 = 1 : mesh.numNodesPerElement
          # Since dxidx is scaled by 1/|J|, we need to get it back,
          # that's why jac is here
          Dx[n1, n2, d] += dxidx[dd, d, n1] * jac[n1] * sbp.Q[n1, n2, dd]
        end
      end
    end

    # Until here Dx stores Qx, we need to left multiply H^(-1)
    for n2=1:mesh.numNodesPerElement
      for n1 = 1:mesh.numNodesPerElement
        Dx[n1, n2, d] /= sbp.w[n1]
      end
    end
  end

  return nothing
end

function calcQx(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                elem::Integer,
                Qx::AbstractArray{Tmsh, 3}) where Tmsh
  @assert(size(Qx, 1) == mesh.numNodesPerElement)
  @assert(size(Qx, 2) == mesh.numNodesPerElement)
  @assert(size(Qx, 3) == size(mesh.dxidx, 1))

  dxidx = sview(mesh.dxidx, :,:,:,elem) # (dim, dim, numNodesPerElement)
  jac = mesh.jac[:, elem]
  fill!(Qx, 0.0)
  dim = size(Qx, 3)

  for d = 1 : dim
    for dd = 1 : dim
      for l=1:mesh.numNodesPerElement
        # Since dxidx is scaled by 1/|J|, we need to get it back,
        # that's why jac is here
        for m = 1:mesh.numNodesPerElement
          Qx[l,m,d] += dxidx[dd,d,l] * jac[l] * sbp.Q[l,m,dd]
        end
      end
    end
  end
end

@doc """

Given variables q at element nodes, compute corresponding gradients

Input:
  sbp      : sbp operator
  dxidx    : derivatives of mapping, i.e., jacobian matrix
  jac      : determinant of jacobian
  q        : element node value
  q_grad   : gradient of q
Output:
  nothing
"""->
function calcGradient(sbp::AbstractOperator{Tsbp},
                      dxidx::AbstractArray{Tmsh, 3},
                      jac::AbstractArray{Tmsh, 1},
                      q::AbstractArray{Tsol, 2},
                      q_grad::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tsbp}
  @assert(size(q, 2) == sbp.numnodes)

  @assert(size(dxidx, 1) == size(dxidx, 2))
  @assert(size(q_grad, 1) == size(dxidx, 1))
  @assert(size(q_grad, 2) == size(q, 1))
  @assert(size(q_grad, 3) == sbp.numnodes)

  numNodes = sbp.numnodes
  numDofs = size(q, 1)
  dim = size(q_grad, 1)

  Dx = Array{Tsbp}(numNodes, numNodes, dim)

  calcDx(sbp, dxidx, jac, Dx)

  for i = 1 : length(q_grad)
    q_grad[i] = 0.0
  end

  for n=1:numNodes
    for iDof=1:numDofs
      for d=1:dim
        for col=1:numNodes
          q_grad[d, iDof, n] += Dx[n,col,d] * q[iDof, col]
        end
      end
    end
  end
  return nothing
end

@doc """

Given variables q at element nodes, compute corresponding gradients

Input:
	mesh:
	sbp:
	q      : element node value
	elem   : index of element
Output :
	q_grad : (in/out) gradient of q
"""->
function calcGradient(mesh::AbstractDGMesh{Tmsh},
                      sbp::AbstractOperator{Tsbp},
                      elem::Integer,
                      q::AbstractArray{Tsol, 2},
                      q_grad::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tsbp}
  @assert(size(q, 2) == mesh.numNodesPerElement)
  @assert(size(q, 1) == mesh.numDofPerNode)

  @assert(size(q_grad, 1) == size(mesh.coords, 1))
  @assert(size(q_grad, 3) == mesh.numNodesPerElement)
  @assert(size(q_grad, 2) == mesh.numDofPerNode)

  numNodes = mesh.numNodesPerElement
  numDofs = mesh.numDofPerNode
  dim = size(q_grad, 1)

  Dx = Array{Tsbp}(numNodes, numNodes, dim)
  # for e=1:numElems
  # First compute Dx for this element
  calcDx(mesh, sbp, elem, Dx)

  q_grad[:,:,:] = 0.0

  for n=1:numNodes
    for iDof=1:numDofs
      for d=1:dim
        for col=1:numNodes
          q_grad[d, iDof, n] += Dx[n,col,d]*q[iDof, col]
        end
      end
    end
  end
  return nothing
end

@doc """
Another(single face) version of interiorfaceintegrate.
Given Q on element L and element R, interpolate Q to interface shared by L and R

Input:
  sbpface     : face SBP operator
  face        : the interface which we are interpolating Q onto
  qvolL       : Q at nodes in element L, qvolL(idof, iNode)
  qvolR       : Q at nodes in element R
Output:
  qface       : Q at nodes on interface, qface(idof, L/R, ifacenode)

function call examples
"""->
function interiorfaceinterpolate(sbpface::AbstractFace{Tsbp},
                                 face::Interface,
                                 qvolL::AbstractArray{Tsol, 2},
                                 qvolR::AbstractArray{Tsol, 2},
                                 qface::AbstractArray{Tsol, 3}) where {Tsbp, Tsol}
  @assert(size(qvolR, 1) == size(qvolL, 1))
  @assert(size(qvolR, 2) == size(qvolL, 2))
  @assert(size(qvolR, 1) == size(qface, 1))
  @assert(size(sbpface.interp, 1) <= size(qvolL, 2))
  @assert(size(sbpface.interp, 2) == size(qface, 3))

  numDofs = size(qvolL, 1)

  for n = 1 : sbpface.numnodes
    iR = sbpface.nbrperm[n, face.orient]

    for dof = 1 : numDofs
      qface[dof, 1, n] = zero(Tsol)
      qface[dof, 2, n] = zero(Tsol)
    end

    for j = 1 : sbpface.stencilsize
      permL = sbpface.perm[j, face.faceL]
      permR = sbpface.perm[j, face.faceR]

      for dof = 1 : numDofs
        qface[dof, 1, n] += sbpface.interp[j, n] * qvolL[dof, permL]
        qface[dof, 2, n] += sbpface.interp[j, iR]* qvolR[dof, permR]
      end
    end
  end

  return nothing
end

@doc """
Another version of boundaryinterpolate.
Given Q on the parent element, interpolate Q to boundary face owned by parent element

Input:
  sbpface     : face SBP operator
  bndface     : the interface which we are interpolating Q onto
  qvol        : Q at nodes in element L, qvolL(idof, iNode)
Output:
  qface       : Q at nodes on interface, qface(idof, L/R, ifacenode)

"""->
function boundaryinterpolate(sbpface::AbstractFace{Tsbp},
                             bndface::Boundary,
                             qvol::AbstractArray{Tsol, 2},
                             qface::AbstractArray{Tsol, 2}) where {Tsbp, Tsol}

  @assert(size(qvol, 1) == size(qface, 1))            # dof
  @assert(size(qface,2) == sbpface.numnodes)
  @assert(size(sbpface.interp, 1) <= size(qvol, 2))    # stencilsize <= numNodesPerElem
  @assert(size(sbpface.interp, 2) == size(qface, 2))    # check number of face nodes

  numDofs = size(qvol, 1)

  for n = 1 : sbpface.numnodes
    for dof = 1 : numDofs
      qface[dof, n] = zero(Tsol)
    end

    for j = 1 : sbpface.stencilsize
      perm = sbpface.perm[j, bndface.face]
      for dof = 1 : numDofs
        qface[dof, n] += sbpface.interp[j, n]*qvol[dof, perm]
      end
    end
  end
end


text:
using OpenSimplex
using Test

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


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
#using Plasm
using ViewerGL
GL = ViewerGL
import Base.show

function show(filename)
	V, EV = Plasm.svg2lar(filename)
	#Plasm.view(V,EV)
	GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],1), GL.GLFrame2 ]);

	#Plasm.view(V,EV)
	return V, EV
end

function show(filename)
	V, EV = Lar.svg2lar(filename)
	GL.VIEW([ GL.GLLines(V,EV) ])
	return V, EV
end

#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/new.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/curved.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/twopaths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/paths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/boundarytest2.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/tile.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/interior.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/holes.svg")

V,EV = show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/Lar.svg")


text:
using Distributed
using SparseArrays
using Printf

import jInv.InverseSolve
import jInv.Mesh
import jInv.LinearSolvers
import jInv.ForwardShare
import jInv.Utils

# @everywhere begin
using jInv.InverseSolve
using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.Utils
using Test

mutable struct LSparam <: ForwardProbType
	A
	Ainv
end

import jInv.ForwardShare.getNumberOfData
function getNumberOfData(pFor::LSparam)
	return size(pFor.A,1)
end
import jInv.ForwardShare.getSensMatSize
function getSensMatSize(pFor::LSparam)
	return size(pFor.A)
end

import jInv.ForwardShare.getData
function getData(m::Vector,pFor::LSparam,doClear::Bool=false)
	  d = pFor.A*m
	  if doClear
	    clear!(pFor)
          end
	  return d,pFor
end

import jInv.ForwardShare.getSensMatVec
function getSensMatVec(v::Vector,m::Vector,pFor::LSparam)
	return pFor.A*v
end

import jInv.ForwardShare.getSensTMatVec
function getSensTMatVec(v::Vector,m::Vector,pFor::LSparam)
	return pFor.A'*v
end

import jInv.Utils.clear!
function clear!(pFor::LSparam)
	pFor.A = SparseMatrixCSC(I, 0, 0);
	pFor.Ainv = [];
end

# end


text:
####################################################################
function gpmodel(x, flux, σ2, idx; maxiter = 10, maxrandom = 100, seed = 1, show_trace = true)
####################################################################

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    JITTER = 1e-6

    N = length(x) # number of observations

    numF  = length(unique(idx)) # number of filters

    @assert(all(sort(unique(idx)) .== collect(1:numF)))

    @assert(N == length(flux) == length(idx) == length(σ2))

    Σdiag = Diagonal(σ2)

    @printf("Running gpmodel with %d observations and %d filters\n", N, numF)

    #--------------------------------------------------------
    # Pre-calculate distances
    #--------------------------------------------------------

    dist2 = pairwise(SqEuclidean(), reshape(x, 1, N), dims=2)

    K     = zeros(N, N)

    #--------------------------------------------------------
    function unpack(param)
    #--------------------------------------------------------

        @assert(length(param) == 2*numF + 1)

        local a     = param[0*numF+1:1*numF]
        local b     = param[1*numF+1:2*numF]
        local logs  = param[2*numF + 1]

        return a, b, logs

    end

    #--------------------------------------------------------
    function objective(param)
    #--------------------------------------------------------

        local a, b, logs = unpack(param)

        -1.0 * logl(a = a, logs = logs, b = b)

    end

    #--------------------------------------------------------
    function blockvector(idx, b)
    #--------------------------------------------------------

        local bv = zeros(eltype(b), N)

        for i in 1:numF
            bv[findall(idx .== i)] .= b[i] # this should be changes, too slow
        end

        return bv

    end

    #--------------------------------------------------------
    function logl(; a = a, b = b, logs = logs)
    #--------------------------------------------------------

        calculatekernelmatrix!(K, dist2, rbf, [0.0; logs])

        local A = Diagonal(blockvector(idx, a))

        local C = (A' * K * A) + Σdiag

        local C = (C.+C')*0.5 + JITTER*I

        local L = try

            cholesky(C).L

        catch err

            if isa(err, PosDefException)

                @warn("Cholesky failed, matrix not posdef. Return -Inf.")

                return -Inf

            else

                throw(err)

            end

        end

        local bvec    = blockvector(idx, b)

        local logdetC = 2*sum(log.(diag(L)))

        local yTCinvy = dot( (flux .- bvec), L'\(L\(flux .- bvec)) ) # this can be improved

        return - 0.5*logdetC - 0.5*yTCinvy - 0.5*N*log(2.0π)

    end

    #--------------------------------------------------------
    # Initialise GP hyperparameters with random search
    # and few iterations
    #--------------------------------------------------------

    initsolutions       = [randn(rg, 2*numF+1)*5 for i=1:maxrandom]
    initialfitness      = @showprogress "\t init random search " map(objective, initsolutions)
    bestinitialsolution = initsolutions[argmin(initialfitness)]

    #--------------------------------------------------------
    # Call optimiser
    #--------------------------------------------------------

    opt    = Optim.Options(iterations = maxiter, show_trace = show_trace, show_every = 100)
    result = Optim.optimize(objective, bestinitialsolution, NelderMead(), opt)

    # opt = NLopt.Opt(:LN_BOBYQA, length(bestinitialsolution))
    #
    # opt.xtol_rel = 1e-8
    #
    # COUNTER, BESTF = 0, Inf
    #
    # function myfunc(x, dummy)
    #
    #     COUNTER +=1
    #
    #     F = objective(x)
    #
    #     BESTF = min(F, BESTF)
    #
    #     mod(COUNTER, 100)==0 ? @printf("Iter %d f=%e\n", COUNTER, BESTF) : nothing
    #
    #     return F
    #
    # end
    #
    # opt.min_objective = myfunc
    #
    # (_minf, minx, _ret) = NLopt.optimize(opt, bestinitialsolution)

    #--------------------------------------------------------
    # Retrieve optimised parameters
    #--------------------------------------------------------

    a, b, logs = unpack(result.minimizer)

    #--------------------------------------------------------
    # Reconstruct curves
    #--------------------------------------------------------
    calculatekernelmatrix!(K, dist2, rbf, [0.0; logs])

    K = K + JITTER*I

    bvec = blockvector(idx, b)
    A    = Diagonal(blockvector(idx, a))
    L    = Σdiag\I
    Λ    = K\I

    Σ = (Λ + A'*L*A)\I

    μ = Σ * (A' * L * (flux .- bvec) )

    figure(101)

    cla()

    clr = ["b","g","r","k","y","m"]

    for (i,l) in enumerate(unique(idx))

        aux = findall(idx .== l)

        plot(t[aux], flux[aux], clr[i]*".", alpha=0.5)

        plot(t, b[i] .+ a[i]*μ, clr[i]*".")

    end

    # xtest = collect(LinRange(minimum(x), maximum(x), 200))
    #
    # # dimensions: N × Ntest
    # k = calculatekernelmatrix(x, xtest, rbf, [0.0;logs])

    # Ntest × 1
    # c = calculatekernelmatrix(xtest, xtest, rbf, uθ)

    # figure(102)
    #
    # cla()
    #
    # μ = k' * ((A'*K*A + σ2*I) \ (flux .- bvec))
    #
    # for (i,l) in enumerate(unique(idx))
    #
    #     plot(xtest, a[i]*μ .+ b[i], clr[i]*"-")
    #
    #     aux = findall(idx .== l)
    #
    #     plot(t[aux], flux[aux], clr[i]*".", alpha=0.5)
    #
    # end

    @show a, b, logs

end


text:
# Implementation of common model utilies

# character based model utilities
cembed(m) = m[1]
wchar(m) = m[2]; bchar(m) = m[3];
wforw(m) = m[4]; bforw(m) = m[5];
wback(m) = m[6]; bback(m) = m[7];
wsoft(m) = m[8]; bsoft(m) = m[9];

# To-load uniquely trained word-based language model
makewmodel1(d)=[ d["cembed"],
                 d["char"][1],
                 d["char"][2],
                 d["forw"][1],
                 d["forw"][2],
                 d["back"][1],
                 d["back"][2],
                 d["soft"][1],
                 d["soft"][2] ]

# In order to load the word-model either as gpu or normal Array
function makewmodel(d)
    d1 = makewmodel1(d)
    if gpu() >= 0
        return map(KnetArray, d1)
    else
        return map(Array, d1)
    end
end

# xavier initialization
function initx(d...; ftype=Float32)
    if gpu() >=0
        KnetArray{ftype}(xavier(d...))
    else
        Array{ftype}(xavier(d...))
    end
end

# random normal initialziation
function initr(d...; ftype=Float32, GPUFEATS=false)
    if GPUFEATS && gpu() >=0
        KnetArray{ftype}(0.1*randn(d...))
    else
        Array{ftype}(0.1*randn(d...))
    end
end

function makepmodel1(d; GPUFEATS=false)
    m = ([d["postagv"],d["deprelv"],d["lcountv"],d["rcountv"],d["distancev"],d["parserv"]],
         [d["postago"],d["deprelo"],d["lcounto"],d["rcounto"],d["distanceo"],d["parsero"]])

    if gpu() >= 0
        if GPUFEATS
            return map2gpu(m)
        else
            m = map2cpu(m)
            m[1][6] = map2gpu(m[1][6])
            m[2][6] = map2gpu(m[2][6])
            return m
        end
    else
        return map2cpu(m)
    end
end

function makepmodel2(o, s; ftype=Float32)
    model = Any[]
    dpostag, ddeprel, dcount = o[:embed]
    for (k,n,d) in ((:postag,17,dpostag),(:deprel,37,ddeprel),(:lcount,10,dcount),(:rcount,10,dcount),(:distance,10,dcount))
        push!(model, [ initr(d) for i=1:n ])
    end
    p = o[:arctype](s)
    f = features([p], o[:feats], model)
    mlpdims = (length(f), o[:hidden]..., p.nmove)
    info("mlpdims=$mlpdims")
    parser = Any[]
    for i=2:length(mlpdims)
        push!(parser, initx(mlpdims[i],mlpdims[i-1]))
        push!(parser, initx(mlpdims[i],1))
    end
    push!(model,parser)
    optim = initoptim(model,o[:optimization])
    return model,optim
end

# Parser model initialization, and parameter selection
makepmodel(d, o, s) = (haskey(d, "parserv") ? makepmodel1(d) : makepmodel2(o, s))
postagv(m)=m[1]; deprelv(m)=m[2]; lcountv(m)=m[3]; rcountv(m)=m[4]; distancev(m)=m[5]; parserv(m)=m[6]

function mlp(w,x; pdrop=(0,0))
    x = dropout(x,pdrop[1])
    for i=1:2:length(w)-2
        x = relu(w[i]*x .+ w[i+1])
        x = dropout(x,pdrop[2])
    end
    return w[end-1]*x .+ w[end]
end

# col-major lstm
function lstm(weight, bias, hidden, cell, input; mask=nothing)
    gates = weight * vcat(input, hidden) .+ bias
    H = size(hidden, 1)
    forget = sigm(gates[1:H, :])
    ingate = sigm(gates[1+H:2H, :])
    outgate = sigm(gates[1+2H:3H, :])
    change = tanh(gates[1+3H:4H, :])
    (mask != nothing) && (mask = reshape(mask, 1, length(mask)))

    cell = cell .* forget + ingate .* change
    hidden = outgate .* tanh(cell)
    if mask != nothing
        hidden = hidden .* mask
        cell = cell .* mask
    end
    return (hidden, cell)
end

# Token-batched words are coming as input to the system, col-major lstm
function charlstm(model, data, mask)
    weight, bias, embeddings = wchar(model), bchar(model), cembed(model)
    T = length(data)
    B = length(data[1])
    H = div(length(bias), 4)

    if isa(weight, KnetArray)
        mask = map(KnetArray, mask)
    end

    czero = fill!(similar(bias, H, B), 0)
    hidden = cell = czero
    for t in 1:T
        (hidden, cell) = lstm(weight, bias, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
    end
    return hidden
end

# col-major bilstm implementation
function wordlstm(model, data, mask, embeddings)
    weight, bias = wforw(model), bforw(model)
    T = length(data)
    B = length(data[1])
    H = div(length(bias), 4)

    if isa(weight, KnetArray)
        mask = map(KnetArray, mask)
    end

    wzero = fill!(similar(bias, H, B), 0)
    # forward lstm
    hidden = cell = wzero
    fhiddens = Array(Any, T-2)
    for t in 1:T-2
        (hidden, cell) = lstm(weight, bias, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
        fhiddens[t] = hidden
    end

    # backward lstm
    weight_b, bias_b = wback(model), bback(model)
    hidden = cell = wzero
    bhiddens = Array(Any, T-2)
    for t in T:-1:3
        (hidden, cell) = lstm(weight_b, bias_b, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
        bhiddens[t-2] = hidden
    end
    return fhiddens, bhiddens
end

# col-major loss function implementation, as explained in paper
function lmloss(model, data, mask, forw, back; result=nothing)
    T = length(data)
    B = length(data[1])
    weight, bias = wsoft(model), bsoft(model)
    idx(t,b,n) = data[t][b] + (b-1)*n

    total = count = 0
    for t in 1:T
        ypred = weight * vcat(forw[t], back[t]) .+ bias
        nrows,ncols = size(ypred)
        index = Int[]
        for b=1:B
            if mask[t][b]==1
                push!(index, idx(t,b,nrows))
            end
        end
        o1 = logp(ypred, 1)
        o2 = o1[index]
        total += sum(o2)
        count += length(o2)
    end

    if result != nothing
        result[1] += AutoGrad.getval(total)
        result[2] += count
    end
    return total
end

max_stars_repo_name: kirnap/KU-dependencyParser

max_issues_repo_name: kirnap/KU-dependencyParser

max_forks_repo_name: kirnap/KU-dependencyParser


text:
# This analysis generates several different common plots used to
# diagnose and explore screen data including: (1) volcano plot of gene
# pvalues and effect sizes (2) precision-recall curve for the hits
# (3) raw count plots
using Gadfly
using ColorBrewer

function gen_plots(filepath; debug=false, quiet=false)
    colors = palette("Accent", 4)

    for screentype in [FacsScreen(), GrowthScreen()]
        for crisprtype in [CRISPRi(), CRISPRn()]
            if typeof(screentype) == GrowthScreen
                max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
                    :inactive => (0.83, Delta(0.0)),
                    :negcontrol => (0.05, Delta(0.0)),
                    :increasing => (0.02, TruncatedNormal(0.1, 0.1, 0.025, 1)),
                    :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
                )
                lib = Library(max_phenotype_dists, crisprtype)
            else
                lib = Library(crisprtype)
            end

            # redefining function each time and abusing lexical scoping
            gen_plots = (bc_counts, genes) -> begin
                curr_dir = dirname(filepath)
                new_filename = joinpath(curr_dir, "$(basename(filepath))_$(typeof(crisprtype))_$(typeof(screentype))_counts.svg")
                nopseudo = bc_counts[(bc_counts[:counts_bin1] .> 0.5) .& (bc_counts[:counts_bin2] .> 0.5), :]
                draw(SVG(new_filename, 10cm, 10cm), plot(nopseudo,
                x=:counts_bin1, y=:counts_bin2, color=:class, Scale.x_log10,
                Scale.y_log10, Scale.color_discrete_manual(colors..., levels=sort(unique(nopseudo[:class]))),
                Theme(highlight_width=0pt),
                Coord.cartesian(fixed=true),
                Guide.title("$(typeof(crisprtype)) $(typeof(screentype))"),
                Guide.xlabel("log counts t0"), Guide.ylabel("log counts t1")))
                new_filename = joinpath(curr_dir, "$(basename(filepath))_$(typeof(crisprtype))_$(typeof(screentype))_volcano.svg")
                draw(SVG(new_filename, 10cm, 10cm), plot(genes,
                x=:mean_bin2_div_bin1, y=:pvalue_bin2_div_bin1, color=:class,
                Scale.color_discrete_manual(colors..., levels=sort(unique(genes[:class]))),
                Theme(highlight_width=0pt),
                Guide.title("$(typeof(crisprtype)) $(typeof(screentype))"),
                Guide.xlabel("mean log2 fold change"), Guide.ylabel("-log10 pvalue")))
                (1)
            end

            run_exp(screentype, lib, gen_plots; run_idx=1)
        end
    end
end


text:
"""
Quadratic objective
    x → 0.5 xᵀQx - bᵀx
"""
struct Quadratic{matT,vecT}
    Q::matT
    b::vecT
end

"Evaluate the quadratic objective."
(q::Quadratic)(x) = 0.5*dot(x,q.Q*x) - dot(q.b,x)

"Evaluate the quaratic gradient."
grad(q::Quadratic, x) = q.Q*x - q.b

struct LBFGSBQuad
    optimizer::L_BFGS_B
    bounds::Matrix{Float64}
    Q::Quadratic
    bl::Vector{Float64}
    bu::Vector{Float64}
    x0::Vector{Float64}
    m::Int64
    function LBFGSBQuad(Q::Quadratic, bl::Vector, bu::Vector; m::Integer=5)
        n = length(bl)
        bounds = Matrix{Float64}(undef, 3, n)
        optimizer = L_BFGS_B(n, m)
        bounds[1,:] = fill(2, 1, length(bl))
        bounds[2,:] = bl'
        bounds[3,:] = bu'
        x0 = (bl+bu)/2
        new(optimizer, bounds, Q, bl, bu, x0, m)
    end
end

""" Solve min_{l ≤ x ≤ u} 0.5x^TQx - b^Tx using LBFGSB."""
function boxquad(lss::LBFGSBQuad;
                 factr=1.0, pgtol=1e-6, iprint=-1,
                 maxfun=15000, maxiter=15000)
    Q, bounds, bl, bu, x0 = lss.Q, lss.bounds, lss.bl, lss.bu, lss.x0
    f, x = lss.optimizer(Q, (z,x) -> z.=grad(Q, x), x0, bounds,
                         m=lss.m, factr=factr, pgtol=pgtol, iprint=iprint,
                         maxfun=maxfun, maxiter=maxiter)
    return x
end

"""
    Solve the QP
    0.5*x'*Q*x + b'*x s.t. x_i = v_i for i in mask.
"""
function _solve_eq_qp(Q::Quadratic, mask, values)
    Q_new = _fix_values(Q, mask, values)
    x_free = Q_new.Q\Q_new.b
    x = deepcopy(values)
    ind_free = 1
    for i in eachindex(x)
        if !mask[i]
            x[i] = x_free[ind_free]
            ind_free+=1
        end
    end
    x
end

function _fix_values(Q_st::Quadratic, mask, values)
     n = length(Q_st.b)
     Q = Q_st.Q
     not_mask = (!).(mask)
     Q_new = Q[not_mask, not_mask]
     b_new = Q_st.b[not_mask]- Q[not_mask,:]*values #?????
     Quadratic(Q_new, b_new)
end

"""Attempts to solve min 0.5 x^TQx +b^Tx s.t. l ≤ x ≤ u by (repeatedly) guessing the active set."""
function active_set_min_bound_constrained_quadratic(Q_st :: Quadratic, l, u, max_iters = 5, tol=1E-10; mask_l = falses(length(l)), mask_u = falses(length(u)))
    Q,b = Q_st.Q, -Q_st.b
    n = length(b)
    eq_mask = falses(n)
    values = zeros(n)
    #if mask_l[i], x_i is constrained to be l_i
    #if mask_u[i], x_i is constrained to be u_i

    for iter in 1:max_iters
        # set up and solve equality constrained QP

        eq_mask .= false
        values .= 0.0
        for i in 1:n
            if mask_u[i]
                eq_mask[i] = true
                values[i] = u[i]
            elseif mask_l[i]
                eq_mask[i] = true
                values[i] = l[i]
            end
        end

        x = _solve_eq_qp(Q_st, eq_mask, values)

        #compute gradient
        g = grad(Q_st, x)
        #check kkt conditions
        if all(((x,l,u),) -> l-tol < x < u+tol, zip(x, l, u)) # primal feasibility
            ### compute dual variables by enforcing stationarity...
            # and check complementary slackness
            comp_slack = true
            for i in eachindex(x)
                if g[i] < -tol
                    if abs(x[i] - u[i]) > tol
                        comp_slack = false
                        mask_u[i] = false
                        mask_l[i] = false
                    end
                elseif g[i] > tol
                    if abs(x[i] - l[i]) > tol
                        comp_slack = false
                        mask_u[i] = false
                        mask_l[i] = false
                    end
                end
            end
            if comp_slack
                return (x, iter)
            end
        else # Not primal feasible, add constraints...
            for i in 1:n
                if x[i] < l[i] - tol
                    mask_l[i] = true
                    mask_u[i] = false
                end
                if x[i] > u[i] + tol
                    mask_u[i] = true
                    mask_l[i] = false
                end
            end
        end
    end
    return (zeros(n), 0)
end

function form_quadratic_from_least_squares(A,b)
    Quadratic(A'*A, A'*b)
end


text:
mutable struct Variable{T}
    value :: T
    dvalue :: T
end

function Variable(x)
    return Variable(x, zero(x))
end

reset!(v::Variable) = (v.dvalue = zero(v.value))
update!(v::Variable) = (v.value += v.dvalue; reset!(v))
update!(v::Variable{<:Array}) = v.value .+= v.dvalue

using Tensors
a = 1.0
b = [1.0, 2.0, 3.0]
c = Tensor{2, 3}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
vara = Variable(a)
varb = Variable(b)
varc = Variable(c)

@info "Initial state: $vara"
vara.dvalue += rand()
@info "After setting dvalue: $vara"
update!(vara)
@info "After update!: $vara"

@info "Initial state: $varb"
varb.dvalue += rand(3)
@info "After setting dvalue: $varb"
update!(varb)
@info "After update!: $varb"

@info "Initial state: $varc"
varc.dvalue += Tensor{2,3}(rand(9))
@info "After setting dvalue: $varc"
update!(varc)
@info "After update!: $varc"

using BenchmarkTools
N = 1000
function bench_float64()
    # Random walk test"
    var = Variable(1.0)
    for i in 1:N
        var.dvalue += randn()
        update!(var)
    end
    return var
end
function bench_array()
    # Random walk test
    var = Variable([1.0, 2.0, 3.0])
    for i in 1:N
        var.dvalue += randn(3)
        update!(var)
    end
    return var
end
function bench_tensor()
    # Random walk test
    var = Variable(Tensor{2, 3}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]))
    for i in 1:N
        var.dvalue += randn(Tensor{2,3})
        update!(var)
    end
end
function bench_symtensor()
    # Random walk test
    var = Variable(SymmetricTensor{2, 3}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]))
    for i in 1:N
        var.dvalue += randn(SymmetricTensor{2,3,Float64})
        update!(var)
    end
end

# println("Benchmark Variable{Float64}")
# @btime bench_float64()
# println("Benchmark Variable{Array{Float64,1}}")
# @btime bench_array()
# println("Benchmark Variable{Tensor{2,3,Float64,9}}")
# @btime bench_tensor()
# println("Benchmark Variable{SymmetricTensor{2,3,Float64,6}}")
# @btime bench_symtensor()

abstract type AbstractVariableState end

mutable struct TestState <: AbstractVariableState
    x :: Variable{Float64}
end

mutable struct VariableState <: AbstractVariableState
   stress::Variable{SymmetricTensor{2,3,Float64,6}}
   strain::Variable{SymmetricTensor{2,3,Float64,6}}
   backstress1::Variable{SymmetricTensor{2,3,Float64,6}}
   backstress2::Variable{SymmetricTensor{2,3,Float64,6}}
   plastic_strain::Variable{SymmetricTensor{2,3,Float64,6}}
   cumeq::Variable{Float64}
   R::Variable{Float64}
end

function update!(state::T) where {T<:AbstractVariableState}
    for fn in fieldnames(T)
        update!(getfield(state, fn))
    end
end

function bench_chaboche_style_variablestate()
    stress = Variable(zero(SymmetricTensor{2,3}))
    strain = Variable(zero(SymmetricTensor{2,3}))
    backstress1 = Variable(zero(SymmetricTensor{2,3}))
    backstress2 = Variable(zero(SymmetricTensor{2,3}))
    plastic_strain = Variable(zero(SymmetricTensor{2,3}))
    cumeq = Variable(0.0)
    R = Variable(0.0)
    state = VariableState(stress, strain, backstress1,
                        backstress2, plastic_strain, cumeq, R)
    for i in 1:N
        state.stress.dvalue = randn(SymmetricTensor{2,3})
        state.strain.dvalue = randn(SymmetricTensor{2,3})
        state.backstress1.dvalue = randn(SymmetricTensor{2,3})
        state.backstress2.dvalue = randn(SymmetricTensor{2,3})
        state.plastic_strain.dvalue = randn(SymmetricTensor{2,3})
        state.cumeq.dvalue = norm(state.plastic_strain.dvalue)
        state.R.dvalue = randn()
        update!(state)
    end
    return state
end

println("Benchmark Chaboche VariableState")
@btime bench_chaboche_style_variablestate()


text:
@deprecate sampleplot(gp::FiniteGP, n::Int; kwargs...) sampleplot(gp; samples=n, kwargs...)
@deprecate sampleplot!(gp::FiniteGP, n::Int; kwargs...) sampleplot!(
    gp; samples=n, kwargs...
)
@deprecate sampleplot!(plt::RecipesBase.AbstractPlot, gp::FiniteGP, n::Int; kwargs...) sampleplot!(
    plt, gp; samples=n, kwargs...
)


text:
# function NSDEBase.cost(problem::AbstractInitialValueProblem, solver::MovingWindow.AbstractMovingWindowSolver, solution::MovingWindow.AbstractMovingWindowSolution)
#     @↓ (t0, tN) ← tspan = problem
#     @↓ τ, Δτ, parallelsolver, budget = solver
#     @↓ finesolver, coarsolver,  P = parallelsolver
#     @↓ hF ← h = finesolver.stepsize
#     @↓ hG ← h = coarsolver.stepsize
#     ξ = hG / hF
#     @↓ B ← budget, a₁, a₂, a₃ = budget
#     M = length(solution)
#     r = solution.restarts
#     α = τ / Δτ
#     return τ / abs(tN - t0) * (M*(M-1)/(2α*ξ) + (P+ξ)/(2P*ξ)*sum(1 + (B+1)*r[m] + length(solution[m].errors) for m=1:M))
# end
# cost(problem, mowa, mowasolution) * Cₛ


text:
using Zygote, Test

function f(x)
  ch = Channel(Inf)
  put!(ch, x^2) + take!(ch)
end

@test gradient(f, 5) == (20,)


text:
module ConditionalLinearGaussianDrivers

end # module


text:
module MicrobiomePlots

using Microbiome
using RecipesBase
using Colors
import StatPlots: GroupedBar

export
    abundanceplot,
    AnnotationBar,
    annotationbar
    # zeroyplot

include("recipes.jl")

end


text:
module Covid19India

using HTTP, JSON
using DataFrames, CSV
using Dates
include("urls.jl")
include("state_ids.jl")
include("helper_functions.jl")
include("timeseriesIndia.jl")
include("statewise.jl")
include("districtwise.jl")
include("countrywise.jl")

end # module