text:
"""
    RectRegion(
        RegID :: AbstractString,
        ParID :: AbstractString,
        name  :: AbstractString,
        bound :: Vector{<:Real},
        savegeo :: Bool = true,
        ST = String,
        FT = Float64
    ) -> RectRegion{ST,FT}

Creates a rectilinear GeoRegion `RegID`.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the GeoRegion.
            If the ID is already in use, then an error will be thrown.

- `ParID` : The ID of the parent GeoRegion where information can be extracted from
- `name`  : A name for the GeoRegion (meta information, can be used in Logging)
- `bound` : The [N,S,E,W] coordinates defining the region
- `savegeo` : Save the GeoRegion into the master list?
"""
function RectRegion(
    RegID :: AbstractString,
    ParID :: AbstractString,
    name  :: AbstractString,
    bound :: Vector{<:Real};
    savegeo :: Bool = true,
    ST = String,
    FT = Float64
)

    if isGeoRegion(RegID,throw=false)
        error("$(modulelog()) - The GeoRegion $(RegID) has already been defined.  Please use another identifier.")
    else
        @info "$(modulelog()) - Adding the GeoRegion $(RegID) to the list."
    end

    if !isGeoRegion(ParID,throw=false)
        error("$(modulelog()) - The GeoRegion $(ParID) was defined to be the parent GeoRegion of $(RegID), but the GeoRegion $(ParID) is not defined.  Please define the GeoRegion $(ParID) and its properties.")
    end

    N,S,E,W = bound; is180,is360 = checkbounds(N,S,E,W)
    geo  = RectRegion{ST,FT}(RegID,ParID,name,N,S,E,W,is180,is360)
    par  = GeoRegion(ParID); isinGeoRegion(geo,par)
    name = replace(name," "=>"-")

    if savegeo
        open(joinpath(DEPOT_PATH[1],"files","GeoRegions","rectlist.txt"),"a") do io
            write(io,"$RegID, $ParID, $N, $W, $S, $E, $name\n")
        end
    end

    return geo

end

"""
    PolyRegion(
        RegID :: AbstractString,
        ParID :: AbstractString,
        name  :: AbstractString,
        lonpt :: Vector{<:Real},
        latpt :: Vector{<:Real};
        savegeo :: Bool = true,
        ST = String,
        FT = Float64
    ) -> PolyRegion{ST,FT}

Creates a rectilinear GeoRegion `RegID`.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the GeoRegion.
            If the ID is already in use, then an error will be thrown.

- `ParID` : The ID of the parent GeoRegion where information can be extracted from
- `name`  : A name for the GeoRegion (meta information, can be used in Logging)
- `lonpt` : A vector containing the longitude points
- `latpt` : A vector containing the latitude points
- `savegeo` : Save the GeoRegion into the master list?

!!! info "Start and End Points"
    The 1st and last elements of `lonpt` and `latpt` must be equal.
"""
function PolyRegion(
    RegID :: AbstractString,
    ParID :: AbstractString,
    name  :: AbstractString,
    lonpt :: Vector{<:Real},
    latpt :: Vector{<:Real},;
    savegeo :: Bool = true,
    ST = String,
    FT = Float64
)

    if isGeoRegion(RegID,throw=false)
        error("$(modulelog()) - The GeoRegion $(RegID) has already been defined.  Please use another identifier.")
    else
        @info "$(modulelog()) - Adding the GeoRegion $(RegID) to the list."
    end

    if !isGeoRegion(ParID)
        error("$(modulelog()) - The GeoRegion $(ParID) was defined to be the parent GeoRegion of $(RegID), but the GeoRegion $(ParID) is not defined.  Please define the GeoRegion $(ParID) and its properties.")
    end

    if lonpt[1] != lonpt[end]
        error("$(modulelog()) - The first and last points of the vector of longitude coordinates must be the same")
    end

    if latpt[1] != latpt[end]
        error("$(modulelog()) - The first and last points of the vector of latitude coordinates must be the same")
    end

    N = maximum(latpt); S = minimum(latpt)
    E = maximum(lonpt); W = minimum(lonpt)
    is180,is360 = checkbounds(N,S,E,W)
    shape = Point2.(lonpt,latpt)
    geo   = PolyRegion{ST,FT}(RegID,ParID,name,N,S,E,W,shape,is180,is360)
    par   = GeoRegion(ParID); isinGeoRegion(geo,par)
    name  = replace(name," "=>"-")

    npt = length(lonpt)
    if savegeo
        open(joinpath(DEPOT_PATH[1],"files","GeoRegions","polylist.txt"),"a") do io
            write(io,"\n  RegID $(RegID), $(ParID), $name\n")
            write(io,"  RegX  $(lonpt[1])")
            for ipt = 2 : npt; write(io,", $(lonpt[ipt])") end
            write(io,"\n")
            write(io,"  RegY  $(latpt[1])")
            for ipt = 2 : npt; write(io,", $(latpt[ipt])") end
            write(io,"\n")
        end
    end

    return geo

end

"""
    removeGeoRegion(RegID::AbstractString)

Creates the GeoRegion associated with the ID `RegID`.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the GeoRegion.
            If the ID is not valid (i.e. not being used), then an error will be thrown.
"""
function removeGeoRegion(
    RegID :: AbstractString
)

    if RegID == "GLB"
        error("$(modulelog()) - The Global GeoRegion \"GLB\" is an integral part of the GeoRegions.jl package and cannot be removed.")
    else
        @info "$(modulelog()) - Removing the GeoRegion $(RegID) ..."
    end

    regvec,filevec,typevec = listGeoRegions(); isgeoregion(RegID,regvec)
    ind = findall(RegID.==regvec)[1]
    geo = getgeoregion(
        RegID,
        joinpath(DEPOT_PATH[1],"files","GeoRegions",filevec[ind]),
        typevec[ind]
    )
    removegeoregion(geo,filevec[ind])

end

function removegeoregion(
    geo  :: PolyRegion,
    fgeo :: AbstractString
)

    rvec = listpolyregions(joinpath(DEPOT_PATH[1],"files","GeoRegions",fgeo))
    ind  = findall(rvec.==geo.regID)[1]
    ind  = (ind) * 4 .+ (0:3)

    flines = readlines(joinpath(DEPOT_PATH[1],"files","GeoRegions",fgeo))
    nlines = length(flines)

    open("tmp.txt","w") do io
        for iline = 1 : nlines
            if !any(iline.==ind)
                write(io,"$(flines[iline])\n")
            end
        end
    end

    mv("tmp.txt",joinpath(DEPOT_PATH[1],"files","GeoRegions",fgeo),force=true)

    return nothing

end

function removegeoregion(
    geo  :: RectRegion,
    fgeo :: AbstractString
)

    flines = readlines(joinpath(DEPOT_PATH[1],"files","GeoRegions",fgeo))
    nlines = length(flines)

    open("tmp.txt","w") do io
        for iline = 1 : nlines
            if !occursin(geo.regID,flines[iline])
                write(io,"$(flines[iline])\n")
            end
        end
    end

    mv("tmp.txt",joinpath(DEPOT_PATH[1],"files","GeoRegions",fgeo),force=true)

    return nothing

end


text:
using OperatingSystems
using Base.Test

@test CURRENT_OS <: AbstractOS

if is_windows()
    @test CURRENT_OS == Windows
else
    @test CURRENT_OS <: Unix
    if is_apple()
        @test CURRENT_OS == MacOS
        @test CURRENT_OS <: BSD
    elseif is_linux()
        @test CURRENT_OS == Linux
    elseif is_bsd()
        @test CURRENT_OS == eval(OperatingSystems, Sys.KERNEL)
        @test CURRENT_OS <: BSD
    end
end

ext: jl


text:
using FFTW, LinearAlgebra, Plots
using BenchmarkTools
gr()

include("src/dataset.jl")
include("src/micmac.jl")
include("src/error.jl")

const epsilon  = 1.0
const xmin     = 0
const xmax     = 2π
const T        = 2π
const nx       = 256
const ntau     = 128
const Tfinal   = 0.25

println(" epsilon : $epsilon ")
println(" nx      : $nx ")
println(" ntau    : $ntau ")

data = CosSin(xmin, xmax, nx, epsilon, T, Tfinal)

dt = 2.0^(-4) * data.Tfinal / 16

solver = MicMac(data, ntau)

@time u, v = solve(solver, dt)
err = compute_error(u,v,data)
println(" error = $err ")


text:
using ArchGDAL
using StaticArrays

struct Geo
    origin::Array{Float64, 1}
    transform::Array{Float64, 2}
end

"""
Find the origin and linear transform from pixel to x-y.
"""
function band_geo(geo)
    Geo([geo[1], geo[4]], [geo[2] geo[3]; geo[5] geo[6]])
end

pixel_corner(geo::Geo, ij) = geo.transform * (ij .- 1) .+ geo.origin
pixel_center(geo::Geo, ij) = geo.transform * (ij .- 0.5) .+ geo.origin
# First column is upper-left (x y). Second column is lower-right (x, y).
raster_rect(geo::Geo, ullr) = pixel_corner(geo, hcat(ullr[:, 1], ullr[:, 2] .+ 1))
partial_pixel_containing(geo::Geo, xy) = (geo.transform \ (xy - geo.origin)) .+ 1
offset_geo(geo::Geo, offset) = Geo(pixel_corner(geo, offset), geo.transform)
pixel_rect(geo::Geo, index) = pixel_corner(geo, [index[1] (index[1] + 1); index[2] (index[2] + 1)])

function cover_rect(rect, geo, bounds)
    # This gives you floating point pixel that contain the rectangle.
    ul_ij = partial_pixel_containing(geo, rect[:, 1])
    lr_ij = partial_pixel_containing(geo, rect[:, 2])
    # Include an upper-left if the center of the pixel is included.
    # 2.2 is 2. 2.7 is 3.
    ul = Int32.(floor.(ul_ij .+ 0.5))
    # For lower-right, 10.2 is 9. 10.7 is 10.
    lr = Int32.(floor.(lr_ij .- 0.5))
    ul = max.(ul, 1)
    lr = min.(lr, bounds)
    if lr[1] < ul[1] || lr[2] < ul[2]
        return CartesianIndices((0, 0))
    else
        return CartesianIndices((ul[1]:lr[1], ul[2]:lr[2]))
    end
end

function over_sweep!(geo, coarse, coarse_missing, fine, fine_missing)
    # fine_rect = raster_rect(geo[2], [1 size(fine, 1); 1 size(fine, 2)])
    for cidx in CartesianIndices(coarse)
        coarse_value = coarse[cidx]
        if coarse_value == coarse_missing
            continue
        end
        @assert coarse_value >= 0
        rect = pixel_rect(geo[1], cidx)
        fine_indices = cover_rect(rect, geo[2], size(fine))
        if length(fine_indices) == 0
            continue
        end
        nonzero = sum([fine[fidx] != fine_missing for fidx in fine_indices])
        if nonzero > 0
            rate = min(100, coarse_value / nonzero)
            for rmidx in fine_indices
                if fine[rmidx] != fine_missing
                    fine[rmidx] = rate
                    # else leave it missing
                end
            end
        end
    end
    return nothing
end

function whole_band(a::ArchGDAL.Dataset, b::ArchGDAL.IDataset)
    band = [ArchGDAL.getband(a, 1), ArchGDAL.getband(b, 1)]
    raw_geo = ArchGDAL.getgeotransform.([a, b])
    geo = band_geo.(raw_geo)
    coarse = ArchGDAL.read(band[1])
    coarse_no_data = ArchGDAL.getnodatavalue(band[1])
    fine = ArchGDAL.read(band[2])
    fine_no_data = ArchGDAL.getnodatavalue(band[2])
    over_sweep!(geo, coarse, coarse_no_data, fine, fine_no_data)
    ArchGDAL.write!(band[2], fine)
end

"""
     assignweights(points_path, weight_path, reweighted_path)

Given the HRSL data, the LandScan data, and a path where you want
a new dataset, this makes a new dataset at `reweighted_path` that has
a point for every HRSL point but is weighted by the LandScan weights.
"""
function assignweights(points_path, weight_path, reweighted_path)
    # ds = dataset
    fine_copy_ds = ArchGDAL.read(expanduser(points_path)) do fine_ds
        reweighted_path = expanduser(reweighted_path)
        ArchGDAL.copy(fine_ds, filename = reweighted_path)
    end
    ArchGDAL.read(expanduser(weight_path)) do coarse_ds
        whole_band(coarse_ds, fine_copy_ds)
    end
    ArchGDAL.destroy(fine_copy_ds)
end


text:
################################################
# Unconstrained
################################################
function fitnessToMass(fitness::Vector{Float64})
    m = minimum(fitness)

    if m < 0
        fitness = 2 * abs(m) .+ fitness
    end

    fitness = 2 * maximum(fitness) .- fitness

    return fitness
end

"""
    getMass(U, max_fitness, epsilon)
    get the mass in a vector `m[i] = f[i] + 2*max_fitness*sum_vio[i] `, tolerance
    in equality constraints are given by `epsilon`
"""
getMass(U::Array{xf_indiv,1}) = fitnessToMass(fvals(U))

function getMass(U::Array{xfgh_indiv,1})
    fitness = fvals(U)
    fitnessToMass(fitness + 2maximum(abs.(fitness))*sum_violations.(U))
end

function center(U, mass)
    d = length(U[1].x)

    c = zeros(Float64, d)

    for i = 1:length(mass)
        c += mass[i] .* U[i].x
    end

    return c / sum(mass)
end

function center(U::Array)
    n = length(U)

    mass = getMass(U)

    return center(U, mass),
    argworst(U), # worst
    argbest(U)  # best
end

function getU(P::Array, K::Int, I::Vector{Int}, i::Int, N::Int)
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    return P[U_ids]
end

function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int)
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    return U_ids
end

function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int, feasible_solutions)
    # at least half of the population is feasible to generate random centers
    if length(feasible_solutions) >= 0.5N || length(feasible_solutions) == 0
        return getU_ids(K, I, i, N)
    end

    # center of mass is generated with at least one feasible solution
    K -= 1
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    push!(U_ids, rand(feasible_solutions))

    return U_ids
end

function crossover(
    x::Vector{Float64},
    y::Vector{Float64},
    p_cr::Vector{Float64},
)
    D = length(x)
    tmp2 = zeros(D)
    for j = 1:D
        if rand() < p_cr[j]
            y[j] = x[j]
            tmp2[j] += 1
        end
    end

    return y, tmp2
end

"""
    ECA_operator(population, K, η_max)

Compute a solution using ECA variation operator, `K` is the number of solutions used to
calculate the center of mass and `η_max` is the maximum stepsize.
"""
function ECA_operator(
        # population::AbstractArray{xf_indiv}, K, η_max;
        population, K, η_max;
        i = rand(1:length(population)),
        U = rand(population, K),
        bounds = nothing
    )

    x = get_position(population[i])

    # generate center of mass
    c, u_worst, u_best = center(U)

    # stepsize
    η = η_max * rand()

    # u: worst element in U
    u = get_position(U[u_worst])

    y = x .+ η .* (c .- u)
    if isnothing(bounds)
        return y
    end

    !isempty(bounds) && evo_boundary_repairer!(y, c, bounds)

    return y
end


text:
using ESDL
using Test

@testset "Axis generation" begin
@test LonAxis(1.0:10.0)==RangeAxis{Float64,:Lon,StepRangeLen{Float64}}(1.0:10.0)
@test LatAxis(1.0:10.0)==RangeAxis{Float64,:Lat,StepRangeLen{Float64}}(1.0:10.0)
end


text:
using Distributed
addprocs(11)

@everywhere using SeisJOLI, GenSPGL, Arpack, MAT, LinearAlgebra, SharedArrays, JLD, Random
@everywhere include("SubFunction/minL2NormSPG.jl")
@everywhere include("SubFunction/minFrobNormspg.jl")
@everywhere include("SubFunction/minnorm.jl")
@everywhere include("SubFunction/RCAMspg.jl")
@everywhere include("SubFunction/l2Norm_project.jl")
@everywhere include("SubFunction/l2Norm_primaldual.jl")

#Set parameters
rank1 = 250;
sigma = 1e-6;
alt = 1;
iter = 1;
weight = 0.75;

obs = load("the observed data in the low-rank domain")
m,n = size(obs);

L = load("load the prior low-rank factor L")

R = load("load the prior low-rank factor R")

#get prior subspace from prior information
U_pre = svd(L);
V_pre = svd(R);

r1 = min(rank(L),rank(R),rank1);
println(r1)
U = U_pre.U[:,1:r1];
V = V_pre.U[:,1:r1];

eta = sigma*norm(obs[:])/(2*m-1);

#interpolate the data with weighted LR reconstruction
time1 = @elapsed L,R = RCAMspg(weight^2*obs, alt, rank1, weight^2*eta, U, V, L, weight, iter);

#obtain the final factors L and R
global L += (1/weight-1) * U * (U' * L);
global R += (1/weight-1) * V * (V' * R);
global L = convert(Array{Complex{Float32},2},L);
global R = convert(Array{Complex{Float32},2},R);


text:
function get_VAR1_rep(V::VAR)
    K = size(V.Σ,1)
    if V.i == true
       v = [V.β[:,1]; zeros(K*(V.p-1),1)]; B = [V.β[:,2:end,]; eye(K*(V.p-1)) zeros(K*(V.p-1),K)]
    else B = [V.β; eye(K*(V.p-1)) zeros(K*(V.p-1),K)]
    end
    #V.i == true ? (return B, v) : (return B)
    return B
end


text:
@reexport module LShapedSolvers

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using StochasticPrograms
using StochasticPrograms: AbstractScenarioProblems, ScenarioProblems, DScenarioProblems
using StochasticPrograms: LQSolver, getsolution, getobjval, getredcosts, getduals, status, SubSolver, get_solver, loadLP, feasibility_problem!
using StochasticPrograms: Execution, Serial, Synchronous, Asynchronous
using StochasticPrograms: PenaltyTerm, Quadratic, Linearized, InfNorm, ManhattanNorm, initialize_penaltyterm!, update_penaltyterm!, solve_penalized!
using MathProgBase
using ProgressMeter
using Clustering

import Base: show, put!, wait, isready, take!, fetch, zero, +, length, size
import StochasticPrograms: StructuredModel, internal_solver, optimize_structured!, fill_solution!, solverstr

const MPB = MathProgBase

export
    LShapedSolver,
    LShaped,
    Crash,
    ncuts,
    niterations,
    StructuredModel,
    add_params!,
    optimsolver,
    hyperoptimal_lshaped,
    optimize_structured!,
    fill_solution!,
    LShaped,
    Serial,
    Synchronous,
    Asynchronous,
    Quadratic,
    Linearized,
    InfNorm,
    ManhattanNorm,
    DontRegularize,
    NoRegularization,
    RegularizedDecomposition,
    WithRegularizedDecomposition,
    RD,
    WithRD,
    TrustRegion,
    WithTrustRegion,
    TR,
    WithTR,
    LevelSet,
    WithLevelSets,
    LV,
    WithLV,
    DontConsolidate,
    Consolidate,
    Consolidation,
    at_tolerance,
    niters,
    tolerance_reached,
    DontAggregate,
    NoAggregation,
    PartialAggregate,
    PartialAggregation,
    Aggregate,
    FullAggregation,
    DynamicAggregate,
    DynamicAggregation,
    ClusterAggregate,
    ClusterAggregation,
    HybridAggregate,
    HybridAggregation,
    SelectUniform,
    SelectDecaying,
    SelectRandom,
    SelectClosest,
    SelectClosestToReference,
    StaticCluster,
    ClusterByReference,
    Kmedoids,
    Hierarchical,
    absolute_distance,
    angular_distance,
    spatioangular_distance

# Include files
include("types/types.jl")
include("consolidators/consolidation.jl")
include("aggregators/aggregation.jl")
include("regularizers/regularization.jl")
include("execution/execution.jl")
include("solver.jl")
include("spinterface.jl")

end # module


text:
# Autogenerated wrapper script for libtree_jll for aarch64-linux-musl
export libtree

JLLWrappers.@generate_wrapper_header("libtree")
JLLWrappers.@declare_executable_product(libtree)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        libtree,
        "bin/libtree",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Finding balance and variability of constraint-based models

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# Here we will use [`flux_balance_analysis`](@ref), [`flux_variability_analysis`](@ref),
# [`parsimonious_flux_balance_analysis`](@ref), and
# [`minimize_metabolic_adjustment_analysis`](@ref), along with the modification functions of
# `COBREXA.jl`, to analyze a toy model of *E. coli*.

# If it is not already present, download the model.

!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")

using COBREXA

#md # !!! tip "Tip: use `?` to get quick help about functions"
#md #       When you are unsure about how a function works, write `?
#md #       function_name` to see the function reference documentation.

model = load_model("e_coli_core.xml")

# ## Optimization solvers in `COBREXA`
#
# To actually perform any optimization based analysis we need to load an
# optimizer. Any [`JuMP.jl`-supported
# optimizers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers)
# will work. Here, we will use [`Tulip.jl`](https://github.com/ds4dm/Tulip.jl)
# to optimize linear programs and
# [`OSQP.jl`](https://osqp.org/docs/get_started/julia.html) to optimize quadratic
# programs.

#md # !!! note "Note: OSQP can be sensitive"
#md #       We recommend reading the docs of `OSQP` before using it, since
#md #       it may give inconsistent results depending on what settings
#md #       you use. Commercial solvers like `Gurobi`, `Mosek`, `CPLEX`, etc.
#md #       require less user engagement.

using Tulip, OSQP, GLPK

# ## Flux balance analysis (FBA)
#
# Most analysis functions come in several variants that produce different types
# of output. All of them usually require a model and `JuMP.jl`-compatible
# optimizer to work in the model.
#
# In the case of FBA, you may choose from these variants (here using the
# `Tulip` optimizer):

vec_soln = flux_balance_analysis_vec(model, Tulip.Optimizer)
#
dict_soln = flux_balance_analysis_dict(model, Tulip.Optimizer)

# ## Extending FBA with modifications

# Often it is desirable to add a slight modification to the problem before
# performing analysis, to see e.g. differences of the model behavior caused by
# the change introduced.
#
# `COBREXA.jl` supports several modifications by default, which include changing objective
# sense, optimizer attributes, flux constraints, optimization objective, reaction and gene
# knockouts, and others. These modifications are applied in the order they are specified. It
# is up to the user to ensure that the changes are sensible.

dict_soln = flux_balance_analysis_dict(
    model,
    OSQP.Optimizer;
    modifications = [ # modifications are applied in order
        ## this changes the objective to maximize the biomass production
        change_objective("R_BIOMASS_Ecoli_core_w_GAM"),

        ## this fixes a specific rate of the glucose exchange
        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),

        ## this knocks out two genes, i.e. constrains their associated reactions to zero.
        knockout(["b0978", "b0734"]), ## the gene IDs are cytochrome oxidase (CYTBD)

        ## ignore the optimizer specified above and change it to Tulip
        change_optimizer(Tulip.Optimizer),

        ## set a custom attribute of the Tulip optimizer (see Tulip docs for more possibilities)
        change_optimizer_attribute("IPM_IterationsLimit", 110),

        ## explicitly tell the optimizer to maximize the new objective
        change_sense(MAX_SENSE),
    ],
)

# This solution can be display using `flux_summary`. Note, this pretty printing only works
# on flux solutions that are represented as dictionaries.
flux_summary(dict_soln)

# ## Flux variability analysis (FVA)

# The default FVA in [`flux_variability_analysis`](@ref) returns maximized and
# minimized reaction fluxes in a matrix. Here we use the dictionary variant in
# flux_variability_analysis_dict, to show how to easily access specific fluxes
# from its results.

fva_mins, fva_maxs = flux_variability_analysis_dict(
    model,
    Tulip.Optimizer;
    bounds = objective_bounds(0.99), # the objective function is allowed to vary by ~1% from the FBA optimum
    modifications = [
        change_optimizer_attribute("IPM_IterationsLimit", 500),
        change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
        change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
    ],
)

#
fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux

# Another option is to display this information using `flux_variability_summary`. This
# pretty printing only works on flux variability analysis results where dictionary keys indicate
# which flux is optimized and the associated value is a flux dictionary.
flux_variability_summary((fva_mins, fva_maxs))

# More sophisticated variants of [`flux_variability_analysis`](@ref) can be used to extract
# specific pieces of information from the solved optimization problems. Here the objective
# value of the minimized flux and the associated biomass growth rate is returned instead
# of every flux.

biomass_idx = first(indexin(["R_BIOMASS_Ecoli_core_w_GAM"], reactions(model))) # index of biomass function
vs = flux_variability_analysis(
    model,
    Tulip.Optimizer;
    bounds = objective_bounds(0.50), # biomass can vary up to 50% less than optimum
    modifications = [
        change_optimizer_attribute("IPM_IterationsLimit", 500),
        change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
        change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
    ],
    ret = m ->
        (COBREXA.JuMP.objective_value(m), COBREXA.JuMP.value(m[:x][biomass_idx])), # m is the model and m[:x] extracts the fluxes from the model
)
#
fva_mins = Dict(rxn => flux for (rxn, flux) in zip(reactions(model), vs[:, 1]))

# ## Parsimonious flux balance analysis (pFBA)

# Parsimonious flux balance analysis (here in
# [`parsimonious_flux_balance_analysis`](@ref) finds a unique flux solution
# that minimizes the squared sum of fluxes of the system subject, while
# maintaining the same objective value as the flux balance analysis solution.
# Since we are optimizing a quadratic objective, we also need to switch to a
# quadratic optimizer. In this case, OSQP will work. We demonstrate it on the
# dictionary-returning variant of pFBA,
# [`parsimonious_flux_balance_analysis_dict`](@ref):

dict_soln = parsimonious_flux_balance_analysis_dict(
    model,
    OSQP.Optimizer;
    modifications = [
        silence, # silence the optimizer (OSQP is very verbose by default)
        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
        change_optimizer_attribute("polish", true),
    ],
)

# The function also has the expectable second variant that returns a vector of
# solutions, in [`parsimonious_flux_balance_analysis_vec`](@ref). Here, we
# utilize it to show how to use different optimizers for finding the optimum
# and for solving the quadratic problem. That may be preferable if the
# optimizer qualities differ for the differing tasks. pFBA allows you to
# specify `qp_modifications` that are applied after the original optimum is
# found, and before the quadratic part of the problem solving begins.

vec_soln = parsimonious_flux_balance_analysis_vec(
    model,
    Tulip.Optimizer; # start with Tulip
    modifications = [
        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
        change_optimizer_attribute("IPM_IterationsLimit", 500), # we may change Tulip-specific attributes here
    ],
    qp_modifications = [
        change_optimizer(OSQP.Optimizer), # now switch to OSQP (Tulip wouldn't be able to finish the computation)
        change_optimizer_attribute("polish", true), # get an accurate solution, see OSQP's documentation
        silence, # and make it quiet.
    ],
)

# ## Minimizing metabolic adjustment analysis (MOMA)

# MOMA is a technique used to find a flux distribution that is closest to some reference
# distribution with respect to the Euclidian norm.

reference_fluxes = parsimonious_flux_balance_analysis_dict( # reference distribution
    model,
    OSQP.Optimizer;
    modifications = [silence, change_optimizer_attribute("polish", true)],
)

moma = minimize_metabolic_adjustment_analysis_dict(
    model,
    reference_fluxes,
    OSQP.Optimizer;
    modifications = [
        silence,
        change_optimizer_attribute("polish", true),
        change_constraint("R_CYTBD"; lb = 0.0, ub = 0.0), # find flux distribution closest to the CYTBD knockout
    ],
)

# ## Composing (more complicated) modifications

# `COBREXA.jl` contains a number of modifications that allow the user to analyze
# non-standard variants of the classic FBA problem. These include thermodynamic
# ([`add_loopless_constraints`](@ref)), capacity ([`add_crowding_constraints`](@ref)), and
# kinetic/capacity ([`add_moment_constraints`](@ref)) modifications. The documentation of
# each modification details what their purpose is. Here, we will demonstrate how these
# modifications can be composed to generate even more interesting analyses.

# Download the json formatted model so that the reaction identifiers correspond to the Escher map.
!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")

model = load_model("e_coli_core.json")

# First find a flux distribution that is thermodyamically loopless and incorporates enzyme
# capacity constraints by composing loopless FBA and FBAwMC.

rid_crowding_weight = Dict(# crowding needs a weight for each flux
    rid => 0.004 for rid in reactions(model) if
    !looks_like_biomass_reaction(rid) && !looks_like_exchange_reaction(rid)
)

loopless_crowding_fluxes = flux_balance_analysis_dict(
    model,
    GLPK.Optimizer;
    modifications = [
        add_crowding_constraints(rid_crowding_weight),
        add_loopless_constraints(),
    ],
)
flux_summary(loopless_crowding_fluxes)

# Next find a flux distribution that satisfies kinetic/capacity constraints using the moment
# algorithm that is closest (using MOMA) to the loopless crowding solution.

ksas = Dict(rid => 1000.0 for rid in reactions(model)) # make up specific activities of the enyzmes
protein_mass_fraction = 0.56

moment_moma = minimize_metabolic_adjustment_analysis_dict(
    model,
    loopless_crowding_fluxes,
    OSQP.Optimizer;
    modifications = [
        silence,
        change_optimizer_attribute("polish", true),
        change_optimizer_attribute("max-iter", 10_000),
        change_constraint("EX_glc__D_e", lb = -1000),
        change_constraint("CYTBD"; lb = 0, ub = 0),
        add_moment_constraints(ksas, protein_mass_fraction;),
    ],
)
flux_summary(moment_moma)

# Finally, plot the results to inspect the flux distributions visually
using CairoMakie, Escher, ColorSchemes

!isfile("e_coli_core_map.json") && download( # download escher map
    "http://bigg.ucsd.edu/escher_map_json/e_coli_core.Core%20metabolism",
    "e_coli_core_map.json",
)

maxflux = maximum(abs.(values(moment_moma)))
minflux = minimum(abs.(values(moment_moma)))

color_interp(x) = begin # Scale color of reaction edges to fluxes (manually binned)
    normed_x = (abs(x) - minflux) / (maxflux - minflux)
    if 0 <= normed_x < 0.01
        ColorSchemes.RdYlBu_4[4]
    elseif 0.01 <= normed_x < 0.25
        ColorSchemes.RdYlBu_4[3]
    elseif 0.25 <= normed_x < 0.5
        ColorSchemes.RdYlBu_4[2]
    else
        ColorSchemes.RdYlBu_4[1]
    end
end
rc = Dict(k => color_interp(v) for (k, v) in moment_moma) # map reaction id to reaction edge color

fig = Figure();
ax = Axis(fig[1, 1]);
escherplot!(ax, "e_coli_core_map.json", reaction_edge_colors = rc)
hidexdecorations!(ax)
hideydecorations!(ax)
fig

#md # !!! tip "Tip: code your own modifications"
#md #       Making custom problem modification functions is really simple due to the
#md #       tight intergration between COBREXA and JuMP. Look at the source code for
#md #       the implemented modifications in `src\analysis\modifications` to get a flavour
#md #       for it.


text:
using FFTW, DSP
function WindowFKSpectrum(dataFreq, datainfo)
	#Input:
	#dataFreq: one shot gather in Frequency domain (nf,nrx,nry)
    #datainfo: a Dictionary contains info for one shot gather
	#Output:
	# one shot gather in Frequency domain after denoising
    nf,nrx,nry = size(dataFreq)
	#fourier transform along Rx
	Data2 = ifft(dataFreq,2);
	Data_FKx = circshift(Data2,[0 ceil(Int, datainfo["nrx"]/2)-1]);

	#fourier transform along Ry
	Data3 = ifft(Data_FKx,3);
	Data_FKxy = circshift(Data3,[0 0 ceil(Int, datainfo["nry"]/2)-1]);

	Rgrid = datainfo["Rgrid"];
	dt = datainfo["dt"];

	Window = zeros(nf,datainfo["nrx"]);
	Center = Int((datainfo["nrx"]+1)/2);
	Data_NFKxy = copy(Data_FKxy);

	for i = 1:nf  #design band pass filter according to the FK spectrum
		if i >= 240
			low = Center-320;
			high = Center+320;
			band = 5;
			Window[i,low+band+1:high-band-1] .= DSP.Windows.tukey(high-band-1-(low+band+1)+1,0.2;padding=0,zerophase=false);
		else
			range1 = round(230/239*i+(90-230/239));
			low = Center-Int(range1);
			high = Center + Int(range1);
			band = 5;
			Window[i,low+band+1:high-band-1] .= DSP.Windows.tukey(high-band-1-(low+band+1)+1,0.2;padding=0,zerophase=false);
		end
	end

	for j = 1:nrx
		Data_NFKxy[:,:,j] = Data_NFKxy[:,:,j] .* Window;
		Data_NFKxy[:,j,:] = Data_NFKxy[:,j,:] .* Window;
	end

	Data3 = circshift(Data_NFKxy,[0 0 -ceil(Int, datainfo["nry"]/2)+1]);
	DataNFKx = fft(Data3,3);

	Data3 = circshift(DataNFKx,[0 -ceil(Int, datainfo["nrx"]/2)+1]);
	DataFreq = fft(Data3,2);

	return DataFreq
end


text:
## Reactive signal update performance: o, h, l, c vectors from TimeArray values-component
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")

s_ohlc = Reactive.Input((ohlc_BA.timestamp[1], vec(ohlc_BA.values[1,:])))
s_close = lift(s -> s[2][4], s_ohlc, typ=Float64)

vsig = [s_close.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in 2:nmax
        push!(s_ohlc, (ohlc_BA.timestamp[i], vec(ohlc_BA.values[i,:])))  ## most time here
        ##tarow = ohlc_BA[i]
        ##push!(s_ohlc, (tarow.timestamp[1], tarow.values))  ## NOT doing this - much longer
        push!(vsig, s_close.value)
    end
end

println("OHLC updates o, h, l, c vector from TimeArray values-component in tuple with its timestamp:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: longer than ohlc_updates_vector_tavalues.jl but less than a factor of two

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
@testset "Finitely generated abelian groups" begin
  include("GrpAb/GrpAbFinGen.jl")
  include("GrpAb/Elem.jl")
  include("GrpAb/SubgroupEnum.jl")
  include("GrpAb/Map.jl")
  include("GrpAb/Lattice.jl")
end


text:
export openbox

"""
    openbox(block::AbstractBlock) -> AbstractBlock

For a black box, like QFTBlock, you can get its white box (faithful simulation) using this function.
"""
openbox(q::QFT{N}) where N = qft_circuit(N)
openbox(q::Daggered{<:QFT, N}) where {N} = qft_circuit(N)'

function openbox(fs::FSimGate)
    if fs.theta ≈ π/2
        return cphase(2,2,1,-fs.phi)*SWAP*rot(kron(Z,Z), -π/2)*put(2,1=>phase(-π/4))
    else
        return cphase(2,2,1,-fs.phi)*rot(SWAP,2*fs.theta)*rot(kron(Z,Z), -fs.theta)*put(2,1=>phase(fs.theta/2))
    end
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
using DataFrames, Plots, MultivariateStats, GLM

# The code is taken from the original
# Bernanke Boivin and Eliasz paper
# FAVAR Bayesian structural factor augmented VAR
# -------------------------------------------------------------------------
# The model is:
#
#   FACTOR (OBSERVATION) EQUATION:
#     _      _           _      _
#    |  X(t)  |         |  F(t)  |
#    |        |  = L  x |        | + e(t)
#    |_ Y(t) _|         |_ Y(t) _|
#
#
#   VAR (STATE) EQUATION:
#     _      _            _        _
#    |  F(t)  |          |  F(t-1)  |
#    |        |  = PHI x |          |  +  u(t)
#    |_ Y(t) _|          |_ Y(t-1) _|
#
#
#  with u(t)~N(0,S_F), and e(t)~N(0,Sigma), with Sigma = diag{(sigma_1)^2 ,
#  ...., (sigma_N)^2}.
# -------------------------------------------------------------------------
#-----------------------------LOAD DATA------------------------------------
# Load quarterly data
# load data used to extract factors
xdata = readtable("/home/lbrugnol/Downloads/FAVAR/xdata.csv", header = true)
# load data on inflation, unemployment and interest rate
ydata = readtable("/home/lbrugnol/Downloads/FAVAR/ydata.csv", header = false)
# load transformation codes (see file transx.m)
# load the slow or fast moving variables codes (see Bernanke, Boivin and
# Eliasz, 2005, QJE)
tcode = readtable("/home/lbrugnol/Downloads/FAVAR/tcode.csv", header=false)
slowcode = readtable("/home/lbrugnol/Downloads/FAVAR/slowcode.csv", header = false)
# load the file with the dates of the data (quarters)
yearlab = readtable("/home/lbrugnol/Downloads/FAVAR/yearlab.csv", header = false)
# load the file with the names of the variables.
namesX = readtable("/home/lbrugnol/Downloads/FAVAR/namesX.csv", header = false)

# namesX.mat contains only the short codes of the data. For a complete
# description see the attached excel file.
stddata = [std(convert(Array,xdata),1) ones(1,3)]

const dRule = Dict{String,Integer}("none" => 1, "first_diff" => 2, "log_yoy"=> 3,
"log" => 4, "log_first_diff" => 5, "yoy" => 6)
const sFile_sd = "/home/lbrugnol/Downloads/FAVAR/xdata.csv"
const bHeader_sd = true
const sFile_tc = "/home/lbrugnol/Downloads/FAVAR/tcode.csv"
const bHeader_tc = false
# const sStart = "1999-01-31"

# Call starting dataset, transformation and create final dataset
sd = StartData(sFile_sd,bHeader_sd)
tc = TransRule(sFile_tc,bHeader_tc,dRule,sd)
fd = FinalData(sd,tc)

# Correct size after stationarity transformation
xdata = convert(Array,fd.dfData[3:end,:])
ydata = convert(Array,ydata[3:end,:])
yearlab = yearlab[3:end,:]

# Demean data (no intercept is estimated as in BBE (2005))
t1 = size(xdata,1)     # time series observations of xdata
t2 = size(ydata,1)     # time series dimension of ydata
stdffr = std(ydata[:,3],1)       # standard deviation of the Fed Funds Rate
xdata = xdata .- mean(xdata,1)
ydata = ydata .- mean(ydata,1)

# Define X and Y matrices
X = xdata    # X contains the 'xdata' which are used to extract factors.
Y = ydata  # Y contains inflation, unemployment and interest rate
# NamesXY has the short codes of X and Y
namesXY = [fd.sd.vNames; "Inflation"; "Unemployment"; "Fed_funds"]

# Number of observations and dimension of X and Y
T,M=size(Y) # T time series observations
N=size(X,2) # N series from which we extract factors
#----------------------------PRELIMINARIES---------------------------------
# Set some Gibbs - related preliminaries
const nrep = 15000  # Number of replications
const nburn = 5000   # Number of burn-in-draws
const nthin = 1   # Consider every thin-th draw (thin value)
const it_print = 100  #Print in the screen every "it_print"-th iteration
const constant = 1  # Set 0:no constant on the FAVAR equation, 1:include constant

# Number of factors & lags:
const K = 2               # Number of Factors
const p = K+M             # p is the dimensionality of [Factors, Y]
const plag = 2            # plag is number of lags in the VAR part
# ==============================| FACTOR EQUATION |========================
# Extract principal components from X (to be used as starting values for the
# MCMC estimation of the factors)
X_st = X./std(X,1)# First standardize data to extract PC
[F0,Lf]=extract(X_st,K)  # F0 are the factors, Lf are the loadings

##### UNDERSTAND WHY DIFFERENT FACTORS WRT MATLAB!!!!!!

function pca(data::Array{Float64}; n::Int64 = 1)
    # Need MultivariateStats
    mX = data.-mean(data,1)
    M = fit(PCA, mX, method = :cov)
    return projection(M)[:,1:n]
  end

F0 = pca(X_st,n = K)

#----If using PC estimates of the factors:
# Now rotate the factor space as in Bernanke, Boivin and Eliasz (2005)
slowindex = find(convert(Array,slowcode).==1)
xslow = X[:,slowindex]
Fslow0 = pca(xslow,n = K)

Ffast0 = Y[:,end]
t1,k1=size(Ffast0[:,:])
ff = ([ones(t1,1) Ffast0 Fslow0]'*[ones(t1,1) Ffast0 Fslow0])
fF = [ones(t1,1) Ffast0 Fslow0]'*F0
b = ff\fF # original paper ols svd
Fr0 = F0 - Ffast0*b[2:k1+1,:]

X=X_st
Y=Y./std(Y)

# #----If using MCMC estimates of the factors:
# # This part is explained in the PhD thesis of Eliasz (2005). Essentially,
# # we rotate the factors in a certain way which gives a loadings matrix L,
# # with the upper KxK block being the identity matrix. This structure of L
# # is the same structure we will have when using MCMC to estimate the
# # factors (because of the identification restrictions imposed on L).
# # regress X on F0 and Y, obtain loadings
# Lfy=olssvd(X(:,K+1:N),[F0 Y])'     # upper KxM block of Ly set to zero
# Lf=[Lf(1:K,:);Lfy(:,1:K)]
# Ly=[zeros(K,M);Lfy(:,K+1:K+M)]
# # transform factors and loadings for LE normalization
# [ql,rl]=qr(Lf')
# Lf=rl  # do not transpose yet, is upper triangular
# F0=F0*ql
# # need identity in the first K columns of Lf, call them A for now
# A=Lf(:,1:K)
# Lf=[eye(K),inv(A)*Lf(:,(K+1):N)]'
# Fr0=F0*A

# Put it all in state-space representation, write obs equ as XY=FY*L+e
XY=[X Y]    # the data X used to extract factors, plus Y (infl,unemp and interest)
FY=[Fr0 Y]  # the extracted factors, plus Y (infl,unemp and interest)

# Obtain L (the loadings matrix)
#----If using PC estimates of the factors:
L = (FY'*FY)\FY'*XY # original paper ols svd decom

# #----If using MCMC estimates of the factors:
# L=[Lf Ly;zeros(M,K),eye(M)]

# Obtain SIGMA (the error varaince in the factor equation)
e = XY - FY*L
SIGMA = e'*e./T
SIGMA = diagm(diag(SIGMA))

# ================================| VAR EQUATION |=========================
# Generate lagged FY matrix.
function mlag2(mX,p)
    T,N = size(mX)
    mXlag = zeros(T,N*p)
    for ii in 1:p
        mXlag[p+1:T,(N*(ii-1)+1):N*ii] = mX[p+1-ii:T-ii,1:N]
    end
    return mXlag
end

FY_lag = mlag2(FY,plag)

FY_lag = FY_lag[plag+1:T,:]
FY_temp = FY[plag+1:T,:]

# Get an initial estimate of PHI (FAVAR autoregressive coefficients) and
# S_F (FAVAR covariance matrix). These values will be used as initial
# values for the Gibbs sampler (but you can use 'zeros' instead, to declare
# 'ignorance' about the initial values
PHI = inv(FY_lag'*FY_lag)*(FY_lag'*FY_temp)
SSE = (FY_temp - FY_lag*PHI)'*(FY_temp - FY_lag*PHI)
S_F = SSE./(T-p)

# If plag is not 1, we need to write the state equation as a VAR(1) model,
# in order to be able to estimate the latent factors using the Kalman
# filter. In this case, PHI and S_F are the parameter matrices of the VAR(p)
# model, while PHI_mat and S_F_mat are the matrices of the VAR(1)
# transformed model (see the time series book of Luthkephol (2005), for
# example). If using PC estimates this is not necessary, but we will need
# the VAR(1) form anyway for the impulse responses, so do not comment this
# part.
S_F_mat = [S_F zeros(p,p*(plag-1));zeros(p*(plag-1),p*plag)]
PHI_mat = [PHI' ; eye(p*(plag-1)) zeros(p*(plag-1),p)]

#-------------------------------- PRIORS ----------------------------------
# #========= INITIAL CONDITION ON LATENT FACTORS:
# # F(0) ~ N(F_0_prmean, F_0_prvar)
# F_0_prmean = zeros(p*plag,1);
# F_0_prvar = 4*eye(p*plag)
#
# # variable indexnM indexes which draws of FY refer only to the factors F
# # (and not Y, since Y is assumed to be observed and known). See also
# # function kfgibbsnv.m
# indexnM = [ones(K,plag);zeros(M,plag)]
# indexnM = find(indexnM==1)

# PRIORS ON FACTOR EQUATION:
# Prior on loadings L_i ~ N(0, I), where i = 1,...,N
Li_prvar = 4*eye(p)
alpha = 0.01
beta = 0.01

# Prior on covariance SIGMA_i ~ iG(a, b), where i = 1,...,N
alpha = 0.01
beta = 0.01

# PRIORS ON VAR EQUATION:

# IMPULSE RESPONSES:
# Note that impulse response and related stuff involves a lot of storage
# and, hence, put istore=0 if you do not want them
istore = 1
if istore == 1
    # Impulse response horizon
    nhor = 21
    shock_init = diagm([zeros(p-1); 1./stdffr]) # in terms of standard deviation, identification is recursive
    imp = zeros(nrep,N+M,nhor)
    bigj = zeros(p,p*plag)
    bigj[1:p,1:p] = eye(p)
end
#-------------------------- END OF PRELIMINARIES --------------------------

# START SAMPLING
println("Number of iterations")

for irep = 1:nrep + nburn    #  GIBBS iterations starts here
    # Print iterations
    # if mod(irep,it_print) == 0
    #     disp(irep); toc
    # end

    #FACTOR (MEASUREMENT) EQUATION:
    # -----------------------------------------------------------------------------------------
    #   STEP I: Sample latent factors F_t using Carter and Kohn (1994) or
    #   use PC estimates
    # -----------------------------------------------------------------------------------------

    # #----If using MCMC estimates of the factors:
    # FY = kfgibbsnv(XY,F_0_prmean,F_0_prvar,L,SIGMA,PHI_mat,S_F_mat,M,indexnM)
    #
    # # Make sure the factors have mean zero
    # FY = FY - repmat(mean(FY),T,1)

    # Commenting the above, means that we are using the PC estimate of FY.

    # -----------------------------------------------------------------------------------------
    #   STEP II: Sample L and SIGMA
    # -----------------------------------------------------------------------------------------
    # Since the covariance matrix of the error (SIGMA) in this equation is
    # diagonal, we can estimate the parameters equation-by-equation
    for i = 1:N
        # Sample L from a Normal distribution. The upper KxK block of L is
        # the identity matrix, so we sample the rest N-K rows only:
        if i > K
            Li_postvar = inv(inv(Li_prvar) + inv(SIGMA[i,i])*FY'*FY)
            Li_postmean = Li_postvar*(inv(SIGMA[i,i])*FY'*X[:,i])
            Lidraw = Li_postmean' + randn(1,p)*chol(Hermitian(Li_postvar))
            L'[i,1:p] = Lidraw
        end

        ed = X[:,i] - FY*L'[i,:]

        # Sample SIGMA(i,i) from iGamma
		S_1 = alpha/2 + T/2
		S_2 = beta/2 + ((X[:,i] - FY*L'[i,:])'*(X[:,i] - FY*L'[i,:]))/2
        Sidraw = inv(rand(Gamma(S_1,1/S_2)))
        SIGMA[i,i] = Sidraw
    end

    #VAR (STATE) EQUATION:
    # Generate lagged FY and then correct sizes (we loose observations
    # because of taking lags, as with the simple VAR model)
    FY_lag = mlag2(FY,plag)
    FY_lag = FY_lag[plag+1:T,:]

    FY_temp = FY[plag+1:T,:]
    # Now we have everything we need to estimate PHI and S_F, as we did
    # with the simple BVAR models. You can use any of the priors we used in
    # the BVAR models, like the Diffuse, the Normal-Wishart, the Independent
    # Normal-Wishart, or the Minnesota (Littermann) prior. Here I use the
    # Diffuse prior, based on the OLS quantities (remember, at each Gibbs
    # iteration we treat the dependent variables (FY_temp) and the lags of
    # the dependent variables (FY_lag) as observed data.

    # -----------------------------------------------------------------------------------------
    #   STEP I: Sample autoregressive coefficients PHI
    # -----------------------------------------------------------------------------------------

    PHI_OLS = inv(FY_lag'*FY_lag)*(FY_lag'*FY_temp)
    phi_OLS = PHI_OLS(:)

    V_post = kron(S_F,inv(FY_lag'*FY_lag))
    phi_vec = phi_OLS + chol(V_post)'*randn(p*p*plag,1)
    PHI = reshape(phi_vec,p*plag,p)'

    # Now create draw of PHI_mat, the parameters in the VAR(1) transformation
    PHI_mat(1:p,:) = PHI

#     # truncate to ensure stationarity
#     while max(abs(eig(PHI_mat)))>0.999
#         phi_vec = phi_OLS + chol(V_post)'*randn(p*p*plag,1)
#         PHI = reshape(phi_vec,p*plag,p)'
#         PHI_mat(1:p,:) = PHI
#     end
    # -----------------------------------------------------------------------------------------
    #   STEP II: Sample covariance matrix S_F
    # -----------------------------------------------------------------------------------------
    SSE = (FY_temp - FY_lag*PHI_OLS)'*(FY_temp - FY_lag*PHI_OLS)
    S_F = inv(wish(inv(SSE),T-p))

    # Now create draw of S_F_mat, the parameter in the VAR(1) transformation
    S_F_mat(1:p,1:p) = S_F

    #--------------------------SAVE AFTER-BURN-IN DRAWS AND IMPULSE RESPONSES -----------------
    if irep > nburn && mod((irep-nburn),nthin)==0

        # Save here draws of the parameters to get means and variances

        #----------------- Impulse response analysis:
        if istore==1
            # Note that Htsd contains the
            # structural error cov matrix
            # Set up things in VAR(1) format as in Lutkepohl page 11
            biga = PHI_mat

            # st dev matrix for structural VAR
            shock = chol(S_F)'
            d = diag(diag(shock))
            shock = inv(d)*shock
            #shock = Hsd2*shock_init
            #shock=eye(p)

            #now get impulse responses for 1 through nhor future periods@
            impresp = zeros(p,p*nhor)
            impresp(1:p,1:p) = shock
            bigai = biga
            for j = 1:nhor-1
                impresp(:,j*p+1:(j+1)*p) = bigj*bigai*bigj'*shock
                bigai = bigai*biga
            end

            imp_m = zeros(p,nhor)
            jj=0
            for ij = 1:nhor
                jj = jj + p
                imp_m(:,ij) = impresp(:,jj)
            end
            imp(irep-nburn,:,:) = L*imp_m
       end #END the impulse response calculation section
    end # END saving after burn-in results
end #END main Gibbs loop (for irep = 1:nrep+nburn)
#=============================GIBBS SAMPLER ENDS HERE==================================

if istore == 1
    scale=repmat(stddata',[1 3 nhor])
    scale=permute(scale,[2 1 3])

    imp = permute(imp,[2 3 1])

    for i=1:N
        if tcode(i)==4
            imp(i,:,:) = exp(imp(i,:,:))-ones(1,nhor,nrep/nthin)
        elseif tcode(i)==5
            imp(i,:,:)=exp(cumsum(imp(i,:,:),2))-ones(1,nhor,nrep/nthin)
        end
    end
    imp = permute(imp,[3 1 2])

    # Set quantiles from the posterior density of the impulse responses
    qus = [.10, .5, .9]
    #if using PC estimates of the factors, use this
    impXY=squeeze(quantile(imp,qus))
#     #if using MCMC estimates of the factors, use this
#     impXY=squeeze(quantile(imp,qus))./scale

    #====================PLOTS
    #---Plot I: impulse responses of inflation, unemployment, interest
    figure
    set(0,'DefaultAxesColorOrder',[0 0 0],...
        'DefaultAxesLineStyleOrder','--|-|--')
    subplot(3,1,1)
    plot(1:nhor,squeeze(impXY(:,end-2,:)))
    hold
    plot(zeros(1,nhor),'-')
    title('Impulse response of inflation')
    xlim([1 nhor])
    set(gca,'XTick',0:3:nhor)
    subplot(3,1,2)
    plot(1:nhor,squeeze(impXY(:,end-1,:)))
    hold
    plot(zeros(1,nhor),'-')
    title('Impulse response of unemployment')
    xlim([1 nhor])
    set(gca,'XTick',0:3:nhor)
    subplot(3,1,3)
    plot(1:nhor,squeeze(impXY(:,end,:)))
    hold
    plot(zeros(1,nhor),'-')
    title('Impulse response of interest rate')
    xlim([1 nhor])
    set(gca,'XTick',0:3:nhor)
    #-------------------------

    #---Plot II: impulse responses of other variables
    figure
    set(0,'DefaultAxesColorOrder',[0 0 0],...
        'DefaultAxesLineStyleOrder','--|-|--')
    # I will plot only 12 out of the 115 variables here
    var_numbers = [2  9 10 28 42 46 77 91 92 108 109 111]
    # These variables have the following short codes:
    var_names = namesXY(var_numbers)

    for i=1:12
        subplot(4,3,i)
        plot(1:nhor,squeeze(impXY(:,var_numbers(i),:)))
        hold
        plot(zeros(1,nhor),'-')
        title(['Impulse response of ' var_names(i)])
        xlim([1 nhor])
        set(gca,'XTick',0:3:nhor)
    end
    #-------------------------

end

clc
toc # Stop timer and print total time


text:
# 25/12/17
# Francisco Romero H.

# Making a video to see how the measurements evolve in time.
# First we plot the error of the method, to see how many successful reconstructions there were.
# Secondly we plot the L1 norm of the video at each time sample, to
# see when there is a constant number of particles, and have an idea of when the particles are reconstructed.
# Thirdly we plot the sum of all the particles in time, so we can see the domain in which they move
# The matrix is flipped so the particles seem to go horizontally instead of vertically.
# Finally we present a screenshot of every time step, where the moving particles can be observed.

# File to include and test
foldername = "../data/2017-11-03T15-27-19-546"

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt
include(foldername * "/ufus_parameters.jl")
errors = np.load(foldername * "/errors.npy")
video = np.load(foldername * "/video.npy")
n_x = 13
n_y = 13
n_im = size(video)[2]

# Plotting the errors.
println(size(errors))
plt.figure()
plt.plot(1:length(errors),errors)
plt.show()

# Plot the L1 norm of the video at each frame, to see how many particles are present at every timestep
norms = zeros(n_im,1);
for i = 1:n_im
    norms[i] = norm(video[:,i],1)
end
plt.figure()
plt.plot(1:n_im, norms)
plt.show()

# Check all the places where the particles have been
plt.figure()
plt.pcolormesh(np.linspace(0, 0.01, n_x), np.linspace(0, 0.01, n_y), reshape(sum(video,2),n_x,n_y))
plt.colorbar()
#plt.savefig("singleframe.pdf")
plt.show()

for i = 1:n_im
plt.figure()
plt.pcolormesh(np.linspace(0, 0.01, n_x), np.linspace(0, 0.01, n_y), reshape(video[:,i],n_x,n_y))
plt.colorbar()
#plt.savefig("singleframe.pdf")
plt.show()
end


text:
# Fire a state transition
function fire(fsm::StateMachine, name::AbstractString, args...)

    # Reserved event names: is, can, finished
    # is: check whether the current state is the same as the requested state
    #     returns true if they are the same, false otherwise
    if name == "is"
        event = args[1]
        fsm.current == event

    # can: ensure a specified event can be triggered from the current state
    #      returns true if it can be triggered, false otherwise
    elseif name == "can"
        event = args[1]
        haskey(fsm.map[event], fsm.current)

    # finished: check if the terminal state has been reached
    #           returns true if current = terminal state, false otherwise
    elseif name == "finished"
        fsm.current == fsm.terminal

    # User-defined events and callbacks
    elseif haskey(fsm.actions, name)
        fsm.actions[name]()

    # Throw an error if there weren't any matching names
    else
        error(name * " not found")
    end
end


text:
using RedPitayaDAQServer
using PyPlot
using ThreadPools

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 64
modulus = 4800
base_frequency = 125000000
periods_per_step = 5
samples_per_period = div(modulus, dec)
periods_per_frame = 50 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency
slow_dac_periods_per_frame = div(50, periods_per_step)

stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
signalTypeDAC(rp, 1 , "SINE")
amplitudeDAC(rp, 1, 1, 0.2)
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "EXTERNAL")

# Sequence
# Global Settings
slowDACStepsPerFrame(rp, slow_dac_periods_per_frame) # This sets PDMClockDivider, but it can also be set directly and with slowDACStepsPerSequence
numSlowDACChan(master(rp), 1)
# Per Sequence settings
lut = collect(range(0,0.7,length=slow_dac_periods_per_frame))
seq = ArbitrarySequence(lut, nothing, slow_dac_periods_per_frame, 2, 0, 0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))

masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)

sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 5)
stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

fig = figure(1)
clf()
plot(vec(uCurrentFrame[:,1,:,:]))
plot(vec(uCurrentFrame[:,2,:,:]))
legend(("Rx1"))
fig

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
module FlowNodeLayersAdditiveCouplingLayerTest

using Test
using ReactiveMP
using ReactiveMP: getf, getflow, getdim, forward, forward!, backward, backward!, jacobian, inv_jacobian
using ReactiveMP: det_jacobian, absdet_jacobian, logdet_jacobian, logabsdet_jacobian

@testset "Additive Coupling Layer" begin

    @testset "Constructor" begin

        # check for standard layer
        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f)
        @test typeof(layer.f)     <: ReactiveMP.PlanarFlowEmpty
        @test typeof(layer)       <: ReactiveMP.AdditiveCouplingLayerPlaceholder
        @test layer.partition_dim == 1
        @test layer.permute       == Val(true)

        f = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false, partition_dim=3)
        @test typeof(layer.f)     <: ReactiveMP.PlanarFlowEmpty
        @test typeof(layer)       <: ReactiveMP.AdditiveCouplingLayerPlaceholder
        @test layer.partition_dim == 3
        @test layer.permute       == Val(false)

    end

    @testset "Prepare-Compile" begin

        layert = AdditiveCouplingLayer(PlanarFlow())
        layerf = AdditiveCouplingLayer(PlanarFlow(); permute = false)
        outt = ReactiveMP._prepare(2, layert)
        outf = ReactiveMP._prepare(2, layerf)

        @test typeof(outt[1]) <: ReactiveMP.AdditiveCouplingLayerEmpty
        @test typeof(outt[2]) <: ReactiveMP.PermutationLayer
        @test typeof(outf)    <: ReactiveMP.AdditiveCouplingLayerEmpty

        @test outt[1].dim     == 2
        @test outt[2].dim     == 2
        @test outf.dim        == 2

        layer_comp  = compile(outf)
        layer_compp = compile(outf, [1.0, 2.0, 3.0])

        @test typeof(layer_comp)        <: AdditiveCouplingLayer
        @test typeof(layer_compp)       <: AdditiveCouplingLayer
        @test typeof(layer_comp.f)      <: Tuple
        @test typeof(layer_comp.f[1])   <: PlanarFlow
        @test typeof(layer_compp.f)     <: Tuple
        @test typeof(layer_compp.f[1])  <: PlanarFlow
        @test layer_compp.f[1].u        == 1.0
        @test layer_compp.f[1].w        == 2.0
        @test layer_compp.f[1].b        == 3.0

        # TODO expend for multiple mappings
        @test nr_params(outf)       == 3
        @test nr_params(layer_comp) == 3

    end

    @testset "Get" begin

        # check get functions for univariate PlanarFlow
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        out = compile(ReactiveMP._prepare(2, layer))
        @test getf(out)             == out.f
        @test getflow(out)          == out.f
        @test getdim(out)           == out.dim
        @test getdim(out)           == 2

    end

    @testset "Base" begin

        # check base functions (univariate)
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        out = compile(ReactiveMP._prepare(2, layer))
        @test eltype(out) == Float64

    end

    @testset "Forward-Backward" begin

        # check forward function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        @test forward(layer, [5.0, 1.5]) == [5.0, 7.4999983369439445]
        @test forward(layer, [4.0, 2.5]) == [4.0, 7.499909204262595]
        @test forward.(layer, [[5.0, 1.5], [4.0, 2.5]]) == [[5.0, 7.4999983369439445], [4.0, 7.499909204262595]]

        # check forward! function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        output = zeros(2)
        forward!(output, layer, [5.0, 1.5])
        @test output == [5.0, 7.4999983369439445]
        forward!(output, layer, [4.0, 2.5])
        @test output == [4.0, 7.499909204262595]

        # check forward function (input > 2)
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(3, layer))
        x = randn(3)
        @test backward(layer, forward(layer, x)) ≈ x

        # check backward function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        @test backward(layer, [5.0, 7.4999983369439445]) == [5.0, 1.5]
        @test backward(layer, [4.0, 7.499909204262595]) == [4.0, 2.5]
        @test backward.(layer, [[5.0, 7.4999983369439445], [4.0, 7.499909204262595]]) == [[5.0, 1.5], [4.0, 2.5]]

        # check backward! function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        output = zeros(2)
        backward!(output, layer, [5.0, 7.4999983369439445])
        @test output == [5.0, 1.5]
        backward!(output, layer, [4.0, 7.499909204262595])
        @test output == [4.0, 2.5]

    end

    @testset "Jacobian" begin

        # check jacobian function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        @test jacobian(layer, [3.0, 1.5]) == [1.0 0.0; 1.0197320743308804 1.0]
        @test jacobian(layer, [2.5, 5.0]) == [1.0 0.0; 1.1413016497063289 1.0]
        @test jacobian.(layer, [[3.0, 1.5], [2.5, 5.0]]) == [[1.0 0.0; 1.0197320743308804 1.0], [1.0 0.0; 1.1413016497063289 1.0]]

        # check jacobian function
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        @test inv_jacobian(layer, [3.0, 1.5]) == [1.0 0.0; -1.0197320743308804 1.0]
        @test inv_jacobian(layer, [2.5, 5.0]) == [1.0 0.0; -1.1413016497063289 1.0]
        @test inv_jacobian.(layer, [[3.0, 1.5], [2.5, 5.0]]) == [[1.0 0.0; -1.0197320743308804 1.0], [1.0 0.0; -1.1413016497063289 1.0]]

        # check for invertibility
        layer = AdditiveCouplingLayer(PlanarFlow(); permute=false)
        x = randn(10)
        layer = compile(ReactiveMP._prepare(10, layer))
        @test inv(jacobian(layer, x)) ≈ inv_jacobian(layer, forward(layer, x))

    end

    @testset "Utility Jacobian" begin

        # check utility functions jacobian (univariate)
        params = [1.0, 2.0, -3.0]
        f  = PlanarFlow()
        layer = AdditiveCouplingLayer(f; permute=false)
        layer = compile(ReactiveMP._prepare(2, layer), params)
        @test det_jacobian(layer, [1.5, 6.9]) == 1.0
        @test det_jacobian(layer, [2.5, 6.4]) == 1.0
        @test absdet_jacobian(layer, [1.5, 6.9]) == 1.0
        @test absdet_jacobian(layer, [2.5, 6.4]) == 1.0
        @test logdet_jacobian(layer, [1.5, 6.9]) == 0.0
        @test logdet_jacobian(layer, [2.5, 6.4]) == 0.0
        @test logabsdet_jacobian(layer, [1.5, 6.9]) == 0.0
        @test logabsdet_jacobian(layer, [2.5, 6.4]) == 0.0

    end

end

end


text:
using Documenter, GPUShowcases

makedocs(
    modules = [GPUShowcases],
    format = :html,
    sitename = "GPU Showcases for Julia",
    pages = [
        "Home" => "index.md",
        "Showcases" => [
            "SmokeSimulation/smokesimulation.md",
            "PDE/pde.md",
            "Poincare/poincare.md",
            "Convolution/convolution.md",
            "ML/ml.md"
        ],
    ]
)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
    repo   = "github.com/JuliaGPU/GPUShowcases.jl.git",
    julia  = "0.6",
    target = "build",
    osname = "linux",
    make = nothing
)


text:
# (1) The gold-standard solution:
function st_solve(k::Int, state::Symbol=:okay)
    r = Tuple{Symbol,Int}[]
    while k > 0
        push!(r, (state, k))
        if state == :okay
            k -= 3
            state = :ping
        elseif state == :ping
            k ÷= 2
            state = :pong
        elseif state == :pong
            k -= 4
            state = :okay
        end
    end
    return (r, k)
end

# (2) How to set up the data:
struct ProblemData
    value::Air.Volatile{Int}
    state::Air.Volatile{Symbol}
    actor::Air.Actor{Air.PVector{Tuple{Symbol,Int}}}
    function ProblemData(k::Int, state::Symbol=:okay)
        return new(Air.Volatile{Int}(k), Air.Volatile(state),
                   Air.Actor(PVector{Tuple{Symbol,Int}}()))
    end
end

# (3) Task functions:
function mt_task(pd::ProblemData, onstate, nextstate, op)
    go = true
    while go
        go = tx() do
            (pd.value[] <= 0) && return false
            yield()
            if pd.state[] == onstate
                let val = pd.value[]
                    Air.send(a -> Air.push(a, (onstate, val)), pd.actor)
                end
                pd.value[] = op(pd.value[])
                pd.state[] = nextstate
            end
            return true
        end
    end
end

# (4) How to solve the problem.
function mt_solve(k::Int, state::Symbol=:okay)
    pd = ProblemData(k, state)
    t1 = Threads.@spawn mt_task($pd, :pong, :okay, x -> x - 4)
    t2 = Threads.@spawn mt_task($pd, :ping, :pong, x -> x ÷ 2)
    t3 = Threads.@spawn mt_task($pd, :okay, :ping, x -> x - 3)
    wait(t1)
    wait(t2)
    wait(t3)
    sleep(0.1)
    return (pd.actor[], pd.value[])
end

function countdown_test(k::Int)
    (st_steps, st_val) = st_solve(k)
    (mt_steps, mt_val) = mt_solve(k)
    @test mt_val == st_val
    @test mt_steps == st_steps
end

@testset "CountDown" begin
    countdown_test(5)
    countdown_test(7)
    countdown_test(15)
    countdown_test(21)
    countdown_test(23)
    countdown_test(351)
    countdown_test(5016)
    countdown_test(18622)
    countdown_test(195122)
end


text:
using BinDeps
using Compat
using Conda
@BinDeps.setup

libsymengine = library_dependency("libsymengine", aliases=["libsymengine", "symengine"])

Conda.add_channel("conda-forge")
Conda.add_channel("symengine")
if is_windows()
    # SymEngine can only be installed to a python 3.5 environment
    Conda.add("python=3.5")
end
provides(Conda.Manager, "symengine==0.2.0", [libsymengine])

@BinDeps.install Dict([(:libsymengine, :libsymengine)])


text:
########################################################################
# Atoms of the 2-norm ball.
########################################################################

"""
    TwoBallAtom(n, z)

Atom in the TwoBall of dimension n.
"""
struct TwoBallAtom <: AbstractAtom
    n::Int64
    z::Vector{Float64}
    function TwoBallAtom(z)
        n = length(z)
        return new(n, z)
    end
end

"""
Multiply a two-norm atom by a linear map.
"""
Base.:(*)(M::AbstractLinearOp, a::TwoBallAtom) = M*a.z
Base.vec(a::TwoBallAtom) = a.z
Base.length(a::TwoBallAtom) = a.n

########################################################################
# Atomic set for the 1-norm.
########################################################################

"""
    TwoBall(n, maxrank=n)

Atomic set defined by 2-norm ball in `n` variables.
The atomic set takes two optional parameters:

`maxrank` is the maximum dimension of the face exposed by a vector.
"""
struct TwoBall <: AbstractAtomicSet
    n::Int64
    maxrank::Int64
    function TwoBall(n, maxrank)
        n ≥ 1 || throw(DomainError(n,"n must be a positive integer"))
        n ≥ maxrank ≥ 1 || throw(DomainError(maxrank,"maxrank must be ≥ 1"))
        return new(n, maxrank)
    end
end
TwoBall(n; maxrank=n) = TwoBall(n, maxrank)

"""
    gauge(A::TwoBall, x::Vector)

Gives the 2-norm of `x`.
"""
gauge(A::TwoBall, x::Vector) = norm(x)
gauge(A::TwoBall, x::Matrix) = gauge(A, vec(x))

"""
    support(A::TwoBall, z::Vector)

Gives the 2-norm of `z`.
"""
support(A::TwoBall, z::Vector) = norm(z)

"""
    expose(A::TwoBall, z::Vector)

A non-overwriting version of [`expose!`](@ref).
"""
expose(A::TwoBall, z; kwargs...) = expose!(A::TwoBall, copy(z); kwargs...)

"""
    expose!(A::TwoBall, z::Vector; tol=1e-12)

Obtain an atom in the face exposed by the vector `z`.
The vector `z` is overwritten. If `norm(z)<tol`, then
`z` is returned untouched.
"""
function expose!(A::TwoBall, z::Vector; tol=1e-12)
    znorm = norm(z)
    if znorm < tol
        return TwoBallAtom(zero(z))
    end
    return TwoBallAtom(z ./= znorm)
end

Base.length(A::TwoBall) = A.n
rank(A::TwoBall) = A.maxrank
atom_name(A::TwoBall) = "2-norm ball"
atom_description(A::TwoBall) = "{ x ∈ ℝⁿ | ||x||_2 ≤ 1 }"
atom_parameters(A::TwoBall) = "n = $(length(A)); maxrank = $(A.maxrank)"

########################################################################
# Face of the 1-norm ball.
########################################################################

struct TwoBallFace <: AbstractFace
    n::Int64
    k::Int64
    S::LinearMap{Float64}
end

TwoBallFace(S::LinearMap) = TwoBallFace(size(S)...,S)

"""
    face(A, z)

Return a face of the atomic set `A` exposed by the vector `z`. The dimension of
the exposed face is limited by `k=maxrank(A)`. If the dimension is being limited
by `k`, then the routine returns a set of rank `k` atoms that define the subset
of the exposed face.
"""
function face(A::TwoBall, z::Vector; rTol=1e-12)
    @assert length(A) == length(z)
    a = expose(A, z, tol=rTol)
    S = LinearMap(reshape(a.z, length(A), 1))
    return TwoBallFace(S)
end

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::TwoBallFace, c::Vector) = F.S*c
Base.:(*)(M::AbstractLinearOp, F::TwoBallFace) = M*F.S
Base.:(*)(λ::Real, F::TwoBallFace) = λ*F.S
Base.length(F::TwoBallFace) = F.n
rank(F::TwoBallFace) = F.k
vec(F::TwoBallFace) = F.S
face_name(F::TwoBallFace) = "2-norm ball"
face_parameters(F::TwoBallFace) = "rank = $(rank(F)); n = $(length(F))"


text:
n = 6 # rand(2:10)

kappa = 5 # rand(2:6)
mode =  2 # rand(1:5)

A = matrixdepot("randsvd", n, kappa, mode)
B = matrixdepot("randsvd", n, kappa)

if mode == 5
    @test cond(A) <= kappa
else
    @test cond(A) ≈ kappa
end

A5 = matrixdepot("randsvd", n, kappa, 5)
A4 = matrixdepot("randsvd", n, kappa, 4)
A3 = matrixdepot("randsvd", n, kappa, 3)
A2 = matrixdepot("randsvd", n, kappa, 2)
A1 = matrixdepot("randsvd", n, kappa, 1)
θ = matrixdepot("randsvd", 1)

@test_throws ArgumentError matrixdepot("randsvd", 5, kappa, 6)

println("'randsvd' passed test...")


text:
# Declarations mirroring the C declarations in xsum.h

struct xsum_small_accumulator
    chunk::NTuple{67,Int64}
    inf::Int64
    nan::Int64
    adds_until_propagate::Cint
end

struct xsum_large_accumulator
    chunk::NTuple{4096,UInt64}
    count::NTuple{4096,Int16}
    chunks_used::NTuple{64,UInt64}
    used_used::UInt64
    sacc::xsum_small_accumulator
end

# The following low-level functions do little checking of their arguments, e.g. they don't check
# that the array strides == 1 or that the dot product is of vectors with equal lengths.
for kind in ("xsum_small_", "xsum_large_")
    T = Symbol(kind, "accumulator")
    @eval begin
        $(Symbol(kind,"init"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"init"))),libxsum), Cvoid, (Ref{$T},), acc)
        $(Symbol(kind,"addv"))(acc, a::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"addv"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Int), acc, a, length(a))
        $(Symbol(kind,"add_sqnorm"))(acc, a::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"add_sqnorm"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Int), acc, a, length(a))
        $(Symbol(kind,"add_dot"))(acc, a::StridedVector{Float64}, b::StridedVector{Float64}) =
            ccall(($(QuoteNode(Symbol(kind,"add_dot"))),libxsum), Cvoid, (Ref{$T},Ptr{Float64},Ptr{Float64},Int), acc, a, b, length(a))
        $(Symbol(kind,"round"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"round"))),libxsum), Cdouble, (Ref{$T},), acc)
        $(Symbol(kind,"chunks_used"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"chunks_used"))),libxsum), Cint, (Ref{$T},), acc)
        $(Symbol(kind,"add_accumulator"))(accdest, accsrc) =
            ccall(($(QuoteNode(Symbol(kind,"add_accumulator"))),libxsum), Cvoid, (Ref{$T},Ref{$T}), accdest, accsrc)
        $(Symbol(kind,"add1"))(acc, x::Real) =
            ccall(($(QuoteNode(Symbol(kind,"add1"))),libxsum), Cvoid, (Ref{$T}, Cdouble), acc, x)
        $(Symbol(kind,"negate"))(acc) =
            ccall(($(QuoteNode(Symbol(kind,"negate"))),libxsum), Cvoid, (Ref{$T},), acc)
    end
end

xsum_large_to_small_accumulator(acc_small, acc_large) =
    ccall((:xsum_large_to_small_accumulator,libxsum), Cvoid, (Ref{xsum_small_accumulator}, Ref{xsum_large_accumulator}), acc_small, acc_large)


text:
# module adapters_postgresql

module Options

for_postgresql = (
    host = "localhost",
    dbname = "postgresqltest",
    user = "postgres",
    password = get(ENV, "PGPASSWORD", ""),
)

end # module Options

# module adapters_postgresql


text:
# refs:
#     "Definition of the Arithmetic Operations and Comparison Relations
#        for an Interval Arithmetic Standard" by Gerd Bohlender and Ulrich Kulisch
#     Reliable Computing 15, 2011 pp 41-42
#
#     Computer Arithmetic and Validity (2nd edition) by Ulrich Kulisch (section 7.4, pp242-43)

(==){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = ((a.lo == b.lo) & (a.hi == b.hi))
(!=){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = ((a.lo != b.lo) | (a.hi != b.hi))
(<=){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = ((a.lo <= b.lo) & (a.hi <= b.hi))
(>=){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = ((a.lo >= b.lo) & (a.hi >= b.hi))
(< ){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = (!(a >= b))
(> ){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = (!(a <= b))
(isless){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = (a < b)
(isequal){S<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,Q,C}) = (a == b)

#=
    if a,b are of differing Sculpt and one boundry value from each of a and b are equal
    then (a <= b or a >= b) and ( !(a<b) and !(a>b) and !(a==b) )

     ⟨1.0, 2.0⟫ <=   ⟨1.0, 2.0⟩ <= ⟪1.0, 2.0 ⟩ ,  ⟪1.0, 2.0⟫
=#
(==){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = false
(!=){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = true
(<=){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = opened(a) < opened(b)
(>=){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = opened(a) > opened(b)
(< ){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = closed(a) < closed(b)
(> ){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = closed(a) > closed(b)
(isequal){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = false
(isless){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C}) = (a < b)

# see foryouruse.jl to access comparison logic when qualia differ


text:
module Painleve

using LinearAlgebra, RiemannHilbert, ApproxFun

export painleve2_6ray, pl2def_no_s2_pos_x, pl2def_no_s2_neg_x1, pl2def_no_s2_neg_x2, pl2def_no_s2_neg_x3

function painleve2_6ray((s1,s2,s3),x; n=600)
    @assert mod(n,6) == 0
	@assert abs(s1 - s2 + s3 + s1*s2*s3) ≤ 100eps()

	Γ1 = Segment(0, 2.5exp(im*π/6))
    Γ2 = Segment(0, 2.5exp(im*π/2))
    Γ3 = Segment(0, 2.5exp(5im*π/6))
    Γ4 = Segment(0, 2.5exp(-5im*π/6))
    Γ5 = Segment(0, 2.5exp(-im*π/2))
    Γ6 = Segment(0, 2.5exp(-im*π/6))
    Γ = Γ1 ∪ Γ2 ∪ Γ3 ∪ Γ4 ∪ Γ5 ∪ Γ6

    Θ(z) = 8/3*z^3+2*x*z

    S1(z) = [1 0; s1*exp(im*Θ(z)) 1]
    S2(z) = [1 s2*exp(-im*Θ(z)); 0 1]
    S3(z) = [1 0; s3*exp(im*Θ(z)) 1]
    S4(z) = [1 -s1*exp(-im*Θ(z)); 0 1]
    S5(z) = [1 0; -s2*exp(im*Θ(z)) 1]
    S6(z) = [1 -s3*exp(-im*Θ(z)); 0 1]

    G = Fun( z ->   if z in component(Γ, 1)   S1(z)
                elseif z in component(Γ, 2)   S2(z)
                elseif z in component(Γ, 3)   S3(z)
                elseif z in component(Γ, 4)   S4(z)
                elseif z in component(Γ, 5)   S5(z)
                elseif z in component(Γ, 6)   S6(z)
                end, Γ);

    Φ = transpose(rhsolve(transpose(G), n));
    z = Fun(ℂ)
    2(z*Φ[1,2])(Inf)
end

function pl2def_no_s2_pos_x((s1,s2,s3),x; n=400)
    @assert mod(n,4) == 0
    @assert abs(s1 - s2 + s3 + s1*s2*s3) ≤ 100eps()
    @assert x > 0
    @assert s2 == 0

    z_0 = (im*sqrt(x))/2

    Γ1 = Segment(z_0, z_0 + 2.5)
    Γ3 = Segment(z_0, z_0 - 2.5)
    Γ4 = Segment(-z_0, -z_0 - 2.5)
    Γ6 = Segment(-z_0, -z_0 + 2.5)
    Γ = Γ1 ∪ Γ3 ∪ Γ4 ∪ Γ6

    Θ(z) = 8/3*z^3+2*x*z

    S1(z) = [1 0; s1*exp(im*Θ(z)) 1]
    S3(z) = [1 0; s3*exp(im*Θ(z)) 1]
    S4(z) = [1 -s1*exp(-im*Θ(z)); 0 1]
    S6(z) = [1 -s3*exp(-im*Θ(z)); 0 1]

    G = Fun( z ->   if z in component(Γ, 1)    S1(z)
                elseif z in component(Γ, 2)    S3(z)
                elseif z in component(Γ, 3)    S4(z)
                elseif z in component(Γ, 4)    S6(z)
                end, Γ);

    Φ = transpose(rhsolve(transpose(G), n));
    z = Fun(ℂ)
    2(z*Φ[1,2])(Inf)
end

function pl2def_no_s2_neg_x1((s1,s2,s3),x; n=500, l=0.7)
    @assert mod(n,5) == 0
    @assert abs(s1 - s2 + s3 + s1*s2*s3) ≤ 100eps()
    @assert s2 == 0
    @assert x < 0
    @assert l*cos(π/4) < 1/2

    Γ6Γ1 = Segment(-0.5, 0.5)
    Γ1   = Segment(0.5, 0.5 + l*exp(im*π/4))
    Γ3   = Segment(-0.5, -0.5 + l*exp(im*3π/4))
    Γ4   = Segment(-0.5, -0.5 + l*exp(-im*3π/4))
    Γ6   = Segment(0.5, 0.5 + l*exp(-im*π/4))
    Γ = Γ6Γ1 ∪ Γ1 ∪ Γ3 ∪ Γ4 ∪ Γ6

    Θ(z) = (2*z/3)*(sqrt(-x)^3)*(4*z^2-3)

    S1(z)   = [1 0; s1*exp(im*Θ(z)) 1]
    S3(z)   = [1 0; s3*exp(im*Θ(z)) 1]
    S4(z)   = [1 -s1*exp(-im*Θ(z)); 0 1]
    S6(z)   = [1 -s3*exp(-im*Θ(z)); 0 1]

    G = Fun( z -> if z in component(Γ, 1)    S6(z)S1(z)
              elseif z in component(Γ, 2)    S1(z)
              elseif z in component(Γ, 3)    S3(z)
              elseif z in component(Γ, 4)    S4(z)
              elseif z in component(Γ, 5)    S6(z)
              end, Γ);

    Φ = transpose(rhsolve(transpose(G), n));
    z = Fun(ℂ)
    2(z*Φ[1,2])(Inf)*sqrt(-x)
end

function pl2def_no_s2_neg_x2((s1,s2,s3),x; n=450, l=0.7)
    @assert mod(n,9) == 0
    @assert abs(s1 - s2 + s3 + s1*s2*s3) ≤ 100eps()
    @assert s2 == 0
    @assert s1*s3 != 1
    @assert x < 0
    @assert l*cos(π/4) < 1/2

    ΓD  = Segment(-0.5, 0.5)
    Γ1  = Segment(0.5, 0.5 + l*exp(im*π/4))
    ΓUi = Segment(0.5, 0.5 + l*exp(im*3π/4))
    ΓU  = Segment(-0.5, -0.5 + l*exp(im*π/4))
    Γ3  = Segment(-0.5, -0.5 + l*exp(im*3π/4))
    Γ4  = Segment(-0.5, -0.5 + l*exp(-im*3π/4))
    ΓL  = Segment(-0.5, -0.5 + l*exp(-im*π/4))
    ΓLi = Segment(0.5, 0.5 + l*exp(-im*3π/4))
    Γ6  = Segment(0.5, 0.5 + l*exp(-im*π/4))
    Γ = ΓD ∪ Γ1 ∪ ΓUi ∪ ΓU ∪ Γ3 ∪ Γ4 ∪ ΓL ∪ ΓLi ∪ Γ6

    Θ(z) = (2*z/3)*(sqrt(-x)^3)*(4*z^2-3)

    D(z)     = [1-s1*s3 0; 0 1/(1-s1*s3)]
    S1(z)    = [1 0; s1*exp(im*Θ(z)) 1]
    Ui(z)    = [1 s3*exp(-im*Θ(z))/(1-s1*s3); 0 1]
    U(z)     = [1 -s3*exp(-im*Θ(z))/(1-s1*s3); 0 1]
    S3(z)    = [1 0; s3*exp(im*Θ(z)) 1]
    S4(z)    = [1 -s1*exp(-im*Θ(z)); 0 1]
    L(z)     = [1 0; s1*exp(im*Θ(z))/(1-s1*s3) 1]
    Li(z)    = [1 0; -s1*exp(im*Θ(z))/(1-s1*s3) 1]
    S6(z)    = [1 -s3*exp(-im*Θ(z)); 0 1]

    G = Fun( z -> if z in component(Γ, 1)    D(z)
              elseif z in component(Γ, 2)    S1(z)
              elseif z in component(Γ, 3)    Ui(z)
              elseif z in component(Γ, 4)    U(z)
              elseif z in component(Γ, 5)    S3(z)
              elseif z in component(Γ, 6)    S4(z)
              elseif z in component(Γ, 7)    L(z)
              elseif z in component(Γ, 8)    Li(z)
              elseif z in component(Γ, 9)    S6(z)
              end, Γ);

    Φ = transpose(rhsolve(transpose(G), n))*sqrt(-x);
    z = Fun(ℂ)
    2(z*Φ[1,2])(Inf)
end

function pl2def_no_s2_neg_x3((s1,s2,s3),x; n=1600, l=0.6, R=0.1)
    @assert mod(n,16) == 0
    @assert abs(s1 - s2 + s3 + s1*s2*s3) ≤ 100eps()
    @assert s2 == 0
    @assert s1*s3 != 1
    @assert x < 0
    @assert (l+R)*cos(π/4) < 0.5

    Γ1  = Segment(0.5 + R*exp(im*π/4), 0.5 + R*exp(im*π/4) + l*exp(im*π/4))
    ΓUi = Segment(0.5 + R*exp(im*3π/4), 0.5 + R*exp(im*3π/4) + l*exp(im*3π/4))
    ΓU  = Segment(-0.5 + R*exp(im*π/4), -0.5 + R*exp(im*π/4) + l*exp(im*π/4))
    Γ3  = Segment(-0.5 + R*exp(im*3π/4), -0.5 + R*exp(im*3π/4) + l*exp(im*3π/4))
    Γ4  = Segment(-0.5 + R*exp(-im*3π/4), -0.5 + R*exp(-im*3π/4) + l*exp(-im*3π/4))
    ΓL  = Segment(-0.5 + R*exp(-im*π/4), -0.5 + R*exp(-im*π/4) + l*exp(-im*π/4))
    ΓLi = Segment(0.5 + R*exp(-im*3π/4), 0.5 + R*exp(-im*3π/4) + l*exp(-im*3π/4))
    Γ6  = Segment(0.5 + R*exp(-im*π/4), 0.5 + R*exp(-im*π/4) + l*exp(-im*π/4))

    C11 = Segment(0.5 + R*exp(im*π/4), 0.5 + R*exp(-im*π/4))
    C12 = Segment(0.5 + R*exp(-im*π/4), 0.5 + R*exp(-im*3π/4))
    C13 = Segment(0.5 + R*exp(-im*3π/4),0.5 + R*exp(im*3π/4))
    C14 = Segment(0.5 + R*exp(3im*π/4), 0.5 + R*exp(im*π/4))

    C21 = Segment(-0.5 + R*exp(im*-3π/4), -0.5 + R*exp(im*-5π/4))
    C22 = Segment(-0.5 + R*exp(im*3π/4), -0.5 + R*exp(im*π/4))
    C23 = Segment(-0.5 + R*exp(im*π/4), -0.5 + R*exp(-im*π/4))
    C24 = Segment(-0.5 + R*exp(-im*π/4), -0.5 + R*exp(-im*3π/4))

    Γ = Γ1 ∪ ΓUi ∪ ΓU ∪ Γ3 ∪ Γ4 ∪ ΓL ∪ ΓLi ∪ Γ6 ∪
        C11 ∪ C12 ∪ C13 ∪ C15 ∪ C21 ∪ C22 ∪ C23 ∪ C25

    Θ(z) = (2*z/3)*(sqrt(-x)^3)*(4*z^2-3)
    f(z) = ((1+2z)/(2z-1))^(im/(2*π))

    D(z)     = [1-s1*s3 0; 0 1/(1-s1*s3)]
    S1(z)    = [1 0; s1*exp(im*Θ(z)) 1]
    Ui(z)    = [1 s3*exp(-im*Θ(z))/(1-s1*s3); 0 1]
    U(z)     = [1 -s3*exp(-im*Θ(z))/(1-s1*s3); 0 1]
    S3(z)    = [1 0; s3*exp(im*Θ(z)) 1]
    S3i(z)   = [1 0; -s3*exp(im*Θ(z)) 1]
    S4(z)    = [1 -s1*exp(-im*Θ(z)); 0 1]
    S4i(z)   = [1 s1*exp(-im*Θ(z)); 0 1]
    L(z)     = [1 0; s1*exp(im*Θ(z))/(1-s1*s3) 1]
    Li(z)    = [1 0; -s1*exp(im*Θ(z))/(1-s1*s3) 1]
    S6(z)    = [1 -s3*exp(-im*Θ(z)); 0 1]

    P(z)     = [f(z)^(log(1-s1*s3)) 0; 0 f(z)^(log(1/(1-s1*s3)))]
    Pi(z)    = [f(z)^(log(1/(1-s1*s3))) 0; 0 f(z)^(log(1-s1*s3))]

    G = Fun( z -> if z in component(Γ, 1)     P(z)S1(z)Pi(z)
              elseif z in component(Γ, 2)     P(z)Ui(z)Pi(z)
              elseif z in component(Γ, 3)     P(z)U(z)Pi(z)
              elseif z in component(Γ, 4)     P(z)S3(z)Pi(z)
              elseif z in component(Γ, 5)     P(z)S4(z)Pi(z)
              elseif z in component(Γ, 6)     P(z)L(z)Pi(z)
              elseif z in component(Γ, 7)     P(z)Li(z)Pi(z)
              elseif z in component(Γ, 8)     P(z)S6(z)Pi(z)

              elseif z in component(Γ, 9)     Li(z)S6(z)Pi(z)
              elseif z in component(Γ, 10)    Li(z)Pi(z)
              elseif z in component(Γ, 11) && imag(z) ≤ 0    Pi(z)
              elseif z in component(Γ, 11) && imag(z) > 0    D(z)Pi(z)
              elseif z in component(Γ, 12)    Li(z)S6(z)S1(z)Pi(z)

              elseif z in component(Γ, 13)    Li(z)S4i(z)Pi(z)
              elseif z in component(Γ, 14)    Li(z)S4i(z)S3i(z)Pi(z)
              elseif z in component(Γ, 15) && imag(z) > 0    D(z)Pi(z)
              elseif z in component(Γ, 15) && imag(z) ≤ 0    Pi(z)
              elseif z in component(Γ, 16)    Li(z)Pi(z)
              end, Γ);

    Φ = transpose(rhsolve(transpose(G), n))
    z = Fun(ℂ)
    2(z*Φ[1,2])(Inf)*sqrt(-x)
end

end # module


text:
import ChainRulesCore.Tangent

mutable struct BatchNorm2d
    gamma::AbstractArray
    beta::AbstractArray
    eps::Real
    momentum::Real
    training::Bool
    running_mean::AbstractArray
    running_var::AbstractArray
    count::Int
end

BatchNorm2d(num_features::Int; eps=1e-5, momentum=0.1, training=true) =
    BatchNorm2d(
        ones(1, 1, num_features, 1),
        zeros(1, 1, num_features, 1),
        eps,
        momentum,
        training,
        zeros(1, 1, num_features, 1),
        zeros(1, 1, num_features, 1),
        0
    )

Base.show(io::IO, m::BatchNorm2d) = print(io, "BatchNorm2d($(length(m.beta)))")

trainmode!(m::BatchNorm2d, train::Bool=true) = (m.training = train; nothing)

function batchnorm_impl(gamma::AbstractArray, beta::AbstractArray, x::AbstractArray,
                        mu::AbstractArray, sigma2::AbstractArray, momentum::Real; eps, training)
    x_hat = (x .- mu) ./ sqrt.(sigma2 .+ eps)
    return gamma .* x_hat .+ beta
end

function ∇batchnorm_impl(gamma::AbstractArray, beta::AbstractArray, x::AbstractArray, dy::AbstractArray,
                         mu::AbstractArray, sigma2::AbstractArray, momentum; eps, training)
    # based on:
    # https://kratzert.github.io/2016/02/12/understanding-the-gradient-flow-through-the-batch-normalization-layer.html
    sum_dims = tuple([i for i=1:ndims(x) if size(x, i) != size(mu, i)]...)   # all dims but channel
    P = prod(size(x, i) for i in sum_dims)   # product of all dims but channel
    typ = eltype(gamma)
    _1 = one(typ)        # 1 of eltype(x)
    _05 = _1 / 2         # 0.5 of eltype(x)
    # part of forward pass - we can optimize it later and take from cache
    xmu = x .- mu
    sqrtvar = sqrt.(sigma2 .+ eps)
    ivar = one(typ) ./ sqrtvar
    xhat = xmu .* ivar
    # reverse pass
    dbeta = Yota.unbroadcast(beta, dy)
    dgamma = Yota.unbroadcast_prod_x(gamma, xhat, dy)
    dxhat = dy .* gamma
    divar = sum(dxhat .* xmu; dims=sum_dims)
    dxmu1 = dxhat .* ivar
    dsqrtvar = -(sqrtvar .^ -2) .* divar
    dvar = _05 ./ sqrtvar .* dsqrtvar
    dsq = one(typ) ./ P .* dvar
    dxmu2 = 2 .* xmu .* dsq
    dx1 = (dxmu1 .+ dxmu2)
    dmu = -sum(dxmu1 .+ dxmu2; dims=sum_dims)
    # dx2 = Yota.∇mean(x, dmu, sum_dims)
    dx2 = Yota.∇mean(dmu, mean, x, sum_dims)[2]
    dx = dx1 .+ dx2
    return dgamma, dbeta, dx
end

batch_mean(x, dims) = mean(x, dims=dims)
function batch_var(x, dims)
    # note that this is equivalent to Statistics.var(x; dims=dims, corrected=false)
    # but also works for CuArrays
    mu = mean(x, dims=dims)
    return mean((x .- mu) .^ 2; dims=dims)
end

# BatchNorm layer is somewhat tricky: on one hand, we want to track running mean and var during
# training, and since these are mutating operations, we need to hide them inside primitives -
# something like `forward(m:BatchNorm, x)`. On other hand, Yota doesn't support diffrules
# w.r.t. struct fields, so `forward()` above cannot be a primitive. To overcome it, we introduce
# a helper function batchnorm2d() which pretends to be a normal primitive but:
#  * takes a @nodiff struct
#  * separately, takes learnable gamma and beta parameters
function batchnorm2d(m::BatchNorm2d, gamma, beta, x)
    if m.training
        dims = (1, 2, 4)
        mu = batch_mean(x, dims)
        sigma2 = batch_var(x, dims)
        # track running mean and variance for each channel separately
        # based on: https://github.com/ptrblck/pytorch_misc/blob/master/batch_norm_manual.py#L50-L69
        # note: running stats aren't used during training, but are required for evaluation
        m.count += 1
        exp_avg_factor = eltype(x)(1.0) / m.momentum
        m.running_mean .= @. exp_avg_factor * mu + (1 - exp_avg_factor) * m.running_mean
        m.running_var .= @. exp_avg_factor * sigma2 * m.count / (m.count - 1) +
            (1 - exp_avg_factor) * m.running_var
        return batchnorm_impl(gamma, beta, x, mu, sigma2, m.momentum;
                         eps=m.eps, training=true)
    else
        return batchnorm_impl(gamma, beta, x, m.running_mean, m.running_var, m.momentum;
                         eps=m.eps, training=false)
    end
end

(m::BatchNorm2d)(x) = batchnorm2d(m, m.gamma, m.beta, x)

function ∇batchnorm2d(dy, m::BatchNorm2d, gamma, beta, x)
    if m.training
        dims = (1, 2, 4)
        return ∇batchnorm_impl(gamma, beta, x, dy, batch_mean(x, dims), batch_var(x, dims), m.momentum;
                          eps=m.eps, training=true)
    else
        return ∇batchnorm_impl(gamma, beta, x, dy, m.running_mean, m.running_var, m.momentum;
                          eps=m.eps, training=false)
    end
end

function ChainRulesCore.rrule(::typeof(batchnorm2d), m::BatchNorm2d, gamma, beta, x)
    y = batchnorm2d(m, gamma, beta, x)
    function batchnorm2d_pullback(dy)
        dgamma, dbeta, dx = ∇batchnorm2d(dy, m, gamma, beta, x)
        return NoTangent(), Tangent{BatchNorm2d}(gamma=dgamma, beta=dbeta), dgamma, dbeta, dx
    end
    return y, batchnorm2d_pullback
end

# function register_batchnorm_derivs()
#     # @nodiff batchnorm2d(_m, _gamma, _beta, x) _m
#     # @diffrule batchnorm2d(_m, _gamma, _beta, x) _gamma getindex(∇batchnorm2d(dy, _m, _gamma, _beta, x), 1)
#     # @diffrule batchnorm2d(_m, _gamma, _beta, x) _beta getindex(∇batchnorm2d(dy, _m, _gamma, _beta, x), 2)
#     # @diffrule batchnorm2d(_m, _gamma, _beta, x) x getindex(∇batchnorm2d(dy, _m, _gamma, _beta, x), 3)
# end

function forward(m::BatchNorm2d, x)
    size(x, ndims(x)-1) == length(BN.β) ||
        error("BatchNorm expected $(length(BN.β)) channels, got $(size(x, ndims(x)-1))")
    dims = length(size(x))
    channels = size(x, dims-1)
    affine_shape = ntuple(i->i == ndims(x) - 1 ? size(x, i) : 1, ndims(x))
    m = div(prod(size(x)), channels)
    γ = reshape(BN.γ, affine_shape...)
    β = reshape(BN.β, affine_shape...)
    if !istraining()
        μ = reshape(BN.μ, affine_shape...)
        σ² = reshape(BN.σ², affine_shape...)
        ϵ = BN.ϵ
    else
        T = eltype(x)
        axes = [1:dims-2; dims] # axes to reduce along (all but channels axis)
        μ = mean(x, dims = axes)
        σ² = sum((x .- μ) .^ 2, dims = axes) ./ m
        ϵ = convert(T, BN.ϵ)
        # update moving mean/std
        mtm = BN.momentum
        S = eltype(BN.μ)
        BN.μ  = (1 - mtm) .* BN.μ .+ mtm .* S.(reshape(μ, :))
        BN.σ² = (1 - mtm) .* BN.σ² .+ (mtm * m / (m - 1)) .* S.(reshape(σ², :))
    end

    let λ = BN.λ
        x̂ = (x .- μ) ./ sqrt.(σ² .+ ϵ)
        λ.(γ .* x̂ .+ β)
    end
end


text:
# the channel example using an activity-based approach
#
using DiscreteEvents, Printf, Random

mutable struct Server
  id::Int64
  name::AbstractString
  input::Channel
  output::Channel
  op     # operation
  token  # current token

  Server(id, name, input, output, op) = new(id, name, input, output, op, nothing)
end

arrive(S::Server) = event!(SF(serve, S), SF(isready, S.input))

function serve(S::Server)
    S.token = take!(S.input)
    @printf("%5.2f: %s %d took token %d\n", tau(), S.name, S.id, S.token)
    event!((SF(put!, S.output, S.op(S.id, S.token)), SF(arrive, S)), after, rand())
end

reset!(𝐶)
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

s = shuffle(1:8)
for i in 1:2:8
    arrive(Server(s[i], "foo", ch1, ch2, +))
    arrive(Server(s[i+1], "bar", ch2, ch1, *))
end

put!(ch1, 1) # put first token into channel 1

run!(𝐶, 10)


text:
import Markdown

render(i::Inline, md::Markdown.MD) = render(i, renderMD(md))

function render(e::Editor, md::Markdown.MD)
  mds = Atom.CodeTools.flatten(md)
  out = length(mds) == 1 ? Text(chomp(sprint(show, MIME"text/markdown"(), md))) :
                           Tree(Text("MD"), [renderMD(md)])
  render(e, out)
end

renderMD(md::Markdown.MD) = renderMD(md.content)

renderMD(md::AbstractString) = renderMD(Markdown.parse(md))

renderMD(md::Vector) = Hiccup.div([renderMD(x) for x in md], class = "markdown")

function renderMD(header::Markdown.Header{T}) where T
  Hiccup.Node(Symbol(:h, T), renderMDinline(header.text)) end

function renderMD(code::Markdown.Code)
  Hiccup.pre(
    Hiccup.code(code.code,
                class = !isempty(code.language) ? code.language : "julia",
                block = true
    )
  )
end

function renderMD(md::Markdown.Paragraph)
  Hiccup.Node(:p, renderMDinline(md.content))
end

function renderMD(md::Markdown.BlockQuote)
  Hiccup.Node(:blockquote, renderMD(md.content))
end

function renderMD(md::Markdown.LaTeX)
  Hiccup.Node(:latex, latex2katex(md.formula), class = "latex block", block = true)
end

function renderMD(f::Markdown.Footnote)
  Hiccup.div([
    Hiccup.Node(:p, f.id, class = "footnote-title"),
    renderMD(f.text)
  ], class = "footnote", id = "footnote-$(f.id)")
end

function renderMD(md::Markdown.Admonition)
  icon = "icon-info"
  if md.category == "warning"
    icon = "icon-alert"
  end
  Hiccup.div([
    Hiccup.Node(:p, md.title, class = "admonition-title $icon"),
    renderMD(md.content)
  ], class = "admonition $(md.category)")
end

function renderMD(md::Markdown.List)
  Hiccup.Node(Markdown.isordered(md) ? :ol : :ul, [Hiccup.li(renderMD(item)) for item in md.items],
              start = md.ordered > 1 ? string(md.ordered) : "")
end

function renderMD(md::Markdown.HorizontalRule)
  Hiccup.Node(:hr)
end

function renderMD(link::Markdown.Link)
  Hiccup.Node(:a, renderMDinline(link.text), href = link.url)
end

function renderMD(md::Markdown.Table)
  Hiccup.table([
    Hiccup.tr([Hiccup.Node(i == 1 ? :th : :td, renderMDinline(c)) for c in row])
    for (i, row) in enumerate(md.rows)
  ])
end

# Inline elements

function renderMDinline(content::Vector)
  [renderMDinline(x) for x in content]
end

function renderMDinline(code::Markdown.Code)
  Hiccup.code(code.code,
              class = !isempty(code.language) ? "language-$(code.language)" : "julia",
              block = false) # htmlesc?
end

function renderMDinline(md::Union{Symbol,AbstractString})
  md # htmlesc?
end

function renderMDinline(md::Markdown.Bold)
  Hiccup.strong(renderMDinline(md.text))
end

function renderMDinline(md::Markdown.Italic)
  Hiccup.Node(:em, renderMDinline(md.text))
end

function renderMDinline(md::Markdown.Image)
  Hiccup.img(src = md.url, alt = md.alt)
end

function renderMDinline(f::Markdown.Footnote)
  Hiccup.Node(:a, Hiccup.span("[$(f.id)]"), href = "#footnote-$(f.id)", class = "footnote")
end

function renderMDinline(link::Markdown.Link)
  Hiccup.Node(:a, renderMDinline(link.text), href = link.url)
end

function renderMDinline(md::Markdown.LaTeX)
  Hiccup.Node(:latex, latex2katex(md.formula), class = "latex inline", block = false)
end

function renderMDinline(br::Markdown.LineBreak)
  Hiccup.Node(:br)
end

function renderMDinline(x::Expr)
  Hiccup.code(x,
              class = "julia",
              block = false) # htmlesc?
end

renderMDinline(x) = x

# katex doesn't support certain latex expressions. Need to transform those to something
# that *is* supported or get rid of them altogether.
function latex2katex(code)
  code = replace(code, "\\operatorname" => "\\mathrm")
  code = replace(code, "\\latex" => "\\katex")
  # TODO: Unicode -> LaTeX sequences.
  # KaTeX unfortunately doesn't handle Unicode very well, so something like π errors. Should
  # be possible to transform that to \pi in many simple cases.
  return code
end


text:
function HFcContainer(::Type{K}, D) where K
    ArrayPartition((zeros(K,(D,D)),zeros(K,D), zeros(K,1)))
end

#=

struct HFcContainer{D,TH,TF,Tc}
    H::TH
    F::TF
    c::Tc

    function HFcContainer(H::TH, F::TF, c::Tc) where {TH,TF,Tc}
        D = length(F)
        new{D,TH,TF,Tc}(H, F, c)
    end
end

function HFcContainer(::Type{T}, D::Integer) where T
    H = zeros(T, (D,D))
    F = zeros(T, D)
    c = zeros(T, 1)
    HFcContainer(H, F, c)
end

Base.size(var::HFcContainer{D}) where D = ((D,D),D,1)
Base.eltype(var::HFcContainer) = eltype(var.H)
Base.length(var::HFcContainer{D}) where D = D*(D+1)+1

function Base.getindex(var::HFcContainer{D}, i::Int) where D
    i <= D*D ? var.H[i] : (i <= D*(D+1) ? var.F[i-D*D] : var.c[i-D*(D+1)])
end

function Base.getindex(var::HFcContainer, i::Int, I::Vararg{Int,N}) where N
    i == 1 ? var.H[I...] : ( i== 2 ? var.F[I...] : var.c[I...] )
end

function Base.getindex(var::HFcContainer, ::Colon)
    var.H[:], var.F[:], var.c[:]
end

function Base.setindex!(var::HFcContainer{D}, v, i::Int) where D
    (
        i <= D*D ?
        (var.H[i]=v) :
        (
            i <= D*(D+1) ?
            (var.F[i-D*D]=v) :
            (var.c[i-D*(D+1)]=v)
        )
    )
end

function Base.setindex!(var::HFcContainer, v, i::Int, I::Vararg{Int,N}) where N
    (
        i == 1 ?
        (var.H[I...] = v) :
        (
            i == 2 ?
            (var.F[I...] = v) :
            (var.c[I...] = v)
        )
    )
end

function Base.setindex!(var::HFcContainer, v, ::Colon)
    var.H[:], var.F[:], var.c[:] .= v
end

function Base.similar(var::HFcContainer{D}) where D
    HFcContainer(similar(var.H), similar(var.F), similar(var.c))
end

function Base.similar(var::HFcContainer{D}, ::Type{K}) where {D,K}
    HFcContainer(similar(var.H, K), similar(var.F, K), similar(var.c, K))
end

function Base.zero(var::HFcContainer{D}) where D
    HFcContainer(zero(var.H), zero(var.F), zero(var.c))
end

OrdinaryDiffEq.recursive_bottom_eltype(a::HFcContainer) = eltype(a.H)
OrdinaryDiffEq.recursive_unitless_bottom_eltype(a::HFcContainer) = eltype(a.H)
=#


text:
function OrientedPoint(arg0::Vector1D, arg1::jboolean, arg2::jdouble)
    return OrientedPoint((Vector1D, jboolean, jdouble), arg0, arg1, arg2)
end

function copy_self(obj::OrientedPoint)
    return jcall(obj, "copySelf", OrientedPoint, ())
end

function empty_hyperplane(obj::OrientedPoint)
    return jcall(obj, "emptyHyperplane", SubOrientedPoint, ())
end

function get_location(obj::OrientedPoint)
    return jcall(obj, "getLocation", Vector1D, ())
end

function get_offset(obj::OrientedPoint, arg0::JVector)
    return jcall(obj, "getOffset", jdouble, (JVector,), arg0)
end

function get_offset(obj::OrientedPoint, arg0::Point)
    return jcall(obj, "getOffset", jdouble, (Point,), arg0)
end

function get_tolerance(obj::OrientedPoint)
    return jcall(obj, "getTolerance", jdouble, ())
end

function is_direct(obj::OrientedPoint)
    return jcall(obj, "isDirect", jboolean, ())
end

function project(obj::OrientedPoint, arg0::Point)
    return jcall(obj, "project", Point, (Point,), arg0)
end

function revert_self(obj::OrientedPoint)
    return jcall(obj, "revertSelf", void, ())
end

function same_orientation_as(obj::OrientedPoint, arg0::Hyperplane)
    return jcall(obj, "sameOrientationAs", jboolean, (Hyperplane,), arg0)
end

function whole_hyperplane(obj::OrientedPoint)
    return jcall(obj, "wholeHyperplane", SubOrientedPoint, ())
end

function whole_space(obj::OrientedPoint)
    return jcall(obj, "wholeSpace", IntervalsSet, ())
end