text:
function solve_worker(m :: Model)
    sol = init_solution(workspan(m), lifespan(m), n_efficiencies(m, 1));
    solve_last_period!(m, sol);
    for t = (lifespan(m) - 1) : -1 : (workspan(m) + 1)
        solve_one_period!(m, sol, t);
    end
    return sol
end

function solve_last_period!(m :: Model, sol :: Solution)
    set_solution!(sol, lifespan(m), solve_last_period(m));
end

# This is copied almost verbatim from our previous model (without shocks).
# Returns functions `k' = G(k, T)` and `c = H(k, T)`.
function solve_last_period(m :: Model)
    T = lifespan(m);
    kGridV = make_k_grid(m, T);
    kPrimeFct = interpolate_kprime(kGridV, zeros(length(kGridV)));
    # kPrime and efficiency are 0
    cV = consumption(m, T, kGridV, 0.0, 0.0);
    cFct = interpolate_kprime(kGridV, cV);
    return WorkerSolution([kPrimeFct], [cFct])
end

function solve_one_period!(m :: Model, sol :: Solution, t :: Integer)
    newSol = solve_one_period(m, get_solution(sol, t+1), t);
    set_solution!(sol, t, newSol);
end

# test this +++++
# The key input here is `kPrimeTomorrow = k'' = G(k', e', t+1)`.
# This is tomorrow's policy function.
function solve_one_period(m :: Model, kPrimeTomorrow :: WorkerSolution, t :: Integer)
    ne = n_efficiencies(m, t);
    kPrimeFctV = Vector{Any}(undef, ne);
    cFctV = Vector{Any}(undef, ne);
    for eIdx = 1 : ne
        kPrimeFctV[eIdx], cFctV[eIdx] =
            solve_one_efficiency(m, t, eIdx, kPrimeTomorrow);
    end
    return WorkerSolution(kPrimeFctV, cFctV)
end

function solve_one_efficiency(m :: Model, t, eIdx, kPrimeTomorrow :: WorkerSolution)
    kGridV = make_k_grid(m, t);
    kPrimeV = zeros(length(kGridV));
    for (ik, k) in enumerate(kGridV)
        kPrimeV[ik] = solve_one_point(m, t, k, eIdx, kPrimeTomorrow);
    end
    kPrimeFct = interpolate_kprime(kGridV, kPrimeV);
    eff = efficiency(m.endow, eIdx);
    cV = consumption(m.budget, t, kGridV, kPrimeV, eff);
    cFct = interpolate_kprime(kGridV, cV);
    return kPrimeFct, cFct
end

# Solve for one point, keeping in mind that the capital grid defines bounds for k'.
# Therefore, we may have corner solutions.
function solve_one_point(m :: Model, t, k, eIdx, kPrimeTomorrow)
    eff = efficiency(m.endow, eIdx);
    kPrimeMin, kPrimeMax = kprime_range(m, t, k, eff);
    @assert (kPrimeMax > kPrimeMin)  "No feasible kPrime for k = $k: $kPrimeMin, $kPrimeMax"

    # Define a closure that can be passed to root finding
    e_dev(kPrime) = euler_dev_one_point(m, t, k, eIdx, kPrime, kPrimeTomorrow);

    # Try corners
    eDev0 = e_dev(kPrimeMin);
    if eDev0 <= 0.0
        kPrime = kPrimeMin;
    else
        eDevMax = e_dev(kPrimeMax);
        if eDevMax >= 0.0
            kPrime = kPrimeMax
        else
            # Now we know the solution is interior
            kPrime = find_zero(e_dev, (kPrimeMin, kPrimeMax), Bisection());
        end
    end
    @assert kPrimeMin <= kPrime <= kprime_max(m.kGrid, t)
    return kPrime
end

# Euler equation deviation.
# `dev > 0` means that `kPrime` implies a `c` today that is too high.
# We use the form `c - U'^{-1}(βR U'(c'))`.
function euler_dev_one_point(m :: Model, t :: Integer, k, eIdx :: Integer,
    kPrime, kPrimeTomorrow)

    cPrimeV, prob_ePrimeV = cprime_distribution(m, t, eIdx,
        kPrime, kPrimeTomorrow);
    # if isretired(m, t+1)
    #     prob_ePrimeV = [1.0];
    # else
    #     prob_ePrimeV = prob_eprime(m.endow, eIdx);
    # end
    ePrimeV = 1 : length(prob_ePrimeV);
    rhs = 0.0;
    # for ePrimeIdx in ePrimeV
        # effPrime = efficiency(m.endow, eIdx);
        # kPrimePrime = kprime(kPrimeTomorrow, kPrime, ePrimeIdx);
        # cPrime = consumption(m.budget, t+1, kPrime,
        #     kPrimePrime; efficiency = effPrime);
    for (j, cPrime) in enumerate(cPrimeV)
        # βR U'(c')
        betaRUprime = betar(m) * marg_utility(m.util, cPrime);
        rhs += prob_ePrimeV[j] * inv_marg_utility(m.util, betaRUprime);
    end

    eff = efficiency(m.endow, eIdx);
    cToday = consumption(m.budget, t, k, kPrime, eff);
    dev = cToday - rhs;
    return dev
end

function cprime_distribution(m :: Model, t :: Integer, eIdx :: Integer,
    kPrime, kPrimeTomorrow :: WorkerSolution)

    if isretired(m, t+1)
        prob_ePrimeV = [1.0];
    else
        prob_ePrimeV = prob_eprime(m.endow, eIdx);
    end
    ePrimeV = 1 : length(prob_ePrimeV);
    cPrimeV = zeros(length(ePrimeV));
    for ePrimeIdx in ePrimeV
        effPrime = efficiency(m.endow, eIdx);
        kPrimePrime = kprime(kPrimeTomorrow, kPrime, ePrimeIdx);
        cPrimeV[ePrimeIdx] = consumption(m.budget, t+1, kPrime,
            kPrimePrime, effPrime);
    end
    return cPrimeV, prob_ePrimeV
end

# E(u'(c'))
function e_uprime_c(m :: Model, t :: Integer, eIdx :: Integer,
    kPrime, kPrimeTomorrow :: WorkerSolution)

    cPrimeV, prob_ePrimeV = cprime_distribution(m, t, eIdx,
        kPrime, kPrimeTomorrow);
    # R U'(c')
    muV = [marg_utility(m.util, cPrime) for cPrime in cPrimeV];
    eMu = sum(muV .* prob_ePrimeV)
    return eMu
end

function interpolate_kprime(kGridV :: AbstractVector{F}, kPrimeV :: AbstractVector{F}) where F
    interp = LinearInterpolation(kGridV, kPrimeV);
    return interp
end

# update for eIdx +++++++++
function check_retired_solution(m :: Model, sol :: Solution)
    bcValid = true;
    eulerValid = true;

    for t = (workspan(m) + 1) : (lifespan(m) - 1)
        kGridV = make_k_grid(m, t);
        kp_fct = kprime_fct(sol, t);
        c_fct = cons_fct(sol, t);
        wage_t = wage(m);
        for (ik, k) in enumerate(kGridV)
            kPrime = kp_fct(k);
            c = consumption(m.budget, t, k, kPrime, eff);
            bcValid = bcValid && isapprox(c, c_fct(k), atol = 1e-4);

            if kprime_min(m.kGrid, t) < kPrime < kprime_max(m.kGrid, t)
                # Interior point
                kpp_fct = kprime_fct(sol, t+1);
                kPrimePrime = kpp_fct(kPrime);
                cPrime = consumption(m.budget, t+1, kPrime, kPrimePrime, effPrime);
                eDev = only(euler_dev(m, [c, cPrime]));
                eulerValid = eulerValid && abs(eDev < 1e-4);
            end
        end
    end

    if !bcValid
        @warn "Budget constraint violated"
    end
    if !eulerValid
        @warn "Euler equation violated"
    end

    return bcValid  &&  eulerValid
end

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


text:
module CircuitViz

using Catlab.WiringDiagrams
using Catlab.Graphs
using Catlab.Graphics: LeftToRight
import Catlab.Graphics.GraphvizGraphs: to_graphviz
using Catlab.Graphs
import Catlab.Graphs: PropertyGraph, add_vertex!,
                      add_edge!, rem_vertex!,
                      rem_edge!
using Catlab.Graphics: Graphviz
using JSON
import Base: show

export Graph, add_vertex!, add_edge!, rem_vertex!, rem_edge!
export to_graphviz, render
export compile!, show!

function Graph()
    g = PropertyGraph{Any}()
    graph_attrs = Dict(
                       :rankdir => "LR",
                       :dpi => "288.0",
                       :nodesep => "0.15",
                       :pack => "true",
                       :ratio => "compress",
                      )
    node_attrs = Dict(
                      :fontname => "Courier",
                      :shape => "triangle",
                      :orientation => "30",
                      :penwidth => "3.0",
                     )
    edge_attrs = Dict(
                      :fontname => "Courier",
                      :labelangle => "25",
                      :labeldistance => "2",
                      :penwidth => "1.2",
                     )
    set_gprop!(g, :graph, graph_attrs)
    set_gprop!(g, :node, node_attrs)
    set_gprop!(g, :edge, edge_attrs)
    return g
end

function add_vertex!(g::PropertyGraph{T}) where T
    return add_vertex!(g, Dict{Symbol, T}())
end

function add_edge!(g::PropertyGraph{T}, src::Int, tg::Int) where T
    return add_edge!(g, src, tg, Dict{Symbol, T}())
end

function restrict_to_edges(g)
    vs = Set([])
    for v in Graphs.vertices(g)
        isempty(Graphs.inneighbors(g, v)) &&
        isempty(Graphs.outneighbors(g, v)) &&
        continue
        push!(vs, v)
    end
    return Graphs.induced_subgraph(g, vs)
end

function render(graph)
    g = to_graphviz(graph)
    return g
end

function compile!(graph;
        base_path = "$(gensym("digraph"))",
        format = :png)
    dot_path = base_path * ".dot"
    output_path = base_path * ".$(format)"
    open(dot_path, "w") do io
        g = render(graph)
        Graphviz.pprint(io, g)
    end
    run(pipeline(`dot -Gdpi=288 -T$(format) $(dot_path)`, stdout=output_path))
end

end # module

ext: jl


text:
# This file is part of Kpax3. License is MIT.

mutable struct KWeight
  c::Vector{Float64}
  w::Matrix{Float64}
  z::Matrix{Float64}
end

mutable struct MCMCSupport
  m::Int
  n::Int

  u::Vector{Int}
  t::Vector{Float64}

  lp::Array{Float64, 3}
  lq::Array{Float64, 3}
  lr::Array{Float64, 3}

  vi::Int
  ni::Vector{Float64}
  ui::Vector{Int}
  wi::KWeight
  lpi::Matrix{Float64}

  vj::Int
  nj::Vector{Float64}
  uj::Vector{Int}
  wj::KWeight
  lpj::Matrix{Float64}

  tmp::Vector{Float64}

  cl::Vector{Int}
  k::Int

  lograR::Float64

  logmlik::Float64
  logmlikcandidate::Float64
end

function MCMCSupport(state::State,
                     priorC::AminoAcidPriorCol)
  n = length(state.R)
  (maxclust, m) = size(state.C)

  u = Int[a for a in 1:n]
  t = zeros(Float64, n)

  lp = zeros(Float64, 4, maxclust, m)
  lq = zeros(Float64, 4, maxclust, m)
  lr = zeros(Float64, 3, maxclust, m)

  g = 0
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b],
                             priorC.B[1, b])
    lp[2, g, b] = logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b],
                             priorC.B[2, b])
    lp[3, g, b] = logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b],
                             priorC.B[3, b])
    lp[4, g, b] = logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b],
                             priorC.B[4, b])
  end

  vi = 0
  ni = zeros(Float64, m)
  ui = zeros(Int, n)
  wi = KWeight(zeros(Float64, m), zeros(Float64, 4, m), zeros(Float64, 4, m))

  lpi = zeros(Float64, 4, m)

  vj = 0
  nj = zeros(Float64, m)
  uj = zeros(Int, n)
  wj = KWeight(zeros(Float64, m), zeros(Float64, 4, m), zeros(Float64, 4, m))

  lpj = zeros(Float64, 4, m)

  tmp = zeros(Float64, 4)

  cl = zeros(Int, n)
  k = 0

  logmlik = logmarglikelihood(state.cl, state.k, lp, priorC)

  MCMCSupport(m, n, u, t, lp, lq, lr, vi, ni, ui, wi, lpi, vj, nj, uj, wj, lpj,
              tmp, cl, k, 0.0, logmlik, 0.0)
end

function resizesupport!(support::MCMCSupport,
                        maxclust::Int)
  if size(support.lp, 2) < maxclust
    # we don't want to allocate new resources too often, so allocate double the
    # previous size. this should guarantee a logarithmic number of allocations
    newsize = min(support.n, max(maxclust, 2 * size(support.lp, 2)))

    lp = zeros(Float64, 4, newsize, support.m)

    for b in 1:support.m, g in 1:size(support.lp, 2)
      lp[1, g, b] = support.lp[1, g, b]
      lp[2, g, b] = support.lp[2, g, b]
      lp[3, g, b] = support.lp[3, g, b]
      lp[4, g, b] = support.lp[4, g, b]
    end

    support.lp = lp
    support.lq = zeros(Float64, 4, newsize, support.m)
    support.lr = zeros(Float64, 3, newsize, support.m)
  end

  nothing
end

mutable struct KOffspring
  R::Vector{Int}
  v::Vector{Int}
end

mutable struct GASupport
  m::Int
  n::Int

  oi::KOffspring
  oj::KOffspring
end

function GASupport(m::Int,
                   n::Int)
  oi = KOffspring(zeros(Int, n), zeros(Int, n))
  oj = KOffspring(zeros(Int, n), zeros(Int, n))
  GASupport(m, n, oi, oj)
end


text:
export ExpWeight

immutable ExpWeight{S,FF,DD} <: WeightSpace{S,RealBasis,DD,1}
    exponent::FF
    space::S
end

function ExpWeight(q,sp)
    @assert domain(q) == domain(sp)
    ExpWeight{typeof(sp),typeof(q),typeof(domain(q))}(q,sp)
end

weight(sp::ExpWeight,x) = exp(sp.exponent(x))
spacescompatible(a::ExpWeight,b::ExpWeight) = a.exponent ≈ b.exponent && spacescompatible(a.space,b.space)
setdomain(sp::ExpWeight,d::Domain) = ExpWeight(setdomain(sp.exponent,d),setdomain(sp.space,d))

for (RULE,TYP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace))
    @eval function $RULE(a::ExpWeight,b::ExpWeight)
        if a.exponent ≈ b.exponent
            sp=$TYP(a.space,b.space)
            if !isa(sp,NoSpace)
                return ExpWeight(a.exponent,sp)
            end
        end
        NoSpace()
    end
end

function Derivative(wsp::ExpWeight)
    q=wsp.exponent
    sp=wsp.space
    D=Derivative(sp) + q'
    DerivativeWrapper(SpaceOperator(D,wsp,ExpWeight(q,rangespace(D))),1)
end

function Derivative(S::ExpWeight,k::Integer)
    if k==1
        Derivative(S)
    else
        D=Derivative(S)
        DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1),D),k)
    end
end

function Multiplication(f::Fun,sp::ExpWeight)
    M=Multiplication(f,sp.space)
    MultiplicationWrapper(f,SpaceOperator(M,sp,ExpWeight(sp.exponent,rangespace(M))))
end

function Conversion(a::ExpWeight,b::ExpWeight)
    @assert a.exponent ≈ b.exponent
    ConversionWrapper(SpaceOperator(Conversion(a.space,b.space),a,b))
end


text:
using SemialgebraicSets
using Test

using MultivariatePolynomials

# Taken from JuMP/test/solvers.jl
function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

if try_import(:DynamicPolynomials)
    Mod = DynamicPolynomials
    include("commutativetests.jl")
end

if try_import(:TypedPolynomials)
    Mod = TypedPolynomials
    include("commutativetests.jl")
end


text:
#Permuation Mutations

function inversion!{T <: Vector}(recombinant::T)
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    l = div(to - from,2)+1
    for i in 0:(l-1)
        swap!(recombinant, from+i, to-i)
    end
    return recombinant
end

function insertion!{T <: Vector}(recombinant::T)
    l = length(recombinant)
    from, to = rand(1:l, 2)
    val = recombinant[from]
    deleteat!(recombinant, from)
    return insert!(recombinant, to, val)
end

function swap2!{T <: Vector}(recombinant::T)
    l = length(recombinant)
    from, to = rand(1:l, 2)
    swap!(recombinant, from, to)
    return recombinant
end

function scramble!{T <: Vector}(recombinant::T)
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    recombinant[from:to] = shuffle!(recombinant[from:to])
    return recombinant
end

function shifting!{T <: Vector}(recombinant::T)
    l = length(recombinant)
    from, to, where = sort(rand(1:l, 3))
    patch = recombinant[from:to]
    diff = where - to
    if diff > 0
        # move values after tail of patch to the patch head position
        for i in 1:diff
            recombinant[from+i-1] = recombinant[to+i]
        end
        # place patch values in order
        start = from + diff
        for i in 1:length(patch)
            recombinant[start+i-1] = patch[i]
        end
    end
    return recombinant
end


text:
push!(LOAD_PATH,"./src")
using Pkg
Pkg.update()
Pkg.add("Revise")
Pkg.add("SearchLight")
Pkg.add("SearchLightPostgreSQL")
Pkg.add("BitemporalPostgres")
Pkg.add("Intervals")
Pkg.add("TimeZones")
Pkg.add("Dates")
using Revise
using SearchLight
using SearchLightPostgreSQL
using BitemporalPostgres
using TimeZones
run(```psql -f sqlsnippets/droptables.sql```)
SearchLight.Configuration.load() |> SearchLight.connect
SearchLight.Migrations.create_migrations_table()
BitemporalPostgres.up()
include("db/migrations/2022020318035833_create_table_bitemporalwebapp.jl")

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
import AbstractAlgebra

function read_MQ_GF(filename)
    @info "Reading MQ $filename"
    polys = []
    apath = abspath("") * "benchmark\\MQ\\$filename"
    open(apath, "r") do fio
        line1, line2, line3 = readline(fio), readline(fio), readline(fio)
        ch = parse(Int, line1[findfirst('(', line1)+1:findfirst(')', line1)-1])
        nv = parse(Int, line2[findfirst(':', line2)+1:end])
        ne = parse(Int, line3[findfirst(':', line3)+1:end])

        R, xs = AbstractAlgebra.PolynomialRing(AbstractAlgebra.GF(ch), nv, ordering=:degrevlex)
        @info "Created ring" R xs

        labels = sort!(union([x*y for x in xs for y in xs], xs, [R(1)]), rev=true)

        rest = readlines(fio)
        rest = rest[5:end]

        for (i, s) in enumerate(rest)
            @info "Reading poly $i/$(length(rest))"
            s = replace(s, " ;" => "")
            poly = sum(map(x -> parse(Int, x), split(s, " ")) .* labels)
            push!(polys, poly)
        end
    end
    @info "Read" length(polys)
    R = parent(first(polys))
    map(R, polys)
end


text:
using Test
using QuantumOptics.sparsematrix
using SparseArrays, LinearAlgebra

# SparseMatrix = quantumoptics.sparsematrix.SparseMatrix
const SparseMatrix = SparseMatrixCSC{ComplexF64, Int}

@testset "sparsematrix" begin

# Set up test matrices
A = rand(ComplexF64, 5, 5)
A_sp = sparse(A)

B = Matrix{ComplexF64}(I, 5, 5)
B_sp = sparse(ComplexF64(1)*I, 5, 5)

C = rand(ComplexF64, 3, 3)
C[2,:] .= 0
C_sp = sparse(C)

R_sp = A_sp + B_sp
R = A + B

# Test arithmetic
@test 0 ≈ norm(Matrix(R_sp) - R)
@test 0 ≈ norm(Matrix(ComplexF64(0.5,0)*A_sp) - 0.5*A)
@test 0 ≈ norm(Matrix(A_sp/2) - A/2)
@test 0 ≈ norm(Matrix(A_sp*B_sp) - A*B)

# Test kronecker product
@test 0 ≈ norm(Matrix(kron(A_sp, C_sp)) - kron(A, C))
@test 0 ≈ norm(Matrix(kron(A_sp, B_sp)) - kron(A, B))

end # testset


text:
export zero_cost_transitions, transition_costs, gold_scores
using Flux

function zero_cost_transitions(state::GoldState)
  filter(trans -> is_valid(state.config, trans.move) && cost(state, trans.move, trans.label) <= 0, state.system.transitions)
end

# I inverse cost& because I need to zero cost transitions have best value after softmax
function transition_costs(state::GoldState)
  map(state.system.transitions) do trans
    is_valid(state.config, trans.move) ? cost(state, trans.move, trans.label) : FORBIDDEN_COST
  end
end

function gold_scores(costs::Vector{Int64})
  map(costs) do cost
    cost <= 0 ? 1 : 0
  end |> vec -> reshape(vec, 1, :)
end

function weight_mask(costs)
  shift_weights = [1, 1, 2]
  other_weights = [[2, 1, 2] for i in 1:(length(costs) - 1)]

  hcat(shift_weights, other_weights...)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using DSGE, JLD2, Dates, ModelConstructors, DataFrames
using Test, FileIO

path = dirname(@__FILE__)

# Initialize models
function make_test_model(year::Int)
    m = Model990()
    vint = Dates.format(Dates.Date(year, 4, 10), "yymmdd")
    m <= Setting(:data_vintage, vint)
    m <= Setting(:cond_vintage, vint)
    m <= Setting(:date_forecast_start, DSGE.quartertodate("$year-Q1"))
    m <= Setting(:date_conditional_end, DSGE.quartertodate("$year-Q1"))
    m <= Setting(:forecast_horizons, 16)
    return m
end

m_new = make_test_model(2016)
m_old = make_test_model(2014)

# Read in data and parameters
@load "$path/../reference/decompose_forecast_args.jld2" df_new df_old params_new params_old

# Read in expected outputs
@load "$path/../reference/decompose_forecast_out.jld2" exp_decomps

cond_types = [(:none, :none), (:none, :full), (:full, :none), (:full, :full)]
for (cond_new, cond_old) in cond_types
    # Test inner method
    decomps = decompose_forecast(m_new, m_old,
                                 cond_new == :full ? df_new : df_new[1:end-1, :],
                                 cond_old == :full ? df_old : df_old[1:end-1, :],
                                 params_new, params_old, cond_new, cond_old,
                                 [:obs, :pseudo]; check = true)

    for k in keys(decomps)
        @test decomps[k] ≈ exp_decomps[(cond_new, cond_old)][k]
    end
end

# Test outer method
m_new = AnSchorfheide(testing = true)
m_new <= Setting(:saveroot, tempdir())
m_new <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m_new <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m_new <= Setting(:use_population_forecast, true)
m_new <= Setting(:forecast_horizons, 12)

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m_new)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "mhsave_test.h5")
m_new <= Setting(:forecast_block_size, 20)

m_old = deepcopy(m_new)
m_old <= Setting(:date_forecast_start, quartertodate("2014-Q4"))
m_old <= Setting(:date_conditional_end, quartertodate("2014-Q4"))

df_new = load_data(m_new)
df_old = df_new[1:end-4, :]

decompose_forecast(m_new, m_old, df_new, df_old, :mode, :none, :none, [:obs, :pseudo];
                         verbose = :none)
decomposition_means(m_new, m_old, :mode, :none, :none, [:obs, :pseudo]; verbose = :none)
decompose_forecast(m_new, m_old, df_new, df_old, :full, :none, :none, [:obs, :pseudo], verbose = :none)
decomposition_means(m_new, m_old, :full, :none, :none, [:obs, :pseudo], verbose = :none)

## Regime switching
custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_presample_start     => Setting(:date_presample_start, Date(1959, 9, 30)),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_conditional_end, DSGE.quartertodate("2016-Q3")),
    :forecast_horizons        => Setting(:forecast_horizons, 16),
    :n_mon_anticipated_shocks => Setting(:n_mon_anticipated_shocks, 6))
m    = Model1002("ss10", testing = true, custom_settings = custom_settings)  # baseline model
m_rs = Model1002("ss51v", testing = true, custom_settings = custom_settings) # pseudo regime switching (identical values for standard deviations)
m_rs <= Setting(:rate_expectations_source, :ois)
m_rs.settings[:regime_switching] = Setting(:regime_switching, true)
m_rs.settings[:n_regimes] = Setting(:n_regimes, 3)
m_rs.settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
m_rs.test_settings[:regime_switching] = Setting(:regime_switching, true)
m_rs.test_settings[:n_regimes] = Setting(:n_regimes, 3)
m_rs.test_settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
m_rs.settings[:regime_dates] = Setting(:regime_dates,
                                       Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31), 3 => Date(2012, 9, 30)))
m_rs.test_settings[:regime_dates] = Setting(:regime_dates,
                                            Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31), 3 => Date(2012, 9, 30)))
df = load("$path/../reference/regime_switch_data.jld2", "none")

for i in 1:3
    adj = (i == 1) ? 1. : .95
    ModelConstructors.set_regime_val!(m_rs[:α], i, adj * m[:α].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_g], i, adj * m[:σ_g].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_b], i, adj * m[:σ_b].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_μ], i, adj * m[:σ_μ].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_ztil], i, adj * m[:σ_ztil].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_λ_f], i, adj * m[:σ_λ_f].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_λ_w], i, adj * m[:σ_λ_w].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_r_m], i, adj * m[:σ_r_m].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_σ_ω], i, adj * m[:σ_σ_ω].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_μ_e], i, adj * m[:σ_μ_e].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_γ], i, adj * m[:σ_γ].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_π_star], i, adj * m[:σ_π_star].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_lr], i, adj * m[:σ_lr].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_z_p], i, adj * m[:σ_z_p].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_tfp], i, adj * m[:σ_tfp].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_gdpdef], i, adj * m[:σ_gdpdef].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_corepce], i, adj * m[:σ_corepce].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_gdp], i, adj * m[:σ_gdp].value; override_bounds = true)
    ModelConstructors.set_regime_val!(m_rs[:σ_gdi], i, adj * m[:σ_gdi].value; override_bounds = true)

    for j = 1:DSGE.n_mon_anticipated_shocks(m_rs)
        ModelConstructors.set_regime_val!(m_rs[Symbol("σ_r_m$(j)")], i, adj * m[Symbol("σ_r_m$(j)")]; override_bounds = true)
    end
end

m_rs_old = deepcopy(m_rs)
m_rs_old <= Setting(:date_forecast_start, quartertodate("2015-Q3"))
m_rs_old <= Setting(:date_conditional_end, quartertodate("2015-Q3"))
m_rs_old.settings[:forecast_horizons] = Setting(:forecast_horizons, 12)
m_rs_old.test_settings[:forecast_horizons] = Setting(:forecast_horizons, 12)
df_old = df[1:end - 4, :]

decompose_forecast(m_rs, m_rs_old, df, df_old, map(x -> x.value, m_rs.parameters),
                   map(x -> x.value, m_rs_old.parameters), :none, :none, [:obs, :pseudo])
@test_broken decomposition_means(m_rs, m_rs_old, :mode, :none, :none, [:obs, :pseudo]; verbose = :none) # TODO: NEED METHOD THAT LOADS THE DRAWS FROM A SAVED FILE B/C THAT ONE WRITES TO STUFF TO A JLD2 FILE

nothing


text:
export SisMcmcOutput, SimMcmcOutput, SpfMcmcOutput

struct SisMcmcOutput
    model::SIS # The model from which the sample was drawn
    sample::Vector{Vector{Path{Int}}}  # The sample
    performance_measures::Dict  # Dictionary of performance measures key => value, e.g. "acceptance probability" => 0.25
end

struct SimMcmcOutput
    model::SIM
    sample::Vector{Vector{Path{Int}}}  # The sample
    performance_measures::Dict  # Dictionary of performance measures key => value, e.g. "acceptance probability" => 0.25
end

struct SpfMcmcOutput{T<:Union{String,Int}}
    model::SPF{T} # The model from which the sample was drawn
    sample::Vector{Path{T}}  # The sample
    a::Real # Acceptance Probability
end

function Base.show(io::IO, output::T) where {T<:SisMcmcOutput}
    title = "MCMC Sample for SIS Model"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    for (key, value) in output.performance_measures
        println(io, key, ": ", value)
    end
end

function Base.show(io::IO, output::T) where {T<:SimMcmcOutput}
    title = "MCMC Sample for SIM Model"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    for (key, value) in output.performance_measures
        println(io, key, ": ", value)
    end
end

function Base.show(io::IO, output::SpfMcmcOutput)
    title = "MCMC Sample for Spherical Path Family (SPF)"
    println(io, title)
    println(io, "-"^length(title))
    println(io, "\nAcceptance probability: $(output.a)")
end


text:
using Documenter
using SAMTools

makedocs(;
    modules  = [SAMTools],
    doctest  = false,
    format   = Documenter.HTML(
        collapselevel = 1,
        prettyurls    = false
    ),
    authors  = "Nathanael Wong <natgeo.wong@outlook.com>",
    sitename = "SAMTools.jl",
    pages    = [
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/natgeo-wong/SAMTools.jl.git",
)


text:
export
    spatialorder,
    coords_spatial,
    size_spatial,
    spacedirections,
    pixelspacing,
    indices_spatial

@deprecate(spatialorder(x::AbstractArray), spatial_order(x))

@deprecate(coords_spatial(x::AbstractArray), spatialdims(x))

@deprecate(size_spatial(x::AbstractArray), spatial_size(x))

@deprecate(indices_spatial(x::AbstractArray), spatial_indices(x))

@deprecate(spacedirections(x::AbstractArray), spatial_directions(x))

@deprecate(pixelspacing(x), pixel_spacing(x))


text:
using LinearAlgebra
κ = 1.0
n = 100
p = 10
x0 = randn(p)
#linear regression
A = randn(n,p)
σ0 = 0.01
y0 = A*x0 .+ randn(n)*σ0

θ0 = rand((-1.0,1.0), n)
T = 100.0
c = 10ones(n)
function ∇ϕ!(x, j, A, y0, σ0, n)
    i = rand(1:n)
    A[i,j]*(dot(A[i,:],x) - y[j])/σ0^2
end

function ∇ϕ!(x, args...)
    i = rand(1:length(x))
    ∇ϕ!(x, i, args...)
return E_tilde, i

#@time trace0, _ = ZigZagBoomerang.spdmp(ϕ, 0.0, x0, θ0, T, c, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)

@time trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, κ, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)
ts0, xs0 = splitpairs(trace0)


text:
using Test

@testset "Broyden ($R)" for R in [Float32, Float64]
    using LinearAlgebra
    using ProximalAlgorithms: Broyden, update!

    H = R[
        0.63287    0.330934   -0.156908   -0.294776    0.10761;
        0.330934   0.673201    0.0459778   0.231011   -0.235265;
       -0.156908   0.0459778   0.635812   -0.232261   -0.388775;
       -0.294776   0.231011   -0.232261    0.726854   -0.0691783;
        0.10761   -0.235265   -0.388775   -0.0691783   0.336262;
    ]
    l = R[1.0, 2.0, 3.0, 4.0, 5.0]

    n = length(l)

    f(x) = dot(x, H * x) / 2 + dot(x, l)

    x_star = -H \ l
    f_star = f(x_star)
    norm_x_star = norm(x_star)

    Lip = opnorm(H)
    gamma = 1 / Lip
    x = zeros(n)
    res = zeros(n)
    x_prev = zeros(n)
    res_prev = zeros(n)

    acc = Broyden(x)

    mul!(res, H, x)
    res .+= l
    res .*= gamma

    d = acc * res

    for it = 1:10
        # store iterate and residual for the operator update later
        res_prev .= res
        x_prev .= x

        # compute accelerated direction
        mul!(d, acc, res)

        # update iterate
        x .-= d

        # compute new residual
        mul!(res, H, x)
        res .+= l
        res .*= gamma

        # update operator
        update!(acc, x - x_prev, res - res_prev)
    end

    @test f(x) <= f_star + (1 + abs(f_star)) * sqrt(eps(R))

end


text:
import BSplineKit
using BSplineKit.BandedTensors

import BSplineKit.BandedTensors: SubMatrix

using BandedMatrices
using LinearAlgebra: dot
using Random
using Test

function test_banded_tensors()
    A = @inferred BandedTensor3D{Float32}(undef, (20, 20, 16), Val(2))
    rand!(A, -99:99)
    BandedTensors.drop_out_of_bands!(A)

    @testset "Slices" begin
        @inferred A[:, :, 4]
        @inferred A[:, :, 3:5]
        S = A[:, :, 4]
        Acut = A[:, :, 3:5]
        @test S isa SubMatrix
        @test Acut isa BandedTensor3D
        @test Acut[:, :, 2] == S  # SubMatrix comparison
        @test S[S.inds, S.inds] == S.data  # slice over non-zero indices
    end

    @testset "Contraction" begin
        b = rand(size(A, 3))
        Y = @inferred A * b
        @test Y isa BandedMatrix
        @test Y[:, 3] ≈ A[:, 3, :] * b
    end

    @testset "Generalised dot product" begin
        for k in axes(A, 3)
            S = A[:, :, k] :: SubMatrix
            u = rand(size(S, 1))
            v = rand(size(S, 2))

            # Efficient implementation
            x = dot(u, S, v)

            # Inefficient computation using full matrix
            y = dot(u, Array(S), v)

            @test x ≈ y
        end
    end

    nothing
end

@testset "BandedTensors" begin
    test_banded_tensors()
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
###############################################################################
#
# Grid data
#
###############################################################################
"""
    grid_RAW!(param::Array)
    grid_RAW!(dt::MOD15A2Hv006LAI{FT},
              params::Array{Any,1}) where {FT<:AbstractFloat}
    grid_RAW!(dt::MOD15A2Hv006LAI{FT},
              params::Array{Array,1},
              nthread::Int) where {FT<:AbstractFloat}

Grid the data and save it to a csv file, given
- `param` Parameter sets to pass to different threads
- `dt` [`AbstractUngriddedData`](@ref) type ungridded data type
- `params` Array of `param`
- `nthread` Number of threads to run the code in parallel
"""
function grid_RAW!(param::Array{Any,1})
    grid_RAW!(param...);

    return nothing
end

function grid_RAW!(
            dt::MOD09A1v006NIRv{FT},
            nlayer::Int,
            path::String,
            prefix::String,
            band::Array{String,1},
            cache::String,
            mm::Array{FT,1}
) where {FT<:AbstractFloat}
    # transform the data to CSV file
    pos_start = findfirst(prefix, path);
    hv_string = path[pos_start[1]:end];
    outp_file = joinpath(cache, hv_string * ".csv");

    # if file exists, skip
    if !isfile(outp_file)
        h,v = parse_HV(hv_string);
        lat_mat,lon_mat = load_MODIS(dt, h, v);
        dat_out = DataFrame();
        try
            red_mat = read_nc(FT, path, band[1]);
            nir_mat = read_nc(FT, path, band[2]);
            dim_lon = size(red_mat,1);
            dim_lat = size(red_mat,2);
            (dat_out).lay = ones(Int, dim_lon * dim_lat) * nlayer;
            (dat_out).lat = ones(dim_lon * dim_lat) .* -9999;
            (dat_out).lon = ones(dim_lon * dim_lat) .* -9999;
            (dat_out).red = ones(dim_lon * dim_lat) .* -9999;
            (dat_out).nir = ones(dim_lon * dim_lat) .* -9999;
            for i in 1:dim_lon, j in 1:dim_lat
                _data = view(red_mat,i,j,:);
                _eata = view(nir_mat,i,j,:);
                _mask = (_data .> mm[1]) .* (_data .< mm[2]) .*
                        (_eata .> mm[1]) .* (_eata .< mm[2]);
                if sum(_mask) > 0
                    (dat_out).lat[dim_lat*(i-1)+j] = lat_mat[i,j];
                    (dat_out).lon[dim_lat*(i-1)+j] = lon_mat[i,j];
                    (dat_out).red[dim_lat*(i-1)+j] = nanmean(_data[_mask]);
                    (dat_out).nir[dim_lat*(i-1)+j] = nanmean(_eata[_mask]);
                end
            end
        catch e
            @warn "An error occurred when reading the dataset, use NaN...";
            @show e;
            (dat_out).lay = ones(Int, 1) * nlayer;
            (dat_out).lat = ones(1) .* -9999;
            (dat_out).lon = ones(1) .* -9999;
            (dat_out).red = ones(1) .* -9999;
            (dat_out).nir = ones(1) .* -9999;
        end

        # save data
        mask = ((dat_out).red .> -1) .* ((dat_out).nir .> -1);
        save_csv!(outp_file, dat_out[mask,:]);
    end

    return nothing
end

function grid_RAW!(
            dt::MOD15A2Hv006LAI{FT},
            nlayer::Int,
            path::String,
            prefix::String,
            band::String,
            cache::String,
            mm::Array{FT,1}
) where {FT<:AbstractFloat}
    # transform the data to CSV file
    pos_start = findfirst(prefix, path);
    hv_string = path[pos_start[1]:end];
    outp_file = joinpath(cache, hv_string * ".csv");

    # if file exists, skip
    if !isfile(outp_file)
        h,v = parse_HV(hv_string);
        lat_mat,lon_mat = load_MODIS(dt, h, v);
        dat_out = DataFrame();
        try
            lai_mat = read_nc(FT, path, band);
            dim_lon = size(lai_mat,1);
            dim_lat = size(lai_mat,2);
            (dat_out).lay = ones(Int, dim_lon * dim_lat) * nlayer;
            (dat_out).lat = ones(dim_lon * dim_lat) .* -9999;
            (dat_out).lon = ones(dim_lon * dim_lat) .* -9999;
            (dat_out).lai = ones(dim_lon * dim_lat) .* -9999;
            for i in 1:dim_lon, j in 1:dim_lat
                _data = view(lai_mat,i,j,:);
                _mask = (_data .> mm[1]) .* (_data .< mm[2]);
                if sum(_mask) > 0
                    (dat_out).lat[dim_lat*(i-1)+j] = lat_mat[i,j];
                    (dat_out).lon[dim_lat*(i-1)+j] = lon_mat[i,j];
                    (dat_out).lai[dim_lat*(i-1)+j] = nanmean(_data[_mask]);
                end
            end
        catch e
            @warn "An error occurred when reading the dataset, use NaN...";
            @show e;
            (dat_out).lay = ones(Int, 1) * nlayer;
            (dat_out).lat = ones(1) .* -9999;
            (dat_out).lon = ones(1) .* -9999;
            (dat_out).lai = ones(1) .* -9999;
        end

        # save data
        mask = ((dat_out).lai .> -1);
        save_csv!(outp_file, dat_out[mask,:]);
    end

    return nothing
end

function grid_RAW!(
            dt::AbstractUngriddedData{FT},
            params::Array{Any,1}
) where {FT<:AbstractFloat}
    # run the gridding process
    @info "Gridding RAW data..."
    @showprogress for param in params
        grid_RAW!(param);
    end

    return nothing
end

function grid_RAW!(
            dt::AbstractUngriddedData{FT},
            params::Array{Any,1},
            nthread::Int
) where {FT<:AbstractFloat}
    # create threads
    dynamic_workers!(nthread);

    # run the gridding process
    @info "Gridding RAW data..."
    @showprogress pmap(grid_RAW!, params);

    return nothing
end


text:
using Optim, BlackBoxOptim, NLopt, PyCall

import BlackBoxOptimizationBenchmarking: minimizer, minimum, optimize
import Base.string

box(D) = fill((-5.5, 5.5), D)
pinit(D) = 10*rand(D).-5

# Define optimize, minimum and minimizer for each optimizer

## NLopt

    mutable struct NLoptOptimMethod
        s::Symbol
    end
    string(opt::NLoptOptimMethod) = string("NLopt.", opt.s)

    function optimize(NLmeth::NLoptOptimMethod, f, D, run_length)
        opt = Opt(NLmeth.s, D)
        min_objective!(opt, (p, g) -> f(p) ) #NLopt expect gradient
        maxeval!(opt, run_length)
        xtol_abs!(opt, 1e-12)
        xtol_rel!(opt, 1e-12)
        ftol_abs!(opt, 1e-12)
        ftol_rel!(opt, 1e-12)
        lower_bounds!(opt, -5.5*ones(D))
        upper_bounds!(opt, +5.5*ones(D))
        minf, minx, ret = NLopt.optimize(opt, pinit(D))
        return NLmeth, minx, minf
    end
    minimum(mfit::Tuple{NLoptOptimMethod, Vector{Float64}, Float64}) = mfit[3]
    minimizer(mfit::Tuple{NLoptOptimMethod, Vector{Float64}, Float64}) = mfit[2]

## chain

    mutable struct Chain{T, K}
        first::T
        second::K
        p::Float64
    end

    function optimize(m::Chain, f, D, run_length)
        rl1 = round(Int, m.p*run_length)
        rl2 = run_length - rl1

        mfit = optimize(m.first, f, D, run_length)
        xinit = minimizer(mfit)
        mfit = optimize(m.second, f, D, run_length, xinit)
    end

    string(m::Chain) = string(string(m.first), " → ", string(m.second))

## python cma

#try

    @pyimport cma

    struct PyCMA
    end

    function optimize(m::PyCMA, f, D, run_length)
        es = cma.CMAEvolutionStrategy(pinit(D), 3, Dict("verb_log"=>0, "verb_disp"=>0, "maxfevals"=>run_length))
        mfit = es.optimize(f).result
        (m, mfit[1], mfit[2])
    end

    string(m::PyCMA) = "PyCMA"
    minimum(mfit::Tuple{PyCMA, Vector{Float64}, Float64}) = mfit[3]
    minimizer(mfit::Tuple{PyCMA, Vector{Float64}, Float64}) = mfit[2]

## scipy

    @pyimport scipy.optimize as scipy_opt

    struct PyMinimize
        method::String
    end

    optimize(m::PyMinimize, f, D, run_length) = (m, scipy_opt.minimize(
        f, pinit(D), method=m.method,
        options = Dict(
            "maxfev"=>run_length, "xatol"=>1e-8, "fatol"=>1e-8,
            "maxiter"=>run_length, "gtol"=>1e-12,
        )
    ))
    minimum(mfit::Tuple{PyMinimize, Dict{Any, Any}}) = mfit[2]["fun"]
    minimizer(mfit::Tuple{PyMinimize, Dict{Any, Any}}) = mfit[2]["x"]

    string(m::PyMinimize) = string("Py.", m.method)

#end

## Optim

    optimize(opt::Optim.AbstractOptimizer, f, D, run_length) =
        Optim.optimize(f, pinit(D), opt, Optim.Options(f_calls_limit=run_length, g_tol=1e-120, iterations=run_length))

    optimize(opt::Optim.AbstractOptimizer, f, D, run_length, xinit) =
        Optim.optimize(f, xinit, opt, Optim.Options(f_calls_limit=run_length, g_tol=1e-120, iterations=run_length))

    optimize(opt::Optim.SAMIN, f, D, run_length) =
        Optim.optimize(
            f, fill(-5.5, D), fill(5.5, D), pinit(D), opt,
            Optim.Options(f_calls_limit=run_length, g_tol=1e-120, iterations=run_length, show_trace=false)
        )

    string(opt::Optim.AbstractOptimizer) = string(typeof(opt).name)

    # Optim with restart
    try
        mutable struct OptimRestart{T}
            opt::T
        end
    catch
    end

    function optimize(opt::OptimRestart, f, D, run_length)
        fits = [optimize(opt.opt, f, D, round(Int, run_length/20)) for i=1:20]
        mins = [minimum(fit) for fit in fits]
        fits[argmin(mins)]
    end

    string(opt::OptimRestart) = string("Restart-", string(opt.opt))

## BlackBoxOptim

    mutable struct BlackBoxOptimMethod
        s::Symbol
    end
    string(opt::BlackBoxOptimMethod) = string("BBO.", opt.s)

    optimize(method::BlackBoxOptimMethod, f, D, run_length) =
        bboptimize(f; SearchRange=box(D), NumDimensions=D, Method=method.s, MaxFuncEvals=run_length, TraceMode=:silent)

    minimum(mfit::BlackBoxOptim.OptimizationResults) = best_fitness(mfit)
    minimizer(mfit::BlackBoxOptim.OptimizationResults) = best_candidate(mfit)


text:
"""
    mol_helmholtz(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_a`, `:a`

Compute the molar helmholtz energy (J/mol), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_helmholtz end

"""
    mass_helmholtz(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_a`

Compute the specific (mass) helmholtz energy (J/kg), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_helmholtz end

"""
    total_helmholtz(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_a`

Compute the total helmholtz energy (J), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_helmholtz end

#gibbs energy

"""
    mol_gibbs(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_g`, `:g`

Compute the molar gibbs energy (J/mol), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_gibbs end

"""
    mass_gibbs(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_g`

Compute the specific (mass) gibbs energy (J/kg), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_gibbs end

"""
    total_gibbs(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_g`

Compute the total gibbs energy (J), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_gibbs end

#enthalpy

"""
    mol_enthalpy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_h`, `:h`

Compute the molar enthalpy (J/mol), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_enthalpy end

"""
    mass_enthalpy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_h`

Compute the specific (mass) enthalpy (J/kg), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_enthalpy end

"""
    total_enthalpy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_h`

Compute the total enthalpy (J), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_enthalpy end

#internal energy

"""
    mol_internal_energy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_u`, `:u`

Compute the molar internal energy (J/mol), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_internal_energy end

"""
    mass_internal_energy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_u`

Compute the specific (mass) internal energy (J/kg), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_internal_energy end

"""
    total_internal_energy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_u`

Compute the total internal energy (J), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_internal_energy end

#entropy

"""
    mol_entropy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_s`, `:s`

Compute the molar entropy (J/(mol K)), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_entropy end

"""
    mass_entropy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_s`

Compute the specific (mass) entropy (J/(kg K)), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_entropy end

"""
    total_entropy(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_s`

Compute the total entropy (J/K), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_entropy end

#temperature an pressure
"""
    temperature(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:t`, `:T`

Compute the temperature (K), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function temperature end

"""
pressure(model,state::ThermodynamicState,unit,[options...])::Real

    Keyword symbols: `:p`, `:P`

Compute the pressure (Pa), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function pressure end

#volume amounts

"""
    mol_volume(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_v`, `:v`

Compute the molar volume (m³/mol), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_volume end

"""
    mass_volume(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_v`

Compute the specific (mass) volume (m³/kg), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_volume end

"""
    total_volume(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:total_v`, `:V`

Compute the total volume (m³), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function total_volume end

"""
    mol_density(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mol_rho`, `:rho`, `:mol_ρ`, `:ρ`

Compute the molar density (mol/m³), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_density end

"""
    mass_density(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass_rho`, `:mass_ρ`

Compute the specific (mass) density (kg/m³), from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_density end

"""
    mol_fraction(model,state::ThermodynamicState,unit,[options...])::AbstractVector

Keyword symbols: `:xn`

Compute the molar fraction of each present compound present in the provided system, from or using the specifications as base.

The unit argument is not used here. but its present for consistency.
"""
function mol_fraction end

"""
    mass_fraction(model,state::ThermodynamicState,unit,[options...])::AbstractVector

Keyword symbols: `:xm`

Compute the mass fraction of each present compound present in the provided system, from or using the specifications as base.

The unit argument is not used here. but its present for consistency.
"""
function mass_fraction end

"""
    mol_number(model,state::ThermodynamicState,unit,[options...])::AbstractVector

Keyword symbols: `:n`

Compute the amount of moles (in mol) of each present compound present in the provided system, from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mol_number end

"""
    mass_number(model,state::ThermodynamicState,unit,[options...])::AbstractVector

Keyword symbols: `:m`

Compute the amount of mass (in kg) of each present compound present in the provided system, from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass_number end

"""
    moles(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:moles`

Compute the total amount of moles (in mol) of the provided system, from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function moles end

"""
    mass(model,state::ThermodynamicState,unit,[options...])::Real

Keyword symbols: `:mass`

Compute the total amount of mass (in kg) of the provided system, from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function mass end

"""
    molar_mass(model,state::ThermodynamicState,unit,[options...])::Real

Compute the amount of mass per unit of moles (in kg/mol) of the provided state, from or using the specifications as base.

If you need unitful units, use the macro `@to_units` before the function call

"""
function molar_mass end

"""
    options(model,state::ThermodynamicState)::NamedTuple

Keyword symbols: `:options`

Returns any options related to the model, stored in the state or the model.

Returns an empty named tuple otherwise.

"""
function options end

"""
    phase(model,state::ThermodynamicState)::NamedTuple

Keyword symbols: `:phase`

Returns a symbol indicating the phase specified in the state or the model.

Returns `:unspecified` otherwise.

"""
function phase end

"""
    quality(model,state::ThermodynamicState)::NamedTuple

Keyword symbols: `:quality`

Returns the vapor quality stored in the state

Returns NaN otherwise.

This function is an alias for `mass_vapor_fraction`.

"""
function quality end

"""
    molecular_weight(model)

Keyword symbols: `:mw`

Returns the molecular weight stored in a model

Returns `nothing` if not defined.

it should return a list of the molecular weights. in g/mol. You can return a single number if the model has only one component.

"""
function molecular_weight end


text:
__precompile__()

module FMIndexes

export
    FMIndex,
    restore,
    count,
    locate,
    locateall

using SuffixArrays
using WaveletMatrices
using IndexableBitVectors
using Humanize
using Mmap

"""
Index for full-text search.

Type parameters:

* `w`: the number of bits required to encode the alphabet
* `T`: the type to represent positions of a sequence
"""
struct FMIndex{w,T}
    bwt::WaveletMatrix{w,UInt8,SucVector}
    sentinel::Int
    samples::Vector{T}
    sampled::SucVector
    count::Vector{Int}
end

function FMIndex(seq, sa, σ, r)
    wm = WaveletMatrix(make_bwt(seq, sa), log2(Int, σ))
    # sample suffix array
    samples, sampled = sample_sa(sa, r)
    sentinel = something(findfirst(isequal(0), sa), 0) + 1
    # count characters
    count = count_bytes(seq, σ)
    count[1] = 1  # sentinel '$' is smaller than any character
    cumsum!(count, count)
    return FMIndex(wm, sentinel, samples, SucVector(sampled), count)
end

"""
    FMIndex(seq, σ=256; r=32, program=:SuffixArrays, mmap::Bool=false, opts...)

Build an FM-Index from a sequence `seq`.
The sequence must support `convert(UInt8, seq[i])` for each character and the
alphabet size should be less than or equal to 256. The second parameter, `σ`, is
the alphabet size. The third parameter, `r`, is the interval of sampling values
from a suffix array. If you set it large, you can save the memory footprint but
it requires more time to locate the position.
"""
function FMIndex(seq, σ=256; r=32, program=:SuffixArrays, mmap::Bool=false, opts...)
    T = index_type(length(seq))
    opts = Dict(opts)
    @assert !haskey(opts, :σ) "σ should be passed as the second argument"
    if program === :SuffixArrays
        @assert 1 ≤ σ ≤ typemax(UInt8) + 1
        sa = make_sa(T, seq, σ, mmap)
    elseif program === :psascan || program === :pSAscan
        @assert 1 ≤ σ ≤ typemax(UInt8)
        psascan = get(opts, :psascan, "psascan")
        workdir = get(opts, :workdir, pwd())
        sa = make_sa_pscan(T, seq, psascan, workdir, mmap)
    else
        error("unknown program name: $program")
    end
    return FMIndex(seq, sa, σ, r)
end

"""
    FMIndex(text::String; opts...)

Build an FM-Index from an ASCII text.
"""
function FMIndex(text::String; opts...)
    return FMIndex(convert(Vector{UInt8}, text), 128; opts...)
end

Base.length(index::FMIndex) = length(index.bwt)

function Base.show(io::IO, fmindex::FMIndex{w,T}) where {w,T}
    println(io, summary(fmindex), ':')
    totalsize = (
        sizeof(fmindex.bwt) +
        sizeof(fmindex.samples) +
        sizeof(fmindex.sampled) +
        sizeof(fmindex.count)
    )
    print("     length: ", length(fmindex), '\n')
    print("  data size: ", Humanize.datasize(totalsize, style=:bin))
end

"""
Restore the original text from the index.
"""
function restore(index::FMIndex)
    n = length(index)
    text = Vector{UInt8}(undef, n)
    p = index.sentinel
    while n > 0
        p = lfmap(index, p)
        text[n] = index.bwt[p ≥ index.sentinel ? p - 1 : p]
        n -= 1
    end
    return text
end

"""
Count the number of occurrences of the given query.
"""
function Base.count(query, index::FMIndex)
    return length(sa_range(query, index))
end

"""
Locate the positions of occurrences of the given query.
This method returns an iterator of positions:

    for pos in locate(query, index)
        # ...
    end
"""
function locate(query, index::FMIndex)
    return LocationIterator(sa_range(query, index), index)
end

"""
Locate the positions of all occurrences of the given query.
"""
function locateall(query, index::FMIndex)
    iter = locate(query, index)
    locs = Vector{Int}(undef, length(iter))
    for (i, loc) in enumerate(iter)
        locs[i] = loc
    end
    return locs
end

include("index.jl")
include("saca.jl")
include("lociter.jl")

end # module


text:
@testset "PB-representation - creation" begin
    ode = @ODEmodel(
        x1'(t) = x3(t),
        x2'(t) = a * x2(t),
        x3'(t) = x1(t),
        y1(t) = x1(t),
        y2(t) = x2(t) + u(t)
    )
    ioeqs = find_ioequations(ode)
    pbr = PBRepresentation(ode, ioeqs)

    R, (y1_0, y1_1, y1_2, y2_0, y2_1, y2_2) = Nemo.PolynomialRing(Nemo.QQ, ["y1_0", "y1_1", "y1_2", "y2_0", "y2_1", "y2_2"])
    @test find_leader([a, x1, x2, u], pbr) == nothing
    @test find_leader([a, x1, y1_0, x2], pbr) == y1_0
    @test find_leader([y1_0, y1_1, y2_0, y2_1], pbr) == y2_1
    @test find_leader([y1_2, y2_1], pbr) == y2_1
    @test find_leader([y1_2, y2_0], pbr) == y1_2
end


text:
using Test

import Pkg
Pkg.add("ToeplitzMatrices")

using ToeplitzMatrices: Toeplitz
using SparseArrays: sparse
using LinearAlgebra: triu, diag, diagm, Diagonal, Bidiagonal, Tridiagonal, det, I

include("../../lib/KTBC.jl")
using .KTBC: CreateKTBC

# Difference matrix
K₄ = Toeplitz([2, -1, 0, 0], [2, -1, 0, 0])
K₁₀ = Toeplitz([2, -1, zeros(10)...], [2, -1, zeros(10)...])

# T - Top boundary condition
RectangularT = Toeplitz([1, -1, 0, 0], [1, 0, 0])

@test RectangularT == [
  1.0   0.0   0.0
 -1.0   1.0   0.0
  0.0  -1.0   1.0
  0.0   0.0  -1.0
]

# Circulant C, periodic
C₄ = Toeplitz([2, -1, 0, -1], [2, -1, 0, -1])

@test C₄ == [
  2.0  -1.0   0.0  -1.0
 -1.0   2.0  -1.0   0.0
  0.0  -1.0   2.0  -1.0
 -1.0   0.0  -1.0   2.0
]

nullspaceVector = C₄ * ones(4)

@test nullspaceVector == [0.0; 0.0; 0.0; 0.0]

# Sparse version of K₄
sparse(K₄)

invK₄ = round.(inv(convert(Array{Float64}, K₄)); digits=1)

@test invK₄ == [
 0.8  0.6  0.4  0.2
 0.6  1.2  0.8  0.4
 0.4  0.8  1.2  0.6
 0.2  0.4  0.6  0.8
]

U = triu(ones(3, 3))

@test U == [
 1.0  1.0  1.0
 0.0  1.0  1.0
 0.0  0.0  1.0
]

function createDiagonal(value=1, dim=4, isMatrix=true)
    DiagonalMatrix = zeros(Int, dim, dim) + I * value
    return isMatrix ? DiagonalMatrix : diag(DiagonalMatrix)
end

# Not memory efficient method
function createKWithMemoryProblem(dim=4)
    DiagonalMatrix = createDiagonal(2, dim)
    E = diagm(1 => -ones(dim - 1))
    K = DiagonalMatrix + E + E'
end

@test createKWithMemoryProblem() == [
  2.0  -1.0   0.0   0.0
 -1.0   2.0  -1.0   0.0
  0.0  -1.0   2.0  -1.0
  0.0   0.0  -1.0   2.0
]

# Very efficient
function createK(dim=4)
    ULDiagonals = ones(Int, dim - 1)
    Tridiagonal(ULDiagonals, createDiagonal(2, dim, false), ULDiagonals)
end

convert(Array{Int}, createK())

@test convert(Array{Int}, createK()) == [
 2  1  0  0
 1  2  1  0
 0  1  2  1
 0  0  1  2
]

# Convert to memory-efficient
Tridiagonal(K₁₀)

K, T, B, C = CreateKTBC(8)

@test [det(K), det(T), det(B), round(det(C))] == [9.0; 1.0; 0.0; 0.0]

# Worked examples
e = ones(8)
u = [1:8...]

@test B * e == [
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
]

f = B * u

@test f' * e == 0.0

# First way to make fixed-free
H = convert(Array{Int}, Toeplitz([2, -1, 0], [2, -1, 0]))
H[3, 3] = 1

@test H == [
  2  -1   0
 -1   2  -1
  0  -1   1
]

# Second way
J = reverse(convert(Array{Int}, I(8)), dims=2)
T = convert(Array{Int}, T)
H = J * T * J

# Inverse ways
invH = round.(Int, inv(H))
invHUsingT = round.(Int, J * inv(T) * J)

@test invH == invHUsingT


text:
using Heuristics, Suppressor

include("testProblems.jl")

# Sphere function test 1
begin
    N = 14
    M = 200
    LB = -5 .* ones(N)
    UB = 5 .* ones(N)
    prob1 = Problem(sphereFunc, LB, UB)
    opts = Options(;display = false, useParallel=true)
    pso1 = MS_PSO(prob1; numParticlesPerSwarm = M)

    res = optimize!(pso1, opts)

    @test res.fbest <= 1.0e-5

    for i in 1:N
        @test res.xbest[i] >= LB[i]
        @test res.xbest[i] <= UB[i]
    end
end

# Sphere function test 2
begin
    N = 2
    M = 1000
    prob2 = Problem(sphereFunc, N)
    opts = Options(;display = true, useParallel = true)
    pso2 = MS_PSO(prob2; numParticlesPerSwarm = M, initMethod = :LogisticsMap)

    @suppress_out begin
        optimize!(pso2, opts)
    end

    @test pso2.results.fbest <= 1.0e-5

    for i in 1:N
        @test pso2.results.xbest[i] >= LB[i]
        @test pso2.results.xbest[i] <= UB[i]
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StorageAccountInternetEndpoints <: SwaggerModel
    blob::Any # spec type: Union{ Nothing, String } # spec name: blob
    file::Any # spec type: Union{ Nothing, String } # spec name: file
    web::Any # spec type: Union{ Nothing, String } # spec name: web
    dfs::Any # spec type: Union{ Nothing, String } # spec name: dfs

    function StorageAccountInternetEndpoints(;blob=nothing, file=nothing, web=nothing, dfs=nothing)
        o = new()
        validate_property(StorageAccountInternetEndpoints, Symbol("blob"), blob)
        setfield!(o, Symbol("blob"), blob)
        validate_property(StorageAccountInternetEndpoints, Symbol("file"), file)
        setfield!(o, Symbol("file"), file)
        validate_property(StorageAccountInternetEndpoints, Symbol("web"), web)
        setfield!(o, Symbol("web"), web)
        validate_property(StorageAccountInternetEndpoints, Symbol("dfs"), dfs)
        setfield!(o, Symbol("dfs"), dfs)
        o
    end
end # type StorageAccountInternetEndpoints

const _property_map_StorageAccountInternetEndpoints = Dict{Symbol,Symbol}(Symbol("blob")=>Symbol("blob"), Symbol("file")=>Symbol("file"), Symbol("web")=>Symbol("web"), Symbol("dfs")=>Symbol("dfs"))
const _property_types_StorageAccountInternetEndpoints = Dict{Symbol,String}(Symbol("blob")=>"String", Symbol("file")=>"String", Symbol("web")=>"String", Symbol("dfs")=>"String")
Base.propertynames(::Type{ StorageAccountInternetEndpoints }) = collect(keys(_property_map_StorageAccountInternetEndpoints))
Swagger.property_type(::Type{ StorageAccountInternetEndpoints }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_StorageAccountInternetEndpoints[name]))}
Swagger.field_name(::Type{ StorageAccountInternetEndpoints }, property_name::Symbol) =  _property_map_StorageAccountInternetEndpoints[property_name]

function check_required(o::StorageAccountInternetEndpoints)
    true
end

function validate_property(::Type{ StorageAccountInternetEndpoints }, name::Symbol, val)
end


text:
export FFT, fft!

function _clp2(x::Int)
    x == 0 && return 1
    x == 1 && return 2
    x = x - 1
    x = x | (x >> 1)
    x = x | (x >> 2)
    x = x | (x >> 4)
    x = x | (x >> 8)
    x = x | (x >> 16)
    return x + 1
end

function _log2(x::Int)
    x == 0 && error("_log2(x): x cannot be zero!")
    r = 0
    x > 4294967295 && (x >>= 32; r += 32)
    x >      65535 && (x >>= 16; r += 16)
    x >        255 && (x >>=  8; r +=  8)
    x >         15 && (x >>=  4; r +=  4)
    x >          3 && (x >>=  2; r +=  2)
    x >          1 && (x >>=  1; r +=  1)
    return r
end

struct FFT
    cxy::Matrix{Float64}
    crθ::Matrix{Float64}
    fac::Matrix{Float64}
    frq::Vector{Float64}

    function FFT(N::Int)
        cxy = Matrix(undef, N, 2)
        crθ = Matrix(undef, N, 2)
        frq = Vector(undef, N)
        return new(cxy, crθ, twiddle(N >> 1), frq)
    end
end

function twiddle(sz::Int)
    ws = Matrix{Float64}(undef, sz, 2)
    hz = sz >> 1
    ix = 1
    rc = 1.0
    rs = 0.0
    cf = cospi(inv(sz))
    sf = sinpi(inv(sz))
    #### initialize: θ = 0
    @inbounds ws[ix,1] = rc
    @inbounds ws[ix,2] = rs
    @inbounds while ix < hz
        ix += 1
        tc  = cf * rc - sf * rs
        ts  = sf * rc + cf * rs
        rc  = tc
        rs  = ts
        ws[ix,1] = rc
        ws[ix,2] = rs
    end
    #### specific case: θ = π / 2
    ix += 1
    @inbounds ws[ix,1] = 0.0
    @inbounds ws[ix,2] = 1.0
    #### remaining cases: θ > π / 2
    jx = ix
    @inbounds while ix < sz
        ix += 1
        jx -= 1
        ws[ix,1] = -ws[jx,1]
        ws[ix,2] =  ws[jx,2]
    end
    return ws
end

function fft!(fft::FFT, sig::VecI, timestep::Real)
    cxy = fft.cxy
    crθ = fft.crθ
    toN = eachindex(sig)

    N = length(toN)
    H = N >> 1

    if isone(_log2(H) & 1)
        @inbounds for i in toN
            cxy[i,1] = sig[i]
            cxy[i,2] = 0.0
        end
        difnn!(cxy, crθ, fft.fac, H)
    else
        @inbounds for i in toN
            crθ[i,1] = sig[i]
            crθ[i,2] = 0.0
        end
        difnn!(crθ, cxy, fft.fac, H)
    end

    fftshift!(cxy)

    #### Phases unwrapping
    halfπ = 0.5 * π
    @inbounds crθ[1,2] = prevθ = atan(cxy[1,2], cxy[1,1])
    @inbounds crθ[1,1] = 0.0
    for i in 2:N
        @inbounds crθ[i,2] = thisθ = atan(cxy[i,2], cxy[i,1])
        diff  = thisθ - prevθ
        θmod  = rem2pi(diff + halfπ, RoundDown) - halfπ
        diff  > 0.0 && θmod == -halfπ && (θmod = halfπ)
        @inbounds crθ[i,1] = abs(diff) < halfπ ? crθ[i-1,1] : θmod - diff + crθ[i-1,1]
        prevθ = thisθ
    end
    @simd for i in toN
        @inbounds crθ[i,2] += crθ[i,1]
    end
    @inbounds refθ = crθ[H+1,2]
    @simd for i in toN
        @inbounds crθ[i,2] -= refθ
    end
    #### Compute amplitudes
    @simd for i in toN
        @inbounds crθ[i,1] = apy2(cxy[i,1], cxy[i,2]) / H
    end
    #### Generate FFT frequency coordinate, frq[H+1] = 0
    fftfreq!(fft.frq, inv(timestep * N))
    return nothing
end

function butterfly!(y::MatO, x::MatI, w::MatI, jx::Int, kx::Int, wx::Int, H::Int, h::Int, S::Int, s::Int)
    @inbounds for _ in eachindex(1:h)
        jy = jx + s

        wr = w[wx, 1]
        wi = w[wx, 2]
        xr = x[kx, 1]
        xi = x[kx, 2]
        Xr = x[kx + H, 1]
        Xi = x[kx + H, 2]

        y[jx, 1] = xr + Xr
        y[jx, 2] = xi + Xi
        y[jy, 1] = wr * (xr - Xr) - wi * (xi - Xi)
        y[jy, 2] = wr * (xi - Xi) + wi * (xr - Xr)

        jx += S
        kx += s
        wx += s
    end
end

function difnn!(a::MatI, b::MatI, w::MatI, H::Int)
    h = H
    s = 1
    S = 2
    r = false

    while h > 0
        for ix in eachindex(1:s)
            r ? butterfly!(a, b, w, ix, ix, 1, H, h, S, s) : butterfly!(b, a, w, ix, ix, 1, H, h, S, s)
        end

        h >>= 1
        s <<= 1
        S <<= 1
        r = !r
    end
    return nothing
end

function fftshift!(x::MatI)
    N = size(x, 1) >> 1
    @inbounds for j in axes(x, 2), i in eachindex(1:N)
        swap!(x, i, j, i+N, j)
    end
end

function fftfreq!(f::VecI, Δf::Real)
    nhalfp1 = length(f) >> 1 + 1
    @simd for i in eachindex(f)
        @inbounds f[i] = Δf * (i - nhalfp1)
    end
    return nothing
end


text:
# Compatible across multiple version of Julia
# (BLAS.get_num_threads() was only introduced in Julia 1.6)
blas_get_num_threads() = Compat.get_num_threads()


text:
defaultchkfn = (ismissing, isnothing, isnan);
# TODO: no test for filtern.jl
"""
`filtern` works based on `filter` that drops rows where the corresponding specific column `col1` contains `NaN`, `missing` or `nothing` by default.

`filtern` is similar to `dropmissing` but you can drop anything by specifying functions using kwarg `chkfun`, e.g. `filtern(...; chkfun=(ismissing, x -> x>1))`.

`filtern(col1, df)` is the shorthand of `filter(cols => fun, df::AbstractDataFrame)` where `fun` is for examing whether the column contains `NaN`, `missing` or `nothing`.
"""
function filtern!(col, df; chkfun=defaultchkfn)
    # see this thread: https://stackoverflow.com/questions/62789334/how-to-remove-drop-rows-of-nothing-and-nan-in-julia-dataframe
    return filter!(col => c1 -> !any(f -> f(c1), chkfun), df)
end

"""
`filtern(col1, col2, df)` automatically drop rows where either of `col1` and `col2`  contains `NaN`, `missing` or `nothing` by default.
"""
function filtern!(col1, col2, df; chkfun=defaultchkfn)
    return filter!([col1, col2] => (c1, c2) -> !any(f -> f(c1) || f(c2), chkfun), df)
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct LogUniform <: AbstractDistDelayed
    low::NestedFloat
    high::NestedFloat
end

function loguniform(low::Float64, high::Float64)::Float64
    return exp(uniform(low, high))
end


text:
using CuArrays

A = [1 2
     3 4
     5 6]
x = [1
     2]

@show A*x

A=cu(A)
x=cu(x)

@show A*x

x = rand(10)

x = x |> gpu
@show x,typeof(x)
x = x |> cpu
@show x,typeof(x)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Spectra
using Test

# Dummy data
x = collect(0.:1.:100.)
y = collect(0.:1.:100.)
x2 = collect(-100.:0.1:100.)
y_poly4 = 1.0 .+ x2 .+ 0.25.*x2.^2 + 0.0001.*x2.^3 + 1.e-6.*x2.^4
y_peak = collect(0.:1.:100.) .+ 100.0.*exp.(-log(2).*((x.-50.0)./5.0).^2) .+ 0.2.*randn(size(x,1))
roi = [0.0 40;60. 100.0]

# testing baselines
y_fit_dierckx,bas_fit_dierckx = baseline(x[:],y,[0. 100.],"unispline",s=0.01)
@test isapprox(y_fit_dierckx,(y-y[:]),atol=1e-4)
@test isapprox(bas_fit_dierckx,y[:],atol=1e-4)

y_fit_linear,bas_fit_linear = baseline(x[:],y,[0. 100.],"poly",polynomial_order=2)
@test isapprox(y_fit_linear,(y-y[:]),atol=1e-4)
@test isapprox(bas_fit_linear,y[:],atol=1e-4)

y_fit_poly,bas_fit_poly = baseline(x2[:],y_poly4[:],[0. 100.],"poly",polynomial_order=4)
@test isapprox(y_fit_poly,(y_poly4[:]-y_poly4[:]),atol=1e-3)
@test isapprox(bas_fit_poly,y_poly4[:],atol=1e-3)

y_fit_als,bas_fit_als = baseline(x[:],y_peak,roi,"als",lam=10^5,p=0.001)
@test isapprox(y_fit_als,(y_peak-y[:]),atol=10.)
@test isapprox(bas_fit_als,y[:],atol=10.)

y_fit_arPLS,bas_fit_arPLS = baseline(x[:],y_peak,roi,"arPLS",lam=10^9,p=0.001)
@test isapprox(y_fit_arPLS,(y_peak-y[:]),atol=10.)
@test isapprox(bas_fit_arPLS,y[:],atol=10.)

#y_fit_whit,bas_fit_whit = baseline(x[:],y_peak[:],roi,"whittaker",lam=10^9)
#@test isapprox(y_fit_whit,(y_peak-y[:]),10.)
#@test isapprox(bas_fit_whit,y[:],10.)


text:
@testset "Regularizer" begin

    # Test Regularizer
    T = Float64
    reg = Regularizer{T}()
    val = 1e-1
    set!(reg, val)

    @test reg.x == val
    @test reg.u == val
    @test reg.λ == val

    mult!(reg, val)
    @test reg.x == val^2
    @test reg.u == val^2
    @test reg.λ == val^2

end


text:
# ---
# title: 73. Set Matrix Zeroes
# id: problem73
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/set-matrix-zeroes/description/>
# hidden: true
# ---
#
# Given an ` _m_ x _n_` matrix. If an element is **0** , set its entire row and
# column to **0**. Do it [**in-place**](https://en.wikipedia.org/wiki/In-
# place_algorithm).
#
# **Follow up:**
#
#   * A straight forward solution using O( _m_ _n_ ) space is probably a bad idea.
#   * A simple improvement uses O( _m_ \+ _n_ ) space, but still not the best solution.
#   * Could you devise a constant space solution?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/17/mat1.jpg)
#
#
#
#     Input: matrix = [[1,1,1],[1,0,1],[1,1,1]]
#     Output: [[1,0,1],[0,0,0],[1,0,1]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/17/mat2.jpg)
#
#
#
#     Input: matrix = [[0,1,2,0],[3,4,5,2],[1,3,1,5]]
#     Output: [[0,0,0,0],[0,4,5,0],[0,3,1,0]]
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[0].length`
#   * `1 <= m, n <= 200`
#   * `-231 <= matrix[i][j] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end