text:
export sparse_coding

import ..SparseCoding

using mlpack._Internal.io

import mlpack_jll
const sparse_codingLibrary = mlpack_jll.libmlpack_julia_sparse_coding

# Call the C binding of the mlpack sparse_coding binding.
function sparse_coding_mlpackMain()
  success = ccall((:sparse_coding, sparse_codingLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module sparse_coding_internal
  import ..sparse_codingLibrary

import ...SparseCoding

# Get the value of a model pointer parameter of type SparseCoding.
function IOGetParamSparseCoding(paramName::String)::SparseCoding
  SparseCoding(ccall((:IO_GetParamSparseCodingPtr, sparse_codingLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type SparseCoding.
function IOSetParamSparseCoding(paramName::String, model::SparseCoding)
  ccall((:IO_SetParamSparseCodingPtr, sparse_codingLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeSparseCoding(stream::IO, model::SparseCoding)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeSparseCodingPtr, sparse_codingLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeSparseCoding(stream::IO)::SparseCoding
  buffer = read(stream)
  SparseCoding(ccall((:DeserializeSparseCodingPtr, sparse_codingLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    sparse_coding(; [atoms, initial_dictionary, input_model, lambda1, lambda2, max_iterations, newton_tolerance, normalize, objective_tolerance, seed, test, training, verbose])

An implementation of Sparse Coding with Dictionary Learning, which achieves
sparsity via an l1-norm regularizer on the codes (LASSO) or an (l1+l2)-norm
regularizer on the codes (the Elastic Net).  Given a dense data matrix X with d
dimensions and n points, sparse coding seeks to find a dense dictionary matrix D
with k atoms in d dimensions, and a sparse coding matrix Z with n points in k
dimensions.

The original data matrix X can then be reconstructed as Z * D.  Therefore, this
program finds a representation of each point in X as a sparse linear combination
of atoms in the dictionary D.

The sparse coding is found with an algorithm which alternates between a
dictionary step, which updates the dictionary D, and a sparse coding step, which
updates the sparse coding matrix.

Once a dictionary D is found, the sparse coding model may be used to encode
other matrices, and saved for future usage.

To run this program, either an input matrix or an already-saved sparse coding
model must be specified.  An input matrix may be specified with the `training`
option, along with the number of atoms in the dictionary (specified with the
`atoms` parameter).  It is also possible to specify an initial dictionary for
the optimization, with the `initial_dictionary` parameter.  An input model may
be specified with the `input_model` parameter.

As an example, to build a sparse coding model on the dataset `data` using 200
atoms and an l1-regularization parameter of 0.1, saving the model into `model`,
use

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> _, _, model = sparse_coding(atoms=200, lambda1=0.1,
            training=data)
```

Then, this model could be used to encode a new matrix, `otherdata`, and save the
output codes to `codes`:

```julia
julia> using CSV
julia> otherdata = CSV.read("otherdata.csv")
julia> codes, _, _ = sparse_coding(input_model=model,
            test=otherdata)
```

# Arguments

 - `atoms::Int`: Number of atoms in the dictionary.  Default value `15`.

 - `initial_dictionary::Array{Float64, 2}`: Optional initial dictionary
      matrix.
 - `input_model::SparseCoding`: File containing input sparse coding
      model.
 - `lambda1::Float64`: Sparse coding l1-norm regularization parameter.
      Default value `0`.

 - `lambda2::Float64`: Sparse coding l2-norm regularization parameter.
      Default value `0`.

 - `max_iterations::Int`: Maximum number of iterations for sparse coding
      (0 indicates no limit).  Default value `0`.

 - `newton_tolerance::Float64`: Tolerance for convergence of Newton
      method.  Default value `1e-06`.

 - `normalize::Bool`: If set, the input data matrix will be normalized
      before coding.  Default value `false`.

 - `objective_tolerance::Float64`: Tolerance for convergence of the
      objective function.  Default value `0.01`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `test::Array{Float64, 2}`: Optional matrix to be encoded by trained
      model.
 - `training::Array{Float64, 2}`: Matrix of training data (X).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `codes::Array{Float64, 2}`: Matrix to save the output sparse codes of
      the test matrix (--test_file) to.
 - `dictionary::Array{Float64, 2}`: Matrix to save the output dictionary
      to.
 - `output_model::SparseCoding`: File to save trained sparse coding model
      to.

"""
function sparse_coding(;
                       atoms::Union{Int, Missing} = missing,
                       initial_dictionary = missing,
                       input_model::Union{SparseCoding, Missing} = missing,
                       lambda1::Union{Float64, Missing} = missing,
                       lambda2::Union{Float64, Missing} = missing,
                       max_iterations::Union{Int, Missing} = missing,
                       newton_tolerance::Union{Float64, Missing} = missing,
                       normalize::Union{Bool, Missing} = missing,
                       objective_tolerance::Union{Float64, Missing} = missing,
                       seed::Union{Int, Missing} = missing,
                       test = missing,
                       training = missing,
                       verbose::Union{Bool, Missing} = missing,
                       points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, sparse_codingLibrary), Nothing, ());

  IORestoreSettings("Sparse Coding")

  # Process each input argument before calling mlpackMain().
  if !ismissing(atoms)
    IOSetParam("atoms", convert(Int, atoms))
  end
  if !ismissing(initial_dictionary)
    IOSetParamMat("initial_dictionary", initial_dictionary, points_are_rows)
  end
  if !ismissing(input_model)
    sparse_coding_internal.IOSetParamSparseCoding("input_model", convert(SparseCoding, input_model))
  end
  if !ismissing(lambda1)
    IOSetParam("lambda1", convert(Float64, lambda1))
  end
  if !ismissing(lambda2)
    IOSetParam("lambda2", convert(Float64, lambda2))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(newton_tolerance)
    IOSetParam("newton_tolerance", convert(Float64, newton_tolerance))
  end
  if !ismissing(normalize)
    IOSetParam("normalize", convert(Bool, normalize))
  end
  if !ismissing(objective_tolerance)
    IOSetParam("objective_tolerance", convert(Float64, objective_tolerance))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("codes")
  IOSetPassed("dictionary")
  IOSetPassed("output_model")
  # Call the program.
  sparse_coding_mlpackMain()

  return IOGetParamMat("codes", points_are_rows),
         IOGetParamMat("dictionary", points_are_rows),
         sparse_coding_internal.IOGetParamSparseCoding("output_model")
end


text:
"""
    Gamma(α,θ)

The *Gamma distribution* with shape parameter `α` and scale `θ` has probability density
function

```math
f(x; \\alpha, \\theta) = \\frac{x^{\\alpha-1} e^{-x/\\theta}}{\\Gamma(\\alpha) \\theta^\\alpha},
\\quad x > 0
```

```julia
Gamma()          # Gamma distribution with unit shape and unit scale, i.e. Gamma(1, 1)
Gamma(α)         # Gamma distribution with shape α and unit scale, i.e. Gamma(α, 1)
Gamma(α, θ)      # Gamma distribution with shape α and scale θ

params(d)        # Get the parameters, i.e. (α, θ)
shape(d)         # Get the shape parameter, i.e. α
scale(d)         # Get the scale parameter, i.e. θ
```

External links

* [Gamma distribution on Wikipedia](http://en.wikipedia.org/wiki/Gamma_distribution)

"""
struct Gamma{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    θ::T
    Gamma{T}(α, θ) where {T} = new{T}(α, θ)
end

function Gamma(α::T, θ::T; check_args::Bool=true) where {T <: Real}
    @check_args Gamma (α, α > zero(α)) (θ, θ > zero(θ))
    return Gamma{T}(α, θ)
end

Gamma(α::Real, θ::Real; check_args::Bool=true) = Gamma(promote(α, θ)...; check_args=check_args)
Gamma(α::Integer, θ::Integer; check_args::Bool=true) = Gamma(float(α), float(θ); check_args=check_args)
Gamma(α::Real; check_args::Bool=true) = Gamma(α, one(α); check_args=check_args)
Gamma() = Gamma{Float64}(1.0, 1.0)

@distr_support Gamma 0.0 Inf

#### Conversions
convert(::Type{Gamma{T}}, α::S, θ::S) where {T <: Real, S <: Real} = Gamma(T(α), T(θ))
Base.convert(::Type{Gamma{T}}, d::Gamma) where {T<:Real} = Gamma{T}(T(d.α), T(d.θ))
Base.convert(::Type{Gamma{T}}, d::Gamma{T}) where {T<:Real} = d

#### Parameters

shape(d::Gamma) = d.α
scale(d::Gamma) = d.θ
rate(d::Gamma) = 1 / d.θ

params(d::Gamma) = (d.α, d.θ)
partype(::Gamma{T}) where {T} = T

#### Statistics

mean(d::Gamma) = d.α * d.θ

var(d::Gamma) = d.α * d.θ^2

skewness(d::Gamma) = 2 / sqrt(d.α)

kurtosis(d::Gamma) = 6 / d.α

function mode(d::Gamma)
    (α, θ) = params(d)
    α >= 1 ? θ * (α - 1) : error("Gamma has no mode when shape < 1")
end

function entropy(d::Gamma)
    (α, θ) = params(d)
    α + loggamma(α) + (1 - α) * digamma(α) + log(θ)
end

mgf(d::Gamma, t::Real) = (1 - t * d.θ)^(-d.α)

cf(d::Gamma, t::Real) = (1 - im * t * d.θ)^(-d.α)

function kldivergence(p::Gamma, q::Gamma)
    # We use the parametrization with the scale θ
    αp, θp = params(p)
    αq, θq = params(q)
    θp_over_θq = θp / θq
    return (αp - αq) * digamma(αp) - loggamma(αp) + loggamma(αq) -
        αq * log(θp_over_θq) + αp * (θp_over_θq - 1)
end

#### Evaluation & Sampling

@_delegate_statsfuns Gamma gamma α θ

gradlogpdf(d::Gamma{T}, x::Real) where {T<:Real} =
    insupport(Gamma, x) ? (d.α - 1) / x - 1 / d.θ : zero(T)

function rand(rng::AbstractRNG, d::Gamma)
    if shape(d) < 1.0
        # TODO: shape(d) = 0.5 : use scaled chisq
        return rand(rng, GammaIPSampler(d))
    elseif shape(d) == 1.0
        return rand(rng, Exponential(d.θ))
    else
        return rand(rng, GammaGDSampler(d))
    end
end

function sampler(d::Gamma)
    if shape(d) < 1.0
        # TODO: shape(d) = 0.5 : use scaled chisq
        return GammaIPSampler(d)
    elseif shape(d) == 1.0
        return sampler(Exponential(d.θ))
    else
        return GammaGDSampler(d)
    end
end

#### Fit model

struct GammaStats <: SufficientStats
    sx::Float64      # (weighted) sum of x
    slogx::Float64   # (weighted) sum of log(x)
    tw::Float64      # total sample weight

    GammaStats(sx::Real, slogx::Real, tw::Real) = new(sx, slogx, tw)
end

function suffstats(::Type{<:Gamma}, x::AbstractArray{T}) where T<:Real
    sx = zero(T)
    slogx = zero(T)
    for xi = x
        sx += xi
        slogx += log(xi)
    end
    GammaStats(sx, slogx, length(x))
end

function suffstats(::Type{<:Gamma}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Real
    n = length(x)
    if length(w) != n
        throw(DimensionMismatch("Inconsistent argument dimensions."))
    end

    sx = zero(T)
    slogx = zero(T)
    tw = zero(T)
    for i = 1:n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        sx += wi * xi
        slogx += wi * log(xi)
        tw += wi
    end
    GammaStats(sx, slogx, tw)
end

function gamma_mle_update(logmx::Float64, mlogx::Float64, a::Float64)
    ia = 1 / a
    z = ia + (mlogx - logmx + log(a) - digamma(a)) / (abs2(a) * (ia - trigamma(a)))
    1 / z
end

function fit_mle(::Type{<:Gamma}, ss::GammaStats;
    alpha0::Float64=NaN, maxiter::Int=1000, tol::Float64=1e-16)

    mx = ss.sx / ss.tw
    logmx = log(mx)
    mlogx = ss.slogx / ss.tw

    a::Float64 = isnan(alpha0) ? (logmx - mlogx)/2 : alpha0
    converged = false

    t = 0
    while !converged && t < maxiter
        t += 1
        a_old = a
        a = gamma_mle_update(logmx, mlogx, a)
        converged = abs(a - a_old) <= tol
    end

    Gamma(a, mx / a)
end

ext: jl


text:
using Random
using DABenchmarks
using LaTeXStrings

Random.seed!(0)

using Gurobi
using ProgressiveHedgingSolvers

function prepare_dbenchmark(nscenarios::Integer, nsamples::Integer; timeout::Integer = 1000)
    ph = ProgressiveHedgingSolver(()->GurobiSolver(OutputFlag=0), penalty = :adaptive, ζ = 100, log = false, τ = 1e-5)
    solvers = [ph]
    solvernames = ["Synchronous progressive-hedging"]
    # Create Day-ahead benchmark
    return Scaling(solvernames, solvers, nscenarios, nsamples; timeout = timeout)
end
@info "Running progressive-hedging benchmarks on 1 worker..."
dab = prepare_dbenchmark(1000, 5); warmup_benchmarks!(dab); run_benchmarks!(dab); save_results!(dab, "ph_1.json")
exit()

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
@testset "basic predicate tests" begin
    @test pred1(1,1)
    @test !pred1(1,2)
    @test !pred2(1,1)
    @test pred2(1,2)
    @test pred3([1,2],[-1,-2])
    @test !pred3([1,2],[2,3])
end

@testset "is predicate equivalence relation" begin
    @test check_for_equivalence_relation(pred1, -11:10)
    @test check_for_equivalence_relation(pred2, -7:10)
    @test check_for_equivalence_relation(pred3, [[i,j] for i in -7:10 for j in -2:8])
    @test !check_for_equivalence_relation(pred4, -7:10)
end

@testset "mappings" begin
    @test m1(1) == [1]
    @test m1(1) == m1.f(1)
end

@testset "ReduceMap" begin
    @test_throws TypeError ReduceMap(1 => 2.0)
    @test ReduceMap(1=>2).map == OrderedDict(1=>2)
    @test ReduceMap(OrderedDict(1=>2)).map == OrderedDict(1=>2)
    @test ReduceMap(OrderedDict(1=>UInt32(2))).map == OrderedDict(1=>UInt32(2))
    @test ReduceMap(1=>2, 3=>4).map == OrderedDict(1=>2,3=>4)
    @test ReduceMap(i => i+1 for i in 1:2) == OrderedDict(1=>2,2=>3)
    @test ReduceMap(zip([1,3],[2,4])).map == OrderedDict(1=>2, 3=>4)
    @test cl2.map == cl2_res_index
    @test toDirectMap(cl2, vl2) == cl2_res_direct
    @test toIndexMap(toDirectMap(cl2, vl2),vl2) == cl2_res_index
    @test minimal_set(cl2) == sort(collect(unique(values(cl2.map))))
    @test minimal_set(cl2, sorted=true) == UInt32.([1,2,3,4,8,9])
end

@testset "ExpandMap" begin
    em1 = ExpandMap(cl1)
    em2 = ExpandMap(cl2)
    em2_2 = ExpandMap(cl2, vl2)
    @test minimal_set(em2) == sort(collect(unique(values(cl2.map))))
    @test minimal_set(em2, sorted=true) == UInt32.([1,2,3,4,8,9])
    @test sort(collect(keys(ExpandMap(cl1, sorted=true).map))) == [1,2,3,4,5]
    @test Int.(sort(collect(keys(em2.map)))) == [1,2,3,4,8,9]
    @test all(sort(collect(values(em2.map))) .== [[-3,3],[-2],[-1,1],[0],[4],[5]])
    @test ExpandMap(toDirectMap(cl2, vl2)) == em2_2
    @test ExpandMap(cl2, vl2) == em2_2
end

@testset "build equivalence list" begin
    indl1 = 1:5
    eqC_1 = EquivalenceClasses(pred1, indl1)
    sol1_cl = [1,2,3,4,5]
    sol1 = OrderedDict{Int64,Array{Int64,1}}(zip(1:5,[[], [], [], [], []]))
    @test typeof(eqC_1) === ReduceMap{Int64,UInt32}
    @test typeof(eqC_1.map) === OrderedDict{Int64,UInt32}
    @test sort(collect(values(eqC_1.map))) == sol1_cl
end

@testset "complete tests" begin
    #pred1 = EquivalenceClassesConstructor.Predicate((x,y) -> (x .== -y || x .== (y .+ 10)))
    #test_med_ind = (i for i in -300:300);
    #test_large_ind = ([el...] for el in Iterators.product(-30:30, -10:30, -20:30))
end


text:
#=
        Float65
=#

type Float65 <: Real
   fp::Float64

   Float65(fp::Float64) = new(project(fp))
end

@inline getstate(x::Float65) = tst_ebit(x.fp)

function clearstate(x::Float65)
   x.fp = clr_ebit(x.fp)
   x
end

function setstate(x::Float65)
   x.fp = set_ebit(x.fp)
   x
end
setstate(x::Float65, s::Bool) = (s ? setstate(x) : x)

clearstate(x::Float64) = convert(Float65,x)
setstate(x::Float64) = setstate(convert(Float65,x))
setstate(x::Float64,s::Bool) = setstate(convert(Float65,x),s)
getstate(x::Float64) = false

@vectorize_1arg Real Float65

@vectorize_1arg Float65 getstate
@vectorize_1arg Float65 setstate
@vectorize_1arg Float65 clearstate

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#    dispatch narrowing
#=
                                     subobject classifiers
   leaf types of shared abstraction
=#
#=
immutable SituatedFloat  <: MetastableFloat
    fp::Float64
end
immutable StatefulFloat  <: MetastableFloat
    fp::Float64
end

convert{F<:MetastableFloat}(::Type{F}, x::SituatedFloat) = x.fp
convert{F<:Float64}(::Type{F}, x::StatefulFloat) = x.fp
=#


text:
#####
##### Helper functions for working with DiffResults
#####
##### Only included when required by AD wrappers.

import .DiffResults

"""
$(SIGNATURES)

Allocate a DiffResults buffer for a gradient, taking the element type of `x` into account
(heuristically).
"""
function _diffresults_buffer(ℓ, x)
    T = eltype(x)
    S = T <: Real ? float(Real) : Float64 # heuristic
    DiffResults.MutableDiffResult(zero(S), (Vector{S}(undef, dimension(ℓ)), ))
end

"""
$(SIGNATURES)

Extract a return value for [`logdensity_and_gradient`](@ref) from a DiffResults buffer,
constructed with [`diffresults_buffer`](@ref). Gradient is not copied as caller created the
vector.
"""
function _diffresults_extract(diffresult::DiffResults.DiffResult)
    DiffResults.value(diffresult)::Real, DiffResults.gradient(diffresult)
end


text:
export SteadyStates
"""
Some steady state properties of common sequences.
"""
module SteadyStates
using ..MRphy
using ..Unitful, ..UnitfulMR

export Signal
"""
Analytical expressions of common steady states sequences signals.
"""
module Signal
using ..MRphy
using ..Unitful, ..UnitfulMR

"""
    bSSFP(α; TR, Δf, T1, T2)
`TE=0` Steady state bSSFP signal. 10.1002/jmri.24163, eq.(4), with `ϕ=2π*Δf*TR`.

*INPUTS*:
- `α::Real` (1,), tip angle in degree;
*KEYWORDS*:
- `TR::T0D` (1,), repetition time;
- `Δf::F0D` (1,), off-resonance in Hz;
- `T1::T0D` (1,), longitudinal relaxation coefficient;
- `T2::T0D` (1,), transverse relaxation coefficient;
*OUTPUTS*:
- `sig::Complex` (1,), steady-state signal.

See also: [`SPGR`](@ref), [`STFR`](@ref).
"""
function bSSFP(α::Real;
               TR::T0D,
               Δf::F0D=0u"Hz",
               T1::T0D=1u"s",
               T2::T0D=(6e-2)u"s")
  E1, E2, ϕ = exp(-TR/T1), exp(-TR/T2), 2π*Δf*TR # jmri.24163 typoed on E₁, E₂.
  C, D = E2*(E1-1)*(1+cosd(α)), (1-E1*cosd(α)) - (E1-cosd(α))*E2^2
  sig = sind(α) * (1-E1)*(1-E2*exp(-1im*ϕ)) / (C*cos(ϕ) + D)
  return sig
end

"""
    SPGR(α; TR, T1)
`TE=0` Steady state SPGR signal. 10.1002/mrm.1910130109, eq.(1), ideal spoiling.

*INPUTS*:
- `α::Real` (1,), tip angle in degree;
*KEYWORDS*:
- `TR::T0D` (1,), repetition time;
- `T1::T0D` (1,), longitudinal relaxation coefficient;
*OUTPUTS*:
- `sig::Real` (1,), steady-state signal.

See also: [`bSSFP`](@ref), [`STFR`](@ref).
"""
SPGR(α::Real; TR::T0D, T1::T0D=1u"s") =
  sind(α) * (1-exp(-TR/T1))/(1-cosd(α)*exp(-TR/T1));

"""
    STFR(α, β; ϕ, Δf, T1, T2, Tg, Tf)
`TE=0` Steady state STFR signal. 10.1002/mrm.25146, eq.(2), ideal spoiling.

*INPUTS*:
- `α::Real` (1,), tip-down angle in degree;
- `β::Real` (1,), tip-up angle in degree;
*KEYWORDS*:
- `ϕ::Real` (1,), phase of the tip-up pulse in radians;
- `Δf::F0D` (1,), off-resonance in Hz;
- `T1::T0D` (1,), longitudinal relaxation coefficient;
- `T2::T0D` (1,), transverse relaxation coefficient;
- `Tg::T0D` (1,), duration of gradient crusher;
- `Tf::T0D` (1,), duration of free precession in each TR;
*OUTPUTS*:
- `sig::Real` (1,), steady-state signal.

See also: [`bSSFP`](@ref), [`SPGR`](@ref).
"""
function STFR(α::Real, β::Real;
              ϕ::Real=0,
              Δf::F0D=0u"Hz",
              T1::T0D=1u"s",
              T2::T0D=(6e-2)u"s",
              Tg::T0D=(2e-3)u"s",
              Tf::T0D=(1e-2)u"s")
  cΔϕ = cos(2π*Δf*Tf - ϕ) # θf ≔ 2π*Δf*Tf
  Eg, Ef1, Ef2 = exp(-Tg/T1), exp(-Tf/T1), exp(-Tf/T2)
  sα, cα, sβ, cβ = sind(α), cosd(α), sind(β), cosd(β)

  sig = sα * (Eg*(1-Ef1)*cβ + (1-Eg)) / (1-Eg*Ef2*sα*sβ*cΔϕ-Eg*Ef1*cα*cβ)
  return sig
end

end # module Signal

export RFSpoiling
"""
Tools for simulating RF spoiling in gradient echo sequences.
"""
module RFSpoiling
using ..MRphy
using ..Unitful, ..UnitfulMR
"""
    QuadPhase(nTR::Integer, C::Real, B::Real, A::Real)
Quadratically cycling phases in (°): Φ(n) = mod.(C⋅n² + B⋅n + A, 360), n=0:nTR-1
"""
QuadPhase(nTR::Integer, C::Real, B::Real=0, A::Real=0) =
  mod.(([0:nTR-1 ...]|>n->C*n.^2 .+B*n .+A), 360)

"""
    FZstates(Φ, α; TR, T1, T2, FZ)
𝐹, 𝑍 from: 10.1002/(SICI)1099-0534(1999)11:5<291::AID-CMR2>3.0.CO;2-J, eq.(7,8).
eq.(7) refined to `÷√(2)`, instead of `÷2`, as in 10.1002/mrm20736: eq.(2).

Assume constant gradient spoiling of m⋅2π dephasing in each TR, m∈ℤ.
In practice, if dephase is a constant but not exactly m⋅2π, the resulting states
can be computed by convolving a sinc with the m⋅2π dephased results.

*INPUTS*:
- `Φ::TypeND(Real,[2])` (nC,nTR), nC: #`C` as `C` in `QuadPhase`. Typically, one
  simulates a range of `C`s, picking a `C` yielding a signal intensity equals to
  that of ideal spgr spoiling.
- `α::Real` (1,), flip-angle.
*KEYWORDS*:
- `TR::T0D` (1,), repetition time;
- `T1::T0D` (1,), longitudinal relaxation coefficient;
- `T2::T0D` (1,), transverse relaxation coefficient;
- `FZ::NamedTuple`, `(Fs,Fcs,Zs)`, simulate from prescribed states if given:\\
  `Fs ::TypeND(Complex,[2])`, transversal dephasing states, 𝐹ₙ;\\
  `Fcs::TypeND(Complex,[2])`, conjugate transversal dephasing states, 𝐹₋ₙ*;\\
  `Zs ::TypeND(Complex,[2])`, longitudinal states, 𝑍ₙ;
*OUTPUTS*:
- `FZ::NamedTuple`, `(Fs,Fcs,Zs)`, simulation results.
"""
function FZstates(Φ ::TypeND(Real,[2]), α::Real;
                  TR::T0D=(50e-3)u"s",
                  T1::T0D=(1.470)u"s",
                  T2::T0D=(71e-3)u"s",
                  FZ::Union{Nothing, NamedTuple}=nothing)
  if isnothing(FZ)
    (Fs, Fcs, Zs) = (zeros(Complex{Float64}, size(Φ)) for _ in 1:3)
    Zs[:, 1] .= 1
    FZ = (Fs=Fs, Fcs=Fcs, Zs=Zs)
  end
  eltype(Φ) <: AbstractFloat || (Φ = float(Φ))
  eltype(Fs)==Complex{Float64}  || (Fs  = Complex{Float64}(Fs))
  eltype(Fcs)==Complex{Float64} || (Fcs = Complex{Float64}(Fcs))
  eltype(Zs)==Complex{Float64}  || (Zs  = Complex{Float64}(Zs))

  size(FZ.Fs)==size(FZ.Fcs)==size(FZ.Zs) || throw(DimensionMismatch)

  E1, E2 = exp(-TR/T1), exp(-TR/T2)

  return FZ = _FZstates(Φ, α; E1=E1, E2=E2, FZ...)
end

function _FZstates(Φ ::TypeND(AbstractFloat,[2]), α::Real;
                   E1::Real,
                   E2::Real,
                   Fs ::TypeND(Complex,[2]),
                   Fcs::TypeND(Complex,[2]),
                   Zs ::TypeND(Complex,[2]))
  Φ .= mod.(Φ, 360)
  Φ .*= π/180

  Fs⁻, Fcs⁻, Zs⁻ = similar(Fs), similar(Fcs), similar(Zs) # pre-rf states
  Fcs⁻[:,end] .= 0

  c½α², s½α², cα, sα = cosd(α/2)^2, sind(α/2)^2, cosd(α), sind(α)
  ieΦsαr2, e2Φs½α² = im*exp.(im*Φ)*sα/√(2), exp.(im*2*Φ)*s½α²

  @inbounds for t in axes(Φ,2)
    # relaxation and updating pre-rf states; eq.(8)
    Fs⁻ .= E2.*[conj(Fcs[:,2]) Fs[:,1:end-1]]
    Fcs⁻[:, 1:end-1] .= E2.*Fcs[:,2:end]
    Zs⁻ .= E1.*Zs
    Zs⁻[:,1] .+= (1-E1)

    # excitation; eq.(7)
    ieϕsαr2, e2ϕs½α² = view(ieΦsαr2,:,t), view(e2Φs½α²,:,t)

    Fs  .=          c½α².*Fs⁻ .+ e2ϕs½α².*Fcs⁻ .-       ieϕsαr2.*Zs⁻
    Fcs .= conj(e2ϕs½α²).*Fs⁻ .+    c½α².*Fcs⁻ .- conj(ieϕsαr2).*Zs⁻
    Zs  .= conj(ieϕsαr2).*Fs⁻ .+ ieϕsαr2.*Fcs⁻ .+            cα.*Zs⁻
  end

  return (Fs=Fs, Fcs=Fcs, Zs=Zs)
end

end # module RFSpoiling

end # module SteadyStates


text:
"""
    Thirty360Excel()

**US (NASD) 30/360** day count convention, as computed via Microsoft Excel `YEARFRAC` with the basis option of `0`.

This differs from [`Thirty360`](@ref) when:
*  if the start date is the last day of February, then
  -  ``d_1`` is 30, and
  - if the end date is also the last day of February ``d_2`` is also 30.

# Reference
- [Microsoft Excel `YEARFRAC` function](https://support.office.com/en-us/article/yearfrac-function-3844141e-c76d-4143-82b6-208454ddc6a8)
- [David A. Wheeler (2008) "YEARFRAC Incompatibilities between Excel 2007 and OOXML (OXML), and the Definitions Actually Used by Excel 2007"](https://dwheeler.com/yearfrac/excel-ooxml-yearfrac.pdf)
"""
struct Thirty360Excel <: DayCount
end
function frac(startdate::Date, enddate::Date, dc::Thirty360Excel)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    end

    y1 = year(startdate)
    y2 = year(enddate)
    m1 = month(startdate)
    m2 = month(enddate)
    d1 = day(startdate)
    d2 = day(enddate)
    if d1 >= 30
        d1 = 30
        if d2 >= 30
            d2 = 30
        end
    elseif m1 == 2 && startdate == lastdayofmonth(startdate)
        d1 = 30
        if m2 == 2 && enddate == lastdayofmonth(enddate)
            d2 = 30
        end
    end
    return thirty360(y2-y1,m2-m1,d2-d1)
end

"""
    ActualActualExcel()

**Actual/Actual** day count convention, as computed via Microsoft Excel `YEARFRAC` with the basis option of `1`.

The year fraction is
```math
\\frac{\\text{# of days}}{D}
```
where:
* if start date and end date are in the same calendar year, then ``D`` is the number of days in that calendar year,
* otherwise if the end date is less than or equal to one year after the start date, then
  - ``D = 366`` if the interval includes February 29 (including both the start and end dates), or
  - ``D = 365`` otherwise.
* otherwise it is the average length of the years included in the interval (this does not depend on where the start or end dates fall within those years).

# Reference
- [Microsoft Excel `YEARFRAC` function](https://support.office.com/en-us/article/yearfrac-function-3844141e-c76d-4143-82b6-208454ddc6a8)
- [David A. Wheeler (2008) "YEARFRAC Incompatibilities between Excel 2007 and OOXML (OXML), and the Definitions Actually Used by Excel 2007"](https://dwheeler.com/yearfrac/excel-ooxml-yearfrac.pdf)
"""
struct ActualActualExcel <: DayCount
end
function frac(startdate::Date, enddate::Date, dc::ActualActualExcel)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    end
    y1 = year(startdate)
    y2 = year(enddate)
    if y1 == y2
        return Dates.value(enddate-startdate) / daysinyear(y1)
    elseif startdate + Year(1) >= enddate
        if (isleapyear(y1) && startdate <= Date(y1,2,29)) ||
            (isleapyear(y2) && Date(y2,2,29) <= enddate)
            return Dates.value(enddate-startdate) / 366
        else
            return Dates.value(enddate-startdate) / 365
        end
    else
        yrange = y1:y2
        return Dates.value(enddate-startdate) / (sum(daysinyear, yrange) / length(yrange))
    end
end


text:
#-----------------------------------------------------------------------# IndexedTable
Tables.istable(::Type{<:IndexedTable}) = true

Tables.materializer(t::IndexedTable) = table

Tables.columnaccess(::Type{<:IndexedTable}) = Tables.columnaccess(StructArray)
Tables.columns(t::IndexedTable) = Tables.columns(columns(t))

Tables.rowaccess(::Type{<:IndexedTable}) = Tables.rowaccess(StructArray)
Tables.rows(t::IndexedTable) = Tables.rows(rows(t))

# table(x; copy=false, kw...) = table(Tables.columntable(x); copy=copy, kw...)


text:
using LoadFlint
using Test
using Libdl

# Just successfully loading LoadFlint is a good first test

@testset "LoadFlint" begin
    @test occursin("libgmp", LoadFlint.libgmp)
    @test occursin("libflint", LoadFlint.libflint)

    l = dllist()
    @test length(filter(x->occursin("libflint", x), l)) == 1
    x = Ptr{Nothing}(0)
    @test (x = ccall((:flint_malloc,LoadFlint.libflint), Ptr{Nothing}, (UInt,), 8)) != C_NULL
    @test ccall((:flint_free,LoadFlint.libflint), Nothing ,(Ptr{Nothing},),x) == nothing

    # TODO maybe test fmpz?
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 06f9a7b0-7a27-11eb-24d2-2398ba0ec9df
using GMT

# ╔═╡ cee6d350-7a23-11eb-1522-71dfaaac4fb9
md"# Mixing UTM and geographic axes annotations

This is a question that comes up regularly such that GMT has a gallery example (famous [example 28](https://docs.generic-mapping-tools.org/dev/gallery/ex28.html#example-28)) to show how to do it. But even with that example it's not a trivial thing to do.

However it's trivial from the Julia wrapper under the condition that the grid or image to be displayed is referenced internally. If it is not, `grdedit` can be used to assign a referencing system via its option `proj`. The other condition is that grid or image to be displayed is already in memory (so that the internal magicks can work). Given that, [example 28](https://docs.generic-mapping-tools.org/dev/gallery/ex28.html#example-28)) boils down to this (with a little bit less fancy details).
"

# ╔═╡ 1e769320-7a28-11eb-3eef-538e1a5734ea
G = gmtread("@Kilauea.utm.nc");

# ╔═╡ 2251f892-7a28-11eb-206e-b1a4702505e1
md"Create a cpt and display the map
"

# ╔═╡ 77b60470-7a28-11eb-368f-631b2b09db3f
C = makecpt(cmap="copper", range=(0,1500));

# ╔═╡ 45b59260-7a28-11eb-2d3b-5d9eed99c801
imshow(G, cmap=C, shade=true, frame=(axes="WS", annot=true),
	coast=(shore=true, ocean=:lightblue, frame=(axes="EN", annot=true, grid=true)))

# ╔═╡ Cell order:
# ╟─cee6d350-7a23-11eb-1522-71dfaaac4fb9
# ╠═06f9a7b0-7a27-11eb-24d2-2398ba0ec9df
# ╠═1e769320-7a28-11eb-3eef-538e1a5734ea
# ╟─2251f892-7a28-11eb-206e-b1a4702505e1
# ╠═77b60470-7a28-11eb-368f-631b2b09db3f
# ╠═45b59260-7a28-11eb-2d3b-5d9eed99c801


text:
module Juliaset

using Images
using Compat
using JuliaWebAPI

export juliaset

function genjuliaset(c1::Float64=0.0, c2::Float64=0.65, w::Int=200, h::Int=200)
    a = Array(UInt8, w, h)

    c = c1 + c2 * im

    idx = 1
    for r_val in linspace(-2.0, 2.0, w), im_val in linspace(-2.0, 2.0, h)
        z = r_val + im_val*im
        n = 255
        while abs(z) < 10 && n >= 5
            z = z*z + c
            n -= 5
        end
        (idx < length(a)) && (a[idx] = n)
        idx += 1
    end
    grayim(a)
#
#    img = grayim(a)
#    imwrite(img, "out.png")
end

function pngbytes(img)
    iob = IOBuffer()
    writemime(iob, MIME"image/png"(), img)
    takebuf_array(iob)
end

juliaset(c1::Float64, c2::Float64) = pngbytes(genjuliaset(c1, c2))
juliaset(s1::AbstractString, s2::AbstractString) = juliaset(parse(Float64, s1), parse(Float64, s2))
juliaset() = juliaset(rand(), rand())

const PNG_HEADER = @compat Dict{AbstractString,AbstractString}("Content-Type" => "image/png")
const JULIABOX_APIS = [
        (juliaset, false, PNG_HEADER)
    ]

serve_juliabox() = process(JULIABOX_APIS)
serve_local() = process(JULIABOX_APIS, "tcp://127.0.0.1:9999"; bind=true)

end # module


text:
# Original Java version Copyright 2014 Google Inc.
# Julia translation Copyright 2022 Ted Dunning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using S2
using Test
using LinearAlgebra

for i in 1:100
    frame = qr(randn(3, 3)).Q
    x, y, z = [frame[:,i] for i in 1:3]

    @test S2.S1Angle_ZERO == toAngle(S1ChordAngle(z, z))
    @test isapprox(π, radians(toAngle(S1ChordAngle(-z, z))), atol=1e-7)
    @test π / 2 ≈ radians(toAngle(S1ChordAngle(x, z)))
    w = normalize(y + z)
    @test π / 4 ≈ radians(toAngle(S1ChordAngle(w, z)))
end

@test 0.0 == degrees(toAngle(S1ChordAngle(0)))
@test 60.0 ≈ degrees(toAngle(S1ChordAngle(1)))
@test 90.0 ≈ degrees(toAngle(S1ChordAngle(2)))
@test 180.0 ≈ degrees(toAngle(S1ChordAngle(4)))
@test 180.0 ≈ degrees(toAngle(S1ChordAngle(5)))

@test S2.S1Angle_ZERO == toAngle(S2.S1ChordAngle_ZERO)

@test S2.isZero(S2.S1ChordAngle_ZERO)
@test !S2.isNegative(S2.S1ChordAngle_ZERO)
@test !S2.isSpecial(S2.S1ChordAngle_ZERO)
@test !S2.isSpecial(S2.S1ChordAngle_STRAIGHT)
@test S2.S2.isNegative(S2.S1ChordAngle_NEGATIVE)
@test S2.isSpecial(S2.S1ChordAngle_NEGATIVE)
@test S2.isInfinity(S2.S1ChordAngle_INFINITY)
@test S2.isSpecial(S2.S1ChordAngle_INFINITY)

@test 0.0 == radians(toAngle(S1ChordAngle(S2.S1Angle_ZERO)))
@test 4.0 == S1ChordAngle(radians(π)).distance2
@test isapprox(π, radians(toAngle(S1ChordAngle(radians(π)))), atol=1e-7)
@test radians(toAngle(S1ChordAngle(radians(-1)))) < 0.0
@test 1.0 ≈ radians(toAngle(S1ChordAngle(radians(1.0))))

let zero = S1ChordAngle(0)
    degree30 = S1ChordAngle(degrees(30))
    degree60 = S1ChordAngle(degrees(60))
    degree90 = S1ChordAngle(degrees(90))
    degree120 = S1ChordAngle(degrees(120))
    degree180 = S1ChordAngle(degrees(180))
    @test 0.0 == degrees(toAngle(zero + zero))
    @test 0.0 == degrees(toAngle(zero - zero))
    @test 0.0 == degrees(toAngle(degree60 - degree60))
    @test 0.0 == degrees(toAngle(degree180 - degree180))
    @test 0.0 == degrees(toAngle(zero - degree60))
    @test 0.0 == degrees(toAngle(degree30 - degree90))
    @test 60.0 ≈ degrees(toAngle(degree60 + zero))
    @test 60.0 ≈ degrees(toAngle(degree60 - zero))
    @test 60.0 ≈ degrees(toAngle(zero + degree60))
    @test 90.0 ≈ degrees(toAngle(degree30 + degree60))
    @test 90.0 ≈ degrees(toAngle(degree60 + degree30))
    @test 60.0 ≈ degrees(toAngle(degree90 - degree30))
    @test 30.0 ≈ degrees(toAngle(degree90 - degree60))
    @test 180.0 == degrees(toAngle(degree180 + zero))
    @test 180.0 == degrees(toAngle(degree180 - zero))
    @test 180.0 == degrees(toAngle(degree90 + degree90))
    @test 180.0 == degrees(toAngle(degree120 + degree90))
    @test 180.0 == degrees(toAngle(degree120 + degree120))
    @test 180.0 == degrees(toAngle(degree30 + degree180))
    @test 180.0 == degrees(toAngle(degree180 + degree180))
end

let N = 20
    for i in 0:N
        θ = π * i / N
        angle = S1ChordAngle(radians(θ))
        @test isapprox(sin(θ), sin(angle), atol=1e-15)
        @test isapprox(cos(θ), cos(angle), atol=1e-15)
        # Since the tan(x) is unbounded near Pi/4, we map the result back to an
        # angle before comparing.  (The assertion is that the result is equal to
        # the tangent of a nearby angle.)
        @test isapprox(atan(tan(θ)), atan(tan(angle)), atol=1e-15)
    end
end

# Unlike S1Angle, S1ChordAngle can represent 90 and 180 degrees exactly.
let
    angle90 = S1ChordAngle(2)
    angle180 = S1ChordAngle(4)
    begin
        @test 1.0 == sin(angle90)
        @test 0.0 == cos(angle90)
        @test Inf == tan(angle90)
        @test 0.0 == sin(angle180)
        @test -1.0 == cos(angle180)
        @test 0.0 == tan(angle180)
    end
end

#  /**
#   * Verifies that the error bound returned by {@link S1ChordAngle#getS2PointConstructorMaxError} is
#   * large enough.
#   */
#  public void testGetS2PointConstructorMaxError() {
#    for (int iter = 0; iter < 10000; ++iter) {
#      rand.setSeed(iter);
#      S2Point x = randomPoint();
#      S2Point y = randomPoint();
#      if (super.oneIn(10)) {
#        // Occasionally test a point pair that is nearly identical or antipodal.
#        S1Angle r = S1Angle.radians(1e-15 * rand.nextDouble());
#        y = S2EdgeUtil.interpolateAtDistance(r, x, y);
#        if (oneIn(2)) {
#          y = S2Point.neg(y);
#        }
#      }
#      S1ChordAngle dist = new S1ChordAngle(x, y);
#      double error = dist.getS2PointConstructorMaxError();
#      String msg = "angle=" + dist + ", iter=" + iter;
#      assertTrue(msg, S2Predicates.compareDistance(x, y, dist.plusError(error).getLength2()) <= 0);
#      assertTrue(msg, S2Predicates.compareDistance(x, y, dist.plusError(-error).getLength2()) >= 0);
#    }
#  }
#
#  public void testS1AngleConsistency() {
#    // This test checks that the error bounds in the S1ChordAngle constructors
#    // are consistent with the maximum error in S1Angle(x, y).
#    double maxS1AngleError = 3.25 * DBL_EPSILON;
#    for (int iter = 0; iter < 10000; ++iter) {
#      S2Point x = randomPoint();
#      S2Point y = randomPoint();
#      S1ChordAngle dist1 = S1ChordAngle.fromS1Angle(new S1Angle(x, y));
#      S1ChordAngle dist2 = new S1ChordAngle(x, y);
#      double maxError =
#          (maxS1AngleError
#              + dist1.getS1AngleConstructorMaxError()
#              + dist2.getS2PointConstructorMaxError());
#      assertTrue(dist1.compareTo(dist2.plusError(maxError)) <= 0);
#      assertTrue(dist1.compareTo(dist2.plusError(-maxError)) >= 0);
#    }
#  }
#
#  /**
#   * Assert that {@code actual} is almost equal to {@code expected}, within floating point error.
#   */
#  private static void assertNearlyEquals(double expected, double actual) {
#    assertEquals(expected, actual, DOUBLE_ERROR);
#  }
#


text:
using VectorEngine, Test

@noinline function f(p)
    @veprintln(p[])
    p[] = 12
end

function pass(p)
    f(p)
    return
end

x = 10
VectorEngine.@sync @veda pass(Ref(x))
@show x

y = VERefValue(100)
VectorEngine.@sync @veda pass(y)
@show y


text:
immutable CachedTarget <: AbstractCached
    funhash::UInt64
    timestamp::DateTime
end

immutable CacheInfo
    fileversion::VersionNumber
    juliaversion::VersionNumber
    hashsize::UInt
end

# Return the data to be cached for `t`.
# Not meant to be used externally.
cached(t::AbstractTarget) = CachedTarget(t.funhash, t.timestamp)

minorversion(v) = VersionNumber(v.major, v.minor, 0, (), ())

"""
```julia
fixcache()
```

Update the `.maker-cache.jld` file with the present version information
for the cache version number and the Julia version. This may be needed
if the file format changes or the Julia version upgrades. Note that this
may not solve all issues.
"""
fixcache() = getjld() do f
    delete!(f, "CACHEINFO")
    write(f, "CACHEINFO", CacheInfo(CACHEVERSION, minorversion(VERSION), sizeof(Int)))
end

compatible(x::CacheInfo) = x.fileversion == CACHEVERSION &&
    x.juliaversion.major == VERSION.major &&
    x.juliaversion.minor == VERSION.minor &&
    x.hashsize == sizeof(Int)

# Run `fun(f)` on the CACHEFILE with `f` as the opened file.
# Not meant to be used externally.
function getjld(fun::Function)
    if !isfile(CACHEFILE)
        f = jldopen(CACHEFILE, "w")
        write(f, "CACHEINFO", CacheInfo(CACHEVERSION, VERSION, sizeof(Int)))
    else
        f = jldopen(CACHEFILE, "r+")
        if !compatible(read(f, "CACHEINFO"))
            error("Cache version in .maker-cache differs from the active version")
        end
    end
    try
        fun(f)
    finally
        close(f)
    end
end

function fixh5name(name)
    if contains(name, "/")
        return string("M", hash(name))
    else
        return name
    end
end

# Store target `t` information in the cache file.
# Not meant to be used externally.
function updatecache!(f::JLD.JldFile, t::AbstractTarget)
    h5name = fixh5name(t.name)
    if h5name in names(f)
        delete!(f, h5name)
    end
    write(f, h5name, cached(t))
end
updatecache!(t::AbstractTarget) = getjld() do f
    updatecache!(f, t)
end

remove_linenodes(x) = x
remove_linenodes(x::Vector{Any}) = map(remove_linenodes, filter(a -> !isa(a, LineNumberNode), x))
remove_linenodes(e::Expr) = Expr(e.head, remove_linenodes(e.args))

# Find the hash of an anonymous or generic function. The global
# is to have the same hash if the function code is the same. This
# is different than the standard `hash(f::Function)`.
# Not meant to be used externally.
function funhash(f::Function)
    if isdefined(f, :code) # handles anonymous functions
        hash(string(remove_linenodes(Base.uncompressed_ast(f.code))))
    else
        hash(code_lowered(f, ())) # handles generic functions
    end
end
funhash(f::Function, x) = hash(funhash(f), hash(convert(Vector{UTF8String}, x)))

"""
```julia
dependencies(t::AbstractTarget)
```

Return the dependencies of a target.
"""
dependencies(t::AbstractTarget) = [resolvedependency(d) for d in t.dependencies]

function has1arg{T<:AbstractTarget}(t::T)
    if isgeneric(t.action)
        !isempty(methods(t.action, (T,)))
    else
        return length(Base.uncompressed_ast(t.action.code).args[1]) == 1
    end
end

function has1arg(f::Function)
    if isgeneric(f)
        !isempty(methods(f, (AbstractTarget,)))
    else
        return length(Base.uncompressed_ast(f.code).args[1]) == 1
    end
end

"""
```julia
execute(t::AbstractTarget)
```

Run the action for target `t`.
"""
function execute(t::AbstractTarget)
    has1arg(t) ? t.action(t) : t.action()
    t.timestamp = Dates.unix2datetime(time())
    t.isstale = false
    updatecache!(t)
end

"""
```julia
isstale(t::AbstractTarget)
```

Return `true` if target `t` is stale and should be updated.
"""
function isstale(t::AbstractTarget)
    if t.isstale
        return true
    end
    ds = dependencies(t)
    isempty(ds) ? false : timestamp(t) < maximum([timestamp(d) for d in ds])
end

"""
```julia
register(t::AbstractTarget)
```

Register target `t`. `make` only looks for registered targets.
"""
function register(t::AbstractTarget)
    t.description = NEXTDOC[1]
    tasks!(t.name, t)
    NEXTDOC[1] = ""
end

"""
```julia
resolve(s::AbstractString)
```

Return the target registered under name `s`.
"""
resolve(s::AbstractString, default) = get(tasks(), utf8(s), default)
function resolve(s::AbstractString = "default")
    t = resolve(s, nothing)
    t === nothing && error("no rule for target '$s'")
    t
end

"""
```julia
timestamp(t::AbstractTarget)
```

Return the DateTime timestamp for a target.
"""
timestamp(t::AbstractTarget) = t.timestamp

originaltimestamp{T<:AbstractTarget}(::Type{T}, name::AbstractString) = DateTime()

"""
    target{T<:AbstractTarget}(::Type{T}, name::AbstractString, action::Function, dependencies::AbstractArray)

Define and register a target of type `T`.
"""
function target{T<:AbstractTarget}(::Type{T}, name::AbstractString, action::Function, dependencies::AbstractArray)
    dependencies = UTF8String[dependencies...]
    t = resolve(name, nothing)
    fh = funhash(action, dependencies)
    datetime = DateTime()
    if t === nothing
        isstale = true
        # check the cache
        getjld() do f
            h5name = fixh5name(name)
            if h5name in names(f)
                x = read(f[h5name])
                if fh == x.funhash
                    datetime = x.timestamp
                    isstale = false
                end
            else
                datetime = originaltimestamp(T, name)
            end
        end
    else
        isstale = fh != t.funhash
    end
    if t === nothing || fh != t.funhash || isstale
        register(T(name, dependencies, "", action, datetime, fh, isstale))
    end
end
target{T<:AbstractTarget}(::Type{T}, name::AbstractString, action::Function, dependencies::AbstractString) =
    target(T, name, action, [dependencies])

"""
```julia
make(name::AbstractString = "default";
     dryrun::Bool = false, verbose::Bool = false)
make(names::AbstractVector;
     dryrun::Bool = false, verbose::Bool = false)
```

Check, and if needed, update target `name` after updating its dependencies. For
multiple targets specified in `names`, run `make` on each item sequentially.

If keyword argument `verbose` is set, the chain of targets and dependencies is
shown.

If keyword argument `dryrun` is set, the chain of targets and actions
is shown but not run.
"""
function make(t::AbstractTarget, level::Int, dryrun::Bool, verbose::Bool)
    if verbose
        println(("-" ^ level) * " \"$(t.name)\"")
    end
    for d in dependencies(t)
        make(d, level + 1, dryrun, verbose)
    end
    if isstale(t)
        if verbose || dryrun
            println(("*" ^ level) * " Execute \"$(t.name)\"")
        end
        if !dryrun
            execute(t)
        end
    end
    nothing
end
make(s::AbstractString = "default"; dryrun = false, verbose = false) =
    make(resolve(s), 1, dryrun, verbose)
function make(x::AbstractVector; dryrun = false, verbose = false)
    if isempty(x)
        make(dryrun = dryrun, verbose = verbose)
    else
        map(s -> make(resolve(s), 1, dryrun, verbose), x)
    end
end


text: fib(n) = n < 2 ? n : fib(n - 1) + fib(n - 2)

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
using Images
using FileIO
using JLD2
using Statistics: std, mean, var

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

# Funcion para crear el dataset directamente, si existe.
function createDataset(dataset_name::String, extraction::Symbol)
    (inputs, targets) = getInputs("datasets"; extr=extraction);
    println("Tamaños en la generación:");
    println(size(inputs));
    println(size(targets));
    write_dataset(dataset_name,inputs,targets);
    while (!isfile("datasets/faces.data"))
        sleep(1);
    end;
end;

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

# Funcion para convertir una foto (en RGB) en un Array de Float64
function imageToColorArray(image::Array{RGB{Normed{UInt8,8}},2})
    matrix = Array{Float64, 3}(undef, size(image,1), size(image,2), 3)
    matrix[:,:,1] = convert(Array{Float64,2}, red.(image));
    matrix[:,:,2] = convert(Array{Float64,2}, green.(image));
    matrix[:,:,3] = convert(Array{Float64,2}, blue.(image));
    return matrix;
end;

# Funcion para convertir una foto (en RGBA) en un Array de Float64
imageToColorArray(image::Array{RGBA{Normed{UInt8,8}},2}) = imageToColorArray(RGB.(image));

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

# Funcion para leer todas las imagenes de una carpeta dada
function loadFolderImages(folderName::String, extr::Symbol)

    # Comprobar que la foto este en formato .JPEG
    isImageExtension(fileName::String) = any(uppercase(fileName[end-4:end]) .== [".JPEG"]);
    images = [];
    # Para cada fichero que detecte en la carpeta dada
    for fileName in readdir(folderName)
        #println(fileName)
        if isImageExtension(fileName)
            # Leemos la imagen
            image = load(string(folderName, "/", fileName));
            # Comprobar que el archivo cargado sea una imagen en color
            if (typeof(image)==Array{RGBX{Normed{UInt8,8}},2})
                image=convert(Array{RGB{Normed{UInt8,8}},2},image)
            end
            @assert(isa(image, Array{RGBA{Normed{UInt8,8}},2}) || isa(image, Array{RGB{Normed{UInt8,8}},2}))
            # Añadimos la foto al array de imagenes
            push!(images, image);
        end;
    end;

    # Devolvemos el array con las imagenes convertidas a Float64
    if extr==:A1
        return imageToColorArray.(images);
    elseif extr==:A21
        return face_features_1.(images);
    elseif extr==:A22
        return face_features_2.(images);
    elseif extr==:A23
        return face_features_3.(images);
    elseif extr==:A31
        return face_features_masc.(images);
    elseif extr==:A32
        return face_features_masc2.(images);
    elseif extr==:A33
        return face_features_masc3.(images);
    elseif extr==:A34
        return face_features_masc4.(images);
    end;

end;

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

# Funcion para cargar todo el dataset (positivos y negativos)
function loadDataset(folderName::String, extr::Symbol)
    caraDataset = loadFolderImages(string(folderName, "/recortes"),extr);
    negativeDataset = loadFolderImages(string(folderName, "/cara_negativo"),extr);
    mascarillaDataset = loadFolderImages(string(folderName, "/cara_mascarilla"),extr);
    return (caraDataset, negativeDataset, mascarillaDataset);
end;

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

# Funcion que obtiene la varianza, desviacion y media de los 3 canales RGB de
# una imagen dada
function getAttributesFromImage(foto)
    inputs = [];
    for canal in 1:3
        mean_image = mean(foto[:,:,canal]);
        std_image = std(foto[:,:,canal]);

        push!(inputs, mean_image)
        push!(inputs, std_image)
    end;
    return inputs;
end

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

# Funcion que obtiene una matriz NxM, donde N es el numero de elementos
# (positivos + negativos) y M es el numero de columnas (varianza, media y
# desviacion tipica) para cada canal RGB y un vector de targets
function getInputs(path::String; extr::Symbol=:A33)
    # Obtenemos todas las fotos clasificadas en positivas y negativas
    (caraDataset, negativeDataset, mascarillaDataset) = loadDataset(path, extr);
    # Generamos la matriz de inputs y targets
    sizeFaceDataset = size(caraDataset,1);
    sizeMaskDataset = size(mascarillaDataset,1);
    sizeNegativeDataset = size(negativeDataset,1);
    rows = sizeFaceDataset + sizeMaskDataset + sizeNegativeDataset;
    if extr==:A1
        cols = 6;
    elseif extr==:A21
        cols = 42;
    elseif (extr==:A22 || extr==:A23)
        cols = 36;
    elseif (extr==:A31)
        cols = 48;
    elseif (extr==:A32)
        cols = 66;
    elseif (extr==:A33)
        cols = 72;
    elseif (extr==:A34)
        cols = 84;
    end;
    inputs = Array{Float64, 2}(undef, rows, cols);
    targets = Array{String, 1}(undef, rows);

    v1 = false;
    if extr==:A1
        v1 = true;
    end;

    # Generamos la primera parte de la matriz de inputs con los elementos
    # que son positivos
    for i in 1:sizeFaceDataset
        foto = caraDataset[i];
        targets[i] = "Face";
        inputs[i,:] = v1 ? getAttributesFromImage(foto) : getAttrFromImgv2(foto);
    end;

    # Generamos la tercera parte de la matriz de inputs con los elementos
    # que son caras con mascarilla
    for i in 1:sizeMaskDataset
        foto = mascarillaDataset[i];
        targets[sizeFaceDataset+i] = "Mask";
        inputs[sizeFaceDataset+i,:] = v1 ? getAttributesFromImage(foto) : getAttrFromImgv2(foto);
    end;

    # Generamos la segunda parte de la matriz de inputs con los elementos
    # que son negativos
    aux = sizeFaceDataset+sizeMaskDataset;
    for i in 1:sizeNegativeDataset
        foto = negativeDataset[i];
        targets[aux+i] = "NotFace";
        inputs[aux+i,:] = v1 ? getAttributesFromImage(foto) : getAttrFromImgv2(foto);
    end;

    return (inputs,targets)
end

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

# Funcion que dado un nombre de archivo, y unos inputs y targets los escribe
# en el archivo dado.
function write_dataset(file_name::String,inputs::Array{Float64, 2},targets::Array{String, 1})
    f = open(file_name, "w")
    for line in 1:size(inputs,1)
        string_line=""
        for item in inputs[line,:]
            string_line=string(string_line,item,",")
        end
        string_line=string(string_line,targets[line])
        print(f,string_line)
        println(f,"")
    end
end;

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

# Función para recolectar características de un array de imágenes (igual que
# hacíamos antes solo que ahora con un array)
function getAttrFromImgv2(array)
    # Para cada imágen, sacamos las características y las juntamos todas en
    # un solo array de atributos
    attr = [];
    for image = array
        for 𝑋 = getAttributesFromImage(image)
            push!(attr, 𝑋);
        end;
    end;
    return attr;
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_1(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_2(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_3(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_masc(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    counc = image[1:(h ÷ 20 * 5), (w ÷ 20 * 8):(w ÷ 20 * 12)];
    push!(array_of_images, counc);
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);
    # visualizamos los recortes
    #visualize_masc(img, h, w);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_masc2(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    eyes_and_eyebrows = image[(1):(h ÷ 3), (1):(w)];
    push!(array_of_images, eyes_and_eyebrows);
    counc = image[1:(h ÷ 20 * 5), (w ÷ 20 * 8):(w ÷ 20 * 12)];
    push!(array_of_images, counc);
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 12):(w ÷ 20 * 18)];
    push!(array_of_images, right_eye);
    left_eye2 = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 2):(w ÷ 20 * 8)];
    push!(array_of_images, left_eye);
    right_eye2 = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);
    # visualizamos los recortes
    #visualize_masc2(img, h, w);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_masc3(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    counc = image[1:(h ÷ 20 * 5), (w ÷ 20 * 8):(w ÷ 20 * 12)];
    push!(array_of_images, counc);
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 12):(w ÷ 20 * 18)];
    push!(array_of_images, right_eye);
    left_eye2 = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 2):(w ÷ 20 * 8)];
    push!(array_of_images, left_eye);
    right_eye2 = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_eyebrow = image[(1):(h ÷ 20 * 2), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eyebrow);
    right_eyebrow = image[(1):(h ÷ 20 * 2), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eyebrow);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);
    # visualizamos los recortes
    # visualize_masc3(img, h, w);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;

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

# Función que realiza una extracción de características dada una imágen
function face_features_masc4(image::Array{RGB{Normed{UInt8,8}},2})::Array{Array{Float64, 3}}

    # devolvemos un array de imágenes
    array_of_images = [image];
    # tamaños de la imagen para recortes
    (h, w) = size(image);
    # empieza la extracción de características
    eyes_and_eyebrows = image[(1):(h ÷ 3), (1):(w)];
    push!(array_of_images, eyes_and_eyebrows);
    counc = image[1:(h ÷ 20 * 5), (w ÷ 20 * 8):(w ÷ 20 * 12)];
    push!(array_of_images, counc);
    left_eye = image[(h ÷ 20):(h ÷ 3), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eye);
    right_eye = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 12):(w ÷ 20 * 18)];
    push!(array_of_images, right_eye);
    left_eye2 = image[(h ÷ 20 * 2):(h ÷ 20 * 5), (w ÷ 20 * 2):(w ÷ 20 * 8)];
    push!(array_of_images, left_eye);
    right_eye2 = image[(h ÷ 20):(h ÷ 3), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eye);
    left_eyebrow = image[(1):(h ÷ 20 * 2), (w ÷ 20):(w ÷ 20 * 9)];
    push!(array_of_images, left_eyebrow);
    right_eyebrow = image[(1):(h ÷ 20 * 2), (w ÷ 20 * 11):(w ÷ 20 * 19)];
    push!(array_of_images, right_eyebrow);
    left_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20):(w ÷ 20 * 7)];
    push!(array_of_images, left_checkb);
    right_checkb = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 13):(w ÷ 20 * 19)];
    push!(array_of_images, right_checkb);
    nose = image[(h ÷ 20 * 6):(h ÷ 20 * 12), (w ÷ 20 * 7):(w ÷ 20 * 13)];
    push!(array_of_images, nose);
    mouth = image[(h ÷ 20 * 10):(h ÷ 20 * 16), (w ÷ 20 * 4):(w ÷ 20 * 16)];
    push!(array_of_images, mouth);
    mouth_and_chin = image[(h ÷ 20 * 8):(h ÷ 20 * 16), (1):(w)];
    push!(array_of_images, mouth_and_chin);
    # visualizamos los recortes
    # visualize_masc3(img, h, w);

    # ya la devolvemos el formato de Float64
    return imageToColorArray.(array_of_images);
end;


text:
Crops = [:wheat,:corn,:beets]
Purchased = [:wheat,:corn]
Sold = [:wheat,:corn,:beets_quota,:beets_extra]

Cost = Dict(:wheat=>150,:corn=>230,:beets=>260)
Required = Dict(:wheat=>200,:corn=>240,:beets=>0)
PurchasePrice = Dict(:wheat=>238,:corn=>210)
SellPrice = Dict(:wheat=>170,:corn=>150,:beets_quota=>36,:beets_extra=>10)
Budget = 500

@everywhere begin
    struct FarmerScenario <: StochasticPrograms.AbstractScenarioData
        Yield::Dict{Symbol,Float64}
    end
    StochasticPrograms.probability(::FarmerScenario) = 1/3
    function StochasticPrograms.expected(scenarios::Vector{FarmerScenario})
        sd = FarmerScenario(Dict(:wheat=>sum([probability(s)*s.Yield[:wheat] for s in scenarios]),
                                 :corn=>sum([probability(s)*s.Yield[:corn] for s in scenarios]),
                                 :beets=>sum([probability(s)*s.Yield[:beets] for s in scenarios])))
    end
end

s1 = FarmerScenario(Dict(:wheat=>3.0,:corn=>3.6,:beets=>24.0))
s2 = FarmerScenario(Dict(:wheat=>2.5,:corn=>3.0,:beets=>20.0))
s3 = FarmerScenario(Dict(:wheat=>2.0,:corn=>2.4,:beets=>16.0))

sp = StochasticProgram((Crops,Cost,Budget),(Required,PurchasePrice,SellPrice),[s1,s2,s3])

@first_stage sp = begin
    (Crops,Cost,Budget) = stage
    @variable(model, x[c = Crops] >= 0)
    @objective(model, Min, sum(Cost[c]*x[c] for c in Crops))
    @constraint(model, sum(x[c] for c in Crops) <= Budget)
end

@second_stage sp = begin
    @decision x
    (Required,PurchasePrice,SellPrice) = stage
    s = scenario
    @variable(model, y[p = Purchased] >= 0)
    @variable(model, w[s = Sold] >= 0)
    @objective(model, Min, sum( PurchasePrice[p] * y[p] for p = Purchased) - sum( SellPrice[s] * w[s] for s in Sold))

    @constraint(model, const_minreq[p=Purchased],
                   s.Yield[p] * x[p] + y[p] - w[p] >= Required[p])
    @constraint(model, const_minreq_beets,
                   s.Yield[:beets] * x[:beets] - w[:beets_quota] - w[:beets_extra] >= Required[:beets])
    @constraint(model, const_aux, w[:beets_quota] <= 6000)
end

push!(problems,(sp,"Farmer"))


text:
############################################################################################
"""
$(TYPEDEF)

Functor to calculate 'ℓfunc' and gradient at unconstrained 'θᵤ', including eventual Jacobian adjustments.

# Fields
$(TYPEDFIELDS)
"""
struct Objective{M<:ModelWrapper, D, T<:Tagged, F<:AbstractFloat} <: BaytesCore.AbstractObjective
    model::M
    data::D
    tagged::T
    temperature::F
    function Objective(
        model::M,
        data::D,
        tagged::T,
        temperature::F = model.info.reconstruct.default.output(1.0)
    ) where {M<:ModelWrapper,D,T<:Tagged,F<:AbstractFloat}
    ArgCheck.@argcheck 0.0 < temperature <= 1.0 "Temperature has to be bounded between 0.0 and 1.0"
        return new{M,D,T,F}(model, data, tagged, temperature)
    end
end
function Objective(
    model::ModelWrapper{M},
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0)
    ) where {M<:ModelName,D,F<:AbstractFloat}
    return Objective(model, data, Tagged(model), temperature)
end
function Objective(
    model::ModelWrapper{M},
    data::D,
    sym::S,
    temperature::F = model.info.reconstruct.default.output(1.0)
) where {M<:ModelName,D,S<:Union{Symbol,NTuple{k,Symbol} where k},F<:AbstractFloat}
    return Objective(model, data, Tagged(model, sym), temperature)
end

############################################################################################
# Basic functions for Model struct
length(objective::Objective) = length(objective.tagged)
paramnames(objective::Objective) = paramnames(objective.tagged)

############################################################################################
# A bunch of functions that can be used/extended for target model in Sampling process
"""
$(SIGNATURES)
Functor to call target function for Model given parameter and data.

# Examples
```julia
```

"""
function (objective::Objective)(θ)
    return 0.0
end

"""
$(SIGNATURES)
Predict new data given model parameter and data.

# Examples
```julia
```

"""
function predict(_rng::Random.AbstractRNG, objective::Objective)
    return nothing
end

"""
$(SIGNATURES)
Generate statistics given model parameter and data.

# Examples
```julia
```

"""
function generate(_rng::Random.AbstractRNG, objective::Objective)
    return nothing
end
function generate(_rng::Random.AbstractRNG, objective::Objective, gen::BaytesCore.UpdateTrue)
    return generate(_rng, objective)
end
function generate(_rng::Random.AbstractRNG, objective::Objective, gen::BaytesCore.UpdateFalse)
    return nothing
end

"""
$(SIGNATURES)
Assign model dynamics for a given `objective`.

# Examples
```julia
```

"""
function dynamics(objective::Objective)
    return nothing
end

############################################################################################
function (objective::Objective)(θᵤ::AbstractVector{T}) where {T<:Real}
    @unpack model, data, tagged, temperature = objective
    ## Convert vector θᵤ back to constrained space as NamedTuple
    θ = constrain(tagged.info.transform, unflattenAD(tagged.info.reconstruct, θᵤ))
    #!NOTE: There are border cases where θᵤ is still finite, but θ no longer after transformation, so have to cover this separately
    _checkfinite(θ) || return -Inf
    ## logabsdet_jac for transformations
    ℓjac = log_abs_det_jac(tagged.info.transform.unconstrain, θ)
    _checkfinite(ℓjac) || return -Inf
    ## Evaluate objective
    ℓℒ = objective(merge(model.val, θ))
    ## Return log posterior
    return temperature * (ℓℒ + ℓjac)
end

############################################################################################
# Export
export Objective, predict, generate, dynamics


text:
f = Open([1, 2], PetriNet(4, (1,3), (2,4)), [3, 4])

g = Open([1,2], PetriNet(3, ((1,2),3)), [3])

h = f ⋅ g

h′ = Open([1,2], PetriNet(5, (1,3), (2,4), ((3,4), 5)), [5])

h_id = h ⋅ id(OpenPetriNetOb(FinSet(1)))

@test dom(f) == OpenPetriNetOb(FinSet(2))
@test codom(f) == OpenPetriNetOb(FinSet(2))
@test dom(g) == codom(f)
@test codom(g) == OpenPetriNetOb(FinSet(1))
@test dom(h) == dom(f)
@test codom(h) == codom(g)

@test h == h′
@test h == h_id

# Test open petri net notations either fully open or by specifying legs

pn = Open(PetriNet(3, ((1,2), 3)), [1], [2], [3])
pn′ = Open(PetriNet(3, ((1,2), 3)))
@test pn == pn′

lpn = Open(LabelledPetriNet([:I, :R], (:rec, :I=>:R)))
lpn′ = Open([:I], LabelledPetriNet([:I, :R], (:rec, :I=>:R)), [:R])
@test lpn == lpn′

rn = Open(ReactionNet{Number,Int}([10, 0], (.25, 1=>2)))
rn′ = Open([1], ReactionNet{Number,Int}([10, 0], (.25, 1=>2)), [2])
@test rn == rn′

lrn = Open(LabelledReactionNet{Number,Int}([:I=>10, :R=>0], ((:rec=>.25), :I=>:R)))
lrn′ = Open([:I], LabelledReactionNet{Number,Int}([:I=>10, :R=>0], ((:rec=>.25), :I=>:R)), [:R])
@test lrn == lrn′


text:
@testset "StructTabular" begin
    @testset "StructSeries" begin
        @test @inferred(Series(2+3im)) isa Series{Complex{Int}}
        # TODO @test_throws ErrorException Series(1)

        s = Series(2+3im)
        @test @inferred(s[l"re"]) === 2
        @test @inferred(s[l"im"]) === 3
        @test_throws Exception @inferred(s[l"bla"])

        @test @inferred(s[:]) == s
        @test @inferred(s[(l"re", l"im")]) === Series(l"re" => 2, l"im" => 3)
    end

    # Also covered in ArrayTabular tests...
    @testset "Tables with StructSeries as rows" begin
        @test @inferred(Table([1+2im, 3+4im])) isa Table{Vector{Complex{Int}}}

        t = Table([1+2im, 3+4im])

        @test @inferred(t[l"re", 1]) === 1
        @test @inferred(t[l"re", 2]) === 3
        @test @inferred(t[l"im", 1]) === 2
        @test @inferred(t[l"im", 2]) === 4

        @test @inferred(t[:, 1]) == Series(1+2im)
        @test @inferred(t[:, 2]) == Series(3+4im)
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t[l"re", :] == Series([1,3])
            #@test_broken @inferred(t[l"re", :]) == Series([1,3])
        else
            @test @inferred(t[l"re", :]) == Series([1,3])
        end
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t[l"im", :] == Series([2,4])
            #@test_broken @inferred(t[l"im", :]) == Series([2,4])
        else
            @test @inferred(t[l"im", :]) == Series([2,4])
        end

        @test @inferred(t[(l"re", l"im"), 1]) === Series(l"re" => 1, l"im" => 2)
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t[(l"re", l"im"), :] == t
            #@test_broken @inferred(t[(l"re", l"im"), :]) == t
        else
            @test @inferred(t[(l"re", l"im"), :]) == t
        end

        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t[:, :] == t
            #@test_broken @inferred(t[:, :]) == t
        else
            @test @inferred(t[:, :]) == t
        end
    end

    @testset "StructTables" begin
        # structs of structs
        @eval struct Foo{A,B}
            a::A
            b::B
        end

        data = Foo(1+2im, 3.0+4.0im)
        @test @inferred(Table(data)) isa Table{<:Foo}
        t = Table(data)

        @test @inferred(t[l"re", l"a"]) === 1
        @test @inferred(t[l"im", l"a"]) === 2
        @test @inferred(t[l"re", l"b"]) === 3.0
        @test @inferred(t[l"im", l"b"]) === 4.0

        @test @inferred(t[:, l"b"]) == Series(3.0+4.0im)
        @test @inferred(t[l"re", :]) == Series(l"a" => 1, l"b" => 3.0)
        @test @inferred(t[:, :]) == t
    end
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# wrap.jl: automatically wrap C API using Clang.jl

# This wrapper uses the Gnuic/refactor branch of Clang.jl to provide proper enumeration
# support. Implementation guided by:
#   https://github.com/JuliaDiffEq/Sundials.jl (general)
#   https://github.com/JuliaAttic/CUDArt.jl/ (error rewriter)

using Clang: init

const incpath = normpath("/usr/include/spinnaker/spinc")

if !isdir(incpath)
  error("Spinnaker SDK C header not found, please install manually.")
end

const outpath = normpath(@__DIR__, "..", "src", "wrapper")
mkpath(outpath)

const spinnaker_headers = readdir(incpath)

const clang_path = "/usr/lib/clang/6.0"
const clang_includes = [joinpath(clang_path, "include")]

# TODO: should we use actually use Clang.jl installation here?
# const CLANG_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include", "clang-c") |> normpath
# const CLANG_HEADERS = [joinpath(CLANG_INCLUDE, cHeader) for cHeader in readdir(CLANG_INCLUDE) if endswith(cHeader, ".h")]

headerlibmap = Dict("SpinnakerC.h"       => "libSpinnaker_C",
                    "QuickSpinC.h"       => "libSpinnaker_C",
                    "SpinnakerGenApiC.h" => "libSpinnaker_C",
                    "SpinVideoC.h"       => "libSpinVideo_C")

function library_file(header::AbstractString)
  header_name = basename(header)
  if(haskey(headerlibmap, header_name))
    return headerlibmap[header_name]
  else
    @warn "$header has unknown library association, using libSpinnaker_C"
    return "libSpinnaker_C"
  end
end

# Test if header should be wrapped
function header_filter(top_hdr::AbstractString, cursor_header::AbstractString)
  return (top_hdr == cursor_header) # Do not wrap nested includes
end

# Test if element should be wrapped
function cursor_filter(name::AbstractString, cursor)
  return true
end
#
# if typeof(cursor) == Clang.LibClang.CXCursor
#     # Wrap only spin, and quickSpin prefaced functions
#     return occursin(r"^(quickSpin|spin)", name)
# else
#     # Else, wrap everything
#     return true
# end

# Add error checking to generated functions
rewriter(arg) = arg
# rewriter(A::Array) = [rewriter(a) for a in A]
# rewriter(s::Symbol) = string(s)

const skip_expr = []
const skip_error_check = []
function rewriter(ex::Expr)

    # Skip sepcified expressions
    if in(ex, skip_expr)
        return :()
    end

    # Only process function calls
    ex.head == :function || return ex

    decl, body = ex.args[1], ex.args[2]
    ccallexpr = body.args[1]
    rettype = ccallexpr.args[3]

    if rettype == :spinError
        fname = decl.args[1]
        if !in(fname, skip_error_check)
            body.args[1] = Expr(:call, :checkerror, deepcopy(ccallexpr))
        end
    end

    return ex
end

# Build wrapping context
const context = init(
  headers = map(x -> joinpath(incpath, x), spinnaker_headers),
  output_file = joinpath(outpath, "spin_api.jl"),
  common_file = joinpath(outpath, "spin_common.jl"),
  clang_args = [
  "-D", "__STDC_LIMIT_MACROS",
  "-D", "__STDC_CONSTANT_MACROS",
  "-v"
  ],
  clang_diagnostics = true,
  clang_includes = [clang_includes; incpath],
  header_library = library_file,
  header_wrapped = header_filter,
  cursor_wrapped = cursor_filter,
  rewriter = rewriter
)

run(context)

# Manual changes
@warn "Manual changes to wrapper are required, see wrap.jl"
# spinStringGetValue Cstring -> Ptr{UInt8}
# spinEnumerationEntryGetSymbolic Cstrong -> Ptr{UInt8}

# Remove unused files
rm(joinpath(outpath, "LibTemplate.jl"))
rm(joinpath(outpath, "ctypes.jl"))

@info "Warpper written to $outpath"


text:
# # Portfolio Optimisation
#
# We consider a single-period Markowitz portfolio optimisation example.
# Assume that we have a portfolio with $n$ assets at the beginning of time period $t$. Given some forecasts on risks and expected returns we try to find the optimal trade vector that rebalances the portfolio to achieve a good balance between expected risk (variance) $x^\top \Sigma x$ and returns $\mu^\top x$. In it's most simple form we want to solve:
# $$
# \begin{array}{ll} \text{maximize} &  \mu^\top x - \gamma (x^\top \Sigma x)\\
# \text{subject to} &  1^\top x = d + 1^\top x^0 \\
#                    &  x  \geq 0,
# \end{array}
# $$
# with variable $x \in \mathbf{R}^n$, $\mu$ forecasted (expected) returns, $\gamma > 0 $ risk aversion parameter. $x^0_i$ represents the initial investment in asset $i$ and $d$ represents the cash reserve. Consequently, the equality constraint tells us that the sum of the new allocation vector $x$ has to equal the initial allocation plus the cash reserve. Furthermore, the covariance matrix of our risk model is given by $\Sigma \in \mathbf{S}_+^n$. Here we assume a factor risk model, i.e. we can write $\Sigma$ as $\Sigma = D + F F^\top$ where $D$ is diagonal and the factor matrix $F$ has a lower rank $k < n$. This approach allows us to reduce the number of nonzeros in the problem. Furthermore, note that we don't consider shortselling in this example. Let's generate some problem data:
#-
using LinearAlgebra, SparseArrays, Random, COSMO, JuMP, Test

## generate the data
rng = Random.MersenneTwister(1)
k = 5; # number of factors
n = k * 10; # number of assets
D = spdiagm(0 => rand(rng, n) .* sqrt(k))
F = sprandn(rng, n, k, 0.5); # factor loading matrix
μ = (3 .+ 9. * rand(rng, n)) / 100. # expected returns between 3% - 12%
γ = 1.0; # risk aversion parameter
d = 1 # we are starting from all cash
x0 = zeros(n);

# We can now write the problem as a QP:
# $$
# \begin{array}{ll} \text{minimize} &  x^\top D x + y^\top y - \gamma^{-1} \mu^\top x \\
# \text{subject to} &  y = F^\top x \\
#                   &  1^\top x = d + 1^\top x^0 \\
#                    &  x  \geq 0.
# \end{array}
# $$
# Before considering other effects, let's create the model in JuMP and solve it using COSMO:
#-

model = JuMP.Model(with_optimizer(COSMO.Optimizer));
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@objective(model, Min,  x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);
JuMP.optimize!(model)

# After solving the problem, we can calculate the expected return and risk $\sigma= \sqrt{x^{* \top} \Sigma x^*}$:
x_opt = JuMP.value.(x);
y_opt = JuMP.value.(y);
expected_return_basic = dot(μ, x_opt)

#-
expected_risk_basic = sqrt(dot(y_opt, y_opt))

# ## Using standard deviation in the model
# It is pointed out in \[1\] that above problem formulation can lead to numerical problems, e.g. if $\Sigma$ is not strictly positive semidefinite. Another option is to formulate the risk constraint in terms of the standard deviation $\|M^\top x \|$ where $M M^\top = D + F F^\top$ and bound it using a second-order cone constraint:
# $$
# \begin{array}{ll} \text{minimize} &  - \mu^\top x \\
# \text{subject to} &    \|M^\top x\| \leq \gamma \\
#                   &  1^\top x = d + 1^\top x^0 \\
#                    &  x  \geq 0.
# \end{array}
# $$

Mt = [D.^0.5; F']
model = JuMP.Model(with_optimizer(COSMO.Optimizer));
@variable(model, x[1:n]);
@objective(model, Min, - μ' * x);
@constraint(model,  [γ; Mt * x] in SecondOrderCone()); # ||M'x|| <= γ
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);
JuMP.optimize!(model)

# Note that the result is different from the example above because $\gamma$ scales the problem in a different way. Here it can be seen as an upper bound on the standard deviation of the portfolio.
x_opt = JuMP.value.(x);
expected_return = dot(μ, x_opt)
# Let us verify that the bound holds:
@test norm(Mt * x_opt) <= γ

# ## Pareto-optimal front
# The above portfolio optimisation approach yields the optimal expected return for a given level of risk. The result is obviously impacted by the risk aversion $\gamma$ parameter. To visualise the trade-off and present the investor with an efficient Pareto optimal portfolio for their risk appetite we can compute the optimal portfolio for many choices of $\gamma$ and plot the corresponding risk-return trade-off curve.

gammas = [ 0.001, 0.01, 0.1,  0.5,  1., 3., 10, 100, 1000]
risks = zeros(length(gammas))
returns = zeros(length(gammas))
model = JuMP.Model(with_optimizer(COSMO.Optimizer, verbose = false));
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@objective(model, Min,  x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);

## solve the same problem for different values of γ
for (k, gamma) in enumerate(gammas)
    coeff = - 1/gamma * μ
    JuMP.set_objective_coefficient.(model, x, coeff)
    JuMP.optimize!(model)
    x_opt = JuMP.value.(x);
    y_opt = JuMP.value.(y);
    returns[k] = dot(μ, x_opt)
    risks[k] = sqrt(dot(y_opt, y_opt))
end
# We can now plot the risk-return trade-off curve:
using Plots
Plots.plot(risks, returns, xlabel = "Standard deviation (risk)", ylabel = "Expected return", title = "Risk-return trade-off for efficient portfolios", legend = false)

# ## Transaction costs
# In the model above we assume that trading the assets is free and does not impact the market. However, this is clearly not the case in reality. To make the example more realistic consider the following cost $c_j$ associated with the trade $δ_j = x_j - x_j^0$:
# $$
# c_j(\delta_j) = a_j |\delta_j| + b_j |\delta_j|^{3/2},
# $$
# where the first term models the bid-ask spread and broker fees for asset $j$. The second term models the impact on the market that our trade has. This is obviously only a factor if the volume of our trade is significant. The constant $b_j$ is a function of the total volume traded in the considered time periode and the price volatility of the asset and has to be estimated by the trader. To make this example simple we consider the same coefficients $a$ and $b$ for every asset. The $|\delta_j|^{3/2}$ term can be easily modeled using a power cone constraint $\mathcal{K}_{pow} = \{(x, y, z) \mid x^\alpha y^{(1-\alpha)} \geq |z|, x \geq 0, y \geq 0, 0 \leq \alpha \leq 1 \}$. In fact this can be used to model any market impact function with exponent greater than 1.
# We can write the total transaction cost $a^\top s + b^\top t$ where $s_j$ bounds the absolute value of $\delta_j$ and $t_{j}$ is used to bound the term $|x_j - x_j^0|^{3/2} \leq t_{j}$ using a power cone formulation: $(t_{j}, 1, x_j - x_j^0) \in \mathcal{K}_{pow}(2/3)$.

#-
a = 1e-3
b = 1e-1
γ = 1.0;
model = JuMP.Model(with_optimizer(COSMO.Optimizer, eps_abs = 1e-5, eps_rel = 1e-5));
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@variable(model, t[1:n]);
@variable(model, s[1:n]);
@objective(model, Min, x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, x .>= 0);

## transaction costs
@constraint(model, sum(x) + a * sum(s) + b * sum(t) == d + sum(x0) );
@constraint(model, [i = 1:n], x[i] - x0[i] <= s[i]); # model the absolute value with slack variable s
@constraint(model, [i = 1:n], x0[i] - x[i] <= s[i]);
@constraint(model, [i = 1:n], [t[i], 1, x[i] - x0[i]] in MOI.PowerCone(2/3));
JuMP.optimize!(model)
# Let's look at the expected return and the total transaction cost:

x_opt = JuMP.value.(x);
y_opt = JuMP.value.(y);
s_opt = JuMP.value.(s);
t_opt = JuMP.value.(t);
expected_return = dot(μ, x_opt)
#-
expected_risk = dot(y_opt, y_opt)
#-
transaction_cost = a * sum(s_opt) + b * sum( t_opt)

# ## References
# [1] [Mosek Case Studies](https://docs.mosek.com/9.2/pythonfusion/case-studies-portfolio.html)


text:
module GenderInference

export
    gendercount,
    proportionfemale,
    proportionmale,
    gender

using DataDeps
using BSON
using SparseArrays
using Logging

import Base: getindex, haskey

include("data_registrations.jl")
include("data_management.jl")
include("inference.jl")

function __init__()
    init_us_census()

    # Define global constant
    @eval const NAMES = parsedataset(datadep"US Census - names", RawDataSet("USCensus"))
end

end # module Gender Inference


text:
using Base: @_noinline_meta, @_inline_meta
using Core: is_top_bit_set
using Core.Intrinsics: bitcast, trunc_int, sext_int, zext_int, sle_int, eq_int, and_int

# onehot erorr
struct OneHotEncodeError <: Exception
  K
  val
  OneHotEncodeError(@nospecialize(K), @nospecialize(val)) = (@_noinline_meta; new(K, val))
end

function Base.showerror(io::IO, e::OneHotEncodeError)
  print(io, "OneHotEncodeError: cannot encode ")
  print(io, e.val)
  print(io, " with OneHot{")
  print(io, e.K)
  print(io, '}')
end

throw_onehotencode_error(K, val) = (@_noinline_meta; throw(OneHotEncodeError(K, val)))

# onehot encode
primitive type OneHot{K} <: AbstractVector{Bool} 32 end

OneHot(k) = OneHot{UInt32(k)}
OneHot{K}(x) where K = convert(OneHot(K), x)
OneHot(k, x) = OneHot{k}(x)

onehotsize(::OneHot{K}) where K = Int(K)

# array interface

Base.size(o::OneHot) = (onehotsize(o),)
function Base.getindex(o::OneHot{K}, i::I) where {K, I<:Integer}
  @boundscheck checkbounds(o, i)
  return convert(I, o) == i
end

# printing

function Base.showarg(io::IO, x::OneHot{K}, toplevel) where K
  toplevel || print(io, "::")
  join(io, ["OneHot{", Int(K), '}'])
end

# convert

Base.UInt32(o::OneHot) = bitcast(UInt32, o)
Base.UInt64(o::OneHot) = zext_int(UInt64, o)
Base.Int32(o::OneHot) = bitcast(Int32, o)
Base.Int64(o::OneHot) = zext_int(Int64, o)

Base.convert(::Type{Any}, o::OneHot) = o
Base.convert(::Type{OneHot{K}}, o::OneHot{K}) where {K} = o
Base.convert(::Type{UInt32},  o::OneHot) = UInt32(o)
Base.convert(::Type{To}, o::OneHot) where {To} = convert(To, convert(UInt32, o))

Base.convert(ot::Type{OneHot{K}}, x::Core.BuiltinInts) where {K} = toOneHot(ot, x)

# zero

Base.zero(o::O) where {O<:OneHot} = toOneHot(O, 0x00000000)
Base.iszero(o::O) where {O<:OneHot} = iszero(convert(UInt32, o))

# bit-op

function check_onehot_top_bit(::Type{OneHot{K}}, x) where {K}
  @_inline_meta
  is_top_bit_set(x) && throw_onehotencode_error(K, x)
  x
end

function check_onehot_encode(ot::Type{OneHot{K}}, x) where {K}
  @_inline_meta
  sle_int(x, K) || throw_onehotencode_error(K, x)
  bitcast(ot, x)
end

function checked_onehot_trunc_sint(ot::Type{OneHot{K}}, x::From) where {K, From}
  @_inline_meta
  y = trunc_int(UInt32, x)
  back = sext_int(From, y)
  eq_int(x, back) || throw_onehotencode_error(K, x)
  check_onehot_encode(ot, y)
end

function checked_onehot_trunc_uint(ot::Type{OneHot{K}}, x::From) where {K, From}
  @_inline_meta
  y = trunc_int(UInt32, x)
  back = zext_int(From, y)
  eq_int(x, back) || throw_onehotencode_error(K, x)
  check_onehot_encode(ot, y)
end

toOneHot(ot::Type{OneHot{K}}, x::Int8) where {K} = check_onehot_encode(ot, sext_int(UInt32, check_onehot_top_bit(ot, x)))
toOneHot(ot::Type{OneHot{K}}, x::Int16) where {K} = check_onehot_encode(ot, sext_int(UInt32, check_onehot_top_bit(ot, x)))
toOneHot(ot::Type{OneHot{K}}, x::Int32) where {K} = check_onehot_encode(ot, bitcast(UInt32, check_onehot_top_bit(ot, x)))
toOneHot(ot::Type{OneHot{K}}, x::Int64) where {K} = checked_onehot_trunc_sint(ot, check_onehot_top_bit(ot, x))
toOneHot(ot::Type{OneHot{K}}, x::Int128) where {K} = checked_onehot_trunc_sint(ot, check_onehot_top_bit(ot, x))
toOneHot(ot::Type{OneHot{K}}, x::UInt8) where {K} = check_onehot_encode(ot, zext_int(UInt32, x))
toOneHot(ot::Type{OneHot{K}}, x::UInt16) where {K} = check_onehot_encode(ot, zext_int(UInt32, x))
toOneHot(ot::Type{OneHot{K}}, x::UInt32) where {K} = check_onehot_encode(ot, x)
toOneHot(ot::Type{OneHot{K}}, x::UInt64) where {K} = checked_onehot_trunc_uint(ot, x)
toOneHot(ot::Type{OneHot{K}}, x::UInt128) where {K} = checked_onehot_trunc_uint(ot, x)
toOneHot(ot::Type{OneHot{K}}, x::Bool) where {K} = and_int(zext_int(ot, x), toOneHot(ot, 0x1))

# onehot array
struct OneHotArray{K, N, var"N+1", A<:AbstractArray{OneHot{K}, N}} <: AbstractArray{Bool, var"N+1"}
  onehots::A
end

OneHotArray(onehots::A) where {K, A<:AbstractArray{OneHot{K}}} = OneHotArray{K, ndims(onehots), ndims(onehots)+1, A}(onehots)
OneHotArray{K}(indices::A) where {K, A<:AbstractArray{<:Integer}} = OneHotArray(K, indices)
OneHotArray(k, xs) = OneHotArray(OneHot{k}.(xs))
OneHotArray(k, o::OneHot) = OneHotArray([o])

onehotsize(::OneHotArray{K}) where K = Int(K)

# array interface
Base.size(oa::OneHotArray{K}) where K = (onehotsize(oa), size(oa.onehots)...)

function Base.getindex(oa::OneHotArray{K, N}, i, is::Vararg{Int, N}) where {K, N}
  @boundscheck checkbounds(oa, i, is...)
  oa.onehots[is...][i]
end

function Base.getindex(oa::OneHotArray{K, N}, i::Colon, is::Vararg{Int, N}) where {K, N}
  @boundscheck checkbounds(oa, i, is...)
  oa.onehots[is...]
end

function Base.getindex(oa::OneHotArray{K}, i::Colon, is...) where {K}
  @boundscheck checkbounds(oa, i, is...)
  OneHotArray(oa.onehots[is...])
end

Base.similar(o::OneHotArray, ::Type{T}, dims::Dims{N}) where {T, N} = similar(o.onehots, T, dims)

# printing

function Base.summary(io::IO, oa::OneHotArray)
  join(io, size(oa), 'x')
  join(io, [" OneHotArray{", onehotsize(oa), ", ", ndims(oa), ", "])
  Base.showarg(io, oa.onehots, true)
  print(io, "}")
end

# cat

Base.vcat(xss::OneHot{K}...) where K = cat(xss...; dims=Val(1))
Base.hcat(xss::OneHot{K}...) where K = cat(xss...; dims=Val(2))

function Base.cat(xss::OneHot{K}...; dims) where K
  isone(::Val{V}) where V = isone(V)
  isone(v) = Base.isone(v)
  if isone(dims)
    @warn "concat OneHot{$K} along dimension 1."
    Base._cat(Val(1), xss...)
  else
    predecessor(::Val{V}) where V = Val(V-1)
    predecessor(v) = Val(v - 1)
    yss = reshape(collect(xss), reverse(Base.rdims(predecessor(dims), axes(xss))))
    OneHotArray(yss)
  end
end

Base.vcat(xss::OneHotArray{K}...) where K = cat(xss...; dims=Val(1))
Base.hcat(xss::OneHotArray{K}...) where K = cat(xss...; dims=Val(2))

function Base.cat(xss::OneHotArray{K}...; dims) where K
  isone(::Val{V}) where V = isone(V)
  isone(v) = Base.isone(v)
  if isone(dims)
    @warn "concat OneHotArray{$K} along dimension 1."
    Base._cat(Val(1), xss...)
  else
    predecessor(::Val{V}) where V = Val(V-1)
    predecessor(v) = v - 1
    sdims = predecessor(dims)
    xidss = map(xs->xs.onehots, xss)
    ret = cat(xidss...; dims=sdims)
    OneHotArray(ret)
  end
end

# reshape

"reshape the onehots"
function ohreshape(parent::OneHotArray{K}, dims) where K
  onehots = parent.onehots
  OneHotArray(reshape(onehots, dims))
end

function Base.reshape(parent::OneHotArray{K}, dims::Dims) where K
  isequal(prod(dims), length(parent)) || throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(parent))"))
  return isequal(K, first(dims)) ?
    ohreshape(parent, Base.tail(dims)) :
    Base._reshape(parent, dims)
end

function Base.reshape(parent::OneHotArray{K}, dims::Tuple{Vararg{Union{Colon, Int64}}}) where K
  rdims = Base._reshape_uncolon(parent, dims)
  return isequal(K, first(rdims)) ?
    ohreshape(parent, Base.tail(rdims)) :
    Base._reshape(parent,
                  rdims)
end

# old utility

onehot2indices(xs::OneHotArray) = convert(AbstractArray{Int}, xs.onehots)

indices2onehot(nums::Int, xs::AbstractArray{Int}) = OneHotArray(nums, xs)

tofloat(::Type{F}, o::OneHotArray{N, <:AbstractArray}) where {F<:AbstractFloat,N} = Array{F}(o)

# gpu
import Adapt: adapt, adapt_structure
adapt_structure(T, oa::OneHotArray) = OneHotArray(adapt(T, oa.onehots))

# AD
using Flux
Flux.@nograd OneHotArray

using ZygoteRules: @adjoint, AContext, pullback
import ZygoteRules: _pullback

import Base: *
*(A::AbstractMatrix{T}, b::OneHotArray{N}) where {T,N} = gather(A, b)

_pullback(cx::AContext, ::typeof(*), A::AbstractMatrix{T}, b::OneHotArray{N}) where {T, N} = _pullback(cx, gather, A, b)


text:
X = rand(100, 20)

# True generating function has intercept, 2 factors and a lag range of 3-5.
genfn(r, X) = 1.56 + 2.4*X[r-5, 3] - 5.6*X[r-3, 13]

y = Float64[genfn(r, X) for r in 6:size(X, 1)]

# Now in reality we don't know the real model se we will have to train it
# from data. Let's do it and try to recover an interpretable model.
using InterpretableModels

# We don't know how many factors so let's aim for 3 and a lag range of 2-10
lm = LaggedFactorLITRM(3, size(X, 2); minlag=2, maxlag=10)

fit!(lm, X, y; MaxTime = 10.0)

println(lm) # The model should have been recovered quite well...


text:
const MAX_WORD_LENGTH = 50

function read_word(stream)
  buffer = IOBuffer(MAX_WORD_LENGTH)
  while !eof(stream)
    ch = read(stream, Uint8)
    if ch == 0x20
      break
    elseif ch != 0x0A
      write(buffer, ch)
    end
  end
  return takebuf_string(buffer)
end

function read_vector(stream, dim)
  res = Array(Float64, dim)
  for i = 1:dim
    res[i] = float64(read(stream, Float32))
  end
  return res / norm(res)
end


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using DifferentialEquations
using Plots

include("SI\\metapopulation.jl")

N = 10
P = 1000
volume = 1
frac_infected = 0.01
ninfected = Int(frac_infected*P)
g = complete_graph(N)
g = random_configuration_model(N,fill(3,N))
g = path_graph(N)
V = fill(volume, N)
β = 0.3
μ = 0.1
tmax = 1.0

mp = Metapopulation(g,V,β,μ)
s0 = fill(P, N)
i0 = fill(ninfected, N)
i0 = zeros(Int, N)
s0[1] = P-ninfected
i0[1] = ninfected
Ptot = sum(s0) + sum(i0)

t, s, i = metapopulation_gillespie(mp,s0,i0, nmax = 1000000, tmax=tmax)

plot(t, sum(i, dims=2)/Ptot,
    label="",
    line=:steppre,
    xlabel="Time",
    ylabel="Fraction of infected individuals"
    )

plot(t,i/P,
    line=:steppre,
    label="",
    xlabel = "Time",
    ylabel = "Normalized #infected individuals"
    )

plot(t,s/Ptot, line=:steppre)

t_mc, s_mc, i_mc = metapopulation_montecarlo(mp,s0,i0,tmax=tmax, nmax=150000, nbins=200)

t_mf, s_mf, i_mf, sol = metapopulation_ode(mp, s0, i0, tmax=tmax, saveat=t_mc)

plot(t_mf, i_mf,
    label="",
    xlabel="Time",
    ylabel="#infected individuals",
    linestyle=:dash
    )
plot!(t, i, line=:steppre, label="")

plot(t_mc, i_mc,
    label="",
    xlabel="Time",
    ylabel="#infected individuals"
    )
plot!(t_mf, i_mf, label="", linestyle=:dash)

plot(t_mc, i_mf-i_mc, label="")

plot(t_mc, sum(i_mc, dims=2)/Ptot,
    label="Average",
    xlabel="Time",
    ylabel="Fraction of infected individuals",
    legend=:bottomright
    )
plot!(t_mf, sum(i_mf, dims=2)/Ptot,
    label="Mean Field")

function transient_time_montecarlo(mp::Metapopulation, s0, i0; tmax=100.0, nmax=1000, nsims = 100)
    totalpop = sum(s0) + sum(i0)
    tts = zeros(nsims)
    for l in 1:nsims
        t, s, i = metapopulation_gillespie(mp, s0, i0, tmax=tmax, nmax=nmax)
        k = 1
        tot = sum(i[k,:])
        while tot < 0.9*totalpop && k < length(t)
            k += 1
            tot += sum(i[k,:])
        end
        tts[l] = t[k]
    end
    return mean(tts)
end

μs = LinRange(0.0,1.0, 200)
ttime = zeros(200)
for k in 1:200
    mp = Metapopulation(g, β, μs[k])
    ttime[k] = transient_time_montecarlo(mp,s0,i0,tmax = 1.0, nmax=5000, nsims=200)
end
scatter(μs, ttime, label="", xlabel="\\mu", ylabel="\\tau")
savefig("transientvsmu.png")

βs = LinRange(0.0,1.0,200)

ttimes = zeros(200,200)

for k in 1:200, l in 1:200
    mp = Metapopulation(g, βs[k], μs[l])
    ttimes[k,l] = transient_time_montecarlo(mp,s0,i0,tmax=1.0,nmax=5000,nsims=200)
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
function furnsh(KernelFile::String)
  ccall((:furnsh_c, sharedLib),Void,(Ptr{Cchar},),KernelFile)
end

function pxform(from, to, et)
  rotMat = zeros(Float64, 3,3)
  ccall((:pxform_c, sharedLib), Void, (Ptr{Cchar}, Ptr{Cchar}, Cdouble, Ptr{Cdouble}),
        from, to, et, rotMat)
  return rotMat
end

function reclat(a::Vector)
  r = Ref{Cdouble}(0.0)
  lon = Ref{Cdouble}(0.0)
  lat = Ref{Cdouble}(0.0)
  ccall((:reclat_c, sharedLib), Void, (Ptr{Cdouble}, Ptr{Cdouble},
                     Ptr{Cdouble}, Ptr{Cdouble}), a, r, lon, lat)
  return r[], lon[], lat[]
end

function str2et(s::AbstractString)
  et = Ref{Cdouble}(0.0)
  ccall((:str2et_c, sharedLib), Float64, (Ptr{Cchar}, Ptr{Cdouble}), s, et)
  return et[]
end

function utc2et(s::AbstractString)
  et = Ref{Cdouble}(0.0)
  ccall((:utc2et_c, sharedLib), Float64, (Ptr{Cchar}, Ptr{Cdouble}), s, et)
  return et[]
end

function scs2e(sc::Int64, sclkch::AbstractString)
  et = Ref{Cdouble}(0.0)
  ccall((:scs2e_c, sharedLib), Void, (Int64, Ptr{Cchar}, Ptr{Cdouble}),
        sc, sclkch, et)
  return et[]
end

function spkpos(target::AbstractString, et::Float64, ref::AbstractString,
                abcorr::AbstractString, obs::AbstractString)
  pos = zeros(Float64, 3)
  lt = Ref{Cdouble}(0.0)
  ccall((:spkpos_c, sharedLib), Void, (Ptr{Cchar}, Float64,
        Ptr{Cchar}, Ptr{Cchar}, Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}),
        target, et, ref, abcorr, obs, pos, lt)
  return pos, lt[]
end

function timout(et::Float64, pictur::AbstractString, lenOutput::Int64)
  s = "............................."
  ccall((:timout_c, sharedLib), Void, (Float64, Ptr{Cchar}, Int64, Ptr{Cchar}),
        et, pictur, lenOutput, s)
  return s
end

function unload(KernelFile::AbstractString)
  ccall((:furnsh_c, sharedLib),Void,(Ptr{Cchar},),KernelFile)
end

function vsep(a::Vector, b::Vector)
  ccall((:vsep_c, sharedLib), Float64, (Ptr{Float64}, Ptr{Float64}), a, b)
end


text:
using BenchmarkTools

# Get the number of data points N from the command line.
N = parse(Int, ARGS[1])

println("-----------------------------------------")
println("Compute FFTs: ", N)
println("-----------------------------------------")
println(" ")

#Take the FFT of the N random data points.
@btime begin
    mat = complex.(rand(N, N), randn(N, N))
    result = fft(mat)
    result = abs.(result)
end