text:
"""
```
decompose_forecast(m_new, m_old, df_new, df_old, input_type, cond_new, cond_old,
    classes; verbose = :low, kwargs...)

decompose_forecast(m_new, m_old, df_new, df_old, params_new, params_old,
    cond_new, cond_old, classes; check = false)
```
explains the differences between an old forecast and a new forecast
by decomposing the differences into three sources:

(1) Data revisions,
(2) News (e.g. new data that has become available since the old forecast),
(3) Re-estimation (i.e. changes in model parameters).

### Inputs

- `m_new::M` and `m_old::M` where `M<:AbstractDSGEModel`
- `df_new::DataFrame` and `df_old::DataFrame`
- `cond_new::Symbol` and `cond_old::Symbol`
- `classes::Vector{Symbol}`: some subset of `[:states, :obs, :pseudo]`

**Method 1 only:**

- `input_type::Symbol`: estimation type to use. Parameters will be loaded using
  `load_draws(m_new, input_type)` and `load_draws(m_old, input_type)` in this
  method

**Method 2 only:**

- `params_new::Vector{Float64}` and `params_old::Vector{Float64}`: single
  parameter draws to use

### Keyword Arguments

- `check::Bool`: whether to check that the individual components add up to the
  correct total difference in forecasts. This roughly doubles the runtime

**Method 1 only:**

- `verbose::Symbol`

### Outputs

The first method returns nothing. The second method returns
`decomp::Dict{Symbol, Matrix{Float64}}`, which has keys of the form
`:decomp<component><class>` and values of size `Ny` x `Nh`, where

- `Ny` is the number of variables in the given class
- `Nh` is the number of common forecast periods, i.e. periods between
  `date_forecast_start(m_new)` and `date_forecast_end(m_old)`
"""
function decompose_forecast(m_new::M, m_old::M, df_new::DataFrame, df_old::DataFrame,
                            input_type::Symbol, cond_new::Symbol, cond_old::Symbol,
                            classes::Vector{Symbol};
                            verbose::Symbol = :low, forecast_string_new = "", forecast_string_old = "", kwargs...) where M<:AbstractDSGEModel
    # Get output file names
    decomp_output_files = get_decomp_output_files(m_new, m_old, input_type, cond_new, cond_old, classes, forecast_string_old = forecast_string_old, forecast_string_new = forecast_string_new)

    info_print(verbose, :low, "Decomposing forecast...")
    println(verbose, :low, "Start time: $(now())")

    # Set up call to lower-level method
    f(params_new::Vector{Float64}, params_old::Vector{Float64}) =
      decompose_forecast(m_new, m_old, df_new, df_old, params_new, params_old,
                         cond_new, cond_old, classes; kwargs...)

    # Single-draw forecasts
    if input_type in [:mode, :mean, :init]

        params_new = load_draws(m_new, input_type, verbose = verbose)
        params_old = load_draws(m_old, input_type, verbose = verbose)
        decomps = f(params_new, params_old)
        write_forecast_decomposition(m_new, m_old, input_type, classes, decomp_output_files, decomps,
                                     forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old,
                                     verbose = verbose)

    # Multiple-draw forecasts
    elseif input_type == :full

        block_inds, block_inds_thin = forecast_block_inds(m_new, input_type)
        nblocks = length(block_inds)
        total_forecast_time = 0.0

        for block = 1:nblocks
            println(verbose, :low)
            info_print(verbose, :low, "Decomposing block $block of $nblocks...")
            begin_time = time_ns()

            # Get to work!
            params_new = load_draws(m_new, input_type, block_inds[block], verbose = verbose)
            params_old = load_draws(m_old, input_type, block_inds[block], verbose = verbose)
            mapfcn = use_parallel_workers(m_new) ? pmap : map
            decomps = mapfcn(f, params_new, params_old)

            # Assemble outputs from this block and write to file
            decomps = convert(Vector{Dict{Symbol, Array{Float64}}}, decomps)
            decomps = assemble_block_outputs(decomps)
            write_forecast_decomposition(m_new, m_old, input_type, classes, decomp_output_files, decomps,
                                         block_number = Nullable(block), block_inds = block_inds_thin[block],
                                         forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old,
                                         verbose = verbose)
            GC.gc()

            # Calculate time to complete this block, average block time, and
            # expected time to completion
            block_time = (time_ns() - begin_time)/1e9
            total_forecast_time += block_time
            total_forecast_time_min     = total_forecast_time/60
            expected_time_remaining     = (total_forecast_time/block)*(nblocks - block)
            expected_time_remaining_min = expected_time_remaining/60

            println(verbose, :low, "\nCompleted $block of $nblocks blocks.")
            println(verbose, :low, "Total time elapsed: $total_forecast_time_min minutes")
            println(verbose, :low, "Expected time remaining: $expected_time_remaining_min minutes")
        end # of loop through blocks

    else
        error("Invalid input_type: $input_type. Must be in [:mode, :mean, :init, :full]")
    end

    combine_raw_forecast_output_and_metadata(m_new, decomp_output_files, verbose = verbose)

    println(verbose, :low, "\nForecast decomposition complete: $(now())")
end

function decompose_forecast(m_new::M, m_old::M, df_new::DataFrame, df_old::DataFrame,
                            params_new::Vector{Float64}, params_old::Vector{Float64},
                            cond_new::Symbol, cond_old::Symbol, classes::Vector{Symbol};
                            check::Bool = false,
                            forecast_string_old = "",
                            forecast_string_new = "") where M<:AbstractDSGEModel
    # Check numbers of periods
    T, k, H = decomposition_periods(m_new, m_old, df_new, df_old, cond_new, cond_old)

    # Forecast
    f(m::AbstractDSGEModel, df::DataFrame, params::Vector{Float64}, cond_type::Symbol; kwargs...) =
        decomposition_forecast(m, df, params, cond_type, T, k, H;
                               kwargs..., check = check)

    out1 = f(m_new, df_new, params_new, cond_new, outputs = [:shockdec])            # new data, new params
    out2 = f(m_old, df_old, params_new, cond_old, outputs = [:shockdec, :forecast]) # old data, new params
    out3 = f(m_old, df_old, params_old, cond_old, outputs = [:forecast])            # old data, old params

    # Initialize output dictionary
    decomp = Dict{Symbol, Array{Float64}}()

    # Decomposition
    for class in classes
        # All elements of out are of size Ny x Nh, where the second dimension
        # ranges from t = T+1:T+H
        forecastvar = Symbol(:histforecast, class) # Z s_{T+h} + D
        trendvar    = Symbol(:trend,        class) # Z D
        dettrendvar = Symbol(:dettrend,     class) # Z T^{T+h} s_0 + D
        shockdecvar = Symbol(:shockdec,     class) # Z \sum_{t=1}^{T+h} T^{T+h-t} R ϵ_t + D
        datavar     = Symbol(:data,         class) # Z \sum_{t=1}^{T-k} T^{T+h-t} R ϵ_t + D
        newsvar     = Symbol(:news,         class) # Z \sum_{t=T-k+1}^{T+h} T^{T+h-t} R ϵ_t + D

        # 1(a). Data revision and news
        data_comp = (out1[dettrendvar] - out2[dettrendvar]) + (out1[datavar] - out2[datavar])
        decomp[Symbol(:decompdata, class)] = data_comp

        news_comp = out1[newsvar] - out2[newsvar]
        decomp[Symbol(:decompnews, class)] = news_comp

        # 1(b). Shock decomposition and deterministic trend
        shockdec_comp = out1[shockdecvar] - out2[shockdecvar] # Ny x Nh x Ne
        decomp[Symbol(:decompshockdec, class)] = shockdec_comp

        dettrend_comp = out1[dettrendvar] - out2[dettrendvar]
        decomp[Symbol(:decompdettrend, class)] = dettrend_comp

        # Check that 1(a) and 1(b) are equal
        check && @assert dettrend_comp + dropdims(sum(shockdec_comp, dims = 3), dims = 3) ≈ data_comp + news_comp

        # 2. Parameter re-estimation
        para_comp = out2[forecastvar] - out3[forecastvar]
        decomp[Symbol(:decomppara, class)] = para_comp

        # 1 + 2. Total difference
        total_diff = para_comp + data_comp + news_comp
        decomp[Symbol(:decomptotal, class)] = total_diff
        check && @assert total_diff ≈ out1[forecastvar] - out3[forecastvar]
    end

    return decomp
end

"""
```
decomposition_periods(m_new, m_old, df_new, df_old, cond_new, cond_old)
```

Returns `T`, `k`, and `H`, where:

- New model has `T` periods of data
- Old model has `T-k` periods of data
- Old and new models both forecast up to `T+H`
"""
function decomposition_periods(m_new::M, m_old::M, df_new::DataFrame, df_old::DataFrame,
                               cond_new::Symbol, cond_old::Symbol) where M<:AbstractDSGEModel
    # Number of presample periods T0 must be the same
    T0 = n_presample_periods(m_new)
    @assert n_presample_periods(m_old) == T0

    # New model has T main-sample periods
    # Old model has T-k main-sample periods
    T = n_mainsample_periods(m_new)
    k = subtract_quarters(date_forecast_start(m_new), date_forecast_start(m_old))
    @assert k >= 0

    # Number of conditional periods T1 may differ
    T1_new = cond_new == :none ? 0 : n_conditional_periods(m_new)
    T1_old = cond_old == :none ? 0 : n_conditional_periods(m_old)

    # Check DataFrame sizes
    @assert size(df_new, 1) == T0 + T + T1_new
    @assert size(df_old, 1) == T0 + T - k + T1_old

    # Old model forecasts up to T+H
    H = subtract_quarters(date_forecast_end(m_old), date_mainsample_end(m_new))

    return T, k, H
end

"""
```
decomposition_forecast(m, df, params, cond_type, keep_startdate, keep_enddate, shockdec_splitdate;
    outputs = [:forecast, :shockdec], check = false)
```

Equivalent of `forecast_one_draw` for forecast decomposition. `keep_startdate =
date_forecast_start(m_new)` corresponds to time T+1, `keep_enddate =
date_forecast_end(m_old)` to time T+H, and `shockdec_splitdate =
date_mainsample_end(m_old)` to time T-k.

Returns `out::Dict{Symbol, Array{Float64}}`, which has keys determined as follows:

- If `:forecast in outputs` or `check = true`:
  - `:forecast<class>`

- If `:shockdec in outputs`:
  - `:trend<class>`
  - `:dettrend<class>`
  - `:data<class>`: like a shockdec, but only applying smoothed shocks up to `shockdec_splitdate`
  - `:news<class>`: like a shockdec, but only applying smoothed shocks after `shockdec_splitdate`
"""
function decomposition_forecast(m::AbstractDSGEModel, df::DataFrame, params::Vector{Float64}, cond_type::Symbol,
                                T::Int, k::Int, H::Int;
                                outputs::Vector{Symbol} = [:forecast, :shockdec], check::Bool = false)
    # Compute state space
    DSGE.update!(m, params)
    system = compute_system(m)

    # Initialize output dictionary
    out = Dict{Symbol, Array{Float64}}()

    # Smooth and forecast
    histstates, histshocks, histpseudo, s_0 = smooth(m, df, system, cond_type = cond_type, draw_states = false)
    histobs = system[:ZZ] * histstates .+ system[:DD]

    if :forecast in outputs || check
        s_T = histstates[:, end]
        _, forecastobs, forecastpseudo, _ =
            forecast(m, system, s_T, cond_type = cond_type, enforce_zlb = false, draw_shocks = false)

        out[:histforecastobs]    = hcat(histobs,    forecastobs)[:, 1:T+H]
        out[:histforecastpseudo] = hcat(histpseudo, forecastpseudo)[:, 1:T+H]
    end

    # Compute trend, dettrend, and shockdecs
    if :shockdec in outputs
        nstates = n_states_augmented(m)
        nshocks = n_shocks_exogenous(m)

        _, out[:trendobs],    out[:trendpseudo]    = trends(system)
        _, out[:dettrendobs], out[:dettrendpseudo] = deterministic_trends(system, s_0, T+H, 1, T+H)

        # Applying all shocks
        _, out[:shockdecobs], out[:shockdecpseudo] =
            shock_decompositions(system, forecast_horizons(m), histshocks, 1, T+H)

        # Applying ϵ_{1:T-k} and ϵ_{T-k+1:end}
        system0 = zero_system_constants(system)

        data_shocks = zeros(nshocks, T+H)
        data_shocks[:, 1:T-k] = histshocks[:, 1:T-k]
        _, out[:dataobs], out[:datapseudo], _ = forecast(system0, zeros(nstates), data_shocks)

        Tstar = size(histshocks, 2) # either T or T+1
        news_shocks = zeros(nshocks, T+H)
        news_shocks[:, T-k+1:Tstar] = histshocks[:, T-k+1:Tstar]
        _, out[:newsobs], out[:newspseudo], _ = forecast(system0, zeros(nstates), news_shocks)
    end

    # Return
    return out
end


text:
using IRTools, Test
using IRTools: Meta, meta

@generated f(x) = :(x+x)

@test meta(Tuple{typeof(gcd),Int,Int}) isa Meta

@test meta(Tuple{typeof(f),Int}) isa Meta

@test @code_ir(map([1,2], [3,4]) do x, y
  x + y
end) isa IR

ext: jl


text:
open("data/dlrp/pre_dlrp.csv", "w") do fp
    open("data/dlrp/dlrp.txt", "r") do fp2
        ligand = true
        pairs = ["",""]
        for line in readlines(fp2)
            col = split(chomp(line), "\t")
            if length(col) != 1
                if ligand == true
                    pairs[1] = col[2]
                    ligand = false
                else
                    pairs[2] = col[2]
                    write(fp, pairs[1] * "," * pairs[2] * ",-," * "DLRP")
                    write(fp, "\n")
                end
            else
                ligand = true
            end
        end
    end
end

max_stars_repo_name: manabuishii/lrbase-workflow

max_issues_repo_name: manabuishii/lrbase-workflow

max_forks_repo_name: manabuishii/lrbase-workflow


text:
struct Iris <: Tabular end
url(::Iris) = "http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"
checksum(::Iris) = "6f608b71a7317216319b4d27b4d9bc84e6abd734eda7872b71a458569e2656c0"
prep(::Iris) = path -> preprocess(path, Iris())
target(::Iris) = 5
size(::Iris) = (150, 0, 0)
headers(::Iris) = ["sepal length", "sepal width", "petal length", "petal width", "class"]
problem(::Iris) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using SlowActors

function fib(D::Dict{Int,BigInt}, n::Int)
    get!(D, n) do
        n == 0 && return big(0)
        n == 1 && return big(1)
        return fib(D, n-1) + fib(D, n-2)
    end
end

myfib = spawn(Func(fib, Dict{Int,BigInt}()))

request!(myfib, 1000)


text:
using LyceumMuJoCoViz
using Test

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


text:
# Coverage submission script on travis-ci.com, delete or keep depending on whether you are using that service.
# Only pushes coverage from the given version and OS.
get(ENV, "TRAVIS_OS_NAME", nothing)       == "linux" || exit(0)
get(ENV, "TRAVIS_JULIA_VERSION", nothing) == "1.5"   || exit(0)

using Coverage

cd(joinpath(@__DIR__, "..", "..")) do
    Codecov.submit(Codecov.process_folder())
end


text:
using StaticArrays
import Base: show

"""Variable for a `DecisionTree{N, C}`."""
struct Variable{N}
    n::Int
    Variable{N}(m) where N = (1 <= m <= N) ? new{N}(m) : error("Illegal construction")
end

"""Decision tree for `N` features and `C` classes, consisting of `Classification`s and `Desicion`s."""
abstract type DecisionTree{N, C} end

struct Classification{N, C} <: DecisionTree{N, C}
    index::Int
    # label
    Classification{N, C}(m) where {N, C} =
        (1 <= m <= C) ? new{N, C}(m) : error("Illegal construction")
end

struct Decision{N, C} <: DecisionTree{N, C}
    conditions::Dict{Variable{N}, Float64}
    threshold::Float64
    iftrue::DecisionTree{N, C}
    iffalse::DecisionTree{N, C}

    Decision{N, C}(c, t, ift::DecisionTree{N, C}, iff::DecisionTree{N, C}) where {N, C} =
        new{N, C}(c, t, ift, iff)
    Decision(c, t, ift::DecisionTree{N, C}, iff::DecisionTree{N, C}) where {N, C} =
        new{N, C}(c, t, ift, iff)
end

function Base.show(io::IO, v::Variable, indent = 0)
    print(io, "x")
    join(io, map(d -> '\u2080' + d, digits(v.n)) |> reverse)
end

function Base.show(io::IO, c::Classification, indent = 0)
    print(io, "{", c.index, "}")
end

function Base.show(io::IO, d::Decision, indent = 0; digits = 2)
    print(io, "if ")
    if !isempty(d.conditions)
        sorted_conditions = sort(collect(d.conditions), by = c -> c.first.n)
        join(io, ("$(round(f, digits)) × $v" for (v, f) in sorted_conditions), " + ")
    else
        print(io, "0.0")
    end
    print(io, " ≤ ", round(d.threshold, digits), "\n")
    print(io, " " ^ (indent + 2), "then ")
    show(io, d.iftrue, indent + 2)
    print(io, "\n", " " ^ (indent + 2), "else ")
    show(io, d.iffalse, indent + 2)
end

treedepth(::Classification) = 1
treedepth(d::Decision) = 1 + max(treedepth(d.iftrue), treedepth(d.iffalse))

treesize(::Classification) = 1
treesize(d::Decision) = 1 + treesize(d.iftrue) + treesize(d.iffalse)

"""
    decide{N}(value::StaticVector, t::DecisionTree)

Evaluate the decision tree `t` on a data point, returning a `Variable` indicating the result.
"""
function decide(datapoint::AbstractVector{Float64}, d::Decision)
    combination = 0.0
    for (v, λ) in d.conditions
        combination += λ * datapoint[v.n]
    end

    if combination <= d.threshold
        return decide(datapoint, d.iftrue)
    else
        return decide(datapoint, d.iffalse)
    end
end

function decide(::AbstractVector{Float64}, c::Classification)
    return c
end

δ{T}(i, j, ::Type{T} = Float64) = ifelse(i == j, one(T), zero(T))

function basisvector_impl(n, i)
    ee = [δ(i, j) for j in 1:n]
    Expr(:call, :SVector, ee...)
end

macro basisvector(n, i)
    @assert n isa Integer
    @assert i isa Integer
    basisvector_impl(n, i)
end

@generated function basisvector{N, I}(::Type{Val{N}}, ::Type{Val{I}})
    basisvector_impl(N, I)
end

onehot{N}(v::Variable{N}) = basisvector(Val{N}, Val{v.n})

"""
    normalize_conditions{N}(conditions::Dict{Variable{N}, Float64})

Convert a condition dictionary into a vector suitable for linear decision.
"""
function normalize_conditions{N}(conditions::Dict{Variable{N}, Float64})
    z = @SVector zeros(N)
    reduce((v, cond) -> v + onehot(cond[1]) * cond[2], z, conditions)
end

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
"""
$(TYPEDEF)

Represents a differential operator.

# Fields
$(FIELDS)

# Examples

```jldoctest
julia> using Symbolics

julia> @variables x y;

julia> D = Differential(x)
(D'~x)

julia> D(y) # Differentiate y wrt. x
(D'~x)(y)

julia> Dx = Differential(x) * Differential(y) # d^2/dxy operator
(D'~x(t)) ∘ (D'~y(t))

julia> D3 = Differential(x)^3 # 3rd order differential operator
(D'~x(t)) ∘ (D'~x(t)) ∘ (D'~x(t))
```
"""
struct Differential <: Function
    """The variable or expression to differentiate with respect to."""
    x
    Differential(x) = new(value(x))
end
(D::Differential)(x) = Term{symtype(x)}(D, [x])
(D::Differential)(x::Num) = Num(D(value(x)))
SymbolicUtils.promote_symtype(::Differential, x) = x

is_derivative(x::Term) = operation(x) isa Differential
is_derivative(x) = false

Base.:*(D1, D2::Differential) = D1 ∘ D2
Base.:*(D1::Differential, D2) = D1 ∘ D2
Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2
Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n)

Base.show(io::IO, D::Differential) = print(io, "Differential(", D.x, ")")

Base.:(==)(D1::Differential, D2::Differential) = isequal(D1.x, D2.x)

_isfalse(occ::Bool) = occ === false
_isfalse(occ::Term) = _isfalse(operation(occ))

function occursin_info(x, expr)
    !istree(expr) && return false
    if isequal(x, expr)
        true
    else
        args = map(a->occursin_info(x, a), arguments(expr))
        if all(_isfalse, args)
            return false
        end
        Term{Real}(true, args)
    end
end
function occursin_info(x, expr::Sym)
    isequal(x, expr)
end

function hasderiv(O)
    istree(O) ? operation(O) isa Differential || any(hasderiv, arguments(O)) : false
end
"""
$(SIGNATURES)

TODO
"""
function expand_derivatives(O::Symbolic, simplify=false; occurances=nothing)
    if istree(O) && isa(operation(O), Differential)
        @assert length(arguments(O)) == 1
        arg = expand_derivatives(arguments(O)[1], false)

        if occurances == nothing
            occurances = occursin_info(operation(O).x, arg)
        end

        _isfalse(occurances) && return 0
        occurances isa Bool && return 1 # means it's a `true`

        D = operation(O)

        if !istree(arg)
            return D(arg) # Cannot expand
        elseif isa(operation(arg), Sym)
            return D(arg) # Cannot expand
        elseif isa(operation(arg), Differential)
            # The recursive expand_derivatives was not able to remove
            # a nested Differential. We can attempt to differentiate the
            # inner expression wrt to the outer iv. And leave the
            # unexpandable Differential outside.
            if isequal(operation(arg).x, D.x)
                return D(arg)
            else
                inner = expand_derivatives(D(arguments(arg)[1]), false)
                # if the inner expression is not expandable either, return
                if istree(inner) && operation(inner) isa Differential
                    return D(arg)
                else
                    return expand_derivatives(operation(arg)(inner), simplify)
                end
            end
        end

        l = length(arguments(arg))
        exprs = []
        c = 0

        for i in 1:l
            t2 = expand_derivatives(D(arguments(arg)[i]),false, occurances=arguments(occurances)[i])

            x = if _iszero(t2)
                t2
            elseif _isone(t2)
                d = derivative_idx(arg, i)
                d isa NoDeriv ? D(arg) : d
            else
                t1 = derivative_idx(arg, i)
                t1 = t1 isa NoDeriv ? D(arg) : t1
                t1 * t2
            end

            if _iszero(x)
                continue
            elseif x isa Symbolic
                push!(exprs, x)
            else
                c += x
            end
        end

        if isempty(exprs)
            return c
        elseif length(exprs) == 1
            term = (simplify ? SymbolicUtils.simplify(exprs[1]) : exprs[1])
            return _iszero(c) ? term : c + term
        else
            x = +((!_iszero(c) ? vcat(c, exprs) : exprs)...)
            return simplify ? SymbolicUtils.simplify(x) : x
        end
    elseif !hasderiv(O)
        return O
    else
        args = map(a->expand_derivatives(a, false), arguments(O))
        O1 = operation(O)(args...)
        return simplify ? SymbolicUtils.simplify(O1) : O1
    end
end

function expand_derivatives(n::Num, simplify=false; occurances=nothing)
    Num(expand_derivatives(value(n), simplify; occurances=occurances))
end

_iszero(x) = false
_isone(x) = false

expand_derivatives(x, simplify=false;occurances=nothing) = x

# Don't specialize on the function here
"""
$(SIGNATURES)

Calculate the derivative of the op `O` with respect to its argument with index
`idx`.

# Examples

```jldoctest label1
julia> using Symbolics

julia> @variables x y;

julia> Symbolics.derivative_idx(Symbolics.value(sin(x)), 1)
cos(x)
```

Note that the function does not recurse into the operation's arguments, i.e., the
chain rule is not applied:

```jldoctest label1
julia> myop = Symbolics.value(sin(x) * y^2)
sin(x)*(y^2)

julia> typeof(Symbolics.operation(myop))  # Op is multiplication function
typeof(*)

julia> Symbolics.derivative_idx(myop, 1)  # wrt. sin(x)
y^2

julia> Symbolics.derivative_idx(myop, 2)  # wrt. y^2
sin(x)
```
"""
derivative_idx(O::Any, ::Any) = 0
function derivative_idx(O::Symbolic, idx)
    istree(O) ? derivative(operation(O), (arguments(O)...,), Val(idx)) : 0
end

# Indicate that no derivative is defined.
struct NoDeriv
end
derivative(f, args, v) = NoDeriv()

# Pre-defined derivatives
import DiffRules
for (modu, fun, arity) ∈ DiffRules.diffrules()
    fun in [:*, :+, :abs, :mod, :rem, :max, :min] && continue # special
    for i ∈ 1:arity

        expr = if arity == 1
            DiffRules.diffrule(modu, fun, :(args[1]))
        else
            DiffRules.diffrule(modu, fun, ntuple(k->:(args[$k]), arity)...)[i]
        end
        @eval derivative(::typeof($modu.$fun), args::NTuple{$arity,Any}, ::Val{$i}) = $expr
    end
end

derivative(::typeof(+), args::NTuple{N,Any}, ::Val) where {N} = 1
derivative(::typeof(*), args::NTuple{N,Any}, ::Val{i}) where {N,i} = *(deleteat!(collect(args), i)...)
derivative(::typeof(one), args::Tuple{<:Any}, ::Val) = 0

function count_order(x)
    @assert !(x isa Symbol) "The variable $x must have an order of differentiation that is greater or equal to 1!"
    n = 1
    while !(x.args[1] isa Symbol)
        n = n+1
        x = x.args[1]
    end
    n, x.args[1]
end

_repeat_apply(f, n) = n == 1 ? f : f ∘ _repeat_apply(f, n-1)
function _differential_macro(x)
    ex = Expr(:block)
    push!(ex.args,  :(Base.depwarn("`@derivatives D'''~x` is deprecated. Use `Differential(x)^3` instead.", Symbol("@derivatives"), force=true)))
    lhss = Symbol[]
    x = x isa Tuple && first(x).head == :tuple ? first(x).args : x # tuple handling
    x = flatten_expr!(x)
    for di in x
        @assert di isa Expr && di.args[1] == :~ "@derivatives expects a form that looks like `@derivatives D''~t E'~t` or `@derivatives (D''~t), (E'~t)`"
        lhs = di.args[2]
        rhs = di.args[3]
        order, lhs = count_order(lhs)
        push!(lhss, lhs)
        expr = :($lhs = $_repeat_apply(Differential($value($rhs)), $order))
        push!(ex.args,  expr)
    end
    push!(ex.args, Expr(:tuple, lhss...))
    ex
end

"""
$(SIGNATURES)

Define one or more differentials.

# Examples

```jldoctest
julia> using Symbolics

julia> @variables x y z;

julia> Dx = Differential(x); Dy = Differential(y);  # Create differentials wrt. x and y

julia> Dx(z)  # Differentiate z wrt. x
Differential(x)(z)

julia> Dy(z)  # Differentiate z wrt. y
Differential(y)(z)
```
"""
macro derivatives(x...)
    esc(_differential_macro(x))
end

### Jacobians & Hessians

"""
$(SIGNATURES)

A helper function for computing the derivative of an expression with respect to
`var`.
"""
function derivative(O, v; simplify=false)
    if O isa AbstractArray
        Num[Num(expand_derivatives(Differential(v)(value(o)), simplify)) for o in O]
    else
        Num(expand_derivatives(Differential(v)(value(O)), simplify))
    end
end

"""
$(SIGNATURES)

A helper function for computing the gradient of an expression with respect to
an array of variable expressions.
"""
function gradient(O, vars::AbstractVector; simplify=false)
    Num[Num(expand_derivatives(Differential(v)(value(O)),simplify)) for v in vars]
end

"""
$(SIGNATURES)

A helper function for computing the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian(ops::AbstractVector, vars::AbstractVector; simplify=false)
    Num[Num(expand_derivatives(Differential(value(v))(value(O)),simplify)) for O in ops, v in vars]
end

"""
$(SIGNATURES)

A helper function for computing the sparse Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function sparsejacobian(ops::AbstractVector, vars::AbstractVector; simplify=false)
    I = Int[]
    J = Int[]
    du = Num[]

    sp = jacobian_sparsity(ops, vars)
    I,J,_ = findnz(sp)

    exprs = Num[]

    for (i,j) in zip(I, J)
        push!(exprs, Num(expand_derivatives(Differential(vars[j])(ops[i]), simplify)))
    end
    sparse(I,J, exprs, length(ops), length(vars))
end

"""
```julia
jacobian_sparsity(ops::AbstractVector, vars::AbstractVector)
```

Return the sparsity pattern of the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian_sparsity(du, u)
    du = map(value, du)
    u = map(value, u)
    dict = Dict(zip(u, 1:length(u)))

    i = Ref(1)
    I = Int[]
    J = Int[]

    # This rewriter notes down which u's appear in a
    # given du (whose index is stored in the `i` Ref)
    r = [@rule ~x::(x->haskey(dict, x)) => begin
        push!(I, i[])
        push!(J, dict[~x])
        nothing
    end] |> Rewriters.Chain |> Rewriters.Postwalk

    for ii = 1:length(du)
        i[] = ii
        r(du[ii])
    end

    sparse(I, J, true, length(du), length(u))
end

"""
    exprs_occur_in(exprs::Vector, expr)

Return an array of booleans `finds` where `finds[i]` is true if `exprs[i]` occurs in `expr`
false otherwise.
"""
function exprs_occur_in(exprs, expr)
    vec(jacobian_sparsity([expr], exprs))
end

"""
$(SIGNATURES)

A helper function for computing the Hessian of an expression with respect to
an array of variable expressions.
"""
function hessian(O, vars::AbstractVector; simplify=false)
    vars = map(value, vars)
    first_derivs = map(value, vec(jacobian([values(O)], vars, simplify=simplify)))
    n = length(vars)
    H = Array{Num, 2}(undef,(n, n))
    fill!(H, 0)
    for i=1:n
        for j=1:i
            H[j, i] = H[i, j] = expand_derivatives(Differential(vars[i])(first_derivs[j]))
        end
    end
    H
end

isidx(x) = x isa TermCombination

"""
    hessian_sparsity(ops::AbstractVector, vars::AbstractVector)

Return the sparsity pattern of the Hessian of an array of expressions with respect to
an array of variable expressions.
"""
function hessian_sparsity end

let
    # we do this in a let block so that Revise works on the list of rules

    _scalar = one(TermCombination)

    simterm(t, f, args) = Term{Any}(f, args)
    linearity_rules = [
          @rule +(~~xs) => reduce(+, filter(isidx, ~~xs), init=_scalar)
          @rule *(~~xs) => reduce(*, filter(isidx, ~~xs), init=_scalar)
          @rule (~f)(~x::(!isidx)) => _scalar
          @rule (~f)(~x::isidx) => if haslinearity_1(~f)
              combine_terms_1(linearity_1(~f), ~x)
          else
              error("Function of unknown linearity used: ", ~f)
          end
          @rule (^)(~x::isidx, ~y) => ~y isa Number && isone(~y) ? ~x : (~x) * (~x)
          @rule (~f)(~x, ~y) => begin
              if haslinearity_2(~f)
                  a = isidx(~x) ? ~x : _scalar
                  b = isidx(~y) ? ~y : _scalar
                  combine_terms_2(linearity_2(~f), a, b)
              else
                  error("Function of unknown linearity used: ", ~f)
              end
          end]
    linearity_propagator = Fixpoint(Postwalk(Chain(linearity_rules); similarterm=simterm))

    global hessian_sparsity

    function hessian_sparsity(f, u)
        @assert !(f isa AbstractArray)
        f = value(f)
        u = map(value, u)
        idx(i) = TermCombination(Set([Dict(i=>1)]))
        dict = Dict(u .=> idx.(1:length(u)))
        f = Rewriters.Prewalk(x->haskey(dict, x) ? dict[x] : x; similarterm=simterm)(f)
        lp = linearity_propagator(f)
        _sparse(lp, length(u))
    end
end

"""
$(SIGNATURES)

Check if an expression is linear with respect to a list of variable expressions.
"""
function islinear(ex, u)
    isempty(hessian_sparsity(ex, u).nzval)
end

"""
$(SIGNATURES)

A helper function for computing the sparse Hessian of an expression with respect to
an array of variable expressions.
"""
function sparsehessian(O, vars::AbstractVector; simplify=false)
    O = value(O)
    vars = map(value, vars)
    S = hessian_sparsity(O, vars)
    I, J, _ = findnz(S)
    exprs = Array{Num}(undef, length(I))
    fill!(exprs, 0)
    prev_j = 0
    d = nothing
    for (k, (i, j)) in enumerate(zip(I, J))
        j > i && continue
        if j != prev_j
            d = expand_derivatives(Differential(vars[j])(O), false)
        end
        expr = expand_derivatives(Differential(vars[i])(d), simplify)
        exprs[k] = expr
        prev_j = j
    end
    H = sparse(I, J, exprs, length(vars), length(vars))
    for (i, j) in zip(I, J)
        j > i && (H[i, j] = H[j, i])
    end
    return H
end

SymbolicUtils.Code.function_to_expr(::Differential, o, st) = toexpr(diff2term(o), st)


text:
#! this file is not used anywhere in the package

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW #!
using CImGui.OpenGLBackend.ModernGL
using Printf

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3)
    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2)
    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE) # 3.2+ only
    glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3)
    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 0)
    # glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE) # 3.2+ only
    # glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
#? error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
#? GLFW.SetErrorCallback(error_callback)

# custom callbacks
#? function CustomUserCallbackMousebutton(window::glfwWindow, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint)
#?     println("Trigger user mouse button callback!")
#?     return nothing
#? end
SetCustomMouseButtonCallback(CustomUserCallbackMousebutton)

# create window
window = glfwCreateWindow(1280, 720, "Demo", C_NULL, C_NULL)
@assert window != C_NULL
glfwMakeContextCurrent(window)
glfwSwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsDark()

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

show_demo_window = true
clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
while !glfwWindowShouldClose(window)
    # oh my global scope
    global show_demo_window
    global clear_color

    glfwPollEvents()
    # start the Dear ImGui frame
    ImGui_ImplOpenGL3_NewFrame()
    ImGui_ImplGlfw_NewFrame()
    CImGui.NewFrame()

    # show the big demo window
    show_demo_window && @c CImGui.ShowDemoWindow(&show_demo_window)

    # rendering
    CImGui.Render()
    glfwMakeContextCurrent(window)
    display_w, display_h = glfwGetFramebufferSize(window)
    glViewport(0, 0, display_w, display_h)
    glClearColor(clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)
    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

    glfwMakeContextCurrent(window)
    glfwSwapBuffers(window)
end

# cleanup
ImGui_ImplOpenGL3_Shutdown()
ImGui_ImplGlfw_Shutdown()
CImGui.DestroyContext(ctx)

glfwDestroyWindow(window)


text:
using Vega

scatterplot(x = rand(100), y = rand(100))

using Vega, Distributions

d1 = MultivariateNormal([0.0, 0.0], [1.0 0.9; 0.9 1.0])
d2 = MultivariateNormal([10.0, 10.0], [4.0 0.5; 0.5 4.0])
points = vcat(rand(d1, 500)', rand(d2, 500)')

x = points[:, 1]
y = points[:, 2]
group = vcat(ones(Int, 500), ones(Int, 500) + 1)

scatterplot(x = x, y = y, group = group)

x = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9]
y = [28, 43, 81, 19, 52, 24, 87, 17, 68, 49, 55, 91, 53, 87, 48, 49, 66, 27, 16, 15]
g = [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]

a = areaplot(x = x, y = y, group = g, stacked = true)

fruit = ["peaches", "plums", "blueberries", "strawberries", "bananas"];
bushels = [100, 32, 180, 46, 21];
piechart(x = fruit, y = bushels, holesize = 125)

x = Array(Int, 900)
y = Array(Int, 900)
color = Array(Float64, 900)
tmp = 0
for counter in 1:30
    for counter2 in 1:30
        tmp += 1
        x[tmp] = counter
        y[tmp] = counter2
        color[tmp] = rand()
    end
end
hm = heatmap(x = x, y = y, color = color)

max_stars_repo_path: Chapter 5/ch05-vega.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 5/ch05-vega.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 5/ch05-vega.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
register_color_library(:colorbrewer, ColorLibrary(Dict(:default => :sequential, :sequential => :Spectral, :divergent => :RdBu)))

register_gradient_colors(:Spectral, [
    RGB(158/255, 1/255, 66/255),
    RGB(213/255, 62/255, 79/255),
    RGB(244/255, 109/255, 67/255),
    RGB(253/255, 174/255, 97/255),
    RGB(254/255, 224/255, 139/255),
    RGB(255/255, 255/255, 191/255),
    RGB(230/255, 245/255, 152/255),
    RGB(171/255, 221/255, 164/255),
    RGB(102/255, 194/255, 165/255),
    RGB(50/255, 136/255, 189/255),
    RGB(94/255, 79/255, 162/255)], :colorbrewer)

register_gradient_colors(:RdYlGn, [
    RGB(165/255, 0/255, 38/255),
    RGB(215/255, 48/255, 39/255),
    RGB(244/255, 109/255, 67/255),
    RGB(253/255, 174/255, 97/255),
    RGB(254/255, 224/255, 139/255),
    RGB(255/255, 255/255, 191/255),
    RGB(217/255, 239/255, 139/255),
    RGB(166/255, 217/255, 106/255),
    RGB(102/255, 189/255, 99/255),
    RGB(26/255, 152/255, 80/255),
    RGB(0/255, 104/255, 55/255)], :colorbrewer)

register_gradient_colors(:RdBu, [
    RGB(103/255, 0/255, 31/255),
    RGB(178/255, 24/255, 43/255),
    RGB(214/255, 96/255, 77/255),
    RGB(244/255, 165/255, 130/255),
    RGB(253/255, 219/255, 199/255),
    RGB(247/255, 247/255, 247/255),
    RGB(209/255, 229/255, 240/255),
    RGB(146/255, 197/255, 222/255),
    RGB(67/255, 147/255, 195/255),
    RGB(33/255, 102/255, 172/255),
    RGB(5/255, 48/255, 97/255)], :colorbrewer)

register_gradient_colors(:PiYG, [
    RGB(142/255, 1/255, 82/255),
    RGB(197/255, 27/255, 125/255),
    RGB(222/255, 119/255, 174/255),
    RGB(241/255, 182/255, 218/255),
    RGB(253/255, 224/255, 239/255),
    RGB(247/255, 247/255, 247/255),
    RGB(230/255, 245/255, 208/255),
    RGB(184/255, 225/255, 134/255),
    RGB(127/255, 188/255, 65/255),
    RGB(77/255, 146/255, 33/255),
    RGB(39/255, 100/255, 25/255)], :colorbrewer)

register_gradient_colors(:PRGn, [
    RGB(64/255, 0/255, 75/255),
    RGB(118/255, 42/255, 131/255),
    RGB(153/255, 112/255, 171/255),
    RGB(194/255, 165/255, 207/255),
    RGB(231/255, 212/255, 232/255),
    RGB(247/255, 247/255, 247/255),
    RGB(217/255, 240/255, 211/255),
    RGB(166/255, 219/255, 160/255),
    RGB(90/255, 174/255, 97/255),
    RGB(27/255, 120/255, 55/255),
    RGB(0/255, 68/255, 27/255)], :colorbrewer)

register_gradient_colors(:RdYlBu, [
    RGB(165/255, 0/255, 38/255),
    RGB(215/255, 48/255, 39/255),
    RGB(244/255, 109/255, 67/255),
    RGB(253/255, 174/255, 97/255),
    RGB(254/255, 224/255, 144/255),
    RGB(255/255, 255/255, 191/255),
    RGB(224/255, 243/255, 248/255),
    RGB(171/255, 217/255, 233/255),
    RGB(116/255, 173/255, 209/255),
    RGB(69/255, 117/255, 180/255),
    RGB(49/255, 54/255, 149/255)], :colorbrewer)

register_gradient_colors(:BrBG, [
    RGB(84/255, 48/255, 5/255),
    RGB(140/255, 81/255, 10/255),
    RGB(191/255, 129/255, 45/255),
    RGB(223/255, 194/255, 125/255),
    RGB(246/255, 232/255, 195/255),
    RGB(245/255, 245/255, 245/255),
    RGB(199/255, 234/255, 229/255),
    RGB(128/255, 205/255, 193/255),
    RGB(53/255, 151/255, 143/255),
    RGB(1/255, 102/255, 94/255),
    RGB(0/255, 60/255, 48/255)], :colorbrewer)

register_gradient_colors(:RdGy, [
    RGB(103/255, 0/255, 31/255),
    RGB(178/255, 24/255, 43/255),
    RGB(214/255, 96/255, 77/255),
    RGB(244/255, 165/255, 130/255),
    RGB(253/255, 219/255, 199/255),
    RGB(255/255, 255/255, 255/255),
    RGB(224/255, 224/255, 224/255),
    RGB(186/255, 186/255, 186/255),
    RGB(135/255, 135/255, 135/255),
    RGB(77/255, 77/255, 77/255),
    RGB(26/255, 26/255, 26/255)], :colorbrewer)

# These are commented out as they constitute non-continuous gradients, and are thus used differently by Plots
# register_gradient_colors(:Set2, [
#     RGB(102/255, 194/255, 165/255),
#     RGB(252/255, 141/255, 98/255),
#     RGB(141/255, 160/255, 203/255),
#     RGB(231/255, 138/255, 195/255),
#     RGB(166/255, 216/255, 84/255),
#     RGB(255/255, 217/255, 47/255),
#     RGB(229/255, 196/255, 148/255),
#     RGB(179/255, 179/255, 179/255)], :colorbrewer)
#
# register_gradient_colors(:Accent, [
#     RGB(127/255, 201/255, 127/255),
#     RGB(190/255, 174/255, 212/255),
#     RGB(253/255, 192/255, 134/255),
#     RGB(255/255, 255/255, 153/255),
#     RGB(56/255, 108/255, 176/255),
#     RGB(240/255, 2/255, 127/255),
#     RGB(191/255, 91/255, 23/255),
#     RGB(102/255, 102/255, 102/255)], :colorbrewer)
#
# register_gradient_colors(:Set1, [
#     RGB(228/255, 26/255, 28/255),
#     RGB(55/255, 126/255, 184/255),
#     RGB(77/255, 175/255, 74/255),
#     RGB(152/255, 78/255, 163/255),
#     RGB(255/255, 127/255, 0/255),
#     RGB(255/255, 255/255, 51/255),
#     RGB(166/255, 86/255, 40/255),
#     RGB(247/255, 129/255, 191/255),
#     RGB(153/255, 153/255, 153/255)], :colorbrewer)
#
# register_gradient_colors(:Set3, [
#     RGB(141/255, 211/255, 199/255),
#     RGB(255/255, 255/255, 179/255),
#     RGB(190/255, 186/255, 218/255),
#     RGB(251/255, 128/255, 114/255),
#     RGB(128/255, 177/255, 211/255),
#     RGB(253/255, 180/255, 98/255),
#     RGB(179/255, 222/255, 105/255),
#     RGB(252/255, 205/255, 229/255),
#     RGB(217/255, 217/255, 217/255),
#     RGB(188/255, 128/255, 189/255),
#     RGB(204/255, 235/255, 197/255),
#     RGB(255/255, 237/255, 111/255)], :colorbrewer)
#
# register_gradient_colors(:Dark2, [
#     RGB(27/255, 158/255, 119/255),
#     RGB(217/255, 95/255, 2/255),
#     RGB(117/255, 112/255, 179/255),
#     RGB(231/255, 41/255, 138/255),
#     RGB(102/255, 166/255, 30/255),
#     RGB(230/255, 171/255, 2/255),
#     RGB(166/255, 118/255, 29/255),
#     RGB(102/255, 102/255, 102/255)], :colorbrewer)
#
# register_gradient_colors(:Paired, [
#     RGB(166/255, 206/255, 227/255),
#     RGB(31/255, 120/255, 180/255),
#     RGB(178/255, 223/255, 138/255),
#     RGB(51/255, 160/255, 44/255),
#     RGB(251/255, 154/255, 153/255),
#     RGB(227/255, 26/255, 28/255),
#     RGB(253/255, 191/255, 111/255),
#     RGB(255/255, 127/255, 0/255),
#     RGB(202/255, 178/255, 214/255),
#     RGB(106/255, 61/255, 154/255),
#     RGB(255/255, 255/255, 153/255),
#     RGB(177/255, 89/255, 40/255)], :colorbrewer)
#
# register_gradient_colors(:Pastel2, [
#     RGB(179/255, 226/255, 205/255),
#     RGB(253/255, 205/255, 172/255),
#     RGB(203/255, 213/255, 232/255),
#     RGB(244/255, 202/255, 228/255),
#     RGB(230/255, 245/255, 201/255),
#     RGB(255/255, 242/255, 174/255),
#     RGB(241/255, 226/255, 204/255),
#     RGB(204/255, 204/255, 204/255)], :colorbrewer)
#
# register_gradient_colors(:Pastel1, [
#     RGB(251/255, 180/255, 174/255),
#     RGB(179/255, 205/255, 227/255),
#     RGB(204/255, 235/255, 197/255),
#     RGB(222/255, 203/255, 228/255),
#     RGB(254/255, 217/255, 166/255),
#     RGB(255/255, 255/255, 204/255),
#     RGB(229/255, 216/255, 189/255),
#     RGB(253/255, 218/255, 236/255),
#     RGB(242/255, 242/255, 242/255)], :colorbrewer)

register_gradient_colors(:PuOr, [
    RGB(127/255, 59/255, 8/255),
    RGB(179/255, 88/255, 6/255),
    RGB(224/255, 130/255, 20/255),
    RGB(253/255, 184/255, 99/255),
    RGB(254/255, 224/255, 182/255),
    RGB(247/255, 247/255, 247/255),
    RGB(216/255, 218/255, 235/255),
    RGB(178/255, 171/255, 210/255),
    RGB(128/255, 115/255, 172/255),
    RGB(84/255, 39/255, 136/255),
    RGB(45/255, 0/255, 75/255)], :colorbrewer)

register_gradient_colors(:OrRd, [
    RGB(255/255, 247/255, 236/255),
    RGB(254/255, 232/255, 200/255),
    RGB(253/255, 212/255, 158/255),
    RGB(253/255, 187/255, 132/255),
    RGB(252/255, 141/255, 89/255),
    RGB(239/255, 101/255, 72/255),
    RGB(215/255, 48/255, 31/255),
    RGB(179/255, 0/255, 0/255),
    RGB(127/255, 0/255, 0/255)], :colorbrewer)

register_gradient_colors(:PuBu, [
    RGB(255/255, 247/255, 251/255),
    RGB(236/255, 231/255, 242/255),
    RGB(208/255, 209/255, 230/255),
    RGB(166/255, 189/255, 219/255),
    RGB(116/255, 169/255, 207/255),
    RGB(54/255, 144/255, 192/255),
    RGB(5/255, 112/255, 176/255),
    RGB(4/255, 90/255, 141/255),
    RGB(2/255, 56/255, 88/255)], :colorbrewer)

register_gradient_colors(:BuPu, [
    RGB(247/255, 252/255, 253/255),
    RGB(224/255, 236/255, 244/255),
    RGB(191/255, 211/255, 230/255),
    RGB(158/255, 188/255, 218/255),
    RGB(140/255, 150/255, 198/255),
    RGB(140/255, 107/255, 177/255),
    RGB(136/255, 65/255, 157/255),
    RGB(129/255, 15/255, 124/255),
    RGB(77/255, 0/255, 75/255)], :colorbrewer)

register_gradient_colors(:Oranges, [
    RGB(255/255, 245/255, 235/255),
    RGB(254/255, 230/255, 206/255),
    RGB(253/255, 208/255, 162/255),
    RGB(253/255, 174/255, 107/255),
    RGB(253/255, 141/255, 60/255),
    RGB(241/255, 105/255, 19/255),
    RGB(217/255, 72/255, 1/255),
    RGB(166/255, 54/255, 3/255),
    RGB(127/255, 39/255, 4/255)], :colorbrewer)

register_gradient_colors(:BuGn, [
    RGB(247/255, 252/255, 253/255),
    RGB(229/255, 245/255, 249/255),
    RGB(204/255, 236/255, 230/255),
    RGB(153/255, 216/255, 201/255),
    RGB(102/255, 194/255, 164/255),
    RGB(65/255, 174/255, 118/255),
    RGB(35/255, 139/255, 69/255),
    RGB(0/255, 109/255, 44/255),
    RGB(0/255, 68/255, 27/255)], :colorbrewer)

register_gradient_colors(:YlOrBr, [
    RGB(255/255, 255/255, 229/255),
    RGB(255/255, 247/255, 188/255),
    RGB(254/255, 227/255, 145/255),
    RGB(254/255, 196/255, 79/255),
    RGB(254/255, 153/255, 41/255),
    RGB(236/255, 112/255, 20/255),
    RGB(204/255, 76/255, 2/255),
    RGB(153/255, 52/255, 4/255),
    RGB(102/255, 37/255, 6/255)], :colorbrewer)

register_gradient_colors(:YlGn, [
    RGB(255/255, 255/255, 229/255),
    RGB(247/255, 252/255, 185/255),
    RGB(217/255, 240/255, 163/255),
    RGB(173/255, 221/255, 142/255),
    RGB(120/255, 198/255, 121/255),
    RGB(65/255, 171/255, 93/255),
    RGB(35/255, 132/255, 67/255),
    RGB(0/255, 104/255, 55/255),
    RGB(0/255, 69/255, 41/255)], :colorbrewer)

register_gradient_colors(:Reds, [
    RGB(255/255, 245/255, 240/255),
    RGB(254/255, 224/255, 210/255),
    RGB(252/255, 187/255, 161/255),
    RGB(252/255, 146/255, 114/255),
    RGB(251/255, 106/255, 74/255),
    RGB(239/255, 59/255, 44/255),
    RGB(203/255, 24/255, 29/255),
    RGB(165/255, 15/255, 21/255),
    RGB(103/255, 0/255, 13/255)], :colorbrewer)

register_gradient_colors(:RdPu, [
    RGB(255/255, 247/255, 243/255),
    RGB(253/255, 224/255, 221/255),
    RGB(252/255, 197/255, 192/255),
    RGB(250/255, 159/255, 181/255),
    RGB(247/255, 104/255, 161/255),
    RGB(221/255, 52/255, 151/255),
    RGB(174/255, 1/255, 126/255),
    RGB(122/255, 1/255, 119/255),
    RGB(73/255, 0/255, 106/255)], :colorbrewer)

register_gradient_colors(:Greens, [
    RGB(247/255, 252/255, 245/255),
    RGB(229/255, 245/255, 224/255),
    RGB(199/255, 233/255, 192/255),
    RGB(161/255, 217/255, 155/255),
    RGB(116/255, 196/255, 118/255),
    RGB(65/255, 171/255, 93/255),
    RGB(35/255, 139/255, 69/255),
    RGB(0/255, 109/255, 44/255),
    RGB(0/255, 68/255, 27/255)], :colorbrewer)

register_gradient_colors(:YlGnBu, [
    RGB(255/255, 255/255, 217/255),
    RGB(237/255, 248/255, 177/255),
    RGB(199/255, 233/255, 180/255),
    RGB(127/255, 205/255, 187/255),
    RGB(65/255, 182/255, 196/255),
    RGB(29/255, 145/255, 192/255),
    RGB(34/255, 94/255, 168/255),
    RGB(37/255, 52/255, 148/255),
    RGB(8/255, 29/255, 88/255)], :colorbrewer)

register_gradient_colors(:Purples, [
    RGB(252/255, 251/255, 253/255),
    RGB(239/255, 237/255, 245/255),
    RGB(218/255, 218/255, 235/255),
    RGB(188/255, 189/255, 220/255),
    RGB(158/255, 154/255, 200/255),
    RGB(128/255, 125/255, 186/255),
    RGB(106/255, 81/255, 163/255),
    RGB(84/255, 39/255, 143/255),
    RGB(63/255, 0/255, 125/255)], :colorbrewer)

register_gradient_colors(:GnBu, [
    RGB(247/255, 252/255, 240/255),
    RGB(224/255, 243/255, 219/255),
    RGB(204/255, 235/255, 197/255),
    RGB(168/255, 221/255, 181/255),
    RGB(123/255, 204/255, 196/255),
    RGB(78/255, 179/255, 211/255),
    RGB(43/255, 140/255, 190/255),
    RGB(8/255, 104/255, 172/255),
    RGB(8/255, 64/255, 129/255)], :colorbrewer)

register_gradient_colors(:Greys, [
    RGB(255/255, 255/255, 255/255),
    RGB(240/255, 240/255, 240/255),
    RGB(217/255, 217/255, 217/255),
    RGB(189/255, 189/255, 189/255),
    RGB(150/255, 150/255, 150/255),
    RGB(115/255, 115/255, 115/255),
    RGB(82/255, 82/255, 82/255),
    RGB(37/255, 37/255, 37/255),
    RGB(0/255, 0/255, 0/255)], :colorbrewer)

register_gradient_colors(:YlOrRd, [
    RGB(255/255, 255/255, 204/255),
    RGB(255/255, 237/255, 160/255),
    RGB(254/255, 217/255, 118/255),
    RGB(254/255, 178/255, 76/255),
    RGB(253/255, 141/255, 60/255),
    RGB(252/255, 78/255, 42/255),
    RGB(227/255, 26/255, 28/255),
    RGB(177/255, 0/255, 38/255)], :colorbrewer)

register_gradient_colors(:PuRd, [
    RGB(247/255, 244/255, 249/255),
    RGB(231/255, 225/255, 239/255),
    RGB(212/255, 185/255, 218/255),
    RGB(201/255, 148/255, 199/255),
    RGB(223/255, 101/255, 176/255),
    RGB(231/255, 41/255, 138/255),
    RGB(206/255, 18/255, 86/255),
    RGB(152/255, 0/255, 67/255),
    RGB(103/255, 0/255, 31/255)], :colorbrewer)

register_gradient_colors(:Blues, [
    RGB(247/255, 251/255, 255/255),
    RGB(222/255, 235/255, 247/255),
    RGB(198/255, 219/255, 239/255),
    RGB(158/255, 202/255, 225/255),
    RGB(107/255, 174/255, 214/255),
    RGB(66/255, 146/255, 198/255),
    RGB(33/255, 113/255, 181/255),
    RGB(8/255, 81/255, 156/255),
    RGB(8/255, 48/255, 107/255)], :colorbrewer)

register_gradient_colors(:PuBuGn, [
    RGB(255/255, 247/255, 251/255),
    RGB(236/255, 226/255, 240/255),
    RGB(208/255, 209/255, 230/255),
    RGB(166/255, 189/255, 219/255),
    RGB(103/255, 169/255, 207/255),
    RGB(54/255, 144/255, 192/255),
    RGB(2/255, 129/255, 138/255),
    RGB(1/255, 108/255, 89/255),
    RGB(1/255, 70/255, 54/255)], :colorbrewer)


text:
"""
    mix(f, ϕ,                ds::DataSet)
    mix(f, ϕ, θ::NamedTuple, ds::DataSet)

Compute the mixed `(f°, ϕ°)` from the unlensed field `f` and lensing potential
`ϕ`, given the definition of the mixing matrices in `ds` evaluated at parameters
`θ` (or at fiducial values if no `θ` provided).
"""
mix(f, ϕ, ds::DataSet) = mix(f,ϕ,NamedTuple(),ds)
function mix(f, ϕ, θ::NamedTuple, ds::DataSet)
    @unpack D,G,L = ds(θ)
    L(ϕ)*D*f, G*ϕ
end

"""
    unmix(f°, ϕ°,                ds::DataSet)
    unmix(f°, ϕ°, θ::NamedTuple, ds::DataSet)

Compute the unmixed/unlensed `(f, ϕ)` from the mixed field `f°` and mixed
lensing potential `ϕ°`, given the definition of the mixing matrices in `ds`
evaluated at parameters `θ` (or at fiducial values if no `θ` provided).
"""
unmix(f°, ϕ°, ds::DataSet) = unmix(f°,ϕ°,NamedTuple(),ds)
function unmix(f°, ϕ°, θ::NamedTuple, ds::DataSet)
    @unpack D,G,L = ds(θ)
    ϕ = G\ϕ°
    D\(L(ϕ)\f°), ϕ
end

@doc doc"""
    lnP(t, fₜ, ϕₜ,                ds::DataSet)
    lnP(t, fₜ, ϕₜ, θ::NamedTuple, ds::DataSet)

Compute the log posterior probability in the joint parameterization as a
function of the field, $f_t$, the lensing potential, $\phi_t$, and possibly some
cosmological parameters, $\theta$. The subscript $t$ can refer to either a
"time", e.g. passing `t=0` corresponds to the unlensed parametrization and `t=1`
to the lensed one, or can be `:mix` correpsonding to the mixed parametrization.
In all cases, the arguments `fₜ` and `ϕₜ` should then be $f$ and $\phi$ in
that particular parametrization.

If any parameters $\theta$ are provided, we also include the determinant terms for
covariances which depend on $\theta$. In the mixed parametrization, we
also include any Jacobian determinant terms that depend on $\theta$.

The argument `ds` should be a `DataSet` and stores the masks, data, etc...
needed to construct the posterior.
"""
lnP(t, fₜ, ϕₜ,                ds::DataSet) = lnP(Val(t), fₜ, ϕₜ, NamedTuple(), ds)
lnP(t, fₜ, ϕₜ, θ::NamedTuple, ds::DataSet) = lnP(Val(t), fₜ, ϕₜ, θ,            ds)

function lnP(::Val{t}, fₜ, ϕ, θ::NamedTuple, ds::DataSet) where {t}

    @unpack Cn,Cf,Cϕ,L,M,B,d = ds
    @unpack T = fieldinfo(d)

    f,f̃ = t==0 ? (fₜ, L(ϕ)*fₜ) : (L(ϕ)\fₜ, fₜ)
    Δ = d - M(θ)*B(θ)*f̃ - nonCMB_data_components(θ,ds)
    (
        -T(1/2) * (
            Δ'*pinv(Cn(θ))*Δ + logdet(Cn,θ) +
            f'*pinv(Cf(θ))*f + logdet(Cf,θ) +
            ϕ'*pinv(Cϕ(θ))*ϕ + logdet(Cϕ,θ)
        )
        + T(lnPriorθ(θ,ds))
    )

end

function lnP(::Val{:mix}, f°, ϕ°, θ::NamedTuple, ds::DataSet)
    lnP(Val(0), unmix(f°,ϕ°,θ,ds)..., θ, ds) - logdet(ds.D,θ) - logdet(ds.G,θ)
end

# can be specialized for specific DataSet types:
nonCMB_data_components(θ, ds::DataSet) = 0
lnPriorθ(θ, ds::DataSet) = 0

### marginal posterior gradients

function δlnP_δϕ(
    ϕ, θ, ds;
    previous_state = nothing,
    use_previous_MF = false,
    Nsims = (previous_state==nothing ? 50 : previous_state.Nsims),
    Nbatch = 1,
    weights = :unlensed,
    return_state = false,
    progress = false,
    conjgrad_kwargs = (tol=1e-1,nsteps=500),
    aggressive_gc = fieldinfo(ϕ).Nside>=512
)

    @unpack d,P,M,B,Cn,Cf,Cf̃,Cϕ,Cn̂,G,L = ds(θ)
    Lϕ = L(ϕ)
    Lϕ_batch = Nbatch==1 ? Lϕ : cache(LenseFlow(identity.(batch(ϕ,Nbatch))),identity.(batch(d,Nbatch)))

    if G!=1
        @warn "δlnP_δϕ does not currently handle the G mixing matrix"
    end

    if (previous_state == nothing)
        f_sims = [simulate(batch(Cf,Nbatch)) for i=1:Nsims÷Nbatch]
        n_sims = [simulate(batch(Cn,Nbatch)) for i=1:Nsims÷Nbatch]
        f_wf_sims_guesses = fill(nothing, Nsims÷Nbatch)
        f_wf_guess = nothing
    else
        @unpack f_sims, n_sims, f_wf_sims_guesses, f_wf_guess = previous_state
    end

    # generate simulated datasets for the current ϕ
    ds_sims = map(f_sims, n_sims) do f,n
        resimulate(ds, f̃=Lϕ_batch*f, n=n).ds
    end

    W = (weights == :unlensed) ? 1 : (Cf̃ * pinv(Cf))

    # gradient of the quadratic (in ϕ) piece of the likelihood
    function get_gQD(Lϕ, ds, f_wf_guess)
        f_wf, hist = argmaxf_lnP(
            Lϕ, θ, ds;
            which = :wf,
            guess = (f_wf_guess==nothing ? 0d : f_wf_guess),
            conjgrad_kwargs = (hist=(:i,:res), conjgrad_kwargs...)
        )
        aggressive_gc && gc()
        v = Lϕ' \ (Cf \ f_wf)
        w = W * f_wf
        g = gradient(ϕ -> v' * (Lϕ(ϕ) * w), getϕ(Lϕ))[1]
        @namedtuple(g, f_wf, hist)
    end

    # gQD for the real data
    gQD = get_gQD(Lϕ, ds, f_wf_guess)

    # gQD for several simulated datasets
    if use_previous_MF
        @unpack gQD_sims, ḡ = previous_state
    else
        gQD_sims = pmap(ds_sims, f_wf_sims_guesses) do ds_sim, f_wf_guess
            get_gQD(Lϕ_batch, ds_sim, f_wf_guess)
        end
        ḡ = mean(mapreduce(vcat, gQD_sims) do gQD
            [batchindex(gQD.g,i) for i=1:batchsize(gQD.g)]
        end)
    end

    # final total posterior gradient, including gradient of the prior
    g = gQD.g - ḡ - Cϕ\ϕ

    if return_state
        g, @namedtuple(g, f_sims, n_sims, ḡ, gQD, gQD_sims, f_wf_guess=gQD.f_wf, f_wf_sims_guesses=getindex.(gQD_sims,:f_wf))
    else
        g
    end

end


text:
# Constraints utilities

##################################################
# Constraint definition
"Add constraint to model."
function KN_add_cons(m::Model, ncons::Integer)
    ptr_cons = zeros(Cint, ncons)
    ret = @kn_ccall(add_cons, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), m.env, ncons, ptr_cons)
    _checkraise(ret)
    return ptr_cons
end

function KN_add_con(m::Model)
    ptr_cons = Cint[0]
    ret = @kn_ccall(add_con, Cint, (Ptr{Cvoid}, Ptr{Cint}), m.env, ptr_cons)
    _checkraise(ret)
    return ptr_cons[1]
end

##################################################
# Constraint bounds
##################################################
#--------------------------------------------------
# Equality constraints
#--------------------------------------------------
function KN_set_con_eqbnd(m::Model, indexCons::Integer, bnds::Cdouble)
    ret = @kn_ccall(set_con_eqbnd, Cint, (Ptr{Cvoid}, Cint, Cdouble), m.env, indexCons, bnds)
    _checkraise(ret)
end
KN_set_con_eqbnds(m::Model, indexCons::Integer, bnds::Cdouble) = KN_set_con_eqbnd(m, indexCons, bnds)

function KN_set_con_eqbnds(m::Model, consIndex::Vector{Cint}, eqBounds::Vector{Cdouble})
    ncons = length(consIndex)
    ret = @kn_ccall(set_con_eqbnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, ncons, consIndex, eqBounds)
    _checkraise(ret)
end

function KN_set_con_eqbnds(m::Model, eqBnds::Vector{Cdouble})
    ret = @kn_ccall(set_con_eqbnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}), m.env, eqBnds)
    _checkraise(ret)
end

#--------------------------------------------------
# Inequality constraints
#--------------------------------------------------
# Upper bounds
function KN_set_con_upbnd(m::Model, indexCons::Integer, bnds::Cdouble)
    ret = @kn_ccall(set_con_upbnd, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, indexCons, bnds)
    _checkraise(ret)
end
KN_set_con_upbnds(m::Model, indexCons::Integer, bnds::Cdouble) = KN_set_con_upbnd(m, indexCons, bnds)

function KN_set_con_upbnds(m::Model, upBnds::Vector{Cdouble})
    ret = @kn_ccall(set_con_upbnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}), m.env, upBnds)
    _checkraise(ret)
end

function KN_set_con_upbnds(m::Model, consIndex::Vector{Cint}, upBounds::Vector{Cdouble})
    ncons = length(consIndex)
    ret = @kn_ccall(set_con_upbnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, ncons, consIndex, upBounds)
    _checkraise(ret)
end

# Lower bounds
function KN_set_con_lobnd(m::Model, indexCons::Integer, bnds::Cdouble)
    ret = @kn_ccall(set_con_lobnd, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, indexCons, bnds)
    _checkraise(ret)
end
KN_set_con_lobnds(m::Model, indexCons::Integer, bnds::Cdouble) = KN_set_con_lobnd(m, indexCons, bnds)

# Lower bounds
function KN_set_con_lobnds(m::Model, loBnds::Vector{Cdouble})
    ret = @kn_ccall(set_con_lobnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}), m.env, loBnds)
    _checkraise(ret)
end

function KN_set_con_lobnds(m::Model, consIndex::Vector{Cint}, loBounds::Vector{Cdouble})
    ncons = length(consIndex)
    ret = @kn_ccall(set_con_lobnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, ncons, consIndex, loBounds)
    _checkraise(ret)
end

##################################################
# Getters
##################################################
if KNITRO_VERSION >= v"12.0"
    @define_getters get_con_lobnds
    @define_getters get_con_upbnds
    @define_getters get_con_eqbnds
end

##################################################
# Dual init values
##################################################
function KN_set_con_dual_init_values(m::Model, nindex::Integer, lambdaInitVal::Cdouble)
    ret = @kn_ccall(set_con_dual_init_value, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, lambdaInitVal)
    _checkraise(ret)
end

function KN_set_con_dual_init_values(m::Model, indexCon::Vector{Cint}, lambdaInitVals::Vector{Cdouble})
    nvar = length(indexCon)
    @assert nvar == length(lambdaInitVals)
    ret = @kn_ccall(set_con_dual_init_values, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, indexCon, lambdaInitVals)
    _checkraise(ret)
end

function KN_set_con_dual_init_values(m::Model, lambdaInitVals::Vector{Cdouble})
    ret = @kn_ccall(set_con_dual_init_values_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, lambdaInitVals)
    _checkraise(ret)
end

##################################################
# Constraint scalings
##################################################
"""
Set an array of constraint scaling values to perform a scaling
  cScaled[i] = cScaleFactors[i] * c[i]
for each constraint. These scaling factors should try to
represent the "typical" values of the inverse of the constraint
values "c" so that the scaled constraints ("cScaled") used
internally by Knitro are close to one.  Scaling factors for
standard constraints can be provided with "cScaleFactors", while
scalings for complementarity constraints can be specified with
"ccScaleFactors".  The values for cScaleFactors/ccScaleFactors
should be positive.  If a non-positive value is specified, that
constraint will use either the standard Knitro scaling
(KN_SCALE_USER_INTERNAL), or no scaling (KN_SCALE_USER_NONE).

"""
function KN_set_con_scalings(m::Model, nindex::Integer, cScaleFactors::Cdouble)
    ret = @kn_ccall(set_con_scaling, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, cScaleFactors)
    _checkraise(ret)
end

function KN_set_con_scalings(m::Model, indexCon::Vector{Cint}, cScaleFactors::Vector{Cdouble})
    nvar = length(indexCon)
    ret = @kn_ccall(set_con_scalings, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, indexCon, cScaleFactors)
    _checkraise(ret)
end

function KN_set_con_scalings(m::Model, cScaleFactors::Vector{Cdouble})
    ret = @kn_ccall(set_con_scalings_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, cScaleFactors)
    _checkraise(ret)
end

##################################################
# Constraints constants
##################################################
function KN_add_con_constants(m::Model, indexCons::Vector{Cint}, constants::Vector{Cdouble})
    nnc = length(constants)
    @assert length(indexCons) == length(constant)
    ret = @kn_ccall(add_con_constants, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnc,
                    indexCons,
                    constants)
    _checkraise(ret)
end

function KN_add_con_constants(m::Model, constants::Vector{Cdouble})
    nnc = length(constants)
    ret = @kn_ccall(add_con_constants_all, Cint,
                    (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, constants)
    _checkraise(ret)
end

function KN_add_con_constant(m::Model, indexCon::Integer, constant::Cdouble)
    ret = @kn_ccall(add_con_constant, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, indexCon, constant)
    _checkraise(ret)
end

##################################################
# Constraint structure
##################################################
#------------------------------
# add structure of linear constraint
#------------------------------
"""
Add linear structure to the constraint unctions.
Each component i of arrays indexCons, indexVars and coefs adds a linear
term:
   coefs[i]*x[indexVars[i]]
to constraint c[indexCons[i]].

"""
function KN_add_con_linear_struct(m::Model,
                                  jacIndexCons::Vector{Cint},
                                  jacIndexVars::Vector{Cint},
                                  jacCoefs::Vector{Cdouble})
    # get number of constraints
    nnz = length(jacIndexCons)
    @assert nnz == length(jacIndexVars) == length(jacCoefs)
    ret = @kn_ccall(add_con_linear_struct,
                    Cint,
                    (Ptr{Cvoid}, KNLONG, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    jacIndexCons,
                    jacIndexVars,
                    jacCoefs)
    _checkraise(ret)
end

function KN_add_con_linear_struct(m::Model,
                                  indexCon::Integer,
                                  indexVar::Vector{Cint},
                                  coefs::Vector{Cdouble})
    # get number of constraints
    nnz = length(indexVar)
    @assert nnz == length(coefs)
    ret = @kn_ccall(add_con_linear_struct_one,
                    Cint,
                    (Ptr{Cvoid}, KNLONG, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    indexCon,
                    indexVar,
                    coefs)
    _checkraise(ret)
end
KN_add_con_linear_struct(m::Model, indexCon::Integer, indexVar::Integer, coef::Cdouble) =
    KN_add_con_linear_struct(m, indexCon, Cint[indexVar], [coef])

#------------------------------
# add constraint quadratic structure
#------------------------------
"""
Add quadratic structure to the constraint functions.
Each component i of arrays indexCons, indexVars1, indexVars2 and coefs adds a
quadratic term:
   coefs[i]*x[indexVars1[i]]*x[indexVars2[i]]
to the constraint c[indexCons[i]].

"""
function KN_add_con_quadratic_struct(m::Model,
                                     indexCons::Vector{Cint},
                                     indexVars1::Vector{Cint},
                                     indexVars2::Vector{Cint},
                                     coefs::Vector{Cdouble})
    # get number of constraints
    nnz = length(indexVars1)
    @assert nnz == length(indexCons) == length(indexVars2) == length(coefs)
    ret = @kn_ccall(add_con_quadratic_struct,
                    Cint,
                    (Ptr{Cvoid}, KNLONG, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    indexCons,
                    indexVars1,
                    indexVars2,
                    coefs)
    _checkraise(ret)
end

function KN_add_con_quadratic_struct(m::Model,
                                     indexCons::Integer,
                                     indexVars1::Vector{Cint},
                                     indexVars2::Vector{Cint},
                                     coefs::Vector{Cdouble})
    # get number of constraints
    nnz = length(indexVars1)
    @assert nnz == length(indexVars2) == length(coefs)
    ret = @kn_ccall(add_con_quadratic_struct_one,
                    Cint,
                    (Ptr{Cvoid}, KNLONG, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
                    m.env,
                    nnz,
                    indexCons,
                    indexVars1,
                    indexVars2,
                    coefs)
    _checkraise(ret)
end
KN_add_con_quadratic_struct(m, indexCons::Integer, indexVar1::Integer, indexVar2::Integer, coef::Cdouble) =
KN_add_con_quadratic_struct(m, indexCons, Cint[indexVar1], Cint[indexVar2], Cdouble[coef])

#------------------------------
# Conic structure
#------------------------------
"""
Add L2 norm structure of the form ||Ax + b||_2 to a constraint.
  indexCon:    The constraint index that the L2 norm term will be added to.
  nCoords:     The number of rows in "A" (or dimension of "b")
  nnz:         The number of sparse non-zero elements in "A"
  indexCoords: The coordinate (row) index for each non-zero element in "A".
  indexVars:   The variable (column) index for each non-zero element in "A"
  coefs:       The coefficient value for each non-zero element in "A"
  constants:   The array "b" - may be set to NULL to ignore "b"

#Note
L2 norm structure can currently only be added to constraints that
otherwise only have linear (or constant) structure.  In this way
they can be used to define conic constraints of the form
||Ax + b|| <= c'x + d.  The "c" coefficients should be added through
"KN_add_con_linear_struct()" and "d" can be set as a constraint bound
or through "KN_add_con_constants()".

#Note
Models with L2 norm structure are currently only handled by the
Interior/Direct (KN_ALG_BAR_DIRECT) algorithm in Knitro.  Any model
with structure defined with KN_add_L2norm() will automatically be
forced to use this algorithm.

"""
function KN_add_con_L2norm(m::Model, indexCon::Integer, nCoords::Integer, nnz::Integer,
                       indexCoords::Vector{Cint}, indexVars::Vector{Cint},
                       coefs::Vector{Cdouble}, constants::Vector{Cdouble})
    @assert length(coefs) == length(indexVars) == length(indexCoords) == nnz
    ret = @kn_ccall(add_con_L2norm,
                    Cint,
                    (Ptr{Cvoid}, Cint, Cint, KNLONG, Ptr{Cint}, Ptr{Cint},
                     Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, indexCon, nCoords, nnz, indexCoords,
                    indexVars, coefs, constants)
    _checkraise(ret)
end

##################################################
# Complementary constraints
##################################################
"""
This function adds complementarity constraints to the problem.
The two lists are of equal length, and contain matching pairs of
variable indices.  Each pair defines a complementarity constraint
between the two variables.  The function can only be called once.
The array "ccTypes" specifies the type of complementarity:
   KN_CCTYPE_VARVAR: two (non-negative) variables
   KN_CCTYPE_VARCON: a variable and a constraint
   KN_CCTYPE_CONCON: two constraints

#Note
Currently only KN_CCTYPE_VARVAR is supported.  The other
"ccTypes" will be added in future releases.

"""
function KN_set_compcons(m::Model,
                         ccTypes::Vector{Cint},
                         indexComps1::Vector{Cint},
                         indexComps2::Vector{Cint})
    # get number of constraints
    nnc = length(ccTypes)
    @assert nnc == length(indexComps1) == length(indexComps2)
    ret = @kn_ccall(set_compcons,
                    Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                    m.env,
                    nnc,
                    ccTypes,
                    indexComps1,
                    indexComps2)
    _checkraise(ret)
end

##################################################
# Complementary constraint scalings
##################################################

function KN_set_compcon_scalings(m::Model, nindex::Integer, cScaleFactors::Cdouble)
    ret = @kn_ccall(set_compcon_scaling, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, cScaleFactors)
    _checkraise(ret)
end

function KN_set_compcon_scalings(m::Model, indexCompCon::Vector{Cint}, cScaleFactors::Vector{Cdouble})
    nvar = length(indexCompCon)
    ret = @kn_ccall(set_compcon_scalings, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, indexCompCon, cScaleFactors)
    _checkraise(ret)
end

function KN_set_compcon_scalings(m::Model, cScaleFactors::Vector{Cdouble})
    ret = @kn_ccall(set_compcon_scalings_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, cScaleFactors)
    _checkraise(ret)
end

##################################################
## Naming constraints
##################################################
function KN_set_con_names(m::Model, nindex::Integer, name::String)
    ret = @kn_ccall(set_con_name, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cchar}),
                    m.env, nindex, name)
    _checkraise(ret)
end

function KN_set_con_names(m::Model, conIndex::Vector{Cint}, names::Vector{String})
    ncon = length(conIndex)
    ret = @kn_ccall(set_con_names, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Ptr{Char}}),
                    m.env, ncon, conIndex, names)
    _checkraise(ret)
end

function KN_set_con_names(m::Model, names::Vector{String})
    ret = @kn_ccall(set_con_names_all, Cint, (Ptr{Cvoid}, Ptr{Ptr{Cchar}}),
                    m.env, names)
    _checkraise(ret)
end

function KN_set_compcon_names(m::Model, nindex::Integer, name::String)
    ret = @kn_ccall(set_compcon_name, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cchar}),
                    m.env, nindex, name)
    _checkraise(ret)
end

function KN_set_compcon_names(m::Model, conIndex::Vector{Cint}, names::Vector{String})
    ncon = length(conIndex)
    ret = @kn_ccall(set_compcon_names, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Ptr{Char}}),
                    m.env, ncon, conIndex, names)
    _checkraise(ret)
end

function KN_set_compcon_names(m::Model, names::Vector{String})
    ret = @kn_ccall(set_compcon_names_all, Cint, (Ptr{Cvoid}, Ptr{Ptr{Cchar}}),
                    m.env, names)
    _checkraise(ret)
end

# Getters
if KNITRO_VERSION >= v"12.0"
    function KN_get_con_names(m::Model, max_length=1024)
        return String[KN_get_con_names(m, Cint(id-1), max_length) for id in 1:KN_get_number_cons(m)]
    end

    function KN_get_con_names(m::Model, index::Vector{Cint}, max_length=1024)
        return String[KN_get_con_names(m, id, max_length) for id in index]
    end

    function KN_get_con_names(m::Model, index::Cint, max_length=1024)
        rawname = zeros(Cchar, max_length)
        ret = @kn_ccall(get_con_name, Cint,
                        (Ptr{Cvoid}, Cint, Cint, Ptr{Cchar}),
                        m.env, index, max_length, rawname)
        _checkraise(ret)
        name = String(strip(String(convert(Vector{UInt8}, rawname)), '\0'))
        return name
    end
end

##################################################
## Feasibility tolerance
##################################################
function KN_set_con_feastols(m::Model, nindex::Integer, cFeasTol::Cdouble)
    ret = @kn_ccall(set_con_feastol, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, cFeasTol)
    _checkraise(ret)
end

function KN_set_con_feastols(m::Model, cIndex::Vector{Cint}, cFeasTols::Vector{Cdouble})
    ncon = length(cIndex)
    @assert length(cFeasTols) == ncon
    ret = @kn_ccall(set_con_feastols, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, ncon, cIndex, cFeasTols)
    _checkraise(ret)
end

function KN_set_con_feastols(m::Model, cFeasTols::Vector{Cdouble})
    ret = @kn_ccall(set_con_feastols_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, cFeasTols)
    _checkraise(ret)
end

function KN_set_compcon_feastols(m::Model, nindex::Integer, cFeasTol::Cdouble)
    ret = @kn_ccall(set_compcon_feastol, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, cFeasTol)
    _checkraise(ret)
end

function KN_set_compcon_feastols(m::Model, cIndex::Vector{Cint}, cFeasTols::Vector{Cdouble})
    ncon = length(cIndex)
    @assert length(cFeasTols) == ncon
    ret = @kn_ccall(set_compcon_feastols, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, ncon, cIndex, cFeasTols)
    _checkraise(ret)
end

function KN_set_compcon_feastols(m::Model, cFeasTols::Vector{Cdouble})
    ret = @kn_ccall(set_compcon_feastols_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, cFeasTols)
    _checkraise(ret)
end

##################################################
## Constraint property
##################################################
"""
Specify some properties of the objective and constraint functions.

# Note
Use bit-wise specification of the features:
bit value   meaning
  0     1   KN_OBJ_CONVEX/KN_CON_CONVEX
  1     2   KN_OBJ_CONCAVE/KN_CON_CONCAVE
  2     4   KN_OBJ_CONTINUOUS/KN_CON_CONTINUOUS
  3     8   KN_OBJ_DIFFERENTIABLE/KN_CON_DIFFERENTIABLE
  4    16   KN_OBJ_TWICE_DIFFERENTIABLE/KN_CON_TWICE_DIFFERENTIABLE
  5    32   KN_OBJ_NOISY/KN_CON_NOISY
  6    64   KN_OBJ_NONDETERMINISTIC/KN_CON_NONDETERMINISTIC
default = 28 (bits 2-4 enabled: e.g. continuous, differentiable, twice-differentiable)

"""
function KN_set_con_properties(m::Model, nindex::Integer, cProperty::Cint)
    ret = @kn_ccall(set_con_property, Cint, (Ptr{Cvoid}, Cint, Cint),
                    m.env, nindex, cProperty)
    _checkraise(ret)
end

function KN_set_con_properties(m::Model, cIndex::Vector{Cint}, cProperties::Vector{Cint})
    ncon = length(cIndex)
    @assert length(cProperties) == ncon
    ret = @kn_ccall(set_con_properties, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}),
                    m.env, ncon, cIndex, cProperties)
    _checkraise(ret)
end

function KN_set_con_properties(m::Model, cProperties::Vector{Cint})
    ret = @kn_ccall(set_con_properties_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, cProperties)
    _checkraise(ret)
end


text:
function perceptron(seed::Integer)
    rng = StableRNG(seed)
    train = gettrain()

    small = filter(e -> e.label == :yes || e.label == :no, train)
    labels = [:yes, :no]
    train, val = splitdata(small, 0.1, rng)
    x_train = convert.(Float32, reduce(hcat, rpad.(getfield.(train, :data), 16000, 0)) ./ 2^15) |> gpu
    y_train = convert.(Float32, onehotbatch(getfield.(train, :label), labels)) |> gpu
    x_val = convert.(Float32, reduce(hcat, rpad.(getfield.(val, :data), 16000, 0)) ./ 2^15) |> gpu
    y_val = convert.(Float32, onehotbatch(getfield.(val, :label), labels)) |> gpu

    model = Chain(
        Dense(16000 => 800, relu, init=glorot_uniform(rng)),
        Dense(800 => 40, relu, init=glorot_uniform(rng)),
        Dense(40 => 2, init=glorot_uniform(rng)),
    ) |> gpu

    @epochs 100 begin
        train_loss = 0
        train_accuracy = 0

        train!(
            t -> begin
                out = model(t.data)
                loss = logitcrossentropy(out, t.label)
                train_loss += loss * size(t.data)[end] / length(train)
                pred = onecold(out, labels)
                corr = onecold(t.label, labels)
                train_accuracy += count(pred .== corr) / length(train)
                loss
            end,
            params(model),
            DataLoader((data=x_train, label=y_train), batchsize=32, shuffle=true, rng=rng),
            ADAM(),
        )

        out = model(x_val)
        val_loss = logitcrossentropy(out, y_val)
        pred = onecold(out, labels)
        corr = onecold(y_val, labels)
        val_accuracy = count(pred .== corr) / length(val)
        @info "Epoch training loss: $(train_loss)"
        @info "Epoch training accuracy: $(train_accuracy)"
        @info "Current validation loss: $(val_loss)"
        @info "Current validation accuracy: $(val_accuracy)"
    end
end

max_stars_repo_name: kpmrozowski/wut-mini-deep-learning

max_issues_repo_name: kpmrozowski/wut-mini-deep-learning

max_forks_repo_name: kpmrozowski/wut-mini-deep-learning


text:
struct RatioMatchBuffer{T3,M,T4}
    ΔE::T3      # q × N × M
    A::T3       # q × N × M
    A₁::T3      # 1 × N × M     (for J & h)
    A₁₂::T3     # 1 × 1 × M     (for W)
    K::M        # qN × qN       (for J)
    B::M        # P × M         (for W)
    B₁::M       # 1 × M         (for W)
    B̃::T4       # P × q × N × M (for W)
    B̃₁::T4      # 1 × q × N × M (for W)
    B̃₂::T4      # P × 1 × N × M (for W)
    B̃₄::T4      # P × q × N × 1 (for W)
    function RatioMatchBuffer(θ::SpinModel, M::Int)
        (;q, N) = θ
        has_h, has_J, has_W = (hasproperty(θ, s) for s ∈ (:h, :J, :W))
        X = has_h ? θ.h : has_J ? θ.J : θ.W
        ΔE, A = ntuple(i -> similar(X, q,N,M), 2)
        A₁ = (has_h | has_J) ? similar(X, 1,N,M) : similar(X, 0,0,0)
        A₁₂= has_W ? similar(X, 1,1,M) : similar(X, 0,0,0)
        K  = has_J ? similar(X, q*N,q*N) : similar(X, 0,0)
        B  = has_W ? similar(X, θ.P,M) : similar(X, 0,0)
        B₁ = has_W ? similar(X, 1,M) : similar(X, 0,0)
        B̃  = has_W ? similar(X, θ.P,q,N,M) : similar(X, 0,0,0,0)
        B̃₁ = has_W ? similar(X, 1,q,N,M) : similar(X, 0,0,0,0)
        B̃₂ = has_W ? similar(X, θ.P,1,N,M) : similar(X, 0,0,0,0)
        B̃₄ = has_W ? similar(X, θ.P,q,N,1) : similar(X, 0,0,0,0)
        return new{typeof(A), typeof(K), typeof(B̃)}(
            ΔE, A, A₁, A₁₂, K, B, B₁, B̃, B̃₁, B̃₂, B̃₄
        )
    end
    RatioMatchBuffer(z::AbstractArray{<:Any,3}, θ::SpinModel) = RatioMatchBuffer(θ, size(z,3))
end

"""
    ratiomatch(z, θ::SpinModel)

Return ratio matching loss (Hyvärinen, Comput Stat Data Anal 2007)

Computed as, for a system of `N` `q`-state spins,

    mean(abs2(1 - sigmoid( E(zᵠᵏ) - E(z) )) for ϕ ∈ 1:q, k ∈ 1:N, z ∈ Data) - 1//4

where `zᵠᵏ` is the same as `z` but with spin `k` in state `ϕ`.
This loss vanishes for unstructured models, in which `E(zᵠᵏ) - E(z) = 0` for all `ϕ, k, x`.
"""
ratiomatch(z, θ::SpinModel; energyscale = 1) = ratiomatch(z, θ, RatioMatchBuffer(z, θ); energyscale = energyscale)
ratiomatch(z, θ::SpinModel, buffer; energyscale = 1) = ratiomatch(z, θ, θ, buffer, false, energyscale)
ratiomatch(z, θ::T, θ̄::T) where T<:SpinModel = ratiomatch(z, θ, θ̄, RatioMatchBuffer(z, θ))
function ratiomatch(z, θ::T, θ̄::T, buffer, gradient = true, energyscale = 1) where T<:SpinModel
    (; ΔE, A, A₁, A₁₂, K, B, B₁, B̃, B̃₁, B̃₂, B̃₄) = buffer
    fill!(ΔE, zero(eltype(ΔE))) # initialize ΔE
    q, N, M = size(z)

    ## First we compute:    ΔE[ϕ,k,m] = E(zᵠᵏ) - E(z)
    if hasproperty(θ, :J)
        J = θ.J
        K .= reshape(J, q*N, q*N)
        K.+= transpose(reshape(J, q*N, q*N))
        mul!(reshape(ΔE, :,M), K, reshape(z, :,M), -1, 1)
        @. A = ΔE * z
        ΔE .-= sum!(A₁, A)
    end

    if hasproperty(θ, :h)
        h = θ.h
        @. A = h * z
        sum!(A₁, A)
        @. ΔE += -(h - A₁)
    end

    if hasproperty(θ, :W) && hasproperty(θ, :b)
        (;W, b) = θ
        mul!(B, reshape(W, :, q*N), reshape(z, :, M))   # B[a,m] = ∑ᵨₖ W[a,ρ,k] * z[ρ,k,m]
        B.+= b                                          # B[a,m] = b[a] + ∑ᵨₖ W[a,ρ,k] * z[ρ,k,m]
        B̃ .= W .* reshape(z, 1,q,N,M)   # B̃[a,ρ,k,m] = W[a,ρ,k] * z[ρ,k,m]
        sum!(B̃₂, B̃)                     # B̃₂[a,.,k,m] = ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        @. B̃ = W - B̃₂                   # B̃[a,σ,k,m] = W[a,σ,k] - ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        B̃.+= reshape(B, :,1,1,M)        # B̃[a,σ,k,m] = B[a,m] + W[a,σ,k] - ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        ΔE.+= θ.rbm_wt .* reshape(sum!(logcosh, B₁, B), 1,1,M)
        ΔE.-= θ.rbm_wt .* reshape(sum!(logcosh, B̃₁, B̃), q,N,M)
    end

    isone(energyscale) || rmul!(ΔE, convert(eltype(ΔE), energyscale))

    ## Loss ∝ ∑ᵩₖₘ (1-sigmoid(E(zᵠᵏ)-E(z)))² - ∑ᵩₖₘ (zᵠᵏ==z) * (1/4)
    ## We subtract off the terms for which zᵠᵏ==z
    f = ΔE                      # non-allocating alias
    @. f = sigmoid(f)           # f[σ,k,m] = sigmoid(ΔE[σ,k,m])
    L = sum(x -> abs2(1-x), f)
    # subtract loss from zᵠᵏ=z (N terms per sample) for which ΔE = 0 and (1-σ(ΔE))² = 1/4
    L-= 1//4 * N * M
    # take the average over M samples and N*(q-1) neighboring states per sample
    L*= 1//(M * N * (q-1))
    # multiply by 4 so that unstructured model (ΔE = 0 for all zᵠᵏ, z) gives L = 1
    L*= 4
    gradient || return L

    #### Gradient calculations start here
    @. f = abs2(1-f) * 2f       # f[σ,k,m] = (1-sigmoid(ΔE[σ,k,m]))² * 2sigmoid(ΔE[σ,k,m])
    A.= z .* sum!(A₁, f)        # A[σ,k,m] = z[σ,k,m] * ∑ᵨ f[ρ,k,m]
    A.= f .- A                  # A[σ,k,m] = f[σ,k,m] - z[σ,k,m] * ∑ᵨ f[ρ,k,m]

    if hasproperty(θ̄, :h)
        h̄ = θ̄.h
        sum!(h̄, A)
        h̄ .*= 4//(M * N * (q-1))
    end

    if hasproperty(θ̄, :J)
        J̄ = θ̄.J
        Ĵ = reshape(J̄, q*N, q*N)
        mul!(Ĵ, reshape(z,:,M), transpose(reshape(A,:,M)))
        Ĵ .+= transpose(Ĵ)
        dropselfcoupling!(J̄)    # donot update self coupling
        J̄ .*= 4//(M * N * (q-1))
    end

    if hasproperty(θ̄, :W) && hasproperty(θ̄, :b)
        W̄, b̄ = θ̄.W, θ̄.b
        @. B = tanh(B)  # B[a,m] = tanh(b[a] + ∑ᵨₖ W[a,ρ,k]*z[ρ,k,m])
        @. B̃ = tanh(B̃)  # B̃[a,σ,k,m] = tanh(B[a,m] + W[a,σ,k] - ∑ᵨ W[a,ρ,k]*z[ρ,k,m])

        sum!(A₁₂, f)
        mul!(b̄, B, vec(A₁₂), -1, 0)
        mul!(b̄, reshape(B̃,:,q*N*M), vec(f), 1, 1)
        b̄ .*= 4//(M * N * (q-1)) .* θ.rbm_wt

        Bf̄ = B  # non-allocating alias
        Bf̄ .*= reshape(A₁₂, :,M)        # Bf̄[a,m] = B[a,m] * ∑ᵨₖ f[ρ,k,m]
        B̃f = B̃  # non-allocating alias
        B̃f .*= reshape(f, 1,q,N,M)      # B̃f[a,σ,k,m] = B̃[a,σ,k,m] * f[σ,k,m]
        sum!(W̄, B̃f)
        mul!(reshape(W̄, :,q*N), Bf̄, transpose(reshape(z, q*N,M)), -1, 1)

        B̃xf = B  # non-allocating alias
        sum!(reshape(B̃xf, :,1,1,M), B̃f) # B̃xf[a,m] = ∑ᵨₖ B̃[a,ρ,k,m] * f[ρ,k,m]
        mul!(reshape(W̄, :,q*N), B̃xf, transpose(reshape(z, q*N,M)), 1, 1)

        B̃ .= .-sum!(B̃₂, B̃f) .* reshape(z,1,q,N,M)

        W̄ .+= sum!(B̃₄, B̃)
        W̄ .*= 4//(M * N * (q-1)) .* θ.rbm_wt
    end

    return L
end


text:
"""
LSIAF{F<:AbstractFloat}(Y::Array{F})

The function `LSIAF` performs the filtering via coordinate descent and takes only following input:

* y : Signal or time-series to be filtered (a vector of reals).
"""
function LSIAF{F<:AbstractFloat}(Y::Array{F})
    ## File Input
    Y0 = Y
    Ym = mean(Y)
    Y  = Y0 - Ym
    n  = length(Y)

    ## Timming
    tic()

    ## Adaptive Least Absolute Shrinkage and Selection Algorithm
    β_lasso, ~, X = adaLASSO(Y)

    ## Post-Lasso
    # Determine Selected Variables
    selection = Int[]
    for i = 1:size(β_lasso,1)
        if abs(β_lasso[i,:])[1] > 1e-5
            push!(selection, i)
        end
    end

    #on error
    # if sum(selection == (1))==0
    #     selection = [selection 1]
    # end
    push!(selection, 1)

    # OLS to remove LASSO BIAS
    β_ols = inv(X[:,selection]' * X[:,selection]) * X[:,selection]' * Y

    ## Output
    OutputMatrix = Array{Number}(size(selection,1), 2)
    OutputMatrix[:, 1] = selection'
    OutputMatrix[:, 2] = β_ols

    Filter = X[:, selection] * β_ols + Ym

    ## Timming
    timeAdaLASSO = toq();
    error = sum( (Filter - Y0) .^ 2 )
    println("Finished LASSO with LMS approximation error of ", error, " and elapsed time of ", timeAdaLASSO, " seconds.")

    return OutputMatrix, Filter, β_lasso
end


text:
using MeteomaticsAPI
using Dates

# We can obtain a time-series of surface temperature and pressure for two coordinates at 1.5 hourly
# resolution over the next 14 days:

df_ts = querytimeseries(
    [(0, 50), (50, 0)], Dates.now(), Dates.now() + Day(14),
    Hour(1)+Minute(30), ["t_2m:C", "msl_pressure:hPa"]
)

# Here a CompoundPeriod was used as the interval (4th argument) to demonstrate how one can in general
# be created in Julia: the connector accepts these as long as constituent parts are of types Day,
# Hour, Minute and/or Second. Start and end dates were provided as Dates.Date objects. Alternatively,
# users can provide well-formed strings.

# We can obtain gridded data for the same variables. Be sure to provide as many variables into which
# to unpack the requests as you have parameters:

df_grid_t, df_grid_p = querygrid(Dates.now(), ["t_2m:C", "msl_pressure:hPa"], 50, 0, 0, 50, 1, 1)

# The login credentials required to run this query are contained in Config.jl. You can overwrite these,
# or provide your own credentials as optional arguments. All additional optional arguments described at
# https://www.meteomatics.com/en/api/request/optional-parameters/ are also available.

# Two types of metadata requests can be made, but require an API subscription to work properly
# (metadata requests don't work for 'mix' models; other models are not available without subscription)
df_tr = querytimeranges(
    ["t_2m:C", "msl_pressure:hPa"], "ecmwf-vareps"
)

df_it = queryinitdatetime(
    Dates.now(), Dates.now() + Day(14),
    Hour(1)+Minute(30), ["t_2m:C", "msl_pressure:hPa"], "ecmwf-vareps"
)

max_stars_repo_name: meteomatics/julia-connector-api

max_issues_repo_name: meteomatics/julia-connector-api

max_forks_repo_name: meteomatics/julia-connector-api


text:
# The ids are appended to each row inside []`. Replace with comma separation.

dataset_path = joinpath(dirname(@__FILE__),
                        "thyroid-disease-thyroid-0387.data.orig")

f = open(dataset_path)
data = readlines(f)
close(f)

f = open("$dataset_path.custom", "w")
for line in data
  line = replace(line, "[", ",")
  line = replace(line, "]", "")
  write(f, line)
end
close(f)

data = readstring(dataset_path)
write("$dataset_path.custom", replace(replace(data, "[", ","), "]", ""))

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
# ---
# title: 1027. Longest Arithmetic Subsequence
# id: problem1027
# author: Indigo
# date: 2022-04-14
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/longest-arithmetic-subsequence/description/>
# hidden: true
# ---
#
# Given an array `A` of integers, return the **length** of the longest
# arithmetic subsequence in `A`.
#
# Recall that a _subsequence_ of `A` is a list `A[i_1], A[i_2], ..., A[i_k]`
# with `0 <= i_1 < i_2 < ... < i_k <= A.length - 1`, and that a sequence `B` is
# _arithmetic_ if `B[i+1] - B[i]` are all the same value (for `0 <= i < B.length
# - 1`).
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [3,6,9,12]
#     Output: 4
#     Explanation:
#     The whole array is an arithmetic sequence with steps of length = 3.
#
#
# **Example 2:**
#
#
#
#     Input: A = [9,4,7,2,10]
#     Output: 3
#     Explanation:
#     The longest arithmetic subsequence is [4,7,10].
#
#
# **Example 3:**
#
#
#
#     Input: A = [20,1,15,3,10,5,8]
#     Output: 4
#     Explanation:
#     The longest arithmetic subsequence is [20,15,10,5].
#
#
#
#
# **Constraints:**
#
#   * `2 <= A.length <= 1000`
#   * `0 <= A[i] <= 500`
#
#
## @lc code=start
using LeetCode

function longest_arith_seq_length(nums::Vector{Int})
    len = length(nums)
    len == 1 && return 0
    dp = fill(0, len, 1001)
    res = 0
    for i in 1:len, j in 1:(i - 1)
        d = nums[i] - nums[j] + 501
        dp[i, d] = max(dp[i, d], dp[j, d] + 1)
        res = max(res, dp[i, d])
    end
    return res + 1
end
## @lc code=end


text:
import ..UncertainValues:
	AbstractUncertainScalarKDE,
	TheoreticalDistributionScalarValue,
	UncertainScalarTheoreticalFit

import Statistics

# The general strategy is to resample

"""
	mean(uv::AbstractUncertainValue, n::Int)

Compute the mean of an uncertain value over an `n`-draw sample of it.
"""
Statistics.mean(uv::AbstractUncertainValue, n::Int) = mean(resample(uv, n))

"""
	mode(uv::AbstractUncertainValue, n::Int)

Compute the mode of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.mode(uv::AbstractUncertainValue, n::Int) = mode(resample(uv, n))

"""
	median(uv::AbstractUncertainValue, n::Int)

Compute the median of an uncertain value over an `n`-draw sample of it.
"""
Statistics.median(uv::AbstractUncertainValue, n::Int) = median(resample(uv, n))

"""
	middle(uv::AbstractUncertainValue, n::Int)

Compute the middle of an uncertain value over an `n`-draw sample of it.
"""
Statistics.middle(uv::AbstractUncertainValue, n::Int) = middle(resample(uv, n))

"""
	quantile(uv::AbstractUncertainValue, q, n::Int)

Compute the quantile(s) `q` of an uncertain value over an `n`-draw sample of it.
"""
Statistics.quantile(uv::AbstractUncertainValue, q, n::Int) = quantile(resample(uv, n), q)

"""
	std(uv::AbstractUncertainValue, n::Int)

Compute the standard deviation of an uncertain value over an `n`-draw sample of it.
"""
Statistics.std(uv::AbstractUncertainValue, n::Int) = std(resample(uv, n))

"""
	variance(uv::AbstractUncertainValue, n::Int)

Compute the variance of an uncertain value over an `n`-draw sample of it.
"""
Statistics.var(uv::AbstractUncertainValue, n::Int) = var(resample(uv, n))

# But for theoretical distributions, we may directly access properties of the furnishing
# distributions

##########################################################################################
# Uncertain values represented by theoretical distributions.
# Base stats can be estimated directly from the distributions, so no need for resampling.
# These methods definitions are just for compatibility with other uncertain values, which
# may not have analytically determined statistics associated with them.
##########################################################################################
StatsBase.mode(uv::TheoreticalDistributionScalarValue, n::Int) = mode(uv)
Statistics.mean(uv::TheoreticalDistributionScalarValue, n::Int) = mean(uv)
Statistics.median(uv::TheoreticalDistributionScalarValue, n::Int) = median(uv)
Statistics.middle(uv::TheoreticalDistributionScalarValue, n::Int) = middle(uv)
Statistics.quantile(uv::TheoreticalDistributionScalarValue, q, n::Int) = quantile(uv, q)
Statistics.std(uv::TheoreticalDistributionScalarValue, n::Int) = std(uv)
Statistics.var(uv::TheoreticalDistributionScalarValue, n::Int) = var(uv)

StatsBase.mode(uv::UncertainScalarTheoreticalFit, n::Int) = mode(uv)
Statistics.mean(uv::UncertainScalarTheoreticalFit, n::Int) = mean(uv)
Statistics.median(uv::UncertainScalarTheoreticalFit, n::Int) = median(uv)
Statistics.middle(uv::UncertainScalarTheoreticalFit, n::Int) = middle(uv)
Statistics.quantile(uv::UncertainScalarTheoreticalFit, q, n::Int) = quantile(uv, q)
Statistics.std(uv::UncertainScalarTheoreticalFit, n::Int) = std(uv)
Statistics.var(uv::UncertainScalarTheoreticalFit, n::Int) = var(uv)


text:
module TestShowEntries

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running getEntries function tests")

######################
## test getEntries ##
######################

## create test case
df1 = DataFrame()
df2 = DataFrame()
df1[:prices1] = @data([100, 120, 140, 170, 200])
df1[:prices2] = @data([110, 120, NA, 130, 150])
df2[:prices1] = @data([110, 120, NA, 130, 150])
df2[:prices2] = @data([100, 120, 140, 170, 200])
dats = [Date(2010, 1, 1):Date(2010, 1, 5)]
testcase1 = TimeData.Timenum(df1, dats)
testcase2 = TimeData.Timenum(df2, dats)

## expected results for dates ordering
idxVals = Date[Date("2010-01-03"),
               Date("2010-01-04"),
               Date("2010-01-05"),
               Date("2010-01-05")]
variables = [:prices1, :prices1, :prices1, :prices2]
values = Any[140, 170, 200, 150]

df = DataFrame(variable = variables, value = values)
expTd = TimeData.Timedata(df, idxVals)

## entries large than 130
actOut = TimeData.chkElw(x -> x.>130,
                     TimeData.eachentry(testcase1)) |>
     x -> TimeData.asArr(x, Bool, false) |>
     x -> TimeData.getEntries(testcase1, x)

@test isequal(expTd, actOut)

## expected results for second test case with variable ordering
idxVals = Date[Date("2010-01-05"),
               Date("2010-01-03"),
               Date("2010-01-04"),
               Date("2010-01-05")]
variables = [:prices1, :prices2, :prices2, :prices2]
values = Any[150, 140, 170, 200]

df = DataFrame(variable = variables, value = values)
expTd2 = TimeData.Timedata(df, idxVals)

actOut = TimeData.chkElw(x -> x.>130,
                         TimeData.eachentry(testcase2)) |>
         x -> TimeData.asArr(x, Bool, false) |>
         x -> TimeData.getEntries(testcase2, x, false)

@test isequal(expTd2, actOut)

####################################
## test getEntries with indexing ##
####################################

idxVals = Date[Date("2010-01-01"),
               Date("2010-01-05"),
               Date("2010-01-01"),
               Date("2010-01-03")]
variables = [:prices1, :prices1, :prices2, :prices2]
values = Any[100, 200, 110, NA]

df = DataFrame(variable = variables, value = values)
expTd = TimeData.Timedata(df, idxVals)

singleInd = [1, 5, 6, 8]
@test isequal(expTd, TimeData.getEntries(testcase1, singleInd))

@test isequal(expTd,
              TimeData.getEntries(testcase1,
                          ind2sub(size(testcase1),
                                  singleInd)...))

end


text:
# # Collecting and visualizing boundary vertices

# In this tutorial, we will learn
#
#    -  How to import an Abaqus file.
#    -  How to find the vertices on the boundary of the mesh.
#    -  How to visualize the vertices on the boundary of the mesh.

# The tutorial will produce files for mesh visualization in the
# [Paraview](https://www.paraview.org/) format (VTK). One can display this
# information by loading the file with `paraview.exe`. When the tutorial is
# executed in `mybinder.org`, the graphics file needs to be downloaded to your
# desktop, and then visualized locally.

# In this tutorial we import an Abaqus mesh. It represents a rectangular block
# with a circular opening, meshed with quadrilaterals.
using MeshSteward: import_ABAQUS
connectivities = import_ABAQUS("block-w-hole.inp");

# The imported incidence relations are stored in an array, `connectivities`. In
# the present case is just one such incidence relation in the file, and we will
# use it to create the mesh of the block.
using MeshSteward: Mesh, attach!
mesh = Mesh()
attach!(mesh, connectivities[1]);

# We can visualize the mesh by exporting the base relation of the mesh as a VTK
# file.
using MeshSteward: vtkwrite
using MeshSteward: baseincrel
vtkwrite("block-w-hole", baseincrel(mesh))

# The boundary of the mesh consists of curve (line) segments, and may be
# extracted as:
using MeshSteward: boundary
bir = boundary(mesh);
vtkwrite("block-w-hole-boundary", bir)

# To collect all vertices connected by the boundary segments, we can do:
using MeshSteward: connectedv
vl = connectedv(bir);

# The result is a list of integer identifiers of the connected vertices. By
# comparing the length of the list of the vertices on the boundary with the total
# number of vertices (the number of shapes on the right of the incidence
# relation), we can see that less than one fourth of the vertices lies on the
# boundary.
using MeshCore: nshapes
@show length(vl)
@show nshapes(baseincrel(mesh).right)

# The shape collection `baseincrel(mesh).right` represents the vertices in the
# mesh. The very same collection may be retrieved from the mesh as:
using MeshSteward: vertices, summary
verts = vertices(mesh);
summary(verts)

# The incidence relation `(0, 0)` that represents just the vertices on the
# boundary can be created as a subset of all the vertices in the mesh.
using MeshCore: ir_subset, summary
ssverts = ir_subset(verts, vl);
summary(ssverts)

# The vertices on the boundary may be exported for visualization into a VTK
# file.
using MeshCore: nshapes
@show nshapes(ssverts.left), nshapes(ssverts.right)

# Here we inspect the incidence relation. We can see that the shape collection on the left is called "shapes".
using MeshSteward: summary
@show summary(ssverts)

# Hence, we instruct the export function to incorporate into the output file the
# shape collection on the left of the incidence relation. We do that by
# specifying the tag "shapes". The information will be written out as cell
# data, and can be visualized with glyphs. The entities represented in the
# file are simply points. In paraview, they may be visualized with "3D glyphs"
# as little balls: choose "Representation" to be "3D glyphs", and then the
# glyph type "Sphere".
vtkwrite("block-w-hole-boundary-vertices", ssverts, [(name = "shapes",)])

# Start "Paraview", load the boundary and boundary-vertex files. Visualize the boundary vertices with glyphs.
@async run(`paraview `)

# The visualization may look like this:
# ![Visualization of the facets and vertices](t6_i1.png)