text:
# cd("/Users/kevintracy/.julia/dev/CPEG")
# Pkg.activate(".")
# cd("/Users/kevintracy/.julia/dev/CPEG/src")
# using LinearAlgebra
# using StaticArrays
# using ForwardDiff
# using SparseArrays
# using SuiteSparse
# using Printf
#
# include("qp_solver.jl")
# include("atmosphere.jl")
# include("scaling.jl")
# include("gravity.jl")
# include("aero_forces.jl")
# include("vehicle.jl")
# include("dynamics.jl")
# include("post_process.jl")

using LinearAlgebra
using StaticArrays
using ForwardDiff
using SparseArrays
using SuiteSparse
using MATLAB

function eg_mpc_quad(
    ev::CPEGWorkspace,
    A::Vector{SMatrix{7, 7, Float64, 49}},
    B::Vector{SMatrix{7, 1, Float64, 7}},
    X::Vector{SVector{7,Float64}})

    # sizes for state and control
    nx = 7
    nu = 1
    N = length(X)

    # indicees for state and control
    idx_x = [(i-1)*(nx+nu) .+ (1:nx) for i = 1:length(X)]
    idx_u = [(i-1)*(nx+nu) .+ nx .+ (1:nu) for i = 1:(length(X)-1)]

    # constraint jacobian (in this form L≦Az≤U)
    nz = (N*nx) + (N-1)*nu
    nc = (N-1)*nx + nx #+ (N)*1
    A_eq = spzeros(nc,nz)
    A_ineq = spzeros(N,nz) # state constraints

    # dynamics constraint (equality)
    idx_c = [(i-1)*(nx) .+ (1:nx) for i = 1:(N-1)]
    for i = 1:(N-1)
        A_eq[idx_c[i],idx_x[i]]   = A[i]
        A_eq[idx_c[i],idx_u[i]]   = B[i]
        A_eq[idx_c[i],idx_x[i+1]] = -I(nx)
    end
    A_eq[(N-1)*nx .+ (1:nx), idx_x[1]] = I(nx)

    # state constraints on δσ (inequality)
    for i = 1:N
        A_ineq[i,idx_x[i][7]] = 1
    end

    # constraint bounds
    dyn_eq = zeros((N)*nx) # N-1 dynamics constraints + N IC constraint
    up_tr = ev.cost.σ_tr*ones(N)
    low_tr = -ev.cost.σ_tr*ones(N)

    # cost function terms
    P = spzeros(nz,nz)
    q = zeros(nz)
    R = 0.01
    for i = 1:(N-1)
        P[idx_u[i],idx_u[i]] = [R]
        q[idx_u[i]] = [R*ev.U[i][1]]
    end

    # projection matrix onto the landing plane
    Qn = ev.cost.γ*landing_plane_proj(ev.cost.rf)

    # cost terms for terminal cost
    P[idx_x[N][1:3],idx_x[N][1:3]] = Qn'*Qn
    q[idx_x[N][1:3]] = -(Qn'*Qn)*(ev.cost.rf/ev.scale.dscale - X[N][1:3])

    # setup standard form QP
    P += 1e-6*I
    G = [A_ineq;-A_ineq]
    h = [up_tr;-low_tr]

    # solve the equality only constrained QP
    sol = lu([P A_eq';A_eq spzeros(N*nx,N*nx)])\[-q;dyn_eq]
    z = sol[1:length(q)]
    qp_iters = 1

    # if this violates the inequality constraints, then we send it to quadprog
    if sum(G*z .> h) != 0
        # solve with QP solver (src/qp_solver.jl)
        z, qp_iters = quadprog(P,q,A_eq,dyn_eq,G,h; verbose = ev.qp_solver_opts.verbose,
                                             atol = ev.qp_solver_opts.atol,
                                        max_iters = ev.qp_solver_opts.max_iters)
    end

    # pull out δx and δu from z
    δx = [z[idx_x[i]] for i = 1:(N)]
    δu = [z[idx_u[i]] for i = 1:(N-1)]

    # update U = U + δu
    ev.U += δu

    return norm(δu), qp_iters
end

function landing_plane_proj(rf)
    rr = normalize(rf)
    I - rr*rr'
end
function update_miss_distance!(ev::CPEGWorkspace,X::Vector{SVector{7,Float64}})
    ev.miss_distance = miss_distance(ev,X)
    return nothing
end
function miss_distance(ev::CPEGWorkspace,X::Vector{SVector{7,Float64}})
    Qn = landing_plane_proj(ev.cost.rf)
    rf_sim = X[end][1:3]*ev.scale.dscale
    norm(Qn*(rf_sim - ev.cost.rf))
end

function update_σ!(ev::CPEGWorkspace,X::Vector{SVector{7,Float64}})
    ev.σ = [X[i][7] for i = 1:length(X)]
end

function main_cpeg(ev, r, v, σ)

    rs, vs = scale_rv(ev.scale,r,v)
    x0_s = SVector{7}([rs;vs;σ])

    old_md = NaN
    new_md = NaN
    ndu = NaN
    qp_iters = NaN

    althist = [zeros(2) for i = 1:ev.max_cpeg_iter]
    drhist = [zeros(2) for i = 1:ev.max_cpeg_iter]
    crhist = [zeros(2) for i = 1:ev.max_cpeg_iter]
    σhist = [zeros(2) for i = 1:ev.max_cpeg_iter]

    if ev.verbose
        @printf "iter    miss (km)    Δmiss (km)     |ΔU|        N   QP iters\n"
        @printf "----------------------------------------------------------\n"
    end

    for i = 1:ev.max_cpeg_iter

        # rollout and check miss distance
        X = rollout(ev, x0_s)
        althist[i], drhist[i], crhist[i] = postprocess_scaled(ev,X,x0_s)
        σhist[i]=[X[i][7] for i = 1:length(X)]
        old_md = new_md
        new_md = miss_distance(ev,X)

        if ev.verbose && i>1
            @printf("%3d     %9.3e    %+9.3e    %9.3e   %3d  %3d\n",
                    i, new_md/1e3,(new_md - old_md)/1e3, ndu, length(X), qp_iters)
        end

        # check the flipping
        # flipping = (max(old_md,new_md)<5e3) & (i>5) & (qp_iters == 1)

        # check termination criteria
        if (abs(new_md - old_md) < ev.miss_distance_tol) | (ndu < ev.ndu_tol) | (i == ev.max_cpeg_iter)#| flipping
            update_σ!(ev,X)
            return althist[1:i], drhist[1:i], crhist[1:i], σhist[1:i]
        end

        # linearize and solve QP
        A,B = get_jacobians(ev,X)
        ndu, qp_iters = eg_mpc_quad(ev,A,B,X)
    end
    error("CPEG failed: reached max iters")
end

function scale_state(ev,r,v,σ)
    rs,vs = scale_rv(ev.scale,r,v)
    SVector{7}([rs;vs;σ])
end
function unscale_state(ev,x)
    rs = x[SA[1,2,3]]
    vs = x[SA[4,5,6]]
    r,v = unscale_rv(ev.scale,rs,vs)
    return r,v,x[7]
end
function tt()

        ev = CPEGWorkspace()

        Rm = ev.params.gravity.R
        r0 = SA[Rm+125e3, 0.0, 0.0] #Atmospheric interface at 125 km altitude
        V0 = 5.845e3 #Mars-relative velocity at interface of 5.845 km/sec
        γ0 = -15.474*(pi/180.0) #Flight path angle at interface
        v0 = V0*SA[sin(γ0), cos(γ0), 0.0]
        σ0 = deg2rad(3)

        Rf = Rm+10.0e3 #Parachute opens at 10 km altitude
        rf = Rf*SA[cos(7.869e3/Rf)*cos(631.979e3/Rf); cos(7.869e3/Rf)*sin(631.979e3/Rf); sin(7.869e3/Rf)]

        # vehicle parameters
        ev.params.aero.Cl = 0.29740410453983374
        ev.params.aero.Cd = 1.5284942035954776
        ev.params.aero.A = 15.904312808798327    # m²
        ev.params.aero.m = 2400.0                # kg

        # qp solver settings
        ev.qp_solver_opts.verbose = false
        ev.qp_solver_opts.atol = 1e-10
        ev.qp_solver_opts.max_iters = 50

        # MPC stuff
        ev.cost.σ_tr = deg2rad(10) # radians
        ev.cost.rf = rf # goal position (meters, MCMF)
        ev.cost.γ = 1e3

        # sim stuff
        ev.dt = 2.0 # seconds

        # CPEG settings
        ev.verbose = true
        ev.ndu_tol = 1e-3
        ev.max_cpeg_iter = 30
        ev.miss_distance_tol = 1e1  # m

        ev.scale.uscale = 1e1

        N = 100
        X = [@SVector zeros(7) for i = 1:N]
        U = [@SVector zeros(1) for i = 1:N-1]
        X[1] = scale_state(ev,r0,v0,σ0)
        # main_cpeg(ev,r0,v0,σ0)
        # main_cpeg(ev,unscale_state(ev,X[1])...)
        # @show unscale_state(ev,X[1])
        # @show r0
        # @show v0
        # @show σ0
        # for i = 1:N-1
            x0s = [3.394493027569473, 0.399591678638693, 0.001009766098441132, -1.456578872683738, 4.90997118202939, 0.11856550572023335, 1.002017922478615]
            ev.U = SVector{1}.([[-1.0446874855627493], [-1.5480207770930021], [-2.0792718023435577], [-2.6165103698995678], [-3.133332980329963], [-3.602272358062922], [-3.998882259324137], [-4.305331038330731], [-4.512508731782757], [-4.6202527616273725], [-4.635962501176615], [-4.5722798946790375], [-4.444568523571758], [-4.268731622216839], [-4.059636636270206], [-3.8301827417699306], [-3.590907860816463], [-3.349978236653403], [-3.1134079271625916], [-2.8853876478909783],
             [-2.6686409267757703], [-2.4647591748467277], [-2.274492048866316], [-2.0979855378265144], [-1.9349693968414852], [-1.7849000190375819], [-1.6470663824705791], [-1.5206666365829429], [-1.4048620459126708], [-1.2988138882311988], [-1.2017077791152193], [-1.1127688912168858], [-1.0312706985091609], [-0.9565392051243061], [-0.8879540970721507], [-0.8249478583856876], [-0.7670035961772956], [-0.7136520995286785], [-0.6644684966304356], [-0.6190687583989716], [-0.5771062134840089],
              [-0.5382681804201349], [-0.5022727811224021], [-0.4688659711138508], [-0.4378188022176291], [-0.40892492032561795], [-0.38199829236131694], [-0.356871151306608], [-0.3333921451660487], [-0.3114246742864395], [-0.29084540102664813], [-0.2715429160242747], [-0.2534165459799516], [-0.23637528880207825], [-0.22033686300468008], [-0.2052268593482769], [-0.1909779838052427], [-0.1775293819838559], [-0.16482603613939062], [-0.1528182268259207], [-0.141461052095098],
              [-0.13071399792837857], [-0.12054055430002916], [-0.11090787191443939], [-0.10178645524833921], [-0.09314988806239005], [-0.08497458803317286], [-0.07723958760150706], [-0.0699263385415113], [-0.06301853813160513], [-0.05650197515788373], [-0.050364394305557204], [-0.0445953777984172], [-0.0391862434320541], [-0.034129958415797114], [-0.029421068692721798], [-0.025055643647920903], [-0.021031236343760105], [-0.017346859638204308], [-0.014002978749707781],
               [-0.011001521030997203], [-0.008345903905917261], [-0.006041082110155127], [-0.004093615560086733], [-0.002511759356214488], [-0.0013055776105475364], [-0.0004870829722087122], [-7.040391306118486e-5]])
            # call cpeg
            # @show X[i]
            # @show ev.U
            # main_cpeg(ev, unscale_state(ev,X[i])...)
            # for i = 1:length(ev.U)
            #     ev.U[i] = SA[0.0]
            # end
            ev.U = [ev.U[i+1] for i = 1:length(ev.U)-1]
            r0, v0, σ0 = unscale_state(ev,x0s)
            althist, drhist, crhist, σhist = main_cpeg(ev, r0,v0,σ0)
            xf_dr, xf_cr = rangedistances(ev,rf,SVector{6}([r0;v0]))

            num2plot = float(length(althist))
            plot_groundtracks(drhist/1e3,crhist/1e3,althist/1e3,σhist,xf_dr/1e3,xf_cr/1e3,num2plot,"quad")
            # U[i] = ev.U[1]

            # dr = 100/ev.scale.dscale * @SVector randn(3)
            # X[i]
            # dv = @SVector randn(3)
            # @show X[i]
            # X[i+1] = rk4(ev,X[i],ev.U[1],ev.dt/ev.scale.tscale) #+ SA[dr[1];dr[2];dr[3];0;0;0;0]

        # end

        # alt, dr, cr = postprocess(ev,X,[r0;v0])

        # # althist, drhist, crhist = main_cpeg(ev,x0_s)
        # # main_cpeg(ev,x0_s)
        # main_cpeg(ev,r0,v0,σ0)
        # r02 = r0 - SA[100,600,4.2e3]
        # main_cpeg(ev,r02,v0,σ0)
        #
        # pp2 = ev.σ
        # pp = [ev.U[i][1] for i = 1:length(ev.U)]
        # #
        # mat"
        # figure
        # hold on
        # title('Control')
        # %plot(rad2deg($pp))
        # plot($pp)
        # hold off"
        #
        # mat"
        # figure
        # title('Bank angle deg')
        # hold on
        # plot(rad2deg($pp2))
        # hold off"

        return nothing
end

function plot_groundtracks(drhist,crhist,althist,σhist,xf_dr,xf_cr,num2plot,id)

    mat"
    figure
    hold on
    rgb1 = [29 38 113]/255;
    rgb2 = 1.3*[195 55 100]/255;
    drgb = rgb2-rgb1;
    cmap = [];
    for i = 1:round($num2plot)
        px = $drhist{i};
        py = $crhist{i};
        plot(px,py,'Color',rgb1 + drgb*(i)/($num2plot),'linewidth',3)
        cmap = [cmap;rgb1 + drgb*(i)/($num2plot)];
        %plot(px(1),py(1),'r.','markersize',20)
    end
    plot($xf_dr,$xf_cr,'g.','markersize',20)
    colormap(cmap);
    pp = colorbar;
    pp.Ticks = 0:(1/$num2plot):1;
    pp.Location = 'northoutside';
    pp.TickLabels = 0:round($num2plot);
    xlabel('downrange (km)')
    ylabel('crossrange (km)')
    %xlim([250,700])
    %ylim([0,20])
    hold off
    fleg = legend('figure()');
    set(fleg,'visible','off')
    %addpath('/Users/kevintracy/devel/WiggleSat/matlab2tikz-master/src')
    %matlab2tikz(strcat('cpeg_examples/bank_angle/tikz/',$id,'_crdr.tex'))
    %close all
    "

    # this one is for plotting
    mat"
    figure
    hold on
    rgb1 = [29 38 113]/255;
    rgb2 = 1.3*[195 55 100]/255;
    drgb = rgb2-rgb1;
    cmap = [];
    for i = 1:round($num2plot)
        px = $drhist{i};
        alt = $althist{i};
        plot(px,alt,'Color',rgb1 + drgb*(i)/($num2plot),'linewidth',3)
        cmap = [cmap;rgb1 + drgb*(i)/($num2plot)];
    end
    %plot($xf_dr,$xf_cr,'g.','markersize',20)
    colormap(cmap);
    pp = colorbar;
    pp.Ticks = 0:(1/$num2plot):1;
    pp.Location = 'northoutside';
    pp.TickLabels = 0:round($num2plot);
    plot([0,800],ones( 2,1)*10,'r' )
    plot($xf_dr,10,'g.','markersize',20)
    %xlim([400,700])
    %ylim([8,30])
    xlabel('downrange (km)')
    ylabel('altitude (km)')
    hold off
    %saveas(gcf,'alt.png')
    fleg = legend('figure()');
    set(fleg,'visible','off')
    %addpath('/Users/kevintracy/devel/WiggleSat/matlab2tikz-master/src')
    %matlab2tikz(strcat('cpeg_examples/bank_angle/tikz/',$id,'_altdr.tex'))
    %matlab2tikz('bankaoa_alt.tex')
    %close all
    "

    mat"
    figure
    hold on
    rgb1 = [29 38 113]/255;
    rgb2 = 1.3*[195 55 100]/255;
    drgb = rgb2-rgb1;
    cmap = [];
    for i = 1:round($num2plot)
        bank = $σhist{i};
        plot(rad2deg(bank),'Color',rgb1 + drgb*(i)/($num2plot),'linewidth',3)
        cmap = [cmap;rgb1 + drgb*(i)/($num2plot)];
    end
    colormap(cmap);
    pp = colorbar;
    pp.Ticks = 0:(1/$num2plot):1;
    pp.Location = 'northoutside';
    pp.TickLabels = 0:round($num2plot);
    %xlim([400,700])
    %ylim([8,30])
    xlabel('step')
    ylabel('bank angle (degrees)')
    hold off
    %saveas(gcf,'alt.png')
    fleg = legend('figure()');
    set(fleg,'visible','off')
    "

    return nothing
end

tt()

# # incorrect
# typeof(ev.params.aero) = AeroParameters
# typeof(ρ) = Float64
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#475#479"{Int64, CPEGWorkspace, Float64}, Float64}, Float64, 7}}
# typeof(v) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#475#479"{Int64, CPEGWorkspace, Float64}, Float64}, Float64, 7}}
#
# # correct
# typeof(ev.params.aero) = AeroParameters
# typeof(ρ) = ForwardDiff.Dual{ForwardDiff.Tag{var"#477#481"{Int64, CPEGWorkspace, Vector{SVector{7, Float64}}, Float64}, Float64}, Float64, 1}
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#477#481"{Int64, CPEGWorkspace, Vector{SVector{7, Float64}}, Float64}, Float64}, Float64, 1}}
# typeof(v) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#477#481"{Int64, CPEGWorkspace, Vector{SVector{7, Float64}}, Float64}, Float64}, Float64, 1}}

# # succeeding
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}
# typeof(h) = ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}
# typeof(ρ) = ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}
# typeof(ev.params.aero) = AeroParameters
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}
# typeof(v) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}
#
# # failing
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}
# typeof(h) = ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}
# typeof(ρ) = Float64
# typeof(ev.params.aero) = AeroParameters
# typeof(r) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}
# typeof(v) = SVector{3, ForwardDiff.Dual{ForwardDiff.Tag{var"#1006#1007"{CPEGWorkspace}, Float64}, Float64, 7}}


text:
"""
A thread member.
More details [here](https://discord.com/developers/docs/resources/channel#thread-member-object).
"""
struct ThreadMember
    id::Optional{Snowflake}
    user_id::Optional{Snowflake}
    join_timestamp::DateTime
    flags::Int
end
@boilerplate ThreadMember :constructors :docs :lower :merge :mock

ext: jl


text:
module FractalGeometry

using Meshing
using FileIO
using Roots
using MeshIO
using ForwardDiff
using GeometryBasics
import GeometryTypes
using DistMesh
using LinearAlgebra
using CUDA
#using LightGraphs
import Makie
using SparseArrays
import GeometryBasics.faces
using StaticArrays
using NearestNeighbors
using Plots
using ForwardDiff
using DifferentialEquations

include("curves.jl")
include("distancefields.jl")
include("geom_example.jl")
include("simplicialsurface.jl")

end

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


text:
module TestIteration
    using Test, DataFrames
    using DataFrames: columns

    df = DataFrame(A = Vector{Union{Int, Missing}}(1:2), B = Vector{Union{Int, Missing}}(2:3))

    @test size(eachrow(df)) == (size(df, 1),)
    @test eachrow(df)[1] == DataFrameRow(df, 1)
    @test collect(eachrow(df)) isa Vector{DataFrameRow{DataFrame}}
    @test eltype(eachrow(df)) == DataFrameRow{DataFrame}
    for row in eachrow(df)
        @test isa(row, DataFrameRow)
        @test (row[:B] - row[:A]) == 1
        # issue #683 (https://github.com/JuliaData/DataFrames.jl/pull/683)
        @test collect(pairs(row)) isa Vector{Pair{Symbol, Int}}
    end

    # TODO - clean up redundant tests after eachcol deprecation
    @test size(eachcol(df)) == (size(df, 2),)
    @test size(eachcol(df, true)) == (size(df, 2),)
    @test size(columns(df)) == (size(df, 2),)
    @test size(eachcol(df, false)) == (size(df, 2),)
    @test length(eachcol(df)) == size(df, 2)
    @test length(eachcol(df, true)) == size(df, 2)
    @test length(columns(df)) == size(df, 2)
    @test length(eachcol(df, false)) == size(df, 2)
    @test eachcol(df)[1] == (:A => df[1]) # this will be df[1] after eachcol deprecation
    @test eachcol(df, true)[1] == (:A => df[1])
    @test columns(df)[1] == df[1]
    @test eachcol(df, false)[1] == df[1]
    @test collect(eachcol(df)) isa Vector{Pair{Symbol, AbstractVector}}
    @test collect(eachcol(df, true)) isa Vector{Pair{Symbol, AbstractVector}}
    @test collect(eachcol(df)) == [:A => [1, 2], :B => [2, 3]]
    @test collect(eachcol(df, true)) == [:A => [1, 2], :B => [2, 3]]
    @test collect(columns(df)) isa Vector{AbstractVector}
    @test collect(eachcol(df, false)) isa Vector{AbstractVector}
    @test collect(columns(df)) == [[1, 2], [2, 3]]
    @test collect(eachcol(df, false)) == [[1, 2], [2, 3]]
    @test eltype(eachcol(df)) == Pair{Symbol, AbstractVector}
    @test eltype(eachcol(df, true)) == Pair{Symbol, AbstractVector}
    @test eltype(columns(df)) == AbstractVector
    @test eltype(eachcol(df, false)) == AbstractVector
    for col in eachcol(df)
        @test typeof(col) <: Pair{Symbol, <:AbstractVector}
    end
    for col in eachcol(df, true)
        @test typeof(col) <: Pair{Symbol, <:AbstractVector}
    end
    for col in columns(df)
        @test isa(col, AbstractVector)
    end
    for col in eachcol(df, false)
        @test isa(col, AbstractVector)
    end

    @test map(x -> minimum(convert(Array, x)), eachrow(df)) == [1,2]
    @test map(Vector, eachrow(df)) == [[1, 2], [2, 3]]
    @test mapcols(minimum, df) == DataFrame(A = [1], B = [2])
    @test map(minimum, eachcol(df)) == DataFrame(A = [1], B = [2]) # this is deprecated
    @test map(minimum, eachcol(df, true)) == DataFrame(A = [1], B = [2]) # this is deprecated
    @test map(minimum, columns(df)) == [1, 2]
    @test map(minimum, eachcol(df, false)) == [1, 2]
    @test eltypes(mapcols(Vector{Float64}, df)) == [Float64, Float64]
    @test eltypes(map(Vector{Float64}, eachcol(df))) == [Float64, Float64] # this is deprecated
    @test eltypes(map(Vector{Float64}, eachcol(df, true))) == [Float64, Float64] # this is deprecated
    @test eltype(map(Vector{Float64}, columns(df))) == Vector{Float64}
    @test eltype(map(Vector{Float64}, eachcol(df, false))) == Vector{Float64}

    # test mapcols corner cases
    # this behavior might change when we rework setindex! to follow standard broadcasting rules
    # notice that now mixing vectors and scalars is allowed in some cases but not in others
    # this is likely to change
    df_mapcols = DataFrame(a=1:10, b=11:20)
    @test mapcols(sum, df_mapcols) == DataFrame(a=55, b=155)
    @test mapcols(x -> x[1] == 1 ? 0 : [0], df_mapcols) == DataFrame(a=0, b=0)
    @test mapcols(x -> x[1] == 1 ? x : 0, df_mapcols) == DataFrame(a=1:10, b=0)
    @test_throws ArgumentError mapcols(x -> x[1] != 1 ? x : 0, df_mapcols)

    row = DataFrameRow(df, 1)

    row[:A] = 100
    @test df[1, :A] == 100

    row[1] = 101
    @test df[1, :A] == 101

    df = DataFrame(A = Vector{Union{Int, Missing}}(1:4), B = Union{String, Missing}["M", "F", "F", "M"])

    s1 = view(df, 1:3, :)
    s1[2,:A] = 4
    @test df[2, :A] == 4
    @test view(s1, 1:2, :) == view(df, 1:2, :)

    s2 = view(df, 1:2:3, :)
    s2[2, :B] = "M"
    @test df[3, :B] == "M"
    @test view(s2, 1:1:2, :) == view(df, [1,3], :)

    # @test_fail for x in df; end # Raises an error
end


text:
"""
    findiff(p,q,r,a,b,alpha,beta,M)
berisi program untuk mencari solusi persamaan differensial
`x''=f(t,x,x')` dengan masalah nilai batas `x(a)=alpha` dan `x(b)=beta` pada interval `[a,b]`.
Program ini secara default berisi 7 masukan, yaitu fungsi `p(t)`, `q(t)` dan `r(t)` dari persamaan
`x''=p(t)x'+q(t)x+r(t)`, titik ujung interval `[a,b]`, nilai batas `[alpha,beta]`, dan
jumlah sub-interval `M`.
# Example
```jl
julia> p(t)= 2*t/(1+t^2);
       q(t)= -2/(1+t^2);
       r(t)= 1 + 0*t;

julia> solusi = findiff(p,q,r,0,4,1.25,-0.95,40)
41×2 Array{Float64,2}:
 0.0   1.25
 0.1   1.29077
 0.2   1.31665
 0.3   1.32791
 ⋮
 3.7  -1.03501
 3.8  -1.01924
 3.9  -0.991063
 4.0  -0.95
```
return solusi masalah nilai batas `sol`.
"""
function findiff(p,q,r,a,b,alpha,beta,M)
  h = (b-a)/M;
  T = a:h:b;
  T = T[2:end-1];
  #% Bangun Matriks B
  B = -h^2*r.(T);
  B[1] = B[1] + (1+h/2*p(T[1]))*alpha;
  B[end] = B[end] + (1-h/2*p(T[end]))*beta;
  # Bangun Matriks A - Bagian Diagonal
  Ad = 2 .+h^2*q.(T);
  #% Bangun Matriks A - Bagian Bawah Diagonal
  Tbawah = T[2:end];
  Abawah = -1 .-h/2*p.(Tbawah);
  #% Bangun Matriks A - Bagian Atas Diagonal
  Tatas = T[1:end-1];
  Aatas = -1 .+h/2*p.(Tatas);
  A = Tridiagonal(Abawah,Ad,Aatas)
  #% Selesaikan AX=B
  X = A\B;
  T = [a; T; b];
  X = [alpha; X ;beta];
  solusi = [T X];
end


text:
using Visualize, GeometryTypes, ModernGL
using Visualize: GLRasterization, Uniforms, VertexCS, get_texture_atlas, Sprite2, Window
using Visualize.GLRasterization: show!, destroy!, swapbuffers!
using GPUArrays

resolution = (1024, 1024)
window = GLFWWindow(Area => resolution)
for event in Visualize.NativeWindowEvents
    add!(window, event)
end
show!(window)
window[Visualize.Open] = true
GLBackend.init(window[Window])

proj = orthographicprojection(SimpleRectangle(0, 0, resolution...), -10_000f0, 10_000f0)

proj
uniforms = Visualize.Uniforms(
    proj,
    Vec4f0(1, 0, 0, 1),
    Vec4f0(1, 0, 1, 1),
    Visualize.circle
)

N = 20
middle = Vec2f0(resolution) / 2f0
radius = min(resolution...) / 2f0
vertices = [VertexCS(
    Vec2f0((sin(2pi * (i / N)) , cos(2pi * (i / N))) .* radius) .+ middle,
    Vec4f0(1, i/N, 0, 1),
    Vec2f0(40, 40)
) for i = 1:N]
uniform_buff = UniformBuffer(uniforms);
vbo = VertexArray(vertices, face_type = Face{1, OffsetInteger{1, GLint}});

draw_particles = GLRasterizer(
    vbo, (uniform_buff,),
    Visualize.vertex_main, Visualize.fragment_main;
    geometryshader = Visualize.geometry_main
)

glClearColor(0, 0, 0, 0)
glDisable(GL_DEPTH_TEST)
GLAbstraction.enabletransparency()

@async begin
    while isopen(window)
        GLFW.PollEvents()
        glViewport(0, 0, widths(window[Area])...)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        draw_particles(vbo, (uniform_buff,))
        swapbuffers!(window)
        sleep(0.001)
        yield()
    end
    destroy!(window)
end


text:
function Psis(ϵ::Real)
    return [1-ϵ ϵ;
            ϵ 1-ϵ]
end

max_stars_repo_name: jpfairbanks/NewsGraphing

max_issues_repo_name: jpfairbanks/NewsGraphing

max_forks_repo_name: jpfairbanks/NewsGraphing


text:
module Knet
using Libdl
# using LinearAlgebra, Statistics, SpecialFunctions, Libdl

# To see debug output, start julia with `JULIA_DEBUG=Knet julia`
# To perform profiling, set ENV["KNET_TIMER"] to "true" and rebuild Knet. (moved this to gpu.jl)
# The @dbg macro below evaluates `ex` only when debugging. The @debug macro prints stuff as documented in Julia.
macro dbg(ex); :(if Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug,:none,Knet)!==nothing; $(esc(ex)); end); end

const libknet8 = Libdl.find_library(["libknet8"], [joinpath(dirname(@__DIR__),"deps")])

using  AutoGrad: @diff, Param, params, grad, gradloss, value, cat1d, @primitive, @zerograd, @primitive1, @zerograd1, forw, back, Value, AutoGrad
export AutoGrad, @diff, Param, params, grad, gradloss, value, cat1d #@primitive, @zerograd, @primitive1, @zerograd1, forw, back, Value, getval

include("gpu.jl");              export gpu
include("uva.jl")
include("kptr.jl");             export knetgc # KnetPtr
include("karray.jl");           export KnetArray
include("gcnode.jl");
include("ops.jl");
include("unary.jl");            export relu, sigm, invx, elu, selu
include("broadcast.jl");        # elementwise broadcasting operations
include("reduction.jl");        # sum, max, mean, etc.
include("linalg.jl");           export mat # matmul, axpy!, transpose, (i)permutedims
include("conv.jl");             export conv4, pool, deconv4, unpool
include("batchnorm.jl");        export batchnorm, bnmoments, bnparams
include("rnn.jl");              export rnnforw, rnninit, rnnparam, rnnparams, RNN # TODO: deprecate old interface
include("data.jl");             export Data, minibatch
include("model.jl");		export param, param0, train!, Train
include("loss.jl");             export logp, logsoftmax, logsumexp, softmax, nll, logistic, bce, accuracy, zeroone # TODO: PR
include("dropout.jl");          export dropout
include("update.jl"); 		export SGD, Sgd, Momentum, Nesterov, Adam, Adagrad, Adadelta, Rmsprop, update!, optimizers
include("distributions.jl"); 	export gaussian, xavier, bilinear
include("random.jl");           export setseed  # TODO: deprecate setseed
include("hyperopt.jl");         export hyperband, goldensection
include("serialize.jl");        export gpucopy,cpucopy
include("jld.jl");              # load, save, @load, @save; not exported use with Knet. prefix.

"""
    Knet.dir(path...)

Construct a path relative to Knet root.

# Example
```julia
julia> Knet.dir("examples","mnist.jl")
"/home/dyuret/.julia/v0.5/Knet/examples/mnist.jl"
```
"""
dir(path...) = joinpath(dirname(@__DIR__),path...)

# See if we have a gpu at initialization:
function __init__()
    try
        r = gpu(true)
        # info(r >= 0 ? "Knet using GPU $r" : "No GPU found, Knet using the CPU")
    catch e
        gpu(false)
        # warn("Knet using the CPU: $e")
    end
end

# @use X,Y,Z calls using on packages installing them if necessary. (WIP)
# 1. still need "using Knet"
# 2. Pkg.insalled gives false for stdlib packages.
# macro use(ps)
#     if isa(ps, Symbol); ps = Expr(:tuple,ps); end
#     a = map(ps.args) do p
#         s=string(p)
#         esc(:(haskey(Pkg.installed(),$s)||Pkg.add($s); using $p))
#     end
#     Expr(:block,:(using Pkg),a...)
# end
# export @use

end # module


text:
include("command_utils.jl")

function handle_update(token::String, message::Dict)::Nothing
    if haskey(message, "message")
        message = message["message"]
        message_text = message["text"]
        sender_id = message["chat"]["id"]

        command_name, command_params = parse_command(message_text)

        if sender_id != MASTER_ID
            send_message(token, sender_id, STRANGER_MESSAGE)
            return
        elseif command_name == ""
            return
        elseif command_name == "start"
            send_message(token, MASTER_ID, START_MESSAGE)
            return
        elseif command_name in COMMANDS
            send_message(token, MASTER_ID, "Got you, master!")
            message = execute_command(command_name, command_params)
        else
            message = "I do not know what to do, master. Your command \"$(command_name[1:end])\" confuses me."
        end

        send_message(token, MASTER_ID, message)
    end
end

max_stars_repo_name: iwsylit/spotify-utils-bot

max_issues_repo_name: iwsylit/spotify-utils-bot

max_forks_repo_name: iwsylit/spotify-utils-bot


text:
module DemoPackageMCF

print_greeting() = printstyled("Hello, world!"; color=:magenta)

end


text:
export rangedomain, Forward, Reverse

_NI(x) = error("not implemented: $x")

rangedomain(op::Operator) = domain(rangespace(op))

containstransfer(M::Operator) = ApproxFun.iswrapper(M::Operator) && containstransfer(M.op)
for optype in (:(ApproxFun.InterlaceOperator),:(ApproxFun.PlusOperator),:(ApproxFun.TimesOperator))
  eval(:(containstransfer(M::$optype) = any(containstransfer,M.ops)))
end

# Overloads

@compat const PureLaurent{D<:Domain} = Hardy{false,D}

@compat const GeneralInterval = Union{Segment,Interval,PeriodicSegment}

# rem()

# InterpolationNode

@compat struct InterpolationNode{S<:Space}<:Number
  sp::S
  k::Int
  n::Int
  function InterpolationNode{S}(sp::S,k::Int,n::Int) where S<:Space
    @assert k ≥ 1
    @assert k ≤ n
    new{typeof(sp)}(sp,k,n)
  end
end
InterpolationNode(sp::Space,k::Int,n::Int) = InterpolationNode{typeof(sp)}(sp,k,n)

Base.convert(::Type{T},p::InterpolationNode) where T<:Number = convert(T,points(p.sp,p.n)[p.k])
Base.convert(::Type{InterpolationNode{S}},p::InterpolationNode{S}) where S<:Space = p
Base.promote_rule(::Type{T},::Type{InterpolationNode{S}}) where {T<:Number,S<:Space} = promote_type(T,eltype(S))
Base.show(p::InterpolationNode) = show(convert(Nu))
ApproxFun.space(p::InterpolationNode) = p.sp
ApproxFun.domain(p::InterpolationNode) = domain(p.sp)
# eltype(p::InterpolationNode) = eltype(p.sp)
ApproxFun.prectype(p::InterpolationNode) = prectype(p.sp)

# Base.zero(p::InterpolationNode) = zero(eltype(p))

# Circle domains

interval_mod(x,a,b) = (b-a)*mod((x-a)/(b-a),1)+a
domain_mod(x,p::PeriodicSegment) = interval_mod(x,leftendpoint(p),rightendpoint(p))

# Interval domains

function approx_in(x,d::GeneralInterval; atol=50eps(arclength(d)), kwargs...)
    x <= rightendpoint(d)+atol && x >= leftendpoint(d) - atol
end
approx_in(x,d) = approx_in(convert(typeof(d),x),d)

function approx_issubset(a,b;kwargs...)
    acupb = a∪b
    atol_isapprox(acupb,b;kwargs...)
end

atol_isapprox(a,b; atol=50eps(max(arclength(a),arclength(b))), kwargs...) = isapprox(a,b;atol=atol,kwargs...)
coveringinterval(dsm::AbstractArray{T}) where {T<:AbstractInterval} = Interval(minimum(leftendpoint(convert(Domain,d)) for d in dsm),maximum(rightendpoint(convert(Domain,d)) for d in dsm))
coveringinterval(ds::AbstractArray) = coveringinterval([convert(Domain,d) for d in ds])
mapinterval(f,d::T) where T<:GeneralInterval = T(extrema((f(leftendpoint(d)),f(rightendpoint(d))))...)
mapinterval(f,d) = mapinterval(f,convert(Domain,d))

# Newton's method

domain_newton(f,df,y::U,D::Domain,x::U=convert(U,ApproxFun.checkpoints(D)[1]),tol=400eps(eltype(U))) where U = basic_newton(f,df,y,x,tol)
domain_guess(x,dom::Domain,ran::Domain) = rand(ran)

domain_guess(x::SVector{N},dom::ProductDomain,ran::ProductDomain) where N = SVector{N}([interval_guess(x[i],dom.domains[i],ran.domains[i]) for i =1:N])

function basic_newton(f,df,y::U,x::U,tol=40eps(U)) where U
  z = x
  rem = f(z)-y
  for i = 1:2log2(-200log(eps(eltype(y))))
    #    while abs(rem) > 30eps(maximum(∂(D)))
    z -= df(z) \ rem
    rem = f(z)-y
    norm(rem) < tol && break
  end

  norm(rem) > tol && error("Newton: failure to converge: y = $y, x estimate = $z, rem = $rem")
  z
end

# function monotonic_newton(f,df,y::U,xl::U,xh::U,tol=40eps(U)) where U ## TODO: REALLY WANT UPPER/LOWER GUESSES OR STH
#   σ = sign(df(x))
#   det(σ) == 0 && error("f has critical point at $x")
#   z = xl
#   σzl = σ*xl; σzu = σ*xh; # actually find x,u
#   rem = f(z)-y
#   σz = σ*z
#   rem < 0 ? (σzu = min(σzu,σz)) : (σzl = max(σzl,σz))
#   println((z,y,rem,σz,σzl,σzu))
#   for i = 1:200#log2(-200log(eps(eltype(y))))
#     #    while abs(rem) > 30eps(maximum(∂(D)))
#     z -= df(z) \ rem
#     rem = f(z)-y
#     norm(rem) < tol && break
#
#     println((z,y,rem,σz,σzl,σzu))
#
#     σz = σ*z
#     rem < 0 ? (σzu = min(σzu,σz)) : (σzl = max(σzl,σz))
#     σzu < σzl && (σzu += 4tol, σzl -= 4tol)
#     (σz > σzu || σz < σzl) && (z = σ*(σzl + (σzu-σzl)*rand()))
#     println((z,y,rem,σz,σzl,σzu))
#   end
#
#   norm(rem) > tol && error("Newton: failure to converge")
#   z
# end

function interval_newton(f,df,y::U,da::T,db::T,x::U=(da+db)/2,tol=40eps(max(abs(da),abs(db)))) where {T,U<:Real}
  rem = f(x)-y
  for i = 1:2log2(-200log(eps(typeof(y))))
    #    while abs(rem) > 30eps(maximum(∂(D)))
    x -= rem / df(x) + eps(x)*(mod(i,11)-5)/4
    x .> db && (x = db)
    x .< da && (x = da)
    rem = f(x)-y
    #  println((x,f(x)-y))
    abs(rem) < tol && break
  end

  abs(rem) > tol && error("Newton: failure to converge: y = $y, x estimate = $x, rem = $rem")
  x
end

function interval_newton(f,df,y::U,da::T,db::T,x::U=(da+db)/2,tol=40eps(max(abs(da),abs(db)))) where {T,U<:Complex}
  rem = f(x)-y
  for i = 1:2log2(-200log(eps(typeof(abs(y)))))
    #    while abs(rem) > 30eps(maximum(∂(D)))
    x -= rem / df(x)
    rem = f(x)-y
    abs(rem) < tol && break
  end

  abs(rem) > tol && error("Newton: failure to converge: y = $y, x estimate = $x, rem = $rem")
  x
end
function domain_newton(f,df,y::U,D::GeneralInterval,x::U=(leftendpoint(D)+rightendpoint(D))/2,
            tol=40eps(max(leftendpoint(D),rightendpoint(D)))) where U<:Union{Real,Complex}
    # x = convert(U,x);
    interval_newton(f,df,y,leftendpoint(D),rightendpoint(D),x,tol)
end

interval_guess(y::Number,dom::Domain,ran::Domain) =
    (leftendpoint(dom)*rightendpoint(ran)-leftendpoint(ran)*rightendpoint(dom)+y*(rightendpoint(dom)-leftendpoint(dom)))/(rightendpoint(ran)-leftendpoint(ran))
domain_guess(y::Number,dom::GeneralInterval,ran::GeneralInterval) =
    interval_guess(y,dom,ran)

function disc_newton(f,df,y::U,rad::T,x::U=y,tol=40eps(rad)) where {T, U}
#   x = convert(typeof(y),rad*rand(typeof(real(y))))
  rem = f(x) - y
  for i = 1:2log2(-200log(eps(typeof(real(y)))))
    x -= rem / df(x)
    rem = f(x) - y
    abs(rem) < tol && break
    abs(x) > rad && (x /= abs(x))
  end
  abs(rem) > tol && error("Newton: failure to converge: y = $y, x estimate = $x, rem = $rem")
  x
end

# ForwardDiff
# function forwarddiff(f)
#   function dfdx(x)
#     ForwardDiff.derivative(f,x)::typeof(x)
#   end
#   dfdx
# end

@compat struct FunctionDerivative{ff<:Function}
  f::ff
end
(fd::FunctionDerivative)(x) = oftype(x,dualpart(fd.f(Dual(x,one(x)))))

@compat struct BasisFun{S<:Space,II<:Integer}
  s::S
  k::II
end
(sk::BasisFun)(x) = getbasisfun(x,sk)
getbasisfun(x,sk::BasisFun) = Fun(sk.s,[zeros(prectype(sk.s),sk.k-1);one(prectype(sk.s))])(x)
getbasisfun_int(x,sk::BasisFun) = cumsum(Fun(sk.s,[zeros(prectype(sk.s),sk.k-1);one(prectype(sk.s))]))(x)

@compat struct BasisFunInt{S<:Space,II<:Integer}
  sk::BasisFun{S,II}
end
cumsum(sk::BasisFun) = BasisFunInt(sk)
(ski::BasisFunInt)(x) = getbasisfun_int(ski.sk,x)

#specialised getbasisfuns

fourierCSk(x,d,k::Integer) = rem(k,2) == 1 ? cos(fld(k,2)*tocanonical(d,x)) : sin(fld(k,2)*tocanonical(d,x))
function fourierCSk_int(x,d,k::Integer)
  k == 1 && return x
  (rem(k,2) == 1 ? -sin(fld(k,2)*tocanonical(d,x)) : cos(fld(k,2)*tocanonical(d,x)))/
    fld(k,2)/tocanonicalD(d,x)
end
getbasisfun(x,sk::BasisFun{Fourier{DD},K}) where {DD, K<:Integer} = fourierCSk(x,domain(sk.s),sk.k)
getbasisfun_int(x,sk::BasisFun{Fourier{DD},K}) where {DD, K<:Integer} = fourierCSk_int(x,domain(sk.s),sk.k)

chebyTk(x,d,k::Integer) = cos((k-1)*acos(tocanonical(d,x))) #roundoff error grows linearly(??) with k may not be bad wrt x too
function chebyTk_int(x,d,k::Integer)
#   k == 1 && return tocanonical(d,x)/tocanonicalD(d,x)
  k == 2 && return tocanonical(d,x)^2/2tocanonicalD(d,x)
  ((k-1)*chebyTk(x,d,k+1) - k*tocanonical(d,x)*chebyTk(x,d,k))/((k-1)^2-1)/tocanonicalD(d,x)
end
getbasisfun(x,sk::BasisFun{F,K}) where {F<:Chebyshev, K<:Integer} = chebyTk(x,domain(sk.s),sk.k)
getbasisfun_int(x,sk::BasisFun{F,K}) where {F<:Chebyshev,K<:Integer} = chebyTk_int(x,domain(sk.s),sk.k)

function getbasisfun(x,sk::BasisFun{F,K}) where {F<:TensorSpace,K<:Integer}
  ks = ApproxFun.tensorizer(sk.s)[sk.k]
  prod(getbasisfun(x[i],BasisFun(sk.s.spaces[i],ks[i])) for i = eachindex(sk.s.spaces))
end
# no getbasisfun_int as you don't have antiderivatives

# # Directions
@compat abstract type Direction; end
@compat struct ForwardDirection <: Direction; end
@compat struct ReverseDirection <: Direction; end
const Forward = ForwardDirection()
const Reverse = ReverseDirection()


text:
using HDF5
using H5WriteBuffer
using Base.Test

mktemp() do p,f
    f2 = h5open(p, "w")
    wb = FileBackedBuffer(f2, "A", Float64, (3,4), 5)
    wb2 = FileBackedBuffer(f2, "B", Float64, (2,), 5)
    @test wb.bsize == 0
    @test size(wb.A) == (3,4,0)
    push!(wb, ones(3,4))
    @test wb.bsize == 1
    @test size(wb.A) == (3,4,0)
    for i=1:3 push!(wb, ones(3,4)) end
    @test wb.bsize == 4
    @test size(wb.A) == (3,4,0)
    push!(wb, ones(3,4))
    @test wb.bsize == 0
    @test size(wb.A) == (3,4,5)

    # Test interleaved writes to different dsets in the same file.
    for i=1:10 push!(wb2, 2*ones(2)) end
    @test d_open(f2, "B")[:,:] == 2*ones(2,10)

    for i=1:5 push!(wb, ones(3,4)) end
    @test wb.bsize == 0
    @test size(wb.A) == (3,4,10)

    @test d_open(f2, "A")[:,:,:] == ones(3,4,10)
end

# Test constructor taking filename.
mktempdir() do d
    cd(d) do
        fname = "somefile.h5"

        # Test new file.
        wb = FileBackedBuffer(fname, "A", Float64, (3,4), 5)
        for i=1:10 push!(wb, ones(3,4)) end
        h5open(fname) do f
            @test d_open(f, "A")[:,:,:] == ones(3,4,10)
        end
        close(wb)

        # Test existing file.
        wb = FileBackedBuffer(fname, "A", Float64, (3,4), 5)
        for i=1:5 push!(wb, ones(3,4)) end
        h5open(fname) do f
            @test d_open(f, "A")[:,:,:] == ones(3,4,15)
        end
        close(wb)

        # Test flush.
        wb = FileBackedBuffer(fname, "A", Float64, (3,4), 5)
        for i=1:4 push!(wb, ones(3,4)) end
        h5open(fname) do f
            @test d_open(f, "A")[:,:,:] == ones(3,4,15)
        end
        flush(wb)
        h5open(fname) do f
            @test d_open(f, "A")[:,:,:] == ones(3,4,19)
        end

        # Test dset function, returning a handle to the underlying data.
        @test read(dset(wb)) == ones(3,4,19)
        close(wb)

        # Test existing empty file.
        rm(fname)
        touch(fname)
        wb = FileBackedBuffer(fname, "A", Float64, (3,4), 5)
        for i=1:5 push!(wb, ones(3,4)) end
        h5open(fname) do f
            @test d_open(f, "A")[:,:,:] == ones(3,4,5)
        end
        close(wb)
    end
end


text:
"""
    c2(array, phase[; len,][directions,] periodic = false)

Calculate `C₂` (cluster) correlation function for one-, two- or
three-dimensional multiphase system.

`C₂(x)` equals to probability that corner elements of a line segment
with the length `x` cut from the array belong to the same
cluster of the specific phase. This implementation calculates C2 for
all `x`es in the range from `1` to `len` which defaults to half of the
minimal dimension of the array.

# Examples
```jldoctest
julia> c2([1,1,1,0,1,1], 1; len = 6)[:x]
6-element Array{Float64,1}:
 0.8333333333333334
 0.6
 0.25
 0.0
 0.0
 0.0
```

For a list of possible dimensions, see also: [`direction1Dp`](@ref),
[`direction2Dp`](@ref), [`direction3Dp`](@ref).
"""
function c2(array      :: AbstractArray,
            phase;
            directions :: Vector{Symbol} = array |> default_directions,
            len        :: Integer = (array |> size |> minimum) ÷ 2,
            periodic   :: Bool = false)
    field = map(x -> x == phase, array)
    s2(label_components(field, periodic ? Torus() : Plane()),
       InseparableIndicator((x, y) -> x == y != 0);
       len        = len,
       directions = directions,
       periodic   = periodic)
end


text:
struct JSONLD_Web <: JSONLD end

function mainpage_url(repo::JSONLD_Web, dataname)
    page=getpage(dataname)
    pattern = sel"script[type=\"application/ld+json\"]"
    jsonld_blocks = eachmatch(pattern, page.root)
    if length(jsonld_blocks)==0
        throw(GeneratorError(repo, "No JSON-LD Linked Data Found"))
    end
    @assert length(jsonld_blocks)==1
    script_block = text_only(first(jsonld_blocks))
    json = JSON.parse(script_block)
    return json, dataname
end


text:
# ---
# title: 420. Strong Password Checker
# id: problem420
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/strong-password-checker/description/>
# hidden: true
# ---
#
# A password is considered strong if the below conditions are all met:
#
#   * It has at least `6` characters and at most `20` characters.
#   * It contains at least **one lowercase** letter, at least **one uppercase** letter, and at least **one digit**.
#   * It does not contain three repeating characters in a row (i.e., `"...aaa..."` is weak, but `"...aa...a..."` is strong, assuming other conditions are met).
#
# Given a string `password`, return _the minimum number of steps required to
# make`password` strong. if `password` is already strong, return `0`._
#
# In one step, you can:
#
#   * Insert one character to `password`,
#   * Delete one character from `password`, or
#   * Replace one character of `password` with another character.
#
#
#
# **Example 1:**
#
#
#
#     Input: password = "a"
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: password = "aA1"
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: password = "1337C0d3"
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= password.length <= 50`
#   * `password` consists of letters, digits, dot `'.'` or exclamation mark `'!'`.
#
#
## @lc code=start
using LeetCode

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


text:
"""
Adds reactions to the model `m`
"""
function add_reactions(
    m::CoreModel,
    s::V1,
    b::V2,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat;
    check_consistency = false,
) where {V1<:VecType,V2<:VecType}
    return add_reactions(
        m,
        sparse(reshape(s, (length(s), 1))),
        sparse(b),
        sparse([c]),
        sparse([xl]),
        sparse([xu]),
        check_consistency = check_consistency,
    )
end

function add_reactions(
    m::CoreModel,
    s::V1,
    b::V2,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat,
    rxn::String,
    mets::K;
    check_consistency = false,
) where {V1<:VecType,V2<:VecType,K<:StringVecType}
    return add_reactions(
        m,
        sparse(reshape(s, (length(s), 1))),
        sparse(b),
        sparse([c]),
        sparse([xl]),
        sparse([xu]),
        [rxn],
        mets,
        check_consistency = check_consistency,
    )
end

function add_reactions(
    m::CoreModel,
    Sp::M,
    b::V,
    c::V,
    xl::V,
    xu::V;
    check_consistency = false,
) where {M<:MatType,V<:VecType}
    rxns = ["r$x" for x = length(m.rxns)+1:length(m.rxns)+length(xu)]
    mets = ["m$x" for x = length(m.mets)+1:length(m.mets)+size(Sp)[1]]
    return add_reactions(
        m,
        Sp,
        b,
        c,
        xl,
        xu,
        rxns,
        mets,
        check_consistency = check_consistency,
    )
end

function add_reactions(m1::CoreModel, m2::CoreModel; check_consistency = false)
    return add_reactions(
        m1,
        m2.S,
        m2.b,
        m2.c,
        m2.xl,
        m2.xu,
        m2.rxns,
        m2.mets,
        check_consistency = check_consistency,
    )
end

function add_reactions(
    m::CoreModel,
    Sp::M,
    b::V,
    c::V,
    xl::V,
    xu::V,
    rxns::K,
    mets::K;
    check_consistency = false,
) where {M<:MatType,V<:VecType,K<:StringVecType}

    Sp = sparse(Sp)
    b = sparse(b)
    c = sparse(c)
    xl = sparse(xl)
    xu = sparse(xu)

    all([length(b), length(mets)] .== size(Sp, 1)) ||
        throw(DimensionMismatch("inconsistent number of metabolites"))
    all(length.([c, xl, xu, rxns]) .== size(Sp, 2)) ||
        throw(DimensionMismatch("inconsistent number of reactions"))

    new_reactions = findall(Bool[!(rxn in m.rxns) for rxn in rxns])
    new_metabolites = findall(Bool[!(met in m.mets) for met in mets])

    if check_consistency
        (newReactions1, newMetabolites1) = verify_consistency(
            m,
            Sp,
            b,
            c,
            xl,
            xu,
            rxns,
            mets,
            new_reactions,
            new_metabolites,
        )
    end

    new_mets = vcat(m.mets, mets[new_metabolites])

    zero_block = spzeros(length(new_metabolites), n_reactions(m))
    ext_s = vcat(sparse(m.S), zero_block)

    mapping = [findfirst(isequal(met), new_mets) for met in mets]
    (I, J, elements) = findnz(sparse(Sp[:, new_reactions]))
    ext_sp = spzeros(length(new_mets), length(new_reactions))
    for (k, i) in enumerate(I)
        new_i = mapping[i]
        ext_sp[new_i, J[k]] = elements[k]
    end

    new_s = hcat(ext_s, ext_sp)
    newb = vcat(m.b, b[new_metabolites])
    #new_c = hcat(m.C, spzeros(size(m.C, 1), length(new_reactions)))
    newc = vcat(m.c, c[new_reactions])
    newxl = vcat(m.xl, xl[new_reactions])
    newxu = vcat(m.xu, xu[new_reactions])
    new_rxns = vcat(m.rxns, rxns[new_reactions])
    #new_lp = CoreModel(new_s, newb, new_c, m.cl, m.cu, newc, newxl, newxu, new_rxns, new_mets)
    new_lp = CoreModel(new_s, newb, newc, newxl, newxu, new_rxns, new_mets)

    if check_consistency
        return (new_lp, new_reactions, new_metabolites)
    else
        return new_lp
    end
end

"""
Verifies the consistency of a given model
"""
function verify_consistency(
    m::CoreModel,
    Sp::M,
    b::V,
    c::V,
    xl::V,
    xu::V,
    names::K,
    mets::K,
    new_reactions,
    new_metabolites,
) where {M<:MatType,V<:VecType,K<:StringVecType}

    if !isempty(new_reactions)
        statuses = Vector{ReactionStatus}(undef, length(names))
        for (i, name) in enumerate(names)
            rxn_index = findfirst(isequal(name), m.rxns)
            reaction = Sp[:, i]
            stoich_index = findfirst(Bool[reaction == m.S[:, j] for j = 1:size(m.S, 2)])
            if isnothing(rxn_index) & isnothing(stoich_index)
                statuses[i] = ReactionStatus(false, 0, "new")
            end

            if !isnothing(rxn_index) & isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same name")
            end

            if isnothing(rxn_index) & !isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same stoichiometry")
            end

            if !isnothing(rxn_index) & !isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same name, same stoichiometry")
            end
        end
    end

    return (new_reactions, new_metabolites)
end

"""
Removes a set of reactions from a CoreModel.
Also removes the metabolites not involved in any reaction.
"""
function remove_reactions(m::CoreModel, rxns::Vector{Int})
    rxns_to_keep = filter(e -> e ∉ rxns, 1:n_reactions(m))
    temp_s = m.S[:, rxns_to_keep]

    (mets_to_keep, J, val) = findnz(temp_s)
    sort!(mets_to_keep)
    unique!(mets_to_keep)
    new_s = m.S[mets_to_keep, rxns_to_keep]
    newb = m.b[mets_to_keep]
    newc = m.c[rxns_to_keep]
    newxl = m.xl[rxns_to_keep]
    newxu = m.xu[rxns_to_keep]
    new_rxns = m.rxns[rxns_to_keep]
    new_mets = m.mets[mets_to_keep]
    new_model = CoreModel(new_s, newb, newc, newxl, newxu, new_rxns, new_mets)
    return new_model
end

function remove_reactions(m::CoreModel, rxn::Integer)
    return remove_reactions(m, [rxn])
end

function remove_reactions(m::CoreModel, rxn::String)
    return remove_reactions(m, [rxn])
end

function remove_reactions(m::CoreModel, rxns::Vector{String})
    rxn_indices = [findfirst(isequal(name), m.rxns) for name in intersect(rxns, m.rxns)]
    if isempty(rxn_indices)
        return m
    else
        return remove_reactions(m, rxn_indices)
    end
end

"""
Returns indices of exchange reactions.
Exchange reactions are identified based on most commonly used prefixes.
"""
function find_exchange_reactions(
    model::CoreModel;
    exclude_biomass = false,
    biomass_str::String = "biomass",
    exc_prefs = ["EX_"; "Exch_"; "Ex_"],
)
    is_exc = falses(n_reactions(model))
    for pref in exc_prefs
        is_exc = is_exc .| startswith.(model.rxns, pref)
    end
    exc_inds = findall(is_exc)
    if exclude_biomass
        biom_inds = findall(x -> occursin(biomass_str, x), model.rxns)
        exc_inds = setdiff(exc_inds, biom_inds)
    end
    return exc_inds
end

"""
Returns indices of exchanged metabolites, ie, the outermost metabolites in the network
In practice returns the metabolites consumed by the reactions given by `find_exchange_reactions`
and if called with the same arguments, the two outputs correspond.
"""
function find_exchange_metabolites(
    model::CoreModel;
    exclude_biomass = false,
    biomass_str::String = "biomass",
    exc_prefs = ["EX_"; "Exch_"; "Ex_"],
)
    exc_rxn_inds = find_exchange_reactions(
        model,
        exclude_biomass = exclude_biomass,
        biomass_str = biomass_str,
        exc_prefs = exc_prefs,
    )
    exc_met_inds = [findfirst(x -> x == -1, model.S[:, j]) for j in exc_rxn_inds]
    return exc_met_inds
end

"""
Change the lower and/or upper bounds ('xl' and 'xu') for given reactions
"""
function change_bounds!(
    model::CoreModel,
    rxns::Vector{Int};
    xl::V = Float64[],
    xu::V = Float64[],
) where {V<:VecType}
    found = [index ∈ 1:n_reactions(model) for index in rxns]
    length(rxns[found]) == length(unique(rxns[found])) ||
        error("`rxns` appears to contain duplicates")

    if !isempty(xl)
        length(rxns) == length(xl) ||
            throw(DimensionMismatch("`rxns` size doesn't match with `xl`"))
        model.xl[rxns[found]] = xl[found]
    end
    if !isempty(xu)
        length(rxns) == length(xu) ||
            throw(DimensionMismatch("`rxns` size doesn't match with `xu`"))
        model.xu[rxns[found]] = xu[found]
    end
end

function change_bounds!(
    model::CoreModel,
    rxns::Vector{String};
    xl::V = Float64[],
    xu::V = Float64[],
) where {V<:VecType}
    found = [name ∈ model.rxns for name in rxns]
    rxn_indices = zeros(Int, length(rxns))
    rxn_indices[found] = [findfirst(isequal(name), model.rxns) for name in rxns[found]]
    change_bounds!(model, rxn_indices, xl = xl, xu = xu)
end


text:
using JuMP
import Gurobi

M = 6
N = 3

# Build model indication, and set the solver.
m = Model(solver = Gurobi.GurobiSolver(Method=2))

# Set Variables
@variable(m, x[1:N] >= 0)
@variable(m, y[1:M] >= 0)

# Set Parameters
Xind = [1 2 3]
yind = [2 3 4 5 6 7]

# define objective
@objective(m, Min,(Xind * x + yind * y)[1])

Xind1 = [2 -1 5]
Yind1 = [9 2 -3 1 -7 12]

@constraint(m, (Xind1 * x + Yind1 * y)[1] <= 6)

Yind2 = [9 8 7 9 7 6]

Xind2 = [3 2 1; 2 5 7; 9 -1 4]
Zind2 = [4 3 2]'

@constraint(m, Xind2 * x .>= Zind2)

status = solve(m)

print(getvalue(x))
print(getvalue(y))


text:
####
#### internal methods, not part of the API
####

####
#### scaled problem
####

"""
Given `Ê`, solves
```math
Ê = \\log(∑ᵢ \\exp(zᵢ + ϵᵢ Ĉ))
```
for `Ĉ`. It is assumed that called does all the sanity checks for parameters.

# Important

The original problem needs to be mapped into this one with

| original problem    | scaled problem |
|---------------------|----------------|
| ``(1 - σ) Ê``       | ``Ê``          |
| ``(1 - σ) Ĉ``       | ``Ĉ``          |
| ``Ω̂ᵢ + (1 - σ) p̂ᵢ`` | ``zᵢ``         |
| ``ϵᵢ``              | ``ϵᵢ`` (same)  |
"""
struct ScaledProblem{N,Tz<:Real,Tϵ<:Real}
    zs::SVector{N,Tz}
    ϵs::SVector{N,Tϵ}
end

Broadcast.broadcastable(problem::ScaledProblem) = Ref(problem)

# workaround for https://github.com/JuliaStats/LogExpFunctions.jl/issues/33
function _logsumexp(x::SVector)
    A = maximum(x)
    log(sum(exp.(x .- A))) + A
end

"""
$(SIGNATURES)

Solve the scaled `problem` using Newton's method. `atol` is the absolute tolerance (should
be positive) between the left and right hand sides, and `M` is the maximum number of
iterations.
"""
function calculate_Ĉ_newton(problem::ScaledProblem, Ê::Real, atol, M = 100)
    @unpack zs, ϵs = problem
    Ĉ = (Ê - _logsumexp(zs)) / mean(ϵs)
    for _ in 1:M
        x = zs .+ ϵs .* Ĉ
        # calculate logsumexp and “weighted softmax” for f′ at the same time
        # it works because all inputs are finite
        A = maximum(x)
        X = exp.(x .- A)        # numerical stability correction
        f = log(sum(X)) + A - Ê
        abs(f) ≤ atol && return Ĉ
        f′ = dot(X, ϵs) / sum(X)
        Δ = f / f′
        Ĉ′ = Ĉ - Δ
        Δ ≠ 0 && Ĉ′ == Ĉ && return Ĉ # no change because abs(Δ) < eps(Ĉ)
        Ĉ = Ĉ′
    end
    error("no solution after $M iterations")
end

###
### consumption aggregator
###

"""
$(SIGNATURES)

Check parameters, throw an error message for invalid ones.
"""
function check_σ_ϵs(σ, ϵs)
    @argcheck σ > 0 DomainError
    @argcheck σ ≠ 1 DomainError
    @argcheck all(ϵs .> 0) DomainError
    nothing
end

"""
Default tolerance for the residual in `Ê = ...` for calculating `Ĉ`.
"""
const DEFAULT_TOL = 1e-10

# `T` is for a separate, conditional dispatch path to catch `ForwardDiff.Dual`.
function calculate_Ĉ(::Type{T}, Ê, σ, Ω̂s, ϵs, p̂s;
                     tol = DEFAULT_TOL,
                     skipcheck::Bool = false) where {T <: Real}
    skipcheck || check_σ_ϵs(σ, ϵs)
    scaled_problem = ScaledProblem(Ω̂s .+ (1 - σ) .* p̂s, ϵs)
    atol = (1 + abs(Ê)) * abs(tol)
    calculate_Ĉ_newton(scaled_problem, (1 - σ) * Ê, atol) / (1 - σ)
end

"""
$(SIGNATURES)

Calculate `Ĉ` from parameters.

# Arguments

- `skipcheck`: skip sanity checks

- `tol`: (relative) tolerance, translated to absolute tolerance `(1 + abs(Ê)) * abs(tol)`.
"""
function calculate_Ĉ(Ê::TÊ, σ::Tσ, Ω̂s::AbstractVector{TΩ̂}, ϵs::AbstractVector{Tϵ},
                     p̂s::AbstractVector{Tp̂};
                     tol = DEFAULT_TOL,
                     skipcheck::Bool = false) where {Tσ,TΩ̂,Tϵ,TÊ,Tp̂}
    T = promote_type(Tσ,TΩ̂,Tϵ,TÊ,Tp̂)
    calculate_Ĉ(T, Ê, σ, Ω̂s, ϵs, p̂s; tol, skipcheck)
end

###
### partial derivatives
###
### Coded here internally so that they can be reused for various AD frameworks.

"""
$(SIGNATURES)

Calculate and return

1. `Zs`, where each element is ``Zᵢ ∝ exp(zᵢ)``, scaled to protect from under- and overflow,
where ``zᵢ = (Ĉ ϵᵢ + pᵢ) (1 - σ) + Ωᵢ``,

2. `sum(Zs .* ϵs)`.

These are useful for partial derivatives.
"""
function calculate_Zs_∑Zϵ(Ĉ, Ê, σ, Ω̂s, ϵs, p̂s)
    zs = (Ĉ .* ϵs .+ p̂s) .* (1 - σ) .+ Ω̂s
    Zs = exp.(zs .- maximum(zs))
    Zs, dot(Zs, ϵs)
end

"""
$(SIGNATURES)

Calculate `∂Ĉ/∂Ê`.
"""
@inline calculate_∂Ê(Zs, ∑Zϵ) = sum(Zs) / ∑Zϵ

"""
$(SIGNATURES)

Calculate `∂Ĉ/∂p̂ᵢ` for each price in `p̂s`.

Also note that the result can be reused, as ``∂Ĉ/∂ϵᵢ=Ĉ⋅∂Ĉ/∂p̂ᵢ`` and
``∂Ĉ/∂Ω̂ᵢ=∂p̂ᵢ/(1-σ)``.
"""
@inline calculate_∂p̂s(Zs, ∑Zϵ) = Zs ./ -∑Zϵ

"""
$(SIGNATURES)

Calculate `∂Ĉ/∂σ`.
"""
@inline function calculate_∂σ(Zs, ∑Zϵ, Ĉ, Ê, σ, ϵs, p̂s)
    (Ĉ + sum(Zs .* (p̂s .- Ê)) / ∑Zϵ) / (1 - σ)
end


text:
############################################################################################
"""
$(SIGNATURES)
Compute ESS as a function of proposed temperature `λ` against current temperature `λₜ₋₁`.
Note that weights are in the original (non-log) space, and should have temperature == 1.
# Examples
```julia
```

"""
function weightedESS(weights::AbstractVector{T}, λ::T, λₜ₋₁::T) where {T<:AbstractFloat}
    ArgCheck.@argcheck 0 < λₜ₋₁ <= λ <= 1.0
    num = sum(weights[iter]^(λ - λₜ₋₁) for iter in eachindex(weights))^2
    denom = sum(weights[iter]^(λ - λₜ₋₁)^2 for iter in eachindex(weights))
    return num/denom
end

############################################################################################
"""
$(SIGNATURES)
Closure to compute new temperature `λ` as a function of weights `weights` and current temperature `λₜ₋₁`.
`ESSTarget` is the benchmark to optimize against.

# Examples
```julia
```

"""
function get_ESSDifference(weights::AbstractVector{T}, λₜ₋₁::T, ESSTarget::T) where {T<:AbstractFloat}
    ArgCheck.@argcheck 0 < ESSTarget
    function ESSDifference(λ)
        return ESSTarget - weightedESS(weights, λ, λₜ₋₁)
    end
end

############################################################################################
"""
$(SIGNATURES)
Temperature adaption that takes into account updating schedule, see Chopin, Papaspiliopoulos (2020).
# Examples
```julia
```

"""
function update(weights::AbstractVector{T}, λₜ₋₁::T, ESSTarget::T) where {T<:AbstractFloat}
    # Assing ESS closure
    ESSDifference = get_ESSDifference(weights, λₜ₋₁, ESSTarget)
    # Check if ESS > ESSTarget; !NOTE: Compute ESS, not ESSDifference!
    if weightedESS(weights, T(1.0), λₜ₋₁) >= ESSTarget
        return T(1.0)
    end
    # Else tune for new stepsize
    λ = bisect(ESSDifference, λₜ₋₁, T(1.0))
    return λ
end

"Inverse CDF of Univariate Normal Mixture Model, where prob can be vectorized as a bivariate array of observations."
function bisect(f, min, max, atol = 10e-4, increasing = sign(f(max)), iter_max = 1000)
    c = Statistics.middle(min, max)
    z = f(c) * increasing
    iter = 0
    if z > 0 #
        maxi = c
        mini = typeof(maxi)(min)
    else
        mini = c
        maxi = typeof(mini)(max)
    end
    @inbounds while abs(mini - maxi) > atol && iter < iter_max
        iter += 1
        c = Statistics.middle(mini,maxi)
        if f(c) * increasing > 0 #
            maxi = c
        else
            mini = c
        end
    end
    #return output
    if iter < iter_max
        return Statistics.middle(mini,maxi)
    else
        return -Inf
    end
end

############################################################################################
"""
$(TYPEDEF)

Tuning container for joint uptdates of temperature.

# Fields
$(TYPEDFIELDS)
"""
struct JointTempering{A<:UpdateBool, T<:AbstractFloat} <: TemperingMethod
    "Checks if temperature will be updated after new iterations"
    adaption::A
    "current temperature."
    val::ValueHolder{T}
    "Target Effective Sample Size."
    ESSTarget::ValueHolder{T}
    "Buffer for weights used to adapt temperature"
    weights::Vector{T}
    function JointTempering(adaption::A, val::ValueHolder{T}, ESSTarget::ValueHolder{T}, weights::Vector{T}
    ) where {A<:UpdateBool,T<:AbstractFloat}
        @argcheck 0.0 < val.current <= 1.0
        @argcheck 0.0 < ESSTarget.current <= length(weights)
        return new{A,T}(adaption, val, ESSTarget, weights)
    end
end

function JointTempering(::Type{T}, adaption::B, val::F, ESSTarget::S, chains::Integer
) where {T<:AbstractFloat, B<:UpdateBool, F<:AbstractFloat, S<:AbstractFloat}
    # Return tuning struct
    return JointTempering(adaption, ValueHolder(T(val)), ValueHolder(T(ESSTarget)), zeros(T, chains))
end

############################################################################################
function update!(tempering::JointTempering, adaption::UpdateTrue, ℓweights::AbstractVector)
    # Update weights
    for iter in eachindex(tempering.weights)
        tempering.weights[iter] = ℓweights[iter]
    end
    #Normalize numerically stable
    tempering.weights .-= logsumexp(tempering.weights)
    tempering.weights .= exp.(tempering.weights)
    # Compute new temperature
    tempering.val.current = update(tempering.weights, tempering.val.current, tempering.ESSTarget.current)
    return tempering.val.current
end

############################################################################################
function update!(tempering::JointTempering, adaption::UpdateFalse, weights::AbstractVector)
    return tempering.val.current
end
function update!(tempering::JointTempering, weights::AbstractVector)
    return update!(tempering, tempering.adaption, weights)
end

############################################################################################
"""
$(SIGNATURES)
Return initial temperature of tuning container.

# Examples
```julia
```

"""
function initial(tempertune::JointTempering)
    return tempertune.val.current
end

############################################################################################
# Export

export
    JointTempering,
    initial,
    update,
    update!


text:
# foreach nonempty Period in a CompoundPeriod, get type of Period
typesof(x::CompoundPeriod) = map(typeof, x.periods)
typesof(x::P) where {P<:Period} = (P,)

# the extremal Periods with nonzero value
max(x::CompoundPeriod) = isempty(x) ? x : x.periods[1]
min(x::CompoundPeriod) = isempty(x) ? x : x.periods[end]
minmax(x::CompoundPeriod) = min(x), max(x)

maxtype(x::CompoundPeriod) = typeof(max(x))
mintype(x::CompoundPeriod) = typeof(min(x))

max(x::ReverseCompoundPeriod) = isempty(x) ? x : x.cperiod.periods[1]
min(x::ReverseCompoundPeriod) = isempty(x) ? x : x.cperiod.periods[end]
minmax(x::ReverseCompoundPeriod) = min(x), max(x)

maxtype(x::ReverseCompoundPeriod) = typeof(max(x))
mintype(x::ReverseCompoundPeriod) = typeof(min(x))

max(x::Period) = x
min(x::Period) = x
minmax(x::Period) = x, x

maxtype(x::T) where {T<:Period} = T
mintype(x::T) where {T<:Period} = T

Base.zero(::Type{CompoundPeriod}) = Nanosecond(0)
Base.zero(::Type{ReverseCompoundPeriod}) = Nanosecond(0)

# CompoundPeriod(::Time) is exported from Dates
CompoundPeriod(dt::Date) = Year(dt)+Month(dt)+Day(dt)
CompoundPeriod(dtm::DateTime) = CompoundPeriod(Date(dtm)) + CompoundPeriod(Time(dtm))

Time(tm::CompoundPeriod) = Time(hour(tm),minute(tm),second(tm),millisecond(tm),microsecond(tm),nanosecond(tm))
function Date(dt::CompoundPeriod)
    yr = year(dt)
    mo = max(1, month(dt))
    dy = max(1, day(dt))
    return Date(yr, mo, dy)
end
DateTime(dtm::CompoundPeriod)  = Date(dtm) + Time(hour(dtm), minute(dtm), second(dtm), millisecond(dtm))

signbit(period::Period) = signbit(period.value)
signbit(cperiod::CompoundPeriod) = signbit(max(canonical(cperiod)))
sign(cperiod::CompoundPeriod) = sign(max(canonical(cperiod)))


text:
immutable Gomadam2008_MaxSINRState
    U::Array{Matrix{Complex128},1}
    W::Array{Hermitian{Complex128},1} # these are only used for rate calculations
    V::Array{Matrix{Complex128},1}
end

function Gomadam2008_MaxSINR(channel, network)
    assignment = get_assignment(network)

    K = get_num_MSs(network)
    Ps = get_transmit_powers(network)
    sigma2s = get_receiver_noise_powers(network)
    ds = get_num_streams(network); max_d = maximum(ds)
    alphas = get_user_priorities(network)
    aux_params = get_aux_precoding_params(network)

    state = Gomadam2008_MaxSINRState(
        initial_receivers(channel, Ps, sigma2s, ds, assignment, aux_params),
        Array(Hermitian{Complex128}, K),
        initial_precoders(channel, Ps, sigma2s, ds, assignment, aux_params))
    objective = Float64[]
    logdet_rates = Array(Float64, K, max_d, aux_params["max_iters"])
    MMSE_rates = Array(Float64, K, max_d, aux_params["max_iters"])
    weighted_logdet_rates = Array(Float64, K, max_d, aux_params["max_iters"])
    weighted_MMSE_rates = Array(Float64, K, max_d, aux_params["max_iters"])
    allocated_power = Array(Float64, K, max_d, aux_params["max_iters"])

    iters = 0; conv_crit = Inf
    while iters < aux_params["max_iters"]
        update_MSs!(state, channel, sigma2s, assignment)
        iters += 1

        # Results after this iteration
        logdet_rates[:,:,iters] = calculate_logdet_rates(state)
        push!(objective, sum(logdet_rates[:,:,iters]))
        MMSE_rates[:,:,iters] = calculate_MMSE_rates(state)
        weighted_logdet_rates[:,:,iters] = calculate_weighted_logdet_rates(state, alphas)
        weighted_MMSE_rates[:,:,iters] = calculate_weighted_MMSE_rates(state, alphas)
        allocated_power[:,:,iters] = calculate_allocated_power(state)

        # Check convergence
        if iters >= 2
            conv_crit = abs(objective[end] - objective[end-1])/abs(objective[end-1])
            if conv_crit < aux_params["stop_crit"]
                Lumberjack.debug("Gomadam2008_MaxSINR converged.",
                    @compat Dict(
                        :num_iters => iters,
                        :final_objective => objective[end],
                        :conv_crit => conv_crit,
                        :stop_crit => aux_params["stop_crit"],
                        :max_iters => aux_params["max_iters"] ))
                break
            end
        end

        # Begin next iteration, unless the loop will end
        if iters < aux_params["max_iters"]
            update_BSs!(state, channel, Ps, sigma2s, assignment, aux_params)
        end
    end
    if iters == aux_params["max_iters"]
        Lumberjack.debug("Gomadam2008_MaxSINR did NOT converge.",
            @compat Dict(
                :num_iters => iters,
                :final_objective => objective[end],
                :conv_crit => conv_crit,
                :stop_crit => aux_params["stop_crit"],
                :max_iters => aux_params["max_iters"] ))
    end

    results = PrecodingResults()
    if aux_params["output_protocol"] == :all_iterations
        results["objective"] = objective
        results["logdet_rates"] = logdet_rates
        results["MMSE_rates"] = MMSE_rates
        results["weighted_logdet_rates"] = weighted_logdet_rates
        results["weighted_MMSE_rates"] = weighted_MMSE_rates
        results["allocated_power"] = allocated_power
    elseif aux_params["output_protocol"] == :final_iteration
        results["objective"] = objective[iters]
        results["logdet_rates"] = logdet_rates[:,:,iters]
        results["MMSE_rates"] = MMSE_rates[:,:,iters]
        results["weighted_logdet_rates"] = weighted_logdet_rates[:,:,iters]
        results["weighted_MMSE_rates"] = weighted_MMSE_rates[:,:,iters]
        results["allocated_power"] = allocated_power[:,:,iters]
    end
    return results
end

function update_MSs!(state::Gomadam2008_MaxSINRState,
    channel::SinglecarrierChannel, sigma2s, assignment)

    ds = [ size(state.V[k], 2) for k = 1:channel.K ]

    for i = 1:channel.I; for k in served_MS_ids(i, assignment)
        Phi = complex(sigma2s[k]*eye(channel.Ns[k]))
        for j = 1:channel.I; for l in served_MS_ids(j, assignment)
            F = channel.H[k,j]*state.V[l]
            Phi += F*F'
        end; end

        # Per-stream receivers
        for n = 1:ds[k]
            f = channel.H[k,i]*state.V[k][:,n]
            u = Hermitian(Phi - f*f')\f
            state.U[k][:,n] = u/norm(u,2)
        end

        # Optimal MSE weights (for rate calculation only)
        F = channel.H[k,i]*state.V[k]
        Ummse = Hermitian(Phi)\F
        Emmse = UniformScaling(1.) - Ummse'*F
        state.W[k] = inv(Hermitian((Emmse + Emmse')/2))
    end; end
end

function update_BSs!(state::Gomadam2008_MaxSINRState,
    channel::SinglecarrierChannel, Ps, sigma2s, assignment, aux_params)

    ds = [ size(state.W[k], 1) for k = 1:channel.K ]

    for i in active_BSs(assignment)
        Gamma = zeros(Complex128, channel.Ms[i], channel.Ms[i])
        for j = 1:channel.I; for l = served_MS_ids(j, assignment)
            G = channel.H[l,i]'*state.U[l]
            Gamma += G*G'
        end; end

        # Per-stream precoders
        served = served_MS_ids(i, assignment)
        Nserved = length(served)
        for k in served
            for n = 1:ds[k]
                g = channel.H[k,i]'*state.U[k][:,n]
                v = Hermitian(Gamma - g*g' + (sigma2s[k]/Ps[i])*eye(channel.Ms[i]))\g
                state.V[k][:,n] = sqrt(Ps[i]/(Nserved*ds[k]))*v/norm(v,2)
            end
        end
    end
end


text:
@snn_kw struct ADEXParameter{FT=Float32}
    a::FT = 4.0
    b::FT = 0.0805
    cm::FT = 0.281
    v_rest::FT = -70.6
    tau_m::FT = 9.3667
    tau_w::FT = 144.0
    v_thresh::FT = -50.4
    delta_T::FT = 2.0
    v_spike::FT = -40.0
    v_reset::FT = -70.6
    spike_delta::FT = 30
end
@snn_kw mutable struct AD{VFT=Vector{Float32},VBT=Vector{Bool}}
    param::ADEXParameter = ADEXParameter()

    N::Int32 = 1
    cnt::Int32 = 1
    v::VFT = fill(param.v_rest, N)
    w::VFT = zeros(N)
    fire::VBT = zeros(Bool, N)
    I::VFT = zeros(N)
    sized::Int32 = 1
    spike_raster::Vector{Int32} = zeros(N)
    records::Dict = Dict()
end

function integrate!(p::AD, param::ADEXParameter, dt::Float32)
    @unpack N, cnt, v, w, fire, I,spike_raster,sized = p
    @unpack a,b,cm,v_rest,tau_m,tau_w,v_thresh,delta_T,v_spike,v_reset,spike_delta = param
    @inbounds for i = 1:N
        if spike_raster[cnt] == 1 || fire[i]
          v[i] = v_reset
          w[i] += b
        end
        dv  = (((v_rest-v[i]) +
                delta_T*exp((v[i] - v_thresh)/delta_T))/tau_m +
                (I[i] - w[i])/cm) *dt
        v[i] += dv
        w[i] += dt * (a*(v[i] - v_rest) - w[i])/tau_w * dt
        fire[i] = v[i] > v_thresh
        if v[i]>v_thresh
            fire[i] = 1
            v[i] = spike_delta
            spike_raster[cnt] = 1

        else
            spike_raster[cnt] = 0
        end
    end
    cnt+=1
end

max_stars_repo_name: russelljjarvis/backupSNN

max_issues_repo_name: russelljjarvis/backupSNN

max_forks_repo_name: russelljjarvis/backupSNN


text:
if Sys.islinux()
    using ImageMagick
end
import APTDecoder
using Test
using RemoteFiles
using SatelliteToolbox

@testset "decoding" begin
    wavname = "gqrx_20190825_182745_137620000.wav"

    if !isfile(wavname)
        download("https://archive.org/download/gqrx_20190825_182745_137620000/gqrx_20190825_182745_137620000.wav",wavname)
    end
    current_tles = APTDecoder.get_tle(:weather)
    tles = SatelliteToolbox.read_tle(joinpath(dirname(pathof(APTDecoder)),"..","examples","weather-20190825.txt"))

    APTDecoder.makeplots(wavname,"NOAA 15", tles = tles)

    @test isfile(replace(wavname,r"\.wav$" => "_raw.png"))

    s = randn(4,5);
    APTDecoder.reshape_signal(s[:],4) == s

    ioffset = 6
    frame_len = 1000
    sync_frame = APTDecoder.gen_sync_frame(3*4160,APTDecoder.sync_frequency)[1];
    y_demod = vcat(ones(Int,ioffset),sync_frame,ones(Int,7));
    sync_frame_index = APTDecoder.find_sync(y_demod,sync_frame,frame_len)
    @test sync_frame_index == [ioffset+1]

    tt = APTDecoder.mark_sync(y_demod,sync_frame,frame_len);
    @test tt[ioffset+1,1] ≈ 1.

    @test_throws ErrorException APTDecoder.starttimename("bogous_filename.wav")

    pngname = joinpath(dirname(@__FILE__),"gqrx_20190825_182745_137620000.png")
    datatime,channel,data = APTDecoder.wxload(pngname)
    @test length(channel) == 2

    satellite_name = "NOAA 15"
    channel = 'a'
    plon,plat,data = APTDecoder.georeference(pngname,satellite_name,channel, tles = tles)

    satellite_name = "NOAA 15"
    channel = 'x' # bogous
    @test_throws ErrorException APTDecoder.georeference(pngname,satellite_name,channel)
end

@testset "GeoMapping" begin
    @test APTDecoder.GeoMapping.azimuth(0,0,1,0) ≈ 0;
    @test APTDecoder.GeoMapping.azimuth(0,0,0,1) ≈ 90;

    @test APTDecoder.GeoMapping.distance(0,0,1,0) ≈ 1;
    @test APTDecoder.GeoMapping.distance(0,0,90,0) ≈ 90;
    @test APTDecoder.GeoMapping.distance(0,0,0,180) ≈ 180;
    @test APTDecoder.GeoMapping.distance(0,0,0,270) ≈ 90;
    @test APTDecoder.GeoMapping.distance(1,2,3,4) ≈ 2.82749366820155;

    @test APTDecoder.GeoMapping.distance(43.016666412353515625,3.3333332538604736328125,43.016666412353515625,3.3333332538604736328125) ≈ 0

    lato,lono = APTDecoder.GeoMapping.reckon(30,-80,20,40);
    @test lato ≈ 44.16661401448592
    @test lono ≈ -62.15251496909770

    lato,lono = APTDecoder.GeoMapping.reckon(-30,80,[5, 10],[40, 45]);
    @test lato ≈ [-26.12155703039504, -22.70996703614572]
    @test lono ≈ [83.57732793979254,  87.64920016442251]

    lato,lono = APTDecoder.GeoMapping.reckon([-30, 31],[80, 81],[5, 10],[40, 45]);
    @test lato ≈ [-26.12155703039504,  37.76782079033356]
    @test lono ≈ [83.57732793979254,  89.93590456974810]
end

max_stars_repo_name: Alexander-Barth/APTDecoder

max_issues_repo_name: Alexander-Barth/APTDecoder

max_forks_repo_name: Alexander-Barth/APTDecoder


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Baseline model
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, RData,  XLSX, LinearAlgebra, Statistics, CSV

dir = "X:/VIVES/1-Personal/Florian/git/IOTables/src/"
include(dir * "model/transform_data.jl") # Script with functions to import and transform raw data
include(dir * "model/price_hat.jl") # Script with function to obtain the price index
include(dir * "model/wage_hat.jl") # Script with function to obtain the wages and gross output
include(dir * "counterfactual/tariffs_function.jl") # Script with functions to create τ_hat_Z, τ_hat_F from tariff data
include(dir * "counterfactual/head_ries_index.jl") # Script with functions to create bilateral Head-Ries index (symmetric bilateral trade costs)
include(dir * "elasticity/elasticities_functions.jl") # Script with functions to compute statistics for estimating sectoral trade elasticities (method of Caliendo and Parro, 2015)

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

dir = "C:/Users/u0148308/Desktop/raw/" # location of raw data

# Data specification

# WIOD rev. 2013
source = "WIOD"
revision = "2013"
year = 1995 # specified year
N = 41 # number of countries
S = 35 # number of industries

# # WIOD rev. 2016
# source = "WIOD"
# revision = "2016"
# year = 2014 # specified year
# N = 44 # number of countries
# S = 56 # number of industries

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare tariff data from WTO
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_tariffs = DataFrame(reporter_lab=String[], year=Int64[], HS2017_code=String[],
    HS2017_lab=String[], duty_lab=String[], avg_tariff=Union{Float64,Missing}[], reporter=String[])

function tariff_data(reporter::String, id::String, year::String)

    col_names_old = ["Reporter", "Year", "Product Code", "Product Description", "Duty Description", "All items excl. NA - Simple Averages"]
    col_names_new = ["reporter_lab", "year", "HS2017_code", "HS2017_lab", "duty_lab", "avg_tariff"]

    path = dir * "MFN/" * reporter * "_" * year * "_TA/" * reporter * "_" * year * "_Duties_TA.txt"

    df = CSV.read(path, DataFrame)
    df = df[:, col_names_old]
    rename!(df, col_names_new)
    df[:, :reporter] .= id

    return df
end

ctry_names_2013 = ["AUS", "AUT", "BEL", "BGR", "BRA", "CAN", "CHN", "CYP", "CZE", "DEU", "DNK", "ESP", "EST", "FIN", "FRA", "GBR",
    "GRC", "HUN", "IDN", "IND", "IRL", "ITA", "JPN", "KOR", "LTU", "LUX", "LVA", "MEX", "MLT", "NLD", "POL", "PRT", "ROM", "RUS",
    "SVK", "SVN", "SWE", "TUR", "TWN", "USA"]

ctry_EU28 = ["AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GBR", "GRC", "HUN", "IRL", "ITA", "LVA",
    "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROM", "SVK", "SVN", "ESP", "SWE"]

for i in ctry_names_2013

    if i in ctry_EU28
        df = tariff_data("EU", i, "2018") # for EU countries
    elseif i in ["RUS", "TUR"]
        df = tariff_data(i, i, "2019") # countries for which tariffs in 2018 were not available
    else
        df = tariff_data(i, i, "2018")
    end

    append!(df_tariffs, df)
end

# Correspondence table HS2017, ISIC rev. 3 and ISIC rev. 4
correspondence_HS_ISIC = DataFrame(XLSX.readtable(dir * "correspondence_HS_ISIC.xlsx", "Sheet1")...)
transform!(correspondence_HS_ISIC, [:ISIC3, :ISIC4, :HS2017] .=> ByRow(x-> lpad(x,2,'0')), renamecols=false)
rename!(correspondence_HS_ISIC, [:HS2017, :ISIC3, :ISIC4] .=> [:HS2017_code, :ISIC3_code, :ISIC4_code])

# join tables
transform!(df_tariffs, :HS2017_code => ByRow(x->replace(x, "'" => "")), renamecols=false) # remove ' from column otherwise we cannot join properly
df = leftjoin(df_tariffs, correspondence_HS_ISIC, on=:HS2017_code)

# -------

# WIOD rev. 2013
gdf = groupby(df, [:reporter, :ISIC3_code]) # group to compute average tariffs according to ISIC rev. 3 classification
df_tariffs = combine(gdf, :avg_tariff => mean, renamecols=false)
sort!(df_tariffs, [:reporter, :ISIC3_code])
dropmissing!(df_tariffs) # drop missing for using average later

# use average of all countries for RoW
gdf = groupby(df_tariffs, :ISIC3_code)
df_RoW = combine(gdf, :avg_tariff => mean, renamecols=false)
df_RoW[:, :reporter] .= "ZoW"

df_all = vcat(df_tariffs, df_RoW)

reporter = [ctry_names_2013; "ZoW"] # use "ZoW" instead of "RoW" so we can sort later on
industries = lpad.(1:S, 2, '0') # left pad with leading zeros so we can sort later on

gdf = groupby(df_all, :reporter)
df_ctry_avg = combine(gdf, :avg_tariff => mean, renamecols=false)
avg_ctry_tariff = df_ctry_avg.avg_tariff # use average country tariffs for sectors which are missing

df = DataFrame([repeat(reporter, inner=S) repeat(industries, outer=N) repeat(avg_ctry_tariff, inner=S)], [:reporter, :ISIC3_code, :ctry_avg])
df = leftjoin(df, df_all, on=[:reporter, :ISIC3_code])
df.avg_tariff = ifelse.(ismissing.(df.avg_tariff) .== true, df.ctry_avg, df.avg_tariff)
rename!(df, :ISIC3_code => :industry)
df_tariffs = df[:, [:reporter, :industry, :avg_tariff]]
sort!(df_tariffs, [:reporter, :industry])

CSV.write(dir * "MFN/df_tariffs.csv", df_tariffs)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# format data to NS×N matrix
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

vec_tariffs = df_tariffs.avg_tariff
matrix_tariffs = repeat(vec_tariffs, 1, N)
countries = [ctry_names_2013; "ZoW"]

τ_Z = ones(N*S,N)

for (i, ctry_i) in enumerate(countries)
    for (j, ctry_j) in enumerate(countries)

        if ctry_i == ctry_j
            τ_Z[1+(i-1)*S:i*S,j] .= 0

        elseif ctry_i in ctry_EU28 && ctry_j in ctry_EU28
            τ_Z[1+(i-1)*S:i*S,j] .= 0

        else
            τ_Z[1+(i-1)*S:i*S,j] = matrix_tariffs[1+(i-1)*S:i*S,j]
        end

    end
end

# export
rows = repeat(countries, inner=S) .* "__" .* repeat(industries, outer=N)
df_tariff_matrix = DataFrame([rows τ_Z], ["reporter"; countries])
XLSX.writetable(dir * "MFN/tariff_matrix.xlsx", df_tariff_matrix, overwrite=true)

max_stars_repo_path: old stuff/tariff_data.jl

max_issues_repo_path: old stuff/tariff_data.jl

max_forks_repo_path: old stuff/tariff_data.jl


text:
# these functions are more useful for interactive usage of the code

"""
  Print all initial condition names for a given physics module

  Inputs:

    mod: the module to print the ICs for
    f: the IO to print to, defaults to STDOUT

  Outputs:
    none
"""
function printICNames(mod::Module, f::IO=STDOUT)

  println(f, "IC names for module ", mod, ":")
  for name in keys(mod.ICDict)
    println(f, "  ", name)
  end

  return nothing
end

"""
  Like printICNames, but for boundary conditions, see that function for details
"""
function printBCNames(mod::Module, f::IO=STDOUT)

  println(f, "BC names for module ", mod, ":")
  for name in keys(mod.BCDict)
    println(f, "  ", name)
  end

  return nothing
end

"""
  Prints the name of all currently registered physics modules and the module
  name itself
"""
function printPhysicsModules(f::IO=STDOUT)

  n = length(keys(PhysicsModDict))
  println(f, "Currently registered physics modules: ", n)
  for (name, val) in PhysicsModDict
    mod = val[1]
    println(f, "  ", name, " => ", mod)
  end

  return nothing
end


text:
###
### Interface Implementations for running full sampling loops
###

using Printf

StatsBase.sample(rng::AbstractRNG, model::AbstractModel, sampler::Nested; kwargs...) =
    mcmcsample(rng, model, sampler, nested_isdone; progressname="Nested Sampling", chain_type=Chains, kwargs...)

StatsBase.sample(model::AbstractModel, sampler::Nested; kwargs...) =
    StatsBase.sample(GLOBAL_RNG, model, sampler; kwargs...)

function nested_isdone(rng, model, smapler, samples, state, i; progress=true, maxiter=Inf, maxcall=Inf, dlogz=0.5, maxlogl=Inf, kwargs...)
    # 1) iterations exceeds maxiter
    done_sampling = state.it > maxiter
    # 2) number of loglike calls has been exceeded
    done_sampling |= state.ncall > maxcall
    # 3) remaining fractional log-evidence below threshold
    logz_remain = maximum(state.logl) + state.logvol
    delta_logz = logaddexp(state.logz, logz_remain) - state.logz
    done_sampling |= delta_logz < dlogz
    # 4) last dead point loglikelihood exceeds threshold
    done_sampling |= state.logl_dead > maxlogl
    # 5) number of effective samples
    # TODO

    if progress
        str = @sprintf "iter=%d\tncall=%d\tΔlogz=%.2g\tlogl=%.2g\tlogz=%.2g" i state.ncall delta_logz state.logl_dead state.logz
        print("\r\33[2K", str)
    end

    return done_sampling
end


text:
# MOI.Indicator
function MOI.supports_constraint(
    o::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{S},
) where {A, S2 <: MOI.AbstractSet, S <: MOI.Indicator{A, S2}}
    # No information about the dimension of the function, only works if 1 + dim(S2).
    return MOI.supports_constraint(o, MOI.VariableIndex, S2) ||
           MOI.supports_constraint(o, MOI.VectorOfVariables, S2)
end

function MOI.supports_constraint(
    o::Optimizer,
    ::Type{MOI.VectorAffineFunction{T}},
    ::Type{S},
) where {
    T <: Int,
    A,
    S2 <: MOI.AbstractSet,
    S <: MOI.Indicator{A, S2},
}
    # No information about the dimension of the function, only works if 1 + dim(S2).
    return MOI.supports_constraint(o, MOI.ScalarAffineFunction{T}, S2) ||
           MOI.supports_constraint(o, MOI.VectorAffineFunction{T}, S2)
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::MOI.Indicator{A, S2},
) where {A, S2 <: MOI.AbstractSet, T <: Int}
    # Split the dimensions in the right parts. Only parse the first component,
    # as the rest will be handled by the reified constraint.
    dim_set = MOI.dimension(s.set)
    f_parsed = _parse(model, f)
    reify_indicator = f_parsed[1]
    reify_set_variables = _slice(f, 2:(dim_set + 1))

    # Ensure that the indicator is a binary variable.
    # TODO: allow this constraint to be deleted at the same time as the reified constraint.
    cpo_java_add(
        model.inner,
        cpo_java_range(model.inner, 0, reify_indicator, 1),
    )

    # Build the constraint.
    indicator_value = if A == MOI.ACTIVATE_ON_ONE
        Int32(true)
    else
        Int32(false)
    end
    indicator = cpo_java_eq(model.inner, reify_indicator, indicator_value)
    set = _build_constraint(model, reify_set_variables, s.set)
    return cpo_java_ifthen(model.inner, indicator, set)
end

# CP.Reification
function MOI.supports_constraint(
    o::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{S},
) where {S2 <: MOI.AbstractSet, S <: CP.Reification{S2}}
    # No information about the dimension of the function, only works if 1 + dim(S2).
    return MOI.supports_constraint(o, MOI.VariableIndex, S2) ||
           MOI.supports_constraint(o, MOI.VectorOfVariables, S2)
end

function MOI.supports_constraint(
    o::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: MOI.VectorAffineFunction{T},
    S2 <: MOI.AbstractSet,
    S <: CP.Reification{S2},
}
    # No information about the dimension of the function, only works if 1 + dim(S2).
    return MOI.supports_constraint(o, MOI.ScalarAffineFunction{T}, S2) ||
           MOI.supports_constraint(o, MOI.VectorAffineFunction{T}, S2)
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.Reification{S2},
) where {S2 <: MOI.AbstractSet, T <: Int}
    # Split the dimensions in the right parts. Only parse the first component,
    # as the rest will be handled by the reified constraint.
    dim_set = MOI.dimension(s.set)
    f_parsed = _parse(model, f)
    reify_indicator = f_parsed[1]
    reify_set_variables = _slice(f, 2:(dim_set + 1))

    # Ensure that the indicator is a binary variable.
    # TODO: allow this constraint to be deleted at the same time as the reified constraint.
    cpo_java_add(
        model.inner,
        cpo_java_range(model.inner, 0, reify_indicator, 1),
    )

    # Build the constraint.
    indicator = cpo_java_eq(model.inner, reify_indicator, Int32(true))
    set = _build_constraint(model, reify_set_variables, s.set)
    return cpo_java_equiv(model.inner, indicator, set)
end

# CP.Equivalence
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S1 <: MOI.AbstractSet,
    S2 <: MOI.AbstractSet,
    S <: CP.Equivalence{S1, S2},
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.Equivalence{S1, S2},
) where {T <: Int, S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet}
    dim_first = MOI.dimension(s.set1)
    dim_second = MOI.dimension(s.set2)
    dim_end = dim_first + dim_second

    first = _slice(f, 1:dim_first)
    second = _slice(f, (1 + dim_first):dim_end)

    return cpo_java_equiv(
        model.inner,
        _build_constraint(model, first, s.set1),
        _build_constraint(model, second, s.set2),
    )
end

# CP.EquivalenceNot
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S1 <: MOI.AbstractSet,
    S2 <: MOI.AbstractSet,
    S <: CP.EquivalenceNot{S1, S2},
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.EquivalenceNot{S1, S2},
) where {T <: Int, S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet}
    dim_first = MOI.dimension(s.set1)
    dim_second = MOI.dimension(s.set2)
    dim_end = dim_first + dim_second

    first = _slice(f, 1:dim_first)
    second = _slice(f, (1 + dim_first):dim_end)

    return cpo_java_neq(
        model.inner,
        _build_constraint(model, first, s.set1),
        _build_constraint(model, second, s.set2),
    )
end

# CP.IfThenElse
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S1 <: MOI.AbstractSet,
    S2 <: MOI.AbstractSet,
    S3 <: MOI.AbstractSet,
    S <: CP.IfThenElse{S1, S2, S3},
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.IfThenElse{S1, S2, S3},
) where {
    T <: Int,
    S1 <: MOI.AbstractSet,
    S2 <: MOI.AbstractSet,
    S3 <: MOI.AbstractSet,
}
    dim_first = MOI.dimension(s.condition)
    dim_second = MOI.dimension(s.true_constraint)
    dim_third = MOI.dimension(s.false_constraint)
    dim_end = dim_first + dim_second + dim_third

    first = _slice(f, 1:dim_first)
    second = _slice(f, (1 + dim_first):(dim_first + dim_second))
    third = _slice(f, (1 + dim_first + dim_second):dim_end)

    return cpo_java_ifthenelse(
        model.inner,
        _build_constraint(model, first, s.condition),
        _build_constraint(model, second, s.true_constraint),
        _build_constraint(model, third, s.false_constraint),
    )
end

# CP.Implication
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S1 <: MOI.AbstractSet,
    S2 <: MOI.AbstractSet,
    S <: CP.Implication{S1, S2},
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.Implication{S1, S2},
) where {T <: Int, S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet}
    dim_first = MOI.dimension(s.antecedent)
    dim_second = MOI.dimension(s.consequent)
    dim_end = dim_first + dim_second

    first = _slice(f, 1:dim_first)
    second = _slice(f, (1 + dim_first):dim_end)

    return cpo_java_imply(
        model.inner,
        _build_constraint(model, first, s.antecedent),
        _build_constraint(model, second, s.consequent),
    )
end

# CP.True
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S <: CP.True,
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.True,
) where {T <: Int}
    return cpo_java_trueconstraint(model.inner)
end

# CP.False
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    T <: Int,
    F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    S <: CP.False,
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
    s::CP.False,
) where {T <: Int}
    return cpo_java_falseconstraint(model.inner)
end


text:
## AlmostBandedMatrix

struct AlmostBandedMatrix{T} <: AbstractMatrix{T}
    bands::BandedMatrix{T}
    fill::LowRankMatrix{T}
    function AlmostBandedMatrix{T}(bands::BandedMatrix{T}, fill::LowRankMatrix{T}) where T
        if size(bands) ≠ size(fill)
            error("Data and fill must be compatible size")
        end
        new{T}(bands,fill)
    end
end

AlmostBandedMatrix(bands::BandedMatrix, fill::LowRankMatrix) =
    AlmostBandedMatrix{promote_type(eltype(bands),eltype(fill))}(bands,fill)

AlmostBandedMatrix{T}(nm::NTuple{2,Integer}, lu::NTuple{2,Integer}, r::Integer) where {T} =
    AlmostBandedMatrix(BandedMatrix{T}(nm,lu), LowRankMatrix{T}(nm,r))

AlmostBandedMatrix{T}(Z::Zeros, lu::NTuple{2,Integer}, r::Integer) where {T} =
    AlmostBandedMatrix(BandedMatrix{T}(Z, lu), LowRankMatrix{T}(Z, r))

AlmostBandedMatrix(Z::AbstractMatrix, lu::NTuple{2,Integer}, r::Integer) =
    AlmostBandedMatrix{eltype(Z)}(Z, lu, r)

for MAT in (:AlmostBandedMatrix, :AbstractMatrix, :AbstractArray)
    @eval convert(::Type{$MAT{T}}, A::AlmostBandedMatrix) where {T} =
        AlmostBandedMatrix(AbstractMatrix{T}(A.bands),AbstractMatrix{T}(A.fill))
end

size(A::AlmostBandedMatrix) = size(A.bands)
Base.IndexStyle(::Type{ABM}) where {ABM<:AlmostBandedMatrix} =
    IndexCartesian()

function getindex(B::AlmostBandedMatrix,k::Integer,j::Integer)
    if j > k + bandwidth(B.bands,2)
        B.fill[k,j]
    else
        B.bands[k,j]
    end
end

# can only change the bands, not the fill
function setindex!(B::AlmostBandedMatrix,v,k::Integer,j::Integer)
        B.bands[k,j] = v
end

function pad!(B::AlmostBandedMatrix,n::Integer,m::Integer)
    pad!(B.bands,n,m)
    pad!(B.fill,n,m)
    B
end