text:
const SupervisedNetwork = Union{DeterministicNetwork,ProbabilisticNetwork}

const GenericNetwork = Union{SupervisedNetwork,UnsupervisedNetwork}

# to suppress inclusion in models():
MMI.is_wrapper(::Type{DeterministicNetwork}) = true
MMI.is_wrapper(::Type{ProbabilisticNetwork}) = true

## FALL-BACKS FOR LEARNING NETWORKS EXPORTED AS MODELS

function update(model::GenericNetwork, verb::Integer, yhat, cache, args...)

    # If any `model` field has been replaced (and not just mutated)
    # then we actually need to fit rather than update (which will
    # force build of a new learning network). If `model` has been
    # created using a learning network export macro, the test used
    # below is perfect. In any other case it is at least conservative:
    network_model_ids = objectid.(models(yhat))
    fields = [getproperty(model, name) for name in fieldnames(typeof(model))]
    submodels    = filter(f->f isa Model, fields)
    submodel_ids = objectid.(submodels)
    if !issubset(submodel_ids, network_model_ids)
        return fit(model, verb, args...)
    end

    is_anonymised = cache isa NamedTuple{(:sources, :data)}

    if is_anonymised
        sources, data = cache.sources, cache.data
        for k in eachindex(sources)
            rebind!(sources[k], data[k])
        end
    end

    fit!(yhat; verbosity=verb)
    if is_anonymised
        for s in sources
            rebind!(s, nothing)
        end
    end

    return yhat, cache, nothing
end

predict(::SupervisedNetwork, fitres, Xnew)     = fitres(Xnew)

transform(::UnsupervisedNetwork, fitres, Xnew) = fitres(Xnew)

function fitted_params(yhat::Node)
    machs  = machines(yhat)
    fitted = [fitted_params(m) for m in machs]
    return (machines=machs, fitted_params=fitted)
end

fitted_params(::GenericNetwork, yhat) = fitted_params(yhat)

## FOR EXPORTING LEARNING NETWORKS BY HAND

"""
    anonymize!(sources...)

Returns a named tuple `(sources=..., data=....)` whose values are the
provided source nodes and their contents respectively, and clears the
contents of those source nodes.

"""
function anonymize!(sources...)
    data = Tuple(s.data for s in sources)
    [MLJBase.rebind!(s, nothing) for s in sources]
    return (sources=sources, data=data)
end

function MLJBase.report(yhat::Node)
    machs = machines(yhat)
    reports = [report(m) for m in machs]
    return (machines=machs, reports=reports)
end

# what should be returned by a fit method for an exported learning
# network:
function fitresults(yhat)
    inputs = sources(yhat, kind=:input)
    targets = sources(yhat, kind=:target)
    weights = sources(yhat, kind=:weights)

    length(inputs) == 1 ||
        error("Improperly exported supervised network does "*
              "not have a unique :input source. ")
    length(targets) < 2 ||
        error("Improperly exported network has multiple :target sources. ")
    length(weights) < 2 ||
        error("Improperly exported network has multiple :weights sources. ")

    cache = anonymize!(vcat(inputs, targets, weights)...)

    r = report(yhat)

    return yhat, cache, r
end

## EXPORTING LEARNING NETWORKS AS MODELS WITH @from_network

"""
    replace(W::Node, a1=>b1, a2=>b2, ...; empty_unspecified_sources=false)

Create a deep copy of a node `W`, and thereby replicate the learning
network terminating at `W`, but replacing any specified sources and
models `a1, a2, ...` of the original network with `b1, b2, ...`.

If `empty_unspecified_sources=ture` then any source nodes not
specified are replaced with empty version of the same kind.

"""
function Base.replace(W::Node, pairs::Pair...; empty_unspecified_sources=false)

    # Note: We construct nodes of the new network as values of a
    # dictionary keyed on the nodes of the old network. Additionally,
    # there are dictionaries of models keyed on old models and
    # machines keyed on old machines. The node and machine
    # dictionaries must be built simultaneously.

    # build model dict:
    model_pairs = filter(collect(pairs)) do pair
        first(pair) isa Model
    end
    models_ = models(W)
    models_to_copy = setdiff(models_, first.(model_pairs))
    model_copy_pairs = [model=>deepcopy(model) for model in models_to_copy]
    newmodel_given_old = IdDict(vcat(model_pairs, model_copy_pairs))

    # build complete source replacement pairs:
    sources_ = sources(W)
    specified_source_pairs = filter(collect(pairs)) do pair
        first(pair) isa Source
    end
    unspecified_sources = setdiff(sources_, first.(specified_source_pairs))
    unspecified_sources_wrapping_something =
        filter(s -> !isempty(s), unspecified_sources)
    if !isempty(unspecified_sources_wrapping_something) &&
        !empty_unspecified_sources
        @warn "No replacement specified for one or more non-empty source "*
        "nodes. Contents will be duplicated. "
    end
    if empty_unspecified_sources
        unspecified_source_pairs = [s => source(kind=MLJBase.kind(s)) for
                                    s in unspecified_sources]
    else
        unspecified_source_pairs = [s => deepcopy(s) for
                                    s in unspecified_sources]
    end

    all_source_pairs = vcat(specified_source_pairs, unspecified_source_pairs)

    # drop source nodes from all nodes of network terminating at W:
    nodes_ = filter(nodes(W)) do N
        !(N isa Source)
    end
    isempty(nodes_) && error("All nodes in network are source nodes. ")
    # instantiate node and machine dictionaries:
    newnode_given_old =
        IdDict{AbstractNode,AbstractNode}(all_source_pairs)
    newmach_given_old = IdDict{NodalMachine,NodalMachine}()

    # build the new network:
    for N in nodes_
       args = [newnode_given_old[arg] for arg in N.args]
         if N.machine === nothing
             newnode_given_old[N] = node(N.operation, args...)
         else
             if N.machine in keys(newmach_given_old)
                 mach = newmach_given_old[N.machine]
             else
                 train_args = [newnode_given_old[arg] for arg in N.machine.args]
                 mach = machine(newmodel_given_old[N.machine.model], train_args...)
                 newmach_given_old[N.machine] = mach
             end
             newnode_given_old[N] = N.operation(mach, args...)
        end
    end

    return newnode_given_old[nodes_[end]]

 end

# closure for later:
function fit_method(network, models...)

    network_Xs = sources(network, kind=:input)[1]

    function _fit(model::M, verb::Integer, args...) where M <: Supervised
        replacement_models = [getproperty(model, fld)
                              for fld in fieldnames(M)]
        model_replacements = [models[j] => replacement_models[j]
                              for j in eachindex(models)]
        network_ys = sources(network, kind=:target)[1]
        Xs = source(args[1])
        ys = source(args[2], kind=:target)
        source_replacements = [network_Xs => Xs, network_ys => ys]
        if length(args) == 3
            ws = source(args[3], kind=:weights)
            network_ws = sources(network, kind=:weights)[1]
            push!(source_replacements, network_ws => ws)
        end
        replacements = vcat(model_replacements, source_replacements)
        yhat = replace(network, replacements...; empty_unspecified_sources=true)

        if length(args) == 2
            issubset([Xs, ys], sources(yhat)) ||
                error("Failed to replace learning network "*
                      ":input or :target source")
        elseif length(args) == 3
            issubset([Xs, ys, ws], sources(yhat)) ||
                error("Failed to replace learning network "*
                      ":input or :target source")
        else
            throw(ArgumentError)
        end

        fit!(yhat, verbosity=verb)

        return fitresults(yhat)
    end

    function _fit(model::M, verb::Integer, X) where M <:Unsupervised
        replacement_models = [getproperty(model, fld)
                              for fld in fieldnames(M)]
        model_replacements = [models[j] => replacement_models[j]
                          for j in eachindex(models)]
        Xs = source(X)
        source_replacements = [network_Xs => Xs,]
        replacements = vcat(model_replacements, source_replacements)
        Xout = replace(network, replacements...)
        Set([Xs]) == Set(sources(Xout)) ||
            error("Failed to replace learning network :input source ")

        fit!(Xout, verbosity=verb)

        return fitresults(Xout)
    end

    return _fit
end

net_alert(message) = throw(ArgumentError("Learning network export error.\n"*
                                     string(message)))
net_alert(k::Int) = throw(ArgumentError("Learning network export error $k. "))

# returns Model supertype - or `missing` if arguments are incompatible
function kind_(is_supervised, is_probabilistic)
    if is_supervised
        if ismissing(is_probabilistic) || !is_probabilistic
            return :DeterministicNetwork
        else
            return :ProbabilisticNetwork
        end
    else
        if ismissing(is_probabilistic) || !is_probabilistic
            return :UnsupervisedNetwork
        else
            return missing
        end
    end
end

function from_network_preprocess(modl, ex,
                                 is_probabilistic::Union{Missing,Bool})

    trait_ex_given_name_ex = Dict{Symbol,Any}()

    ex isa Expr || net_alert(1)
    ex.head == :call || net_alert(2)
    ex.args[1] == :(<=) || net_alert(3)
    ex.args[2] isa Expr || net_alert(4)
    ex.args[2].head == :call || net_alert(5)
    modeltype_ex = ex.args[2].args[1]
    modeltype_ex isa Symbol || net_alert(6)
    if length(ex.args[2].args) == 1
        kw_exs = []
    else
        kw_exs = ex.args[2].args[2:end]
    end
    fieldname_exs = []
    model_exs = []
    for ex in kw_exs
        ex isa Expr || net_alert(7)
        ex.head == :kw || net_alert(8)
        variable_ex = ex.args[1]
        value_ex = ex.args[2]
        variable_ex isa Symbol || net_alert(9)
        push!(fieldname_exs, variable_ex)
        value = modl.eval(value_ex)
        value isa Model ||
            net_alert("Got $value but expected something of type `Model`.")
        push!(model_exs, value_ex)
    end
    N_ex = ex.args[3]
    N = modl.eval(N_ex)
    N isa AbstractNode ||
        net_alert("Got $N but expected something of type `AbstractNode`. ")

    inputs = sources(N, kind=:input)
    targets = sources(N, kind=:target)
    weights = sources(N, kind=:weights)

    length(inputs) == 0 &&
        net_alert("Network has no source with `kind=:input`.")
    length(inputs) > 1  &&
        net_alert("Network has multiple sources with `kind=:input`.")
    length(targets) > 1 &&
        net_alert("Network has multiple sources with `kind=:target`.")
    length(weights) > 1 &&
        net_alert("Network has multiple sources with `kind=:weights`.")

    if length(weights) == 1
        trait_ex_given_name_ex[:supports_weights] = true
    end

    is_supervised = length(targets) == 1

    kind = kind_(is_supervised, is_probabilistic)
    ismissing(kind) &&
        net_alert("Network appears unsupervised (has no source with "*
                  "`kind=:target`) and so `is_probabilistic=true` "*
                  "declaration is not allowed. ")

    models_ = [modl.eval(e) for e in model_exs]
    issubset(models_, models(N)) ||
        net_alert("One or more specified models are not in the learning network "*
                  "terminating at $N_ex.\n Use models($N_ex) to inspect models. ")

    nodes_  = nodes(N)

    return modeltype_ex, fieldname_exs, model_exs, N_ex,
           kind, trait_ex_given_name_ex

end

from_network_preprocess(modl, ex) = from_network_preprocess(modl, ex, missing)

function from_network_preprocess(modl, ex, kw_ex)
    kw_ex isa Expr || net_alert(10)
    kw_ex.head == :(=) || net_alert(11)
    if kw_ex.args[1] == :is_probabilistic
        value = kw_ex.args[2]
        if value isa Bool
            return from_network_preprocess(modl, ex, value)
        else
            net_alert("`is_probabilistic` can only be `true` or `false`.")
        end
    elseif kw_ex.args[1] == :prediction_type
        value = kw_ex.args[2].value
        if value == :probabilistic
            return from_network_preprocess(modl, ex, true)
        elseif value == :deterministic
            return from_network_preprocess(modl, ex, false)
        else
            net_alert("`prediction_type` can only be `:probabilistic` "*
                      "or `:deterministic`.")
        end
    else
        net_alert("Unrecognized keywork `$(kw_ex.args[1])`.")
    end
end

function from_network_(modl, modeltype_ex, fieldname_exs, model_exs,
                          N_ex, kind, trait_value_ex_given_name_ex)

    args = gensym(:args)

    # code defining the composite model struct and fit method:
    program1 = quote

        import MLJBase
        import MLJModelInterface

        mutable struct $modeltype_ex <: MLJBase.$kind
            $(fieldname_exs...)
        end

        MLJModelInterface.fit(model::$modeltype_ex, verb::Integer, $args...) =
            MLJBase.fit_method($N_ex, $(model_exs...))(model, verb, $args...)

    end

    program2 = quote
        # defined keyword constructor for composite model:
        defaults =
            MLJBase.@set_defaults $modeltype_ex deepcopy.([$(model_exs...)])
    end

    modl.eval(program1)

    # define composite model traits:
    for (name_ex, value_ex) in trait_value_ex_given_name_ex
        program = quote
            MLJBase.$name_ex(::Type{<:$modeltype_ex}) = $value_ex
        end
        modl.eval(program)
    end

    modl.eval(program2)

end

"""
    @from_network(NewCompositeModel(fld1=model1, fld2=model2, ...) <= N
    @from_network(NewCompositeModel(fld1=model1, fld2=model2, ...) <= N is_probabilistic=false

Create a new stand-alone model type called `NewCompositeModel`, using
a learning network as a blueprint. Here `N` refers to the terminal
node of the learning network (from which final predictions or
transformations are fetched).

*Important.* If the learning network is supervised (has a source with
`kind=:target`) and makes probabilistic predictions, then one must
declare `is_probabilistic=true`. In the deterministic case the keyword
argument can be omitted.

The model type `NewCompositeModel` is equipped with fields named
`:fld1`, `:fld2`, ..., which correspond to component models `model1`,
`model2`, ...,  appearing in the network (which must therefore be elements of
`models(N)`).  Deep copies of the specified component models are used
as default values in an automatically generated keyword constructor
for `NewCompositeModel`.

### Return value

 A new `NewCompositeModel` instance, with default field values.

For details and examples refer to the "Learning Networks" section of
the documentation.

"""
macro from_network(exs...)

    args = from_network_preprocess(__module__, exs...)
    modeltype_ex = args[1]

    from_network_(__module__, args...)

    esc(quote
        $modeltype_ex()
        end)

end


text:
#------------------------------------------------------------------------------------
import Pkg
Pkg.instantiate()
Pkg.update()

# high level CSV API
# https://juliadata.github.io/CSV.jl/stable/index.html
# https://www.geeksforgeeks.org/working-with-csv-files-in-julia/
Pkg.add("CSV")
import CSV
using CSV
#------------------------------------------------------------------------------------

file_path_csv = "../../../../../../externals/Core.Math.Data/data/HolisticWare.Core.Sample.Data/CSV/Pejcic_318.csv"

data = CSV.read(file_path_csv, copycols = true)

println(size(data),"\n")
println(names(data),"\n")
println(first(data, 6),"\n")
# println(describe(data),"\n")

Pkg.add("FileIO")
Pkg.add("CSVFiles")
using FileIO, CSVFiles, DataFrames

data_frame_csv_api_1 = load(file_path_csv, spacedelim=true, header_exists=false) |> DataFrame

println(size(data_frame_csv_api_1),"\n")
println(names(data_frame_csv_api_1),"\n")
println(first(data_frame_csv_api_1, 6),"\n")
# println(describe(data_frame_csv_api),"\n")

data_frame_csv_api_2 = load(file_path_csv) |> DataFrame

println(size(data_frame_csv_api_2),"\n")
println(names(data_frame_csv_api_2),"\n")
println(first(data_frame_csv_api_2, 6),"\n")
# println(describe(data_frame_csv_api),"\n")

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

ext: jl


text:
# Right division
@inline /(tvec::TransposedVector, mat::AbstractMatrix) = transpose(transpose(mat) \ transpose(tvec))

# Base only has specializations of A_rdiv_Bc and A_rdiv_Bt
@inline A_rdiv_Bt(tvec::TransposedVector, mat::AbstractMatrix) = transpose(mat \ transpose(tvec))

@inline A_rdiv_Bc(tvec::TransposedVector, mat::AbstractMatrix) = ctranspose(mat  \ ctranspose(tvec))

# These others fall back to applying the transpose, which is perfect for us:
#   Ac_rdiv_B
#   Ac_rdiv_Bc
#   At_rdiv_B
#   At_rdiv_Bt


text:
@testset "io" begin
    folder = joinpath(pkgdir(SnowyOwl), "test", "5e5a9f34-734a-4501-920c-3ecf5e2e2f36.mtx")
    mtx = read_mtx(joinpath(folder, "matrix.mtx.gz"))
    features = read_genes(joinpath(folder, "genes.tsv.gz"))
    cells = read_cells(joinpath(folder, "cells.tsv.gz"))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Test

@testset "Check that non-floats are correctly promoted" begin
    @test [1 0 0; 0 1 0]\[1,1] ≈ [1;1;0]
end

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

@testset "Matrix condition number" begin
    ainit = rand(n,n)
    @testset "for $elty" for elty in (Float32, Float64, Complex64, Complex128)
        ainit = convert(Matrix{elty}, ainit)
        for arraytype in ("Array", "SubArray")
            if arraytype == "Array"
                a = ainit
            else
                a = view(ainit, 1:n, 1:n)
            end
            @test cond(a,1) ≈ 4.837320054554436e+02 atol=0.01
            @test cond(a,2) ≈ 1.960057871514615e+02 atol=0.01
            @test cond(a,Inf) ≈ 3.757017682707787e+02 atol=0.01
            @test cond(a[:,1:5]) ≈ 10.233059337453463 atol=0.01
            @test_throws ArgumentError cond(a,3)
        end
    end
end

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

@testset "For A containing $eltya" for eltya in (Float32, Float64, Complex64, Complex128, Int)
    ainit = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    ainit2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    ε = εa = eps(abs(float(one(eltya))))

    apd  = ainit'*ainit # symmetric positive-definite
    @testset "Positive definiteness" begin
        @test !isposdef(ainit)
        @test isposdef(apd)
        if eltya != Int # cannot perform cholfact! for Matrix{Int}
            @test !isposdef!(copy(ainit))
            @test isposdef!(copy(apd))
        end
    end
    @testset "For b containing $eltyb" for eltyb in (Float32, Float64, Complex64, Complex128, Int)
        binit = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        for arraytype in ("Array", "SubArray")
            if arraytype == "Array"
                a = ainit
                b = binit
            else
                a = view(ainit, 1:n, 1:n)
                b = view(binit, 1:n, 1:2)
            end

            @testset "Solve square general system of equations" begin
                κ = cond(a,1)
                x = a \ b
                @test_throws DimensionMismatch b'\b
                @test_throws DimensionMismatch b\b'
                @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit!
                @test zeros(eltya,n)\ones(eltya,n) ≈ (zeros(eltya,n,1)\ones(eltya,n,1))[1,1]
            end

            @testset "Test nullspace" begin
                a15null = nullspace(a[:,1:n1]')
                @test rank([a[:,1:n1] a15null]) == 10
                @test norm(a[:,1:n1]'a15null,Inf) ≈ zero(eltya) atol=300ε
                @test norm(a15null'a[:,1:n1],Inf) ≈ zero(eltya) atol=400ε
                @test size(nullspace(b), 2) == 0
                @test nullspace(zeros(eltya,n)) == eye(eltya,1)
            end
        end
    end # for eltyb

    for arraytype in ("Array", "SubArray")
        if arraytype == "Array"
            a = ainit
            a2 = ainit2
        else
            a = view(ainit, 1:n, 1:n)
            a2 = view(ainit2, 1:n, 1:n)
        end

        @testset "Test pinv" begin
            pinva15 = pinv(a[:,1:n1])
            @test a[:,1:n1]*pinva15*a[:,1:n1] ≈ a[:,1:n1]
            @test pinva15*a[:,1:n1]*pinva15 ≈ pinva15

            @test size(pinv(ones(eltya,0,0))) == (0,0)
        end

        @testset "Lyapunov/Sylvester" begin
            x = lyap(a, a2)
            @test -a2 ≈ a*x + x*a'
            x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n])
            @test -a2[1:3, 4:n] ≈ a[1:3, 1:3]*x2 + x2*a[4:n, 4:n]
        end

        @testset "Matrix square root" begin
            asq = sqrtm(a)
            @test asq*asq ≈ a
            asym = a'+a # symmetric indefinite
            asymsq = sqrtm(asym)
            @test asymsq*asymsq ≈ asym
        end

        @testset "Powers" begin
            if eltya <: AbstractFloat
                z = zero(eltya)
                t = convert(eltya,2)
                r = convert(eltya,2.5)
                @test a^z ≈ eye(a)
                @test a^t ≈ a^2
                @test eye(eltya,n,n)^r ≈ eye(a)
            end
        end
    end # end for loop over arraytype

    @testset "Numbers" begin
        α = rand(eltya)
        A = zeros(eltya,1,1)
        A[1,1] = α
        @test diagm(α) == A # Test behavior of `diagm` when passed a scalar
        @test expm(α) == exp(α) # `expm` should behave like `exp` with scalar argument
    end

    @testset "Factorize" begin
        d = rand(eltya,n)
        e = rand(eltya,n-1)
        e2 = rand(eltya,n-1)
        f = rand(eltya,n-2)
        A = diagm(d)
        @test factorize(A) == Diagonal(d)
        A += diagm(e,-1)
        @test factorize(A) == Bidiagonal(d,e,:L)
        A += diagm(f,-2)
        @test factorize(A) == LowerTriangular(A)
        A = diagm(d) + diagm(e,1)
        @test factorize(A) == Bidiagonal(d,e,:U)
        if eltya <: Real
            A = diagm(d) + diagm(e,1) + diagm(e,-1)
            @test full(factorize(A)) ≈ full(factorize(SymTridiagonal(d,e)))
            A = diagm(d) + diagm(e,1) + diagm(e,-1) + diagm(f,2) + diagm(f,-2)
            @test inv(factorize(A)) ≈ inv(factorize(Symmetric(A)))
        end
        A = diagm(d) + diagm(e,1) + diagm(e2,-1)
        @test full(factorize(A)) ≈ full(factorize(Tridiagonal(e2,d,e)))
        A = diagm(d) + diagm(e,1) + diagm(f,2)
        @test factorize(A) == UpperTriangular(A)
    end
end # for eltya

@testset "test triu/tril bounds checking" begin
    ainit = rand(5,7)
    for arraytype in ("Array", "SubArray")
        if arraytype == "Array"
            a = ainit
        else
            a = view(ainit, 1:size(ainit, 1), 1:size(ainit, 2))
        end
        @test_throws(ArgumentError,triu(a,8))
        @test_throws(ArgumentError,triu(a,-6))
        @test_throws(ArgumentError,tril(a,8))
        @test_throws(ArgumentError,tril(a,-6))
    end
end

@testset "Test gradient for $elty" for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128)
    if elty <: Real
        x = convert(Vector{elty}, [1:3;])
        g = ones(elty, 3)
    else
        x = convert(Vector{elty}, complex.([1:3;], [1:3;]))
        g = convert(Vector{elty}, complex.(ones(3), ones(3)))
    end
    xsub = view(x, 1:size(x, 1))
    @test gradient(x) ≈ g
    @test gradient(xsub) ≈ g # Test gradient on SubArray
    @test gradient(ones(elty,1)) == zeros(elty,1)
end

@testset "Tests norms" begin
    nnorm = 10
    mmat = 10
    nmat = 8
    @testset "For $elty" for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
        x = ones(elty,10)
        @testset "Vector" begin
            xs = view(x,1:2:10)
            @test norm(x, -Inf) ≈ 1
            @test norm(x, -1) ≈ 1/10
            @test norm(x, 0) ≈ 10
            @test norm(x, 1) ≈ 10
            @test norm(x, 2) ≈ sqrt(10)
            @test norm(x, 3) ≈ cbrt(10)
            @test norm(x, Inf) ≈ 1
            if elty <: Base.LinAlg.BlasFloat
                @test norm(x, 1:4) ≈ 2
                @test_throws BoundsError norm(x,-1:4)
                @test_throws BoundsError norm(x,1:11)
            end
            @test norm(xs, -Inf) ≈ 1
            @test norm(xs, -1) ≈ 1/5
            @test norm(xs, 0) ≈ 5
            @test norm(xs, 1) ≈ 5
            @test norm(xs, 2) ≈ sqrt(5)
            @test norm(xs, 3) ≈ cbrt(5)
            @test norm(xs, Inf) ≈ 1
        end

        @testset "Issue #12552:" begin
            if real(elty) <: AbstractFloat
                for p in [-Inf,-1,1,2,3,Inf]
                    @test isnan(norm(elty[0,NaN],p))
                    @test isnan(norm(elty[NaN,0],p))
                end
            end
        end

        @testset "Number" begin
            norm(x[1:1]) === norm(x[1], -Inf)
            norm(x[1:1]) === norm(x[1], 0)
            norm(x[1:1]) === norm(x[1], 1)
            norm(x[1:1]) === norm(x[1], 2)
            norm(x[1:1]) === norm(x[1], Inf)
        end

        @testset "Absolute homogeneity, triangle inequality, & vectorized versions" begin
            for i = 1:10
                xinit = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) :
                        elty <: Complex ? convert(Vector{elty}, complex.(randn(nnorm), randn(nnorm))) :
                        convert(Vector{elty}, randn(nnorm))
                yinit = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) :
                        elty <: Complex ? convert(Vector{elty}, complex.(randn(nnorm), randn(nnorm))) :
                        convert(Vector{elty}, randn(nnorm))
                α = elty <: Integer ? randn() :
                    elty <: Complex ? convert(elty, complex(randn(),randn())) :
                    convert(elty, randn())
                for arraytype in ("Array", "SubArray")
                    if arraytype == "Array"
                        x = xinit
                        y = yinit
                    else
                        x = view(xinit,1:2:nnorm)
                        y = view(yinit,1:2:nnorm)
                    end
                    # Absolute homogeneity
                    @test norm(α*x,-Inf) ≈ abs(α)*norm(x,-Inf)
                    @test norm(α*x,-1) ≈ abs(α)*norm(x,-1)
                    @test norm(α*x,1) ≈ abs(α)*norm(x,1)
                    @test norm(α*x) ≈ abs(α)*norm(x) # two is default
                    @test norm(α*x,3) ≈ abs(α)*norm(x,3)
                    @test norm(α*x,Inf) ≈ abs(α)*norm(x,Inf)

                    # Triangle inequality
                    @test norm(x + y,1) <= norm(x,1) + norm(y,1)
                    @test norm(x + y) <= norm(x) + norm(y) # two is default
                    @test norm(x + y,3) <= norm(x,3) + norm(y,3)
                    @test norm(x + y,Inf) <= norm(x,Inf) + norm(y,Inf)

                    # Against vectorized versions
                    @test norm(x,-Inf) ≈ minimum(abs.(x))
                    @test norm(x,-1) ≈ inv(sum(1 ./ abs.(x)))
                    @test norm(x,0) ≈ sum(x .!= 0)
                    @test norm(x,1) ≈ sum(abs.(x))
                    @test norm(x) ≈ sqrt(sum(abs2.(x)))
                    @test norm(x,3) ≈ cbrt(sum(abs.(x).^3.))
                    @test norm(x,Inf) ≈ maximum(abs.(x))
                end
            end
        end

        @testset "Matrix (Operator)" begin
            A = ones(elty,10,10)
            As = view(A,1:5,1:5)
            @test norm(A, 1) ≈ 10
            elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A, 2) ≈ 10
            @test norm(A, Inf) ≈ 10
            @test norm(As, 1) ≈ 5
            elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(As, 2) ≈ 5
            @test norm(As, Inf) ≈ 5
        end

        @testset "Absolute homogeneity, triangle inequality, & vecnorm" begin
            for i = 1:10
                Ainit = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) :
                        elty <: Complex ? convert(Matrix{elty}, complex.(randn(mmat, nmat), randn(mmat, nmat))) :
                        convert(Matrix{elty}, randn(mmat, nmat))
                Binit = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) :
                        elty <: Complex ? convert(Matrix{elty}, complex.(randn(mmat, nmat), randn(mmat, nmat))) :
                        convert(Matrix{elty}, randn(mmat, nmat))
                α = elty <: Integer ? randn() :
                    elty <: Complex ? convert(elty, complex(randn(),randn())) :
                    convert(elty, randn())

                for arraytype in ("Array", "SubArray")
                    if arraytype == "Array"
                        A = Ainit
                        B = Binit
                    else
                        A = view(Ainit,1:nmat,1:nmat)
                        B = view(Binit,1:nmat,1:nmat)
                    end

                    # Absolute homogeneity
                    @test norm(α*A,1) ≈ abs(α)*norm(A,1)
                    elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(α*A) ≈ abs(α)*norm(A) # two is default
                    @test norm(α*A,Inf) ≈ abs(α)*norm(A,Inf)

                    # Triangle inequality
                    @test norm(A + B,1) <= norm(A,1) + norm(B,1)
                    elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A + B) <= norm(A) + norm(B) # two is default
                    @test norm(A + B,Inf) <= norm(A,Inf) + norm(B,Inf)

                    # vecnorm:
                    for p = -2:3
                        @test norm(reshape(A, length(A)), p) == vecnorm(A, p)
                    end
                end
            end

            @testset "issue #10234" begin
                if elty <: AbstractFloat || elty <: Complex
                    z = zeros(elty, 100)
                    z[1] = -Inf
                    for p in [-2,-1.5,-1,-0.5,0.5,1,1.5,2,Inf]
                        @test norm(z, p) == (p < 0 ? 0 : Inf)
                        @test norm(elty[Inf],p) == Inf
                    end
                end
            end
        end
    end

    @testset "issue #10234" begin
        @test norm(Any[Inf],-2) == norm(Any[Inf],-1) == norm(Any[Inf],1) == norm(Any[Inf],1.5) == norm(Any[Inf],2) == norm(Any[Inf],Inf) == Inf
    end

    @testset "overflow/underflow in norms" begin
        @test norm(Float64[1e-300, 1], -3)*1e300 ≈ 1
        @test norm(Float64[1e300, 1], 3)*1e-300 ≈ 1
    end
end

## Issue related tests
@testset "issue #1447" begin
    A = [1.0+0.0im 0; 0 1]
    B = pinv(A)
    for i = 1:4
        @test A[i] ≈ B[i]
    end
end

@testset "issue #2246" begin
    A = [1 2 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]
    Asq = sqrtm(A)
    @test Asq*Asq ≈ A
    A2 = view(A, 1:2, 1:2)
    A2sq = sqrtm(A2)
    @test A2sq*A2sq ≈ A2

    N = 3
    @test log(det(eye(N))) ≈ logdet(eye(N))
end

@testset "issue #2637" begin
    a = [1, 2, 3]
    b = [4, 5, 6]
    @test kron(eye(2),eye(2)) == eye(4)
    @test kron(a,b) == [4,5,6,8,10,12,12,15,18]
    @test kron(a',b') == [4 5 6 8 10 12 12 15 18]
    @test kron(a,b')  == [4 5 6; 8 10 12; 12 15 18]
    @test kron(a',b)  == [4 8 12; 5 10 15; 6 12 18]
    @test kron(a,eye(2)) == [1 0; 0 1; 2 0; 0 2; 3 0; 0 3]
    @test kron(eye(2),a) == [ 1 0; 2 0; 3 0; 0 1; 0 2; 0 3]
    @test kron(eye(2),2) == 2*eye(2)
    @test kron(3,eye(3)) == 3*eye(3)
    @test kron(a,2) == [2, 4, 6]
    @test kron(b',2) == [8 10 12]
end

@testset "issue #4796" begin
    dim=2
    S=zeros(Complex,dim,dim)
    T=zeros(Complex,dim,dim)
    T[:] = 1
    z = 2.5 + 1.5im
    S[1] = z
    @test S*T == [z z; 0 0]

    # similar issue for Array{Real}
    @test Real[1 2] * Real[1.5; 2.0] == Real[5.5]
end

@testset "Matrix exponential" begin
    @testset "Tests for $elty" for elty in (Float32, Float64, Complex64, Complex128)
        A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
        eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
                                     183.765138646367 183.765138646366  163.679601723179;
                                     71.797032399996  91.8825693231832 111.968106246371]')
        @test expm(A1) ≈ eA1

        A2  = convert(Matrix{elty},
                      [29.87942128909879    0.7815750847907159 -2.289519314033932;
                       0.7815750847907159 25.72656945571064    8.680737820540137;
                       -2.289519314033932   8.680737820540137  34.39400925519054])
        eA2 = convert(Matrix{elty},
                      [  5496313853692458.0 -18231880972009236.0 -30475770808580460.0;
                       -18231880972009252.0  60605228702221920.0 101291842930249760.0;
                       -30475770808580480.0 101291842930249728.0 169294411240851968.0])
        @test expm(A2) ≈ eA2

        A3  = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52])
        eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812  -4.934938326092;
                                     0.367879439109187 1.47151775849686  1.10363831732856;
                                     0.135335281175235 0.406005843524598 0.541341126763207]')
        @test expm(A3) ≈ eA3

        A4 = convert(Matrix{elty}, [0.25 0.25; 0 0])
        eA4 = convert(Matrix{elty}, [1.2840254166877416 0.2840254166877415; 0 1])
        @test expm(A4) ≈ eA4

        A5 = convert(Matrix{elty}, [0 0.02; 0 0])
        eA5 = convert(Matrix{elty}, [1 0.02; 0 1])
        @test expm(A5) ≈ eA5

        # Hessenberg
        @test hessfact(A1)[:H] ≈ convert(Matrix{elty},
                                                 [4.000000000000000  -1.414213562373094  -1.414213562373095
                                                  -1.414213562373095   4.999999999999996  -0.000000000000000
                                                  0  -0.000000000000002   3.000000000000000])
    end

    @testset "Additional tests for $elty" for elty in (Float64, Complex{Float64})
        A4  = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
                                     1/3 1/4 1/5 1/6;
                                     1/4 1/5 1/6 1/7;
                                     1/5 1/6 1/7 1/8])
        @test expm(logm(A4)) ≈ A4

        A5  = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
        @test expm(logm(A5)) ≈ A5

        A6  = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
        @test expm(logm(A6)) ≈ A6

        A7  = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
        @test expm(logm(A7)) ≈ A7
    end

    A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]
    @test expm(logm(A8)) ≈ A8
end

@testset "issue 5116" begin
    A9  = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
    eA9 = [-0.999786072879326  -0.065407069689389   0.0   0.0
           0.006540706968939  -0.999786072879326   0.0   0.0
           0.0                 0.0                 1.0   0.0
           0.013081413937878  -3.999572145758650   0.0   1.0]
    @test expm(A9) ≈ eA9

    A10  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
    eA10 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
            0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
            0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
            0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
    @test expm(A10) ≈ eA10
end

@testset "Additional matrix logarithm tests" for elty in (Float64, Complex{Float64})
    A11 = convert(Matrix{elty}, [3 2; -5 -3])
    @test expm(logm(A11)) ≈ A11

    A12 = convert(Matrix{elty}, [1 -1; 1 -1])
    @test typeof(logm(A12)) == Array{Complex{Float64}, 2}

    A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
    logmA1 = convert(Matrix{elty}, [1.329661349 0.5302876358 -0.06818951543;
                                    0.2310490602 1.295566591 0.2651438179;
                                    0.2310490602 0.1969543025 1.363756107])
    @test logm(A1) ≈ logmA1
    @test expm(logm(A1)) ≈ A1

    A4  = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
                                 1/3 1/4 1/5 1/6;
                                 1/4 1/5 1/6 1/7;
                                 1/5 1/6 1/7 1/8])
    logmA4 = convert(Matrix{elty}, [-1.73297159 1.857349738 0.4462766564 0.2414170219;
                                    1.857349738 -5.335033737 2.994142974 0.5865285289;
                                    0.4462766564 2.994142974 -7.351095988 3.318413247;
                                    0.2414170219 0.5865285289 3.318413247 -5.444632124])
    @test logm(A4) ≈ logmA4
    @test expm(logm(A4)) ≈ A4
end

@testset "issue #7181" begin
    A = [ 1  5  9
          2  6 10
          3  7 11
          4  8 12 ]
    @test_throws ArgumentError diag(A, -5)
    @test diag(A,-4) == []
    @test diag(A,-3) == [4]
    @test diag(A,-2) == [3,8]
    @test diag(A,-1) == [2,7,12]
    @test diag(A, 0) == [1,6,11]
    @test diag(A, 1) == [5,10]
    @test diag(A, 2) == [9]
    @test diag(A, 3) == []
    @test_throws ArgumentError diag(A, 4)

    @test diag(zeros(0,0)) == []
    @test_throws ArgumentError diag(zeros(0,0),1)
    @test_throws ArgumentError diag(zeros(0,0),-1)

    @test diag(zeros(1,0)) == []
    @test diag(zeros(1,0),-1) == []
    @test_throws ArgumentError diag(zeros(1,0),1)
    @test_throws ArgumentError diag(zeros(1,0),-2)

    @test diag(zeros(0,1)) == []
    @test diag(zeros(0,1),1) == []
    @test_throws ArgumentError diag(zeros(0,1),-1)
    @test_throws ArgumentError diag(zeros(0,1),2)
end

@testset "Matrix to real power" for elty in (Float64, Complex{Float64})
# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
    #Aa : only positive real eigenvalues
    Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2])

    #Ab : both positive and negative real eigenvalues
    Ab = convert(Matrix{elty}, [1 2 3; 4 7 1; 2 1 4])

    #Ac : complex eigenvalues
    Ac = convert(Matrix{elty}, [5 4 2 1;0 1 -1 -1;-1 -1 3 6;1 1 -1 5])

    #Ad : defective Matrix
    Ad = convert(Matrix{elty}, [3 1; 0 3])

    #Ah : Hermitian Matrix
    Ah = convert(Matrix{elty}, [3 1; 1 3])
    if elty <: Base.LinAlg.BlasComplex
        Ah += [0 im; -im 0]
    end

    #ADi : Diagonal Matrix
    ADi = convert(Matrix{elty}, [3 0; 0 3])
    if elty <: Base.LinAlg.BlasComplex
        ADi += [im 0; 0 im]
    end

    for A in (Aa, Ab, Ac, Ad, Ah, ADi)
        @test A^(1/2) ≈ sqrtm(A)
        @test A^(-1/2) ≈ inv(sqrtm(A))
        @test A^(3/4) ≈ sqrtm(A) * sqrtm(sqrtm(A))
        @test A^(-3/4) ≈ inv(A) * sqrtm(sqrtm(A))
        @test A^(17/8) ≈ A^2 * sqrtm(sqrtm(sqrtm(A)))
        @test A^(-17/8) ≈ inv(A^2 * sqrtm(sqrtm(sqrtm(A))))
        @test (A^0.2)^5 ≈ A
        @test (A^(2/3))*(A^(1/3)) ≈ A
        @test (A^im)^(-im) ≈ A
    end
end

@testset "issue #23366 (Int Matrix to Int power)" begin
    @testset "Tests for $elty" for elty in (Int128, Int16, Int32, Int64, Int8,
                                            UInt128, UInt16, UInt32, UInt64, UInt8,
                                            BigInt)
        @test_throws DomainError elty[1 1;1 0]^-2
        @test (@inferred elty[1 1;1 0]^2) == elty[2 1;1 1]
        I_ = elty[1 0;0 1]
        @test I_^-1 == I_
        if !(elty<:Unsigned)
            @test (@inferred (-I_)^-1) == -I_
            @test (@inferred (-I_)^-2) == I_
        end
        # make sure that type promotion for ^(::Matrix{<:Integer}, ::Integer)
        # is analogous to type promotion for ^(::Integer, ::Integer)
        # e.g. [1 1;1 0]^big(10000) should return Matrix{BigInt}, the same
        # way as 2^big(10000) returns BigInt
        for elty2 = (Int64, BigInt)
            TT = Base.promote_op(^, elty, elty2)
            @test (@inferred elty[1 1;1 0]^elty2(1))::Matrix{TT} == [1 1;1 0]
        end
    end
end

@testset "Least squares solutions" begin
    a = [ones(20) 1:20 1:20]
    b = reshape(eye(8, 5), 20, 2)
    @testset "Tests for type $elty" for elty in (Float32, Float64, Complex64, Complex128)
        a = convert(Matrix{elty}, a)
        b = convert(Matrix{elty}, b)

        # Vector rhs
        x = a[:,1:2]\b[:,1]
        @test ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] ≈ convert(elty, 2.546616541353384)

        # Matrix rhs
        x = a[:,1:2]\b
        @test det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) ≈ convert(elty, 4.437969924812031)

        # Rank deficient
        x = a\b
        @test det((a*x-b)'*(a*x-b)) ≈ convert(elty, 4.437969924812031)

        # Underdetermined minimum norm
        x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1])
        @test x ≈ convert(Vector{elty}, [1, 0.5, -0.5])

        # symmetric, positive definite
        @test inv(convert(Matrix{elty}, [6. 2; 2 1])) ≈ convert(Matrix{elty}, [0.5 -1; -1 3])

        # symmetric, indefinite
        @test inv(convert(Matrix{elty}, [1. 2; 2 1])) ≈ convert(Matrix{elty}, [-1. 2; 2 -1]/3)
    end
end

function test_rdiv_pinv_consistency(a, b)
    @test (a*b)/b ≈ a*(b/b) ≈ (a*b)*pinv(b) ≈ a*(b*pinv(b))
    @test typeof((a*b)/b) == typeof(a*(b/b)) == typeof((a*b)*pinv(b)) == typeof(a*(b*pinv(b)))
end
function test_ldiv_pinv_consistency(a, b)
    @test a\(a*b) ≈ (a\a)*b ≈ (pinv(a)*a)*b ≈ pinv(a)*(a*b)
    @test typeof(a\(a*b)) == typeof((a\a)*b) == typeof((pinv(a)*a)*b) == typeof(pinv(a)*(a*b))
end
function test_div_pinv_consistency(a, b)
    test_rdiv_pinv_consistency(a, b)
    test_ldiv_pinv_consistency(a, b)
end

@testset "/ and \\ consistency with pinv for vectors" begin
    @testset "Tests for type $elty" for elty in (Float32, Float64, Complex64, Complex128)
        c = rand(elty, 5)
        r = rand(elty, 5)'
        cm = rand(elty, 5, 1)
        rm = rand(elty, 1, 5)
        @testset "inner prodcuts" begin
            test_div_pinv_consistency(r, c)
            test_div_pinv_consistency(rm, c)
            test_div_pinv_consistency(r, cm)
            test_div_pinv_consistency(rm, cm)
        end
        @testset "outer prodcuts" begin
            test_div_pinv_consistency(c, r)
            test_div_pinv_consistency(cm, rm)
        end
        @testset "matrix/vector" begin
            m = rand(5, 5)
            test_ldiv_pinv_consistency(m, c)
            test_rdiv_pinv_consistency(r, m)
        end
    end
end

@testset "test ops on Numbers for $elty" for elty in [Float32,Float64,Complex64,Complex128]
    a = rand(elty)
    @test expm(a) == exp(a)
    @test isposdef(one(elty))
    @test sqrtm(a) == sqrt(a)
    @test logm(a) ≈ log(a)
    @test lyap(one(elty),a) == -a/2
end

@testset "stride1" begin
    a = rand(10)
    b = view(a,2:2:10)
    @test Base.LinAlg.stride1(b) == 2
end

max_stars_repo_name: baggepinnen/julia

max_issues_repo_name: baggepinnen/julia

max_forks_repo_name: afniedermayer/julia


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

using Markdown
using InteractiveUtils

# ╔═╡ bf3615b1-3a51-4ac3-bfe1-47f56518a4f3
begin
	using CSV
	using DataFrames
	using Downloads
end

# ╔═╡ cd60cf5b-1e5d-4e3d-8171-ba87b50ad71d
using Pipe: @pipe

# ╔═╡ 62e8b182-417e-11ec-1b21-ab0009235240
md"""
# Data Dojo 6 - Lemurs

We are analyzing the data set presented by [Dr. Ludmilla Figueiredo](https://ludmillafigueiredo.github.io/) in this month's [Data Crunch Magic](https://ddojo.github.io/dcm): [Lemurs](https://github.com/rfordatascience/tidytuesday/blob/master/data/2021/2021-08-24/readme.md) from the Duke Lemur Center.

For inspiration and interesting lessons on learning R, look at [her presentation](https://github.com/ludmillafigueiredo/data_crunch_magic_wue).

Today we are using `julia` instead of `R` and we'll try to answer a couple of questions:
- Which lemur had the most offspring
- What is the mean (adult) weight per species
- What is the age distribution by species (separately for male/female)
- How do the weight curves (weight by age) differ per species?
- Which combinations of species have hybrids (and which don’t)?
- What is the name of the oldest (living) lemur?
"""

# ╔═╡ 9511af04-6f40-44ba-a253-5d439d7be550
md"""
## Loading the Data
"""

# ╔═╡ 38427d47-3e1c-4088-b835-d95327a8c526
lemurs = CSV.File(Downloads.download("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-08-24/lemur_data.csv")) |> DataFrame

# ╔═╡ 2c3491af-87ad-4bcb-b2f3-df246c761b3f
summary(lemurs)

# ╔═╡ 27a58762-55c9-429f-888b-86347e4df570
describe(lemurs)

# ╔═╡ 9b50510e-501c-4ba4-b434-e66b330fe988
md"""
## 1. Which lemur had the most offspring?

There is a column `n_known_offspring`, so we'll select the lemur (Id and Name) and that column, make entries unique and sort by it:
"""

# ╔═╡ b541f6f8-6a9b-4092-bec3-ef1cb2f220e1
to_number(x) = parse.(Int, x)

# ╔═╡ 21feee31-c860-4142-8d79-d96652041b54
	lemurs_with_offspring = @pipe lemurs |>
		select(_, :dlc_id, :name, :sex, :n_known_offspring) |>
		unique |>
		filter(:n_known_offspring => !=("NA"), _) |>
		transform(_, :n_known_offspring => (x -> parse.(Int,x)) => :n_known_offspring) |>
		sort(_, :n_known_offspring, rev=true)

# ╔═╡ 32334e5f-4b1b-457b-bb05-3074ef5dfcfd
md"
👉 The female lemur with most offsprings is named 'Thurber' and the male lemur with most offspring is named 'Dias'. Both have 36 children 😱.

🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒🐒
"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"

[compat]
CSV = "~0.9.10"
DataFrames = "~1.2.2"
Pipe = "~1.3.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "74147e877531d7c172f70b492995bc2b5ca3a843"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.10"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "5440c1d26aa29ca9ea848559216e5ee5f16a8627"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.14"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "19cb49649f8c41de7fea32d089d37de917b553da"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.0.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pipe]]
git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d"
uuid = "b98c9c47-44ae-5843-9183-064241ee97a0"
version = "1.3.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "f45b34656397a1f6e729901dc9ef679610bd12b5"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.8"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─62e8b182-417e-11ec-1b21-ab0009235240
# ╟─9511af04-6f40-44ba-a253-5d439d7be550
# ╠═bf3615b1-3a51-4ac3-bfe1-47f56518a4f3
# ╠═38427d47-3e1c-4088-b835-d95327a8c526
# ╠═2c3491af-87ad-4bcb-b2f3-df246c761b3f
# ╠═27a58762-55c9-429f-888b-86347e4df570
# ╟─9b50510e-501c-4ba4-b434-e66b330fe988
# ╠═cd60cf5b-1e5d-4e3d-8171-ba87b50ad71d
# ╠═b541f6f8-6a9b-4092-bec3-ef1cb2f220e1
# ╠═21feee31-c860-4142-8d79-d96652041b54
# ╟─32334e5f-4b1b-457b-bb05-3074ef5dfcfd
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module LinearIsotropicMod

import Base.copy

# Function sto extend
import FEM: create_matstat, stiffness, stress

# Types needed
import FEM: AbstractMaterialStatus, AbstractMaterial, GaussPoint2

export LinearIsotropicMS, LinearIsotropic

immutable LinearIsotropicMS <:AbstractMaterialStatus
    strain::Vector{Float64}
    stress::Vector{Float64}
end

LinearIsotropicMS() = LinearIsotropicMS(zeros(6), zeros(6))

copy(matstat::LinearIsotropicMS) = LinearIsotropicMS(copy(matstat.strain), copy(matstat.stress))

type LinearIsotropic <: AbstractMaterial
    E::Float64
    ν::Float64
    G::Float64
    k::Matrix{Float64}
    σ::Vector{Float64}
end

function LinearIsotropic(E, ν)
    G = E / (2 * (1 + ν))
    f = E / ((1.0 + ν) * (1.0 - 2.0 * ν))

    k = zeros(4, 4)

    k[1, 1] = k[2, 2] = k[3, 3] = f * (1.0 - ν)
    k[4, 4]  = G

    # Off diagonals
    k[1, 2] = k[1, 3] = k[2, 1] = k[2, 3] = k[3, 1] = k[3, 2] = f * ν

    LinearIsotropic(E, ν, G, k, zeros(4))
end

stiffness(mat::LinearIsotropic, ::GaussPoint2) = mat.k

create_matstat(::Type{LinearIsotropic}) = LinearIsotropicMS()

function stress(mat::LinearIsotropic, ɛ::Vector{Float64}, gp::GaussPoint2)
    D = stiffness(mat, gp)
    # TODO: Unroll this matrix multiplication
    A_mul_B!(mat.σ, D, ɛ)
    return mat.σ
end

end


text:
#haar.jl
#Algorithm 2.4 of Random Eigenvalues by Alan Edelman

#Experiment:    Generate random orthogonal/unitary matrices
#Plot:          Histogram eigenvalues
#Theory:        Eigenvalues are on unit circle

## Parameters
t = 5000        # trials
dx = 0.05       # binsize
n = 10          # matrix size

function haar_experiment(t,dx,n)
    v = Complex{Float64}[] # eigenvalue samples
    for i = 1:t
        # Sample random orthogonal matrix
        # X = QR(randn(n,n))[:Q]
        #  If you have non-uniformly sampled eigenvalues, you may need this fix:
        #X = X*diagm(sign(randn(n,1)))
        #
        # Sample random unitary matrix
        X = QR(randn(n,n)+im*randn(n,n))[:Q]
        #  If you have non-uniformly sampled eigenvalues, you may need this fix:
        #X = X*diagm(sign(randn(n)+im*randn(n)))
        append!(v, eigvals(full(X)))
    end
    x = ((-1+dx/2):dx:(1+dx/2))*pi
    x, v
end
x, v = haar_experiment(t,dx,n)
grid, count = hist(angle(v), x)

## Theory
h2 = t*n*dx/2*x.^0

## Plot
using Winston
p = FramedPlot()
h = Histogram(count, step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, h2, "color", "blue", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "haar.png")
end


text:
function testnetwork(rn, rn2)
    rxids = 1:numreactions(rn)
    rxids2 = 1:numreactions(rn2)
    @test params(rn) == params(rn2)
    @test paramsmap(rn) == paramsmap(rn2)
    @test species(rn) == species(rn2)
    @test speciesmap(rn) == speciesmap(rn2)
    @test dependents.(rn,rxids) == dependents.(rn2,rxids2)
    @test substrates.(rn,rxids) == substrates.(rn2,rxids2)
    @test products.(rn,rxids) == products.(rn2,rxids2)
    @test rateexpr.(rn,rxids) == rateexpr.(rn2,rxids2)
    @test oderatelawexpr.(rn,rxids) == oderatelawexpr.(rn2,rxids2)
    @test ssaratelawexpr.(rn,rxids) == ssaratelawexpr.(rn2,rxids2)
    @test ismassaction.(rn,rxids) == ismassaction.(rn2,rxids2)
end

# baseline network to test against
repressilator = @min_reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₃
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    10., m₃ --> m₃ + P₃
    2, P₁ --> ∅
    β*P₂, 2m₁ + 3P₁ --> P₂
end α K n δ γ β μ;
numrxs = numreactions(repressilator)

# empty network
rmin = @empty_reaction_network
foreach(param -> addparam!(rmin, param), params(repressilator))
foreach(spec -> addspecies!(rmin, spec), species(repressilator))
addreaction!(rmin, :(hillr(P₃,α,K,n)), :(∅ --> m₁) )
addreaction!(rmin, :(hillr(P₁,α,K,n)), :(∅ --> m₃) )
addreaction!(rmin, :((δ,γ)), :(m₁ ↔ ∅) )
addreaction!(rmin, :((δ,γ)), :(m₃ ↔ ∅) )
addreaction!(rmin, :β, :(m₁ --> m₁ + P₁) )
addreaction!(rmin, 10., :(m₃ --> m₃ + P₃) )
addreaction!(rmin, 2, :(P₁ --> ∅) )
addreaction!(rmin, :(β*P₂), :(2m₁ + 3P₁ --> P₂) )
testnetwork(rmin, repressilator)

# empty network stoichiometry interface
rmin = @empty_reaction_network
foreach(param -> addparam!(rmin, param), params(repressilator))
foreach(spec -> addspecies!(rmin, spec), species(repressilator))
addreaction!(rmin, :(hillr(P₃,α,K,n)), (), (:m₁=>1,))
addreaction!(rmin, :(hillr(P₁,α,K,n)), (), (:m₃=>1,) )
addreaction!(rmin, :δ, (:m₁ => 1,), ())
addreaction!(rmin, :γ, (), (:m₁ => 1,))
addreaction!(rmin, :δ, (:m₃ => 1,), ())
addreaction!(rmin, :γ, (), (:m₃ => 1,))
addreaction!(rmin, :β, (:m₁=>1,), (:m₁=>1, :P₁=>1) )
addreaction!(rmin, 10., (:m₃=>1,), (:m₃=>1, :P₃=>1) )
addreaction!(rmin, 2, (:P₁=>1,),() )
addreaction!(rmin, :(β*P₂), (:m₁=>2,:P₁=>3), (:P₂=>1,) )
testnetwork(rmin, repressilator)

# partial min_reaction_network
rmin2 = @min_reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₃
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    10., m₃ --> m₃ + P₃
end α K n δ γ β;
foreach(param -> addparam!(rmin2, param), params(repressilator))
foreach(spec -> addspecies!(rmin2, spec), species(repressilator))
addreaction!(rmin2, 2, :(P₁ --> ∅) )
addreaction!(rmin2, :(β*P₂), :(2m₁ + 3P₁ --> P₂) )
testnetwork(rmin2, repressilator)

# partial min_reaction_network, stoichiometry interface
rmin2 = @min_reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₃
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    10., m₃ --> m₃ + P₃
end α K n δ γ β;
foreach(param -> addparam!(rmin2, param), params(repressilator))
foreach(spec -> addspecies!(rmin2, spec), species(repressilator))
addreaction!(rmin2, 2, (:P₁=>1,), ())
addreaction!(rmin2, :(β*P₂), (:m₁=>2,:P₁=>3), (:P₂=>1,) )
testnetwork(rmin2, repressilator)

# simple example to test numerical agreement of solutions
fullrn = @reaction_network begin
    (k*X*(X-1)*(X-2),Y * X), X + 2X ⟺ Y+ X
end k
p = (1/6,)
u0 = [1000.,1000.]
tspan = (0.,5e-4)
foprob = ODEProblem(fullrn, u0, tspan, p)
fsol = solve(foprob,Tsit5(), abstol=1e-10, reltol=1e-10)

emptyrn = @empty_reaction_network
addspecies!(emptyrn, :X)
addspecies!(emptyrn, "Y")
addparam!(emptyrn, :k)
addreaction!(emptyrn, :((k*X*(X-1)*(X-2),Y * X)), :(X + 2X ⟺ Y+ X))
addodes!(emptyrn)
eoprob = ODEProblem(emptyrn, u0, tspan, p)
esol = solve(eoprob,Tsit5(), abstol=1e-10, reltol=1e-10)
@test norm(esol .- fsol, Inf) < 1e-7

minrn = @min_reaction_network begin
    X*Y, X + Y ⟾ X + 2X
end
addparam!(minrn, "k")
addreaction!(minrn, :(k*X*(X-1)*(X-2)), :(X + 2X ⟾ Y + X))
addodes!(minrn)
moprob = ODEProblem(minrn, u0, tspan, p)
msol = solve(moprob,Tsit5(), abstol=1e-10, reltol=1e-10)
@test norm(msol .- fsol, Inf) < 1e-7

# another simple example to test numerical agreement of solutions
# using stoichiometry interface
fullrn = @reaction_network begin
    (X^2/6,1.0), X ↔ Y + X
end
u0 = [1000.,1000.]
tspan = (0.,5e-4)
foprob = ODEProblem(fullrn, u0, tspan)
fsol = solve(foprob,Tsit5(), abstol=1e-10, reltol=1e-10)

emptyrn = @empty_reaction_network
addspecies!(emptyrn, :X)
addspecies!(emptyrn, "Y")
addreaction!(emptyrn, :(X^2/6), (:X=>1,), (:X=>1, :Y=>1))
addreaction!(emptyrn, 1.0, (:X=>1,:Y=>1), (:X=>1,))
addodes!(emptyrn)
eoprob = ODEProblem(emptyrn, u0, tspan)
esol = solve(eoprob,Tsit5(), abstol=1e-10, reltol=1e-10)
@test norm(esol .- fsol, Inf) < 1e-7


text:
"""

```
wrap_cell(nc::AbstractVector, i::Int, j::Int, k::Int)
```

Given the indexes of a cell, return the periodic cell which correspondst to
it, if the cell is outside the main box.

"""
function wrap_cell(nc::AbstractVector, i::Int, j::Int, k::Int)

  if i < 1
    i = nc[1] + i
  elseif i > nc[1]
    i = i - nc[1]
  end

  if j < 1
    j = nc[2] + j
  elseif j > nc[2]
    j = j - nc[2]
  end

  if k < 1
    k = nc[3] + k
  elseif k > nc[3]
    k = k - nc[3]
  end

  return i, j, k

end


text:
using BinaryBuilder

# Collection of sources required to build ZMQ
sources = [
    "https://github.com/zeromq/libzmq.git" =>
    "2cb1240db64ce1ea299e00474c646a2453a8435b", # v4.3.1
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/libzmq

if [[ ${target} == x86_64-linux-musl ]] && [[ $(gcc --version | head -1 | awk '{ print $3 }') == "4.8.5" ]]; then

cat > mm_malloc.patch << 'END'
--- mm_malloc.h
+++ mm_malloc.h
@@ -31,7 +31,7 @@
 #ifndef __cplusplus
 extern int posix_memalign (void **, size_t, size_t);
 #else
-extern "C" int posix_memalign (void **, size_t, size_t) throw ();
+extern "C" int posix_memalign (void **, size_t, size_t);
 #endif

 static __inline void *
END

cd /opt/x86_64-linux-musl/lib/gcc/x86_64-linux-musl/4.8.5/include/

patch < "${WORKSPACE}/srcdir/libzmq/mm_malloc.patch"

cd $WORKSPACE/srcdir/libzmq

fi

sh autogen.sh
./configure --prefix=$prefix --host=${target} \
    --without-docs --disable-libunwind --disable-perf --disable-Werror \
    --disable-eventfd --without-gcov --disable-curve-keygen \
    CXXFLAGS="-g -O2 -fms-extensions"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms() # build on all supported platforms
platforms = expand_gcc_versions(platforms)

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

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "ZMQ", v"4.3.1", sources, script, platforms, products, dependencies)

max_stars_repo_name: kdheepak/ZMQBuilder

max_issues_repo_name: kdheepak/ZMQBuilder

max_forks_repo_name: kdheepak/ZMQBuilder


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

"Approximation parameters"
N = 2 # The order of approximation
K1D = 12
CFL = 2 # CFL goes up to 3 OK...
T = 1.0 # endtime

"Mesh related variables"
Kx = convert(Int,4/3*K1D)
Ky = K1D
VX, VY, EToV = uniform_quad_mesh(Kx, Ky)
@. VX = 15*(1+VX)/2
@. VY = 5*VY

# initialize ref element and mesh
# specify lobatto nodes to automatically get DG-SEM mass lumping
rd = init_reference_quad(N,gauss_quad(0,0,N))
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

# construct hybridized SBP operators
@unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd
Qr = Pq'*M*Dr*Pq
Qs = Pq'*M*Ds*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Bs = diagm(wf.*nsJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
         -Br*Ef  Br]
Qsh = .5*[Qs-Qs' Ef'*Bs;
         -Bs*Ef  Bs]

# make sparse skew symmetric versions of the operators"
# precompute union of sparse ids for Qr, Qs
Qrhskew = .5*(Qrh-transpose(Qrh))
Qshskew = .5*(Qsh-transpose(Qsh))
Qrh_sparse = droptol!(sparse(Qrhskew),1e-12)
Qsh_sparse = droptol!(sparse(Qshskew),1e-12)
Qrsids = [unique([Qrh_sparse[i,:].nzind; Qsh_sparse[i,:].nzind]) for i = 1:size(Qrhskew,1)]

# Make node maps periodic
@unpack Nfaces = rd
@unpack x,y,xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY))
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

# convert operators to a quadrature-node basis
@unpack wq = rd
Vh = droptol!(sparse([eye(length(wq)); Ef]),1e-12)
Ph = droptol!(sparse(diagm(@. 1/wq)*transpose(Vh)),1e-12)
Lf = droptol!(sparse(diagm(@. 1/wq)*(transpose(Ef)*diagm(wf))),1e-12)

# define initial conditions at quadrature nodes
@unpack xq,yq = md
rho,u,v,p = vortex(xq,yq,0)
Q = primitive_to_conservative(rho,u,v,p)

# interpolate geofacs to both vol/surf nodes
@unpack rxJ, sxJ, ryJ, syJ = md
rxJ, sxJ, ryJ, syJ = (x->Vh*x).((rxJ, sxJ, ryJ, syJ)) # interp to hybridized points
@pack! md = rxJ, sxJ, ryJ, syJ

# pack SBP operators into tuple
ops = (Qrhskew,Qshskew,Qrh_sparse,Qsh_sparse,Qrsids,Ph,Lf)

"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)/2  # estimated trace constant for CFL
h  = 2/K1D
dt = CFL * h / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"sparse version - precompute sparse row ids for speed"
function sparse_hadamard_sum(Qhe,ops,vgeo,flux_fun)

    (Qr,Qs,Qnzids) = ops
    (rxJ,sxJ,ryJ,syJ) = vgeo
    nrows = size(Qr,1)
    nfields = length(Qhe)

    # precompute logs for logmean
    (rho,u,v,beta) = Qhe
    Qlog = (log.(rho), log.(beta))

    rhsQe = ntuple(x->zeros(nrows),nfields)
    rhsi = zeros(nfields) # prealloc a small array
    for i = 1:nrows
        Qi = (x->x[i]).(Qhe)
        Qlogi = (x->x[i]).(Qlog)

        fill!(rhsi,0) # reset rhsi before accumulation
        for j = Qnzids[i] # nonzero row entries
            Qj = (x->x[j]).(Qhe)
            Qlogj = (x->x[j]).(Qlog)

            # assumes affine meshes here!
            Fx,Fy = flux_fun(Qi,Qj,Qlogi,Qlogj)
            Fr = @. rxJ*Fx + ryJ*Fy
            Fs = @. sxJ*Fx + syJ*Fy

            @. rhsi += Qr[i,j]*Fr + Qs[i,j]*Fs
        end

        # faster than one-line fixes (no return args)
        for fld in eachindex(rhsQe)
            rhsQe[fld][i] = rhsi[fld]
        end
    end
    return rhsQe
end

"Qh = (rho,u,v,beta), while Uh = conservative vars"
function rhs(Q,md::MeshData,ops,flux_fun,compute_rhstest=false)

    @unpack rxJ,sxJ,ryJ,syJ,J,wJq = md
    @unpack nxJ,nyJ,sJ,mapP,mapB,K = md
    Qrh,Qsh,Qrh_sparse,Qsh_sparse,Qrsids,Ph,Lf = ops
    Nq,Nh = size(Ph)

    # entropy var projection
    VU = v_ufun(Q...)
    Uf = u_vfun((x->Ef*x).(VU)...)
    (rho,rhou,rhov,E) = vcat.(Q,Uf) # assume volume nodes are collocated

    # convert to rho,u,v,beta vars
    beta = betafun(rho,rhou,rhov,E)
    Qh = (rho, rhou./rho, rhov./rho, beta) # redefine Q = (rho,u,v,β)

    # compute face values
    QM = (x->x[Nq+1:end,:]).(Qh)
    QP = (x->x[mapP]).(QM)

    # simple lax friedrichs dissipation
    (rhoM,rhouM,rhovM,EM) = Uf
    rhoUM_n = @. (rhouM*nxJ + rhovM*nyJ)/sJ
    lam = abs.(wavespeed(rhoM,rhoUM_n,EM))
    LFc = .5*max.(lam,lam[mapP]).*sJ

    fSx,fSy = flux_fun(QM,QP)
    normal_flux(fx,fy,u) = fx.*nxJ + fy.*nyJ - LFc.*(u[mapP]-u)
    flux = normal_flux.(fSx,fSy,Uf)
    rhsQ = (x->Lf*x).(flux)

    # compute volume contributions using flux differencing
    for e = 1:K
        Qhe = tuple((x->x[:,e]).(Qh)...) # force tuples for fast splatting
        vgeo_local = (x->x[1,e]).((rxJ,sxJ,ryJ,syJ)) # assumes affine elements for now

        Qops = (Qrh_sparse,Qsh_sparse,Qrsids)
        QFe = sparse_hadamard_sum(Qhe,Qops,vgeo_local,flux_fun)

        mxm_accum!(X,x) = X[:,e] += 2*Ph*x
        mxm_accum!.(rhsQ,QFe)
    end

    rhsQ = (x -> -x./J).(rhsQ)

    rhstest = 0
    if compute_rhstest
        for fld in eachindex(rhsQ)
            rhstest += sum(wJq.*VU[fld].*rhsQ[fld])
        end
    end

    return rhsQ,rhstest # scale by Jacobian
end

Q = collect(Q) # make Q,resQ arrays of arrays for mutability
resQ = [zeros(size(x)) for i in eachindex(Q)]
for i = 1:Nsteps

    rhstest = 0
    for INTRK = 1:5
        compute_rhstest = INTRK==5
        rhsQ,rhstest = rhs(Q,md,ops,euler_fluxes,compute_rhstest)

        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. Q += rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Time step: $i out of $Nsteps with rhstest = $rhstest")
    end
end

"project solution back to GLL nodes"
Q = (x->Pq*x).(Q)

# use a higher degree quadrature for error evaluation
@unpack VDM = rd
@unpack J = md
rq2,sq2,wq2 = quad_nodes_2D(N+2)
Vq2 = vandermonde_2D(N,rq2,sq2)/VDM
wJq2 = diagm(wq2)*(Vq2*J)
xq2,yq2 = (x->Vq2*x).((x,y))

Qq = (x->Vq2*x).(Q)
Qex = primitive_to_conservative(vortex(xq2,yq2,T)...)
L2err = 0.0
for fld in eachindex(Q)
    global L2err
    L2err += sum(@. wJq2*(Qq[fld]-Qex[fld])^2)
end
L2err = sqrt(L2err)
println("L2err at final time T = $T is $L2err\n")

#plotting nodes
@unpack Vp = rd
gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2)

vv = Vp*Q[1]
scatter(Vp*x,Vp*y,vv,zcolor=vv,camera=(0,90))

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
include("tabular_deterministic_policy.jl")
include("tabular_random_policy.jl")
include("exploring_start_policy.jl")
include("reinforce_policy.jl")


text:
#################################################################
### This test was lifted from the vingette for the R gmm package.
#################################################################
using Distributions
using Linearalgebra
srand(42)
x = rand(Normal(4, 2), 1000)

function h_norm(θ)
    m1 = [θ[1]] .- x
    m2 = [θ[2]].^2 .- (x .- [θ[1]]).^2
    m3 = x.^3 .- [θ[1]].*([θ[1]].^2 .+ 3*[θ[2]].^2)
    return [m1 m2 m3]
end

step_1 = GMMEstimator(h_norm, [1.0, 1.0], initialW = Matrix(I, 3, 3), mgr = OneStepGMM())
@time estimate!(step_1);

ch(θ) = [1 -1]*θ
hlb = [0.]
hub = [0.]

cstep_1 = GMMEstimator(h_norm, [1.0, 1.0], initialW = Matrix(I, 3, 3), mgr = OneStepGMM(), constraints = Constrained(ch, hlb, hub, 1))

step_2_hac = GMMEstimator(h_norm, coef(step_1), initialW = MomentBasedEstimators.optimal_W(step_1, QuadraticSpectralKernel(0.91469)));

estimate!(step_2_hac);

step_2_iid = GMMEstimator(h_norm, coef(step_1), initialW = optimal_W(step_1, HC0()));
estimate!(step_2_iid);

gmm_qs_mgr = GMMEstimator(h_norm, [1.,1.], initialW = Matrix(I, 3, 3), mgr = TwoStepGMM(QuadraticSpectralKernel(0.91469)));
estimate!(gmm_qs_mgr);

gmm_iid_mgr = GMMEstimator(h_norm, [1.,1.], initialW = Matrix(I, 3, 3), mgr = TwoStepGMM(HC0()));
estimate!(gmm_iid_mgr);

### MDE
kl_iid = MDEstimator(h_norm, [1.0, 1.0])
estimate!(kl_iid)

kl_smt = MDEstimator(h_norm, [1.0, 1.0], kernel = MomentBasedEstimators.TruncatedSmoother(2))
estimate!(kl_smt)

el_iid = MDEstimator(h_norm, [1.0, 1.0], div = Divergences.ReverseKullbackLeibler())
cu_iid = MDEstimator(h_norm, [1.0, 1.0], div = Divergences.ChiSquared())


text:
    ####################################################################
#######                Main Muffin function                  #######
####################################################################

function muffin(;folder="",dataobj="2gauss",datapsf="",nitermax = 500, rhop = 1,
                rhot = 5, rhov = 2, rhos = 1, μt = 5e-1, μv = 1e-0, mueps = 1e-3,
                bw = 5, ws="",parallel="")

println("")
println("MUFFIN initialisation")

                 ##################################
    ################### data initialisation #################
                 ##################################

    if typeof(dataobj) == ASCIIString
        if dataobj == "m31"
            psf = "data/meerkat_m30_25pix.psf.fits"
            obj = "data/M31.fits"
            tmp = string(Pkg.dir("Muffin"))
            psf = string(tmp,tmp[1],psf)
            obj = string(tmp,tmp[1],obj)
        elseif dataobj == "andro"
            psf = "data/meerkat_m30_25pix.psf.fits"
            obj = "data/andro.fits"
            tmp = string(Pkg.dir("Muffin"))
            psf = string(tmp,tmp[1],psf)
            obj = string(tmp,tmp[1],obj)
        elseif dataobj == "2gauss"
            psf = "data/meerkat_m30_25pix.psf.fits"
            obj = "data/2gauss.fits"
            tmp = string(Pkg.dir("Muffin"))
            psf = string(tmp,tmp[1],psf)
            obj = string(tmp,tmp[1],obj)
        elseif dataobj == "chiara"
            psf = "/home/deguignet/Julia/example_sim_psf.fits"
            obj = "/home/deguignet/Julia/example_sim_dirty.fits"
        elseif isempty(folder)
            tmp = pwd()
            psf = string(tmp,tmp[1],datapsf)
            obj = string(tmp,tmp[1],dataobj)
        elseif typeof(folder) == ASCIIString
                 psf = string(folder,folder[1],datapsf)
                 obj = string(folder,folder[1],dataobj)
        else
                 error("data folder is not correct")
        end

    elseif isempty(dataobj)
        psf = "data/meerkat_m30_25pix.psf.fits"
        obj = "data/M31.fits"
    end

println("psf something :"," ",psf)
println("obj :"," ",obj)

                ##################################
    ################# Structure initialisation #################
                ##################################

if dataobj == "chiara"
    ##################################
    println("loading psf chiara...")
    psfst = loadpsf_dirty(psf)
    println("loading sky...")
    skyst = loadsky_dirty(obj,psfst.mypsf,psfst.nu)
    ##################################
else
    ##################################
    println("loading psf other...")
    psfst = loadpsf(psf,bw)
    println("loading sky...")
    skyst = loadsky(obj,psfst.mypsf,psfst.nu)
    ##################################
end

    ##################################
    spatialwlt  = [WT.db1,WT.db2,WT.db3,WT.db4,WT.db5,WT.db6,WT.db7,WT.db8]
    # spatialwlt  = [WT.db1,WT.db2,WT.db3,WT.db4,WT.db5,WT.db6,WT.db7,WT.db8,WT.haar]

    const nspat = length(spatialwlt)
    const nfreq = size(psfst.mypsf)[3]
    const nspec = 1
    const nxy = size(skyst.mydata)[1]
    niter = 0
    lastiter = 0
    #################################

    #################################
    println("loading param...")
    algost = loadparam(nspat,nfreq,nspec,nxy,niter,lastiter,nitermax)

    println("loading arrays...")
    if ws == "true"
        file = string("/home/deguignet/resultats_400ite_woshar.jld")
        # file = string("/Users/deguignet/test.jld")
        admmst = loadarray_ws(rhop,rhot,rhov,rhos,μt,μv,mueps,nspat,nfreq,nxy,
                            skyst.mydata,psfst.mypsfadj,file)
    algost.niter = load(file,"algost.lastiter")
    toolst = loadtools(nitermax,nfreq,nxy)
    toolst.tol1 = load(file,"toolst.tol1")
    toolst.tol2 = load(file,"toolst.tol2")
    toolst.tol3 = load(file,"toolst.tol3")
    toolst.tol4 = load(file,"toolst.tol4")
    toolst.tol5 = load(file,"toolst.tol5")
    skyst.mydata = load(file,"skyst.mydata")

    else
        admmst = loadarray(rhop,rhot,rhov,rhos,μt,μv,mueps,nspat,nfreq,nxy,
                            skyst.mydata,psfst.mypsfadj,parallel)
        toolst = loadtools(nitermax,nfreq,nxy)
    end

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

                 ##################################
    #####################  Main Admm Loop  #####################
                 ##################################
    println("Starting ADMM")
    #################################
    psfst, skyst, algost, admmst, toolst = muffinadmm(psfst, skyst, algost, admmst, toolst)
    #################################

    return psfst, skyst, algost, admmst, toolst
end

####################################################################
#######                   Main Admm Loop                     #######
####################################################################

function muffinadmm(psfst, skyst, algost, admmst, toolst)

    const rhop = admmst.rhop
    const rhot = admmst.rhot
    const rhov = admmst.rhov
    const rhos = admmst.rhos
    const μt = admmst.μt
    const μv = admmst.μv
    const mueps = admmst.mueps
    const tt = admmst.tt
    const mu = admmst.mu
    const nspat = algost.nspat
    const nfreq = algost.nfreq
    const nspec = algost.nspec
    const nxy = algost.nxy
    const fty = admmst.fty
    const nitermax = algost.nitermax

    spatialwlt  = [WT.db1,WT.db2,WT.db3,WT.db4,WT.db5,WT.db6,WT.db7,WT.db8,WT.haar]

    niter = algost.niter

    loop = true

    println("")
    println("                                =======================                                ")
    println("                               | MUFFIN RECONSTRUCTION |                               ")

    println("========================================================================================")
    println("ADMM"," ","|"," ","||x - xm||"," ","|"," ","||x - xp||"," ","|"," ",
            "||Hx - t||"," ","|"," ","||x - s||"," "," |"," ","||sh - v||"," ","|"," ","     time      |")
    println("========================================================================================")
    tic()
        while loop
            tic()
            niter +=1
            println("    "," ","|"," ","          "," ","|"," ","          "," ","|"," ",
                    "          "," ","|"," ","         "," "," |"," ","          "," ","|"," ")
            # println("ADMM iteration: $niter")

    ######################################
    ######################################
            # tic()
            # @parallel for i in 2:nfreq+1
            #     z = i-1
            #     admmst.wlt[:,:,z],admmst.x[:,:,z],admmst.t[:,:,z,:],admmst.taut[:,:,z,:],admmst.p[:,:,z],admmst.taup[:,:,z] =
            #                                 @fetchfrom(i,parallelmuffin(admmst.wlt[:,:,z], admmst.taut[:,:,z,:], admmst.t[:,:,z,:], rhot, admmst.x[:,:,z],
            #                                 psfst.mypsf[:,:,z], admmst.p[:,:,z], admmst.taup[:,:,z],
            #                                 fty[:,:,z], rhop, admmst.taus[:,:,z], admmst.s[:,:,z], rhos, admmst.mu, spatialwlt,
            #                                 μt, nspat))
########################

            @sync @parallel for z in 1:nfreq
                admmst.wlt[:,:,z],admmst.x[:,:,z],admmst.t[:,:,z,:],admmst.taut[:,:,z,:],admmst.p[:,:,z],admmst.taup[:,:,z] =
                                            parallelmuffin(admmst.wlt[:,:,z], admmst.taut[:,:,z,:], admmst.t[:,:,z,:], rhot, admmst.x[:,:,z],
                                            psfst.mypsf[:,:,z], admmst.p[:,:,z], admmst.taup[:,:,z],
                                            fty[:,:,z], rhop, admmst.taus[:,:,z], admmst.s[:,:,z], rhos, admmst.mu, spatialwlt,
                                            μt, nspat)

            end

            ##############################
            ######### prox spec ##########

            tmp = admmst.tauv + rhov*admmst.v

            admmst.s, admmst.sh = estime_ssh(admmst.s,admmst.sh,tmp,nxy,nspec,admmst.spectralwlt,
                                              admmst.x,admmst.taus,rhov,rhos)

            tmp = admmst.sh - admmst.tauv/rhov
            admmst.v = prox_u(tmp,μv/rhov)

            ########################################
            #### update of Lagrange multipliers ####

            admmst.tauv = admmst.tauv + rhov*(admmst.v-admmst.sh)
            admmst.taus = admmst.taus + rhos*(admmst.s-admmst.x)

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

            #
            # ##############################
            # ########## update x ##########
            # for z in 1:nfreq
            #     admmst.wlt[:,:,z] = myidwt((admmst.wlt)[:,:,z], nspat, (admmst.taut)[:,:,z,:], rhot,
            #                         (admmst.t)[:,:,z,:], spatialwlt)
            # end
            #
            # b = fty + admmst.taup + rhop*admmst.p + admmst.taus + rhos*admmst.s
            #
            # admmst.x = estime_x_par(admmst.x,psfst.mypsf,admmst.wlt + b,mu,nfreq)
            # ##############################
            # ######### prox spat ##########
            # tmp1 = 0.0
            # tmp2 = zeros(Float64,nxy,nxy)
            # for z in 1:nfreq
            #     for b in 1:nspat
            #             hx = dwt(admmst.x[:,:,z],wavelet(spatialwlt[b]))
            #             tmp = hx - admmst.taut[:,:,z,b]/rhot
            #             admmst.t[:,:,z,b] = prox_u(tmp,μt/rhot)
            #             admmst.taut[:,:,z,b] = admmst.taut[:,:,z,b] + rhot*(admmst.t[:,:,z,b]-hx)
            #             tmp1 = vecnorm([tmp2 (hx-(admmst.t)[:,:,z,b])],2)
            #             tmp2 = (hx-(admmst.t)[:,:,z,b])
            #     end
            # end
            # tmp2[:] = 0
            # ##############################
            # ###### prox positivity #######
            #
            # tmp = admmst.x-admmst.taup/rhop
            #
            # admmst.p = max(0,tmp)
            #
            # ##############################
            # ######### prox spec ##########
            #
            # tmp = admmst.tauv + rhov*admmst.v
            #
            # admmst.s, admmst.sh = estime_ssh(admmst.s,admmst.sh,tmp,nxy,nspec,admmst.spectralwlt,
            #                                   admmst.x,admmst.taus,rhov,rhos)
            #
            # tmp = admmst.sh - admmst.tauv/rhov
            # admmst.v = prox_u(tmp,μv/rhov)
            #
            #
            # ########################################
            # #### update of Lagrange multipliers ####
            #
            # admmst.taup = admmst.taup + rhop*(admmst.p-admmst.x)
            # admmst.tauv = admmst.tauv + rhov*(admmst.v-admmst.sh)
            # admmst.taus = admmst.taus + rhos*(admmst.s-admmst.x)

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

            ##############################
            ##### computer residues ######

            push!(toolst.tol1,vecnorm(admmst.x - admmst.xmm, 2)^2)
            push!(toolst.tol2,vecnorm(admmst.x - admmst.p, 2)^2)
            # push!(toolst.tol3,tmp1^2)
            push!(toolst.tol4,vecnorm(admmst.x - admmst.s, 2)^2)
            push!(toolst.tol5,vecnorm(admmst.sh - admmst.v, 2)^2)

            # ##############################
            # ####### stopping rule ########

            if (niter >= nitermax) #|| ((toolst.tol1[niter] < 1E-6) && (toolst.tol2[niter] < 1E-4))
                loop = false
                algost.lastiter = niter
            end

            admmst.xmm[:] = admmst.x

            #
            # println("ADMM"," ","|"," ","||x - xm||"," ","|"," ","||x - xp||"," ","|"," ",
            #         "||Hx - t||"," ","|"," ","||x - s||"," ","|"," ","||sh - v||"," ","|"," ","time")
            # println(niter," ","|"," ",toolst.tol1[niter]," ","|"," ",toolst.tol2[niter]," ","|"," ",
            #         toolst.tol3[niter]," ","|"," ",toolst.tol4[niter]," ","|"," ",toolst.tol5[niter]," ","|"," ",toq())

            @printf("%03d  | ", niter)
            @printf("%02.04e | ", toolst.tol1[niter])
            @printf("%02.04e | ", toolst.tol2[niter])
            # @printf("%02.04e | ", toolst.tol3[niter])
            @printf("%02.04e | ", 0)
            @printf("%02.04e | ", toolst.tol4[niter])
            @printf("%02.04e | ", toolst.tol5[niter])
            @printf("%f seconds  \n", toq())

        end
    println("")
    @printf("time for ADMM : %f seconds \n", toq())
    ####################################################################
    ####################################################################
    ####################################################################

    return psfst, skyst, algost, admmst, toolst

end

####################################################################
#######                  Data functions                      #######
####################################################################

##################################
function cubefilter{T<:FloatingPoint}(imagecube::Array{T,3},psfcube::Array{T,3})
    nximg, nyimg, nfreq = size(imagecube)
    nxpsf, nypsf, nfreq = size(psfcube)
    rescube = Array(Float64,nximg,nyimg,nfreq)
    for k in 1:nfreq
        rescube[:,:,k] = imfilter_fft(imagecube[:,:,k],psfcube[:,:,k],"circular")
    end
    return rescube
end

function cubefilter{T<:FloatingPoint}(imagecube::SharedArray{T,3},psfcube::Array{T,3})
    nximg, nyimg, nfreq = size(imagecube)
    nxpsf, nypsf, nfreq = size(psfcube)
    rescube = Array(Float64,nximg,nyimg,nfreq)
    for k in 1:nfreq
        rescube[:,:,k] = imfilter_fft(imagecube[:,:,k],psfcube[:,:,k],"circular")
    end
    return rescube
end

function cubefilter{T<:FloatingPoint}(imagecube::Array{T,2},psfcube::Array{T,3})
    nximg, nyimg = size(imagecube)
    nxpsf, nypsf, nfreq = size(psfcube)
    rescube = Array(Float64,nximg,nyimg,nfreq)
    for k in 1:nfreq
        rescube[:,:,k] = imfilter_fft(imagecube,psfcube[:,:,k],"circular")
    end
    return rescube
end
##################################

##################################
function cubeaverage{T<:FloatingPoint}(imagecube::Array{T,3},M::Int)
    nxpsf, nypsf, nfreq = size(imagecube)
    if nfreq < M
        error("Not enough channels to average!")
    end
    nfreqavg = itrunc(nfreq/M)
    rescube = Array(Float64, nxpsf, nypsf, nfreqavg)
    for k in 1:nfreqavg
        rescube[:,:,k] = sum(imagecube[:,:,(k-1)*M+1:k*M], 3)/M
    end

    return rescube
end
##################################

##################################
function cubefreq(psf::ASCIIString,imagecube::Array,M::Int)
    nxpsf, nypsf, nfreq = size(imagecube)
    if nfreq < M
        error("Not enough channels to average!")
    end
    nfreqavg = itrunc(nfreq/M)
    file = FITS(psf)
    header = read_header(file[1])
    nustart = header["CRVAL4"]
    nustep = header["CDELT4"]
    nu0 = header["RESTFRQ"]
    nu = zeros(Float64,nfreqavg)
    for i in 1:nfreqavg
        nu[i] = nustart + (M-1)*nustep + (i-1)*M*nustep
    end
    close(file)
    return nu,nu0
end
##################################

##################################
function cropcubexy{T<:FloatingPoint}(imagecube::Array{T,3},M::Int)
    nxpsf, nypsf, nfreq = size(imagecube)
    if ~(nxpsf == nypsf)
        error("Image must be square!")
    end
    if nxpsf < M
        println("no crop")
        return imagecube
    end

    rescube = Array(Float64, M, M, nfreq)
    if (iseven(nxpsf) && iseven(M)) | (isodd(nxpsf) && isodd(M))
        nstart = (nxpsf-M)/2 +1
        nend = (nxpsf+M)/2
        rescube = imagecube[nstart:nend,nstart:nend,:]
    else
        error("size of image and crop must have same parity!")
    end
    return rescube
end
##################################

##################################
function sky2cube{T<:FloatingPoint}(sky::Array{T,2},nu::Array{T,1})
    nx, ny = size(sky)
    nbands = length(nu)
    corl = minimum([nx,ny])/5.0
    rho(x,y) = exp(-norm([x, y])/corl)
    tx = linspace(0,nx-1,nx)
    ty = linspace(0,ny-1,ny)
    field = genghf(tx, ty, rho)
    skycube = Array(Float64, nx, ny, nbands)

    sm, sM = extrema(sky)
    fm, fM = extrema(field)

    alpha = 0.8 + 0.4*(sky-sm)/(sM-sm) + 0.4*(field-fm)/(fM-fm)

    nu0 = (nu[end]+nu[1])/2
    for k in 1:nbands
        skycube[:,:,k] = sky.* (nu[k]/nu0) .^(-alpha)
    end
    return skycube,alpha
end
##################################

##################################
function lecture(directory::ASCIIString)
    file = FITS(directory)
    data = float64(read(file[1]))
    close(file)
    data = squeeze(data,find(([size(data)...].==1)))
    data = data[:,:,:]

    return data
end
##################################

####################################################################
#######                  Admm functions                      #######
####################################################################

#################################
######### x estimation ##########
# function estime_x_par(x::Array{Float64,3},mypsf::Array{Float64,3},
#                         wlt_b::Array{Float64,3},mu::Float64,nfreq::Int64)
#
#     nxy = (size(x))[1]
#     nxypsf = (size(mypsf))[1]
#     fftpsf = zeros(Complex64,nxy,nxy,nfreq)
#     psfcbe = zeros(Complex64,nxy,nxy,nfreq)
#     psfpad = zeros(Float64,nxy,nxy,nfreq)
#
#     for z in 1:nfreq
#         psfpad[1:nxypsf,1:nxypsf,z] = fftshift(mypsf[:,:,z])
#         psfcbe[:,:,z] = 1./(abs(fft(psfpad[:,:,z])).^2+mu)
#         x[:,:,z] = real(ifft(psfcbe[:,:,z].*fft(wlt_b[:,:,z])))
#     end
#
#     return x
# end
#
# function estime_x_par(x::SharedArray{Float64,3},mypsf::Array{Float64,3},
#                         wlt_b::SharedArray{Float64,3},mu::Float64,nfreq::Int64)
#
#     nxy = (size(x))[1]
#     nxypsf = (size(mypsf))[1]
#     fftpsf = zeros(Complex64,nxy,nxy,nfreq)
#     psfcbe = zeros(Complex64,nxy,nxy,nfreq)
#     psfpad = zeros(Float64,nxy,nxy,nfreq)
#
#     for z in 1:nfreq
#         psfpad[1:nxypsf,1:nxypsf,z] = fftshift(mypsf[:,:,z])
#         psfcbe[:,:,z] = 1./(abs(fft(psfpad[:,:,z])).^2+mu)
#         x[:,:,z] = real(ifft(psfcbe[:,:,z].*fft(wlt_b[:,:,z])))
#     end
#
#     return x
# end

#################################
###### proximity operators ######
function prox_u(u::Array,μ::Float64)
    return max(0, 1-μ./abs(u)).*u
end

function prox_u(u::SharedArray,μ::Float64)
    return max(0, 1-μ./abs(u)).*u
end

# function prox_u(u::Array,μ::Float64)
# 	res=zeros(size(u));
# 	x=u;
# 	ind = (abs(x).>(.75*μ^(2/3))) ;
# 	res[ind]= 2/3*x[ind].*(1+cos(2*pi/3-2/3*acos(μ/8*(abs(x[ind])/3).^(-3/2)))) ;
#    return res
# end

#################################
####### s / sh estimation #######
function estime_ssh(s::Array{Float64,3},sh::Array{Float64,3},tmp::Array{Float64,3},
                    nxy::Int64,nspec::Int64,spectralwlt::Array{Float64,3},
                    x::Array{Float64,3},taus::Array{Float64,3},rhov::Float64,rhos::Float64)

    spectralwlt = idct(tmp,3)
    s = (spectralwlt + rhos*x - taus)/(rhov*nspec + rhos)
    sh = dct(s,3)

    return s,sh
end

function estime_ssh(s::Array{Float64,3},sh::Array{Float64,3},tmp::Array{Float64,3},
                    nxy::Int64,nspec::Int64,spectralwlt::Array{Float64,3},
                    x::SharedArray{Float64,3},taus::Array{Float64,3},rhov::Float64,rhos::Float64)

    spectralwlt = idct(tmp,3)
    s = (spectralwlt + rhos*x - taus)/(rhov*nspec + rhos)
    sh = dct(s,3)

    return s,sh
end

function myidwt(wlt::Array{Float64,2},nspat::Int,taut::Array{Float64,4},rhot::Float64,t::Array{Float64,4},spatialwlt)
        wlt = idwt(taut[:,:,1,1] + rhot*t[:,:,1,1],wavelet(spatialwlt[1]))
            for b in 2:nspat
                wlt = wlt + idwt(taut[:,:,1,b] + rhot*t[:,:,1,b],wavelet(spatialwlt[b]))
            end
        return wlt
end
function myidwt(wlt,nspat::Int,taut,rhot::Float64,t,spatialwlt)
        wlt = idwt(taut[:,:,1,1] + rhot*t[:,:,1,1],wavelet(spatialwlt[1]))
            for b in 2:nspat
                wlt = wlt + idwt(taut[:,:,1,b] + rhot*t[:,:,1,b],wavelet(spatialwlt[b]))
            end
        return wlt
end

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

##################################
function cubefreqchiara(nfrequencies::Int)
    nfreq = nfrequencies
    nustart = 9.85e8
    nustep = 2e6

    nu = zeros(Float64,nfreq)
    for i in 1:nfreq
        nu[i] = nustart + (i-1)*nustep
    end
    nu0 = (nu[1] + nu[nfreq])/2

    return nu,nu0
end

function parallelmuffin(wlt::Array{Float64,2},taut::Array{Float64,4},t::Array{Float64,4},rhot::Float64,
                        x::Array{Float64,2},psf::Array{Float64,2},p::Array{Float64,2},taup::Array{Float64,2},
                        fty::Array{Float64,2},rhop::Float64,taus::Array{Float64,2},s::Array{Float64,2},rhos::Float64,
                        mu::Float64,spatialwlt,μt::Float64,nspat::Int)

        wlt = myidwt(wlt, nspat, taut[:,:,1,:], rhot, t[:,:,1,:], spatialwlt)
        b = fty + taup + rhop*p + taus + rhos*s
        wlt_b = wlt + b

        nxy = (size(x))[1]
        nxypsf = (size(psf))[1]
        psfcbe = zeros(Complex64,nxy,nxy)
        psfpad = zeros(Float64,nxy,nxy)
        psfpad[1:nxypsf,1:nxypsf] = psf[:,:]
        psfcbe = 1./(abs(fft(psfpad)).^2+mu)
        x = real(ifft(psfcbe.*fft(wlt_b)))

        # tmp1 = 0.0
        # tmp2 = zeros(Float64,nxy,nxy)
        for b in 1:nspat
                    hx = dwt(x,wavelet(spatialwlt[b]))
                    tmp = hx - taut[:,:,1,b]/rhot
                    t[:,:,1,b] = prox_u(tmp,μt/rhot)
                    taut[:,:,1,b] = taut[:,:,1,b] + rhot*(t[:,:,1,b]-hx)
                    # tmp1 = vecnorm([tmp2 (hx-(t)[:,:,z,b])],2)
                    # tmp2 = (hx-(t)[:,:,z,b])
        end
        # tmp2[:] = 0

        tmp = x-taup/rhop
        p = max(0,tmp)
        taup = taup + rhop*(p-x)

    return wlt,x,t,taut,p,taup

end

function parallelmuffin(wlt::SharedArray{Float64,2},taut::SharedArray{Float64,4},t::SharedArray{Float64,4},rhot::Float64,
                        x::SharedArray{Float64,2},psf::Array{Float64,2},p::SharedArray{Float64,2},taup::SharedArray{Float64,2},
                        fty::Array{Float64,2},rhop::Float64,taus::Array{Float64,2},s::Array{Float64,2},rhos::Float64,
                        mu::Float64,spatialwlt,μt::Float64,nspat::Int)

        wlt = myidwt(wlt, nspat, taut[:,:,1,:], rhot, t[:,:,1,:], spatialwlt)
        b = fty + taup + rhop*p + taus + rhos*s
        wlt_b = wlt + b

        nxy = (size(x))[1]
        nxypsf = (size(psf))[1]
        psfcbe = zeros(Complex64,nxy,nxy)
        psfpad = zeros(Float64,nxy,nxy)
        psfpad[1:nxypsf,1:nxypsf] = psf[:,:]
        psfcbe = 1./(abs(fft(psfpad)).^2+mu)
        x = real(ifft(psfcbe.*fft(wlt_b)))

        # tmp1 = 0.0
        # tmp2 = zeros(Float64,nxy,nxy)
        for b in 1:nspat
                    hx = dwt(x,wavelet(spatialwlt[b]))
                    tmp = hx - taut[:,:,1,b]/rhot
                    t[:,:,1,b] = prox_u(tmp,μt/rhot)
                    taut[:,:,1,b] = taut[:,:,1,b] + rhot*(t[:,:,1,b]-hx)
                    # tmp1 = vecnorm([tmp2 (hx-(t)[:,:,z,b])],2)
                    # tmp2 = (hx-(t)[:,:,z,b])
        end
        # tmp2[:] = 0

        tmp = x-taup/rhop
        p = max(0,tmp)
        taup = taup + rhop*(p-x)

    return wlt,x,t,taut,p,taup

end


text:
module PerfTuning

using InteractiveUtils, Libdl, Random
if VERSION ≥ v"0.7.0"
    using Profile
end
using BenchmarkTools, PolrModels

Random.seed!(123)
n, p, J, q = 1_000_000, 5, 7, 5
Xtrue = randn(n, p + q)
X = Xtrue[:, 1:p] # null model matrix
Z = Xtrue[:, p+1:p+q] # covariates to be tested
βtrue = [ones(p); 0.01ones(q)]
θ = collect(1.0:J-1)
α = [θ[1]; [log(θ[i] - θ[i-1]) for i in 2:J-1]]
link = LogitLink() # LogitLink(), ProbitLink(), CauchitLink(), CloglogLink()

@info "generate responses by rpolr"
Y = rpolr(Xtrue, βtrue, θ, link)
# @code_warntype rpolyr(X[1, :], β, θ, :logit)
# @benchmark rpolyr(X, β, θ, :logit)
# Profile.clear_malloc_data()
# Profile.clear()
# @profile rpolyr(X, β, θ, :logit)
# Profile.print(format=:flat)

m = PolrModel(X, Y, link)
println()
@info "loglikelihood calculated by `polrfun!`"
@show polrfun!(m, false, false)
# println()
# @info "type inference on `polrfun`"
# @code_warntype polrfun!(m, true, true)
println()
@info "benchmark `polrfun`"
display(@benchmark polrfun!($m, true, false))
println(); println()
@info "profiling `polrfun`"
Profile.clear()
@profile polrfun!(m, true, true)
Profile.print(format=:flat)

@info "optimization"
# Gradient based: LD_LBFGS, :LD_MMA, :LD_SLSQP, :LD_CCSAQ, :LD_TNEWTON_PRECOND_RESTART, :LD_TNEWTON_PRECOND, :LD_TNEWTON_RESTART, :LD_VAR2, :LD_VAR1
# Gradient free: :LN_COBYLA
# solver = NLoptSolver(algorithm=:LD_SLSQP)
solver = IpoptSolver() # more stable but take a lot more iterations
# solver = IpoptSolver(mehrotra_algorithm="yes")
dd = polr(X, Y, link, solver)
display(coeftable(dd))

@info "score test"
ts = PolrScoreTest(dd, Z)
@show polrtest(ts)
@info "type inference"
@code_warntype polrtest(ts)
println()
@info "benchmark `polrtest`"
display(@benchmark polrtest($ts))
println(); println()
@info "profiling `polrtest`"
Profile.clear()
@profile polrtest(ts)
Profile.print(format=:flat)

end