text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

using QuadGK
import Optim

const ϵ_reg = 1e-6

struct ParametricNoiseModel
    m_hat::Vector{Float64}
    σ_hat::Vector{Float64}
    model::Function
    model_init::Function
    params
    params_init
end

struct localEstimation{T<:AbstractFloat}
    y::T
    σ::T
    κ::T
    n::Int
end

function foi_noiseestimation(z::AbstractArray{Float32}; τ=0.2, maxnum_pairs::Int, verbose::Bool=false)
    # τ = 0.2 is good to reject a lot of outliers.

    @assert maximum(z) <= 1.0
    @assert minimum(z) >= 0.0

    # Compute the local sample mean and standard deviation over a smooth region
    println("Computing local noise variance.")
    num_imgs = size(z,3)
    z_stack = (z[:,:,i] for i in 1:num_imgs)
    est_pairs_array = map(x->localnoiseestimation(x, τ)[1], z_stack)
    est_pairs = vcat(est_pairs_array...)

    println("Initializing parameters by least-squares.")
    ab_init, y_hat, σ_hat = initialize_parameters(est_pairs)

    num_pairs = length(est_pairs)
    if num_pairs > maxnum_pairs
        idx = rand(1:num_pairs, maxnum_pairs)
        est_pairs = est_pairs[idx]
    end

    println("Initialization done.")
    println("Starting likelihood maximization.")

    likelihood0(x::Vector{Float64}) = nonclipped_negloglikelihood(x, est_pairs)
    result = Optim.optimize(likelihood0, ab_init, Optim.NelderMead(), Optim.Options(show_trace = verbose))

    println("Finished the maximization.")

    ab_hat = Optim.minimizer(result)
    @assert Optim.converged(result)

    return ab_hat
end

function localnoiseestimation(z::AbstractArray{Float32,2}, τ)
    # Wavelet and scaling functions used in the original paper are the followings:
    # ψ = Array{Float32}([0.035, 0.085, -0.135, -0.460, 0.807, -0.333])
    # ϕ = Array{Float32}([0.025, -0.060, -0.095, 0.325, 0.571, 0.235])
    ϕ = [0.035226291882100656f0, -0.08544127388224149f0, -0.13501102001039084f0, 0.4598775021193313f0, 0.8068915093133388f0, 0.3326705529509569f0]
    ϕ ./= sum(ϕ)

    ψ = [-0.3326705529509569f0, 0.8068915093133388f0, -0.4598775021193313f0, -0.13501102001039084f0, 0.08544127388224149f0, 0.035226291882100656f0]
    ψ ./= norm(ψ)

    σ_gauss = 1.2f0
    gauss = [exp(-x^2 / (2.0f0 * σ_gauss^2)) for x in -10.f0:10.f0]
    gauss ./= sum(gauss)

    z_wdet = circconv(z, ψ)[1:2:end, 1:2:end]
    z_wapp = circconv(z, ϕ)[1:2:end, 1:2:end]

    ω = ones(Float32, 7) ./ 7.0f0
    z_smo = circconv(z_wapp, ω, ω)
    s = sqrt(0.5 * π) .* circconv(abs.(z_wdet), ω, ω)

    g = [-0.5f0, 0.0f0, 0.5f0]
    smoothed_zwapp = circconv(z_wapp, gauss, gauss)
    dx_wapp = circconv(smoothed_zwapp, [1.0f0], g)
    dy_wapp = circconv(smoothed_zwapp, g, [1.0f0])

    x_smo = sqrt.(dx_wapp.^2 .+ dy_wapp.^2) .< τ .* s
    N = length(z_smo)
    num_bins = 300
    histogram_zwapp = [Vector{Float32}() for _ in 1:num_bins]
    histogram_zwdet = [Vector{Float32}() for _ in 1:num_bins]
    min_wapp = sum(ϕ[ϕ .<  0.f0])
    max_wapp = sum(ϕ[ϕ .>= 0.f0])
    Δ = (max_wapp - min_wapp) / num_bins

    for i in 1:N
        if x_smo[i]
            idx = floor(Int, (z_smo[i] - min_wapp) / Δ)
            push!(histogram_zwapp[idx], z_wapp[i])
            push!(histogram_zwdet[idx], z_wdet[i])
        end
    end

    est_pairs = Vector{localEstimation{Float64}}()
    for i in 1:num_bins
        if histogram_zwdet[i] != []
            n = length(histogram_zwdet[i])
            κ = madBiasFactor(n)
            σ = sqrt(max(Float64(mad(histogram_zwdet[i], κ))^2, .0))
            y = Float64(mean(histogram_zwapp[i]))
            push!(est_pairs, localEstimation(y, σ, κ, n))
        end
    end

    return est_pairs, x_smo
end

function initialize_parameters(est_pairs::Vector{localEstimation{T}}) where T
    num_pairs = length(est_pairs)

    y_hat = zeros(num_pairs)
    σ_hat = zeros(num_pairs)

    Φ_tmp = Vector{Float64}()
    v = Vector{Float64}()

    for i = 1:num_pairs
        tmp = est_pairs[i]
        y_hat[i] = tmp.y
        σ_hat[i] = tmp.σ
        push!(Φ_tmp, tmp.y)
        push!(v, tmp.σ^2)
    end
    Φ = hcat(Φ_tmp, ones(length(Φ_tmp)))
    ab0 =  Φ \ v
    return ab0, y_hat, σ_hat
end

function nonclipped_negloglikelihood(ab::Vector{Float64},  est_pairs::Vector{localEstimation{T}}) where T

    Δ = .1
    total_val = .0

    for l in est_pairs
        c_i = 1.0 / l.n
        d_i = 1.35 / (l.n + 1.5)
        y_i = l.y
        σ_i = l.σ

        # This integrand is sometimes a very sharp peak-like function like a Dirac funciton.
        # Therefore, the direct numerical integration over [0.0, 1.0] is realatively difficult
        # because the algorithm may not evaluate the integrand around the peak.
        # To avoid this issue, the integration interval is separated into multiple intervals.
        # Since the peak is known to be close to y_i, the function value at y_i should be enough large than zero.
        integrand(y::Float64)::Float64 = 1.0 / σsq_reg(y, ab) *
            exp(-1.0 / (2.0 * σsq_reg(y, ab)) * ( (y_i - y)^2 / c_i + (σ_i - sqrt(σsq_reg(y, ab)) )^2 / d_i))

        val = .0
        if y_i - Δ > .0
            val += quadgk(integrand, .0, y_i - Δ)[1]
        end
        val += quadgk(integrand, max(.0, y_i - Δ), y_i)[1]
        val += quadgk(integrand, y_i, min(1.0, y_i + Δ))[1]
        if y_i + Δ < 1.0
            val += quadgk(integrand, y_i + Δ, 1.0)[1]
        end

        total_val -= log(1 / (2π * sqrt(c_i * d_i) ) * abs(val) + ϵ_reg)
    end

    return total_val
end

σsq_reg(y::Float64, ab::Vector{Float64}) = max(ϵ_reg^2, σsq(y, ab))
σsq(y::Float64, ab::Vector{Float64}) = ab[1] * y + ab[2]


text:
"""
    pdesolve(f, grid, yend [, τs]; is_algebraic = OrderedDict(k => false for k in keys(yend)), bc = nothing)

    solves a system of pdes with an arbitrary number of functions and up to state variable. Denote  F the number of functions to solve for and S the number of state variables.

    ### Arguments
    *  f: is a function (state, y) -> out
    where:
       state is a tuple of real numbers of  length S (for the state values)
       y is a tuple of real numbers of length F * 3 with the functions as well as their first and second derivatives at the state
       out is a named tuple of length F with the value of the time derivatives of the functions at the state
    * grid is an OrderedDict that, for each state, associates an AbstractVector (for the grid)
    * yend is an OrderedDict that gives the terminal value of each function in the system of PDEs
    * τs (optional) is a time grid on which to solve the pde on. In this case, yend corresponds to the solution at time τs[end]
"""
function pdesolve(apm, @nospecialize(grid), @nospecialize(yend), τs::Union{Nothing, AbstractVector} = nothing; is_algebraic = OrderedDict(k => false for k in keys(yend)), bc = nothing, verbose = true, kwargs...)
    stategrid = StateGrid(NamedTuple(grid))
    S = size(stategrid)
    all(size(v) == S for v in values(yend)) || throw(ArgumentError("The length of initial guess (e.g. terminal value) does not equal the length of the state space"))
    all(keys(is_algebraic) .== keys(yend)) || throw(ArgumentError("the terminal guess yend and the is_algebric keyword argument must have the same names"))
    is_algebraic = OrderedDict(first(p) => fill(last(p), S) for p in pairs(is_algebraic))
    if τs isa AbstractVector
        issorted(τs) || throw(ArgumentError("The set of times must be increasing."))
        ys = [OrderedDict(first(p) => collect(last(p)) for p in pairs(yend)) for t in τs]
    else
        y = OrderedDict(first(p) => collect(last(p)) for p in pairs(yend))
    end
    # convert to Matrix
    yend_M = catlast(values(yend))
    is_algebraic_M = catlast(values(is_algebraic))
    bc_M = _Array_bc(bc, yend, grid)

    # create sparsity
    J0c = sparsity_jac(stategrid, yend)
    Tsolution = Type{tuple(keys(yend)...)}

    # iterate on time
    if τs isa AbstractVector
        apm_onestep = hasmethod(apm, Tuple{NamedTuple, NamedTuple, Number}) ? (state, grid) -> apm(state, grid, τs[end]) : apm
        a = get_a(apm_onestep, stategrid, Tsolution, yend_M, bc_M)
        as = (a === nothing) ? nothing : [deepcopy(a) for τ in τs]
        y_M = yend_M
        residual_norms = zeros(length(τs))
        if verbose
            @printf "    Time Residual\n"
            @printf "-------- --------\n"
        end
        for iτ in length(τs):(-1):1
            _setindex!(ys[iτ], y_M)
            apm_onestep = hasmethod(apm, Tuple{NamedTuple, NamedTuple, Number}) ? (state, grid) -> apm(state, grid, τs[iτ]) : apm
            if a !== nothing
                _setindex!(as[iτ], apm_onestep, stategrid, Tsolution, y_M, bc_M)
                as[iτ] = merge(ys[iτ], as[iτ])
            end
            if iτ > 1
                y_M, residual_norms[iτ] = implicit_timestep((ydot, y) -> hjb!(apm_onestep, stategrid, Tsolution, ydot, y, bc_M, size(yend_M)), vec(y_M), τs[iτ] - τs[iτ-1]; is_algebraic = vec(is_algebraic_M), verbose = false, J0c = J0c, kwargs...)
                if verbose
                    @printf "%8g   %8.4e\n" τs[iτ-1] residual_norms[iτ]
                end
                y_M = reshape(y_M, size(yend_M)...)
            end
        end
        return EconPDEResult(ys, residual_norms, as)
    else
        a = get_a(apm, stategrid, Tsolution, yend_M, bc_M)
        y_M, residual_norm = finiteschemesolve((ydot, y) -> hjb!(apm, stategrid, Tsolution, ydot, y, bc_M, size(yend_M)), vec(yend_M); is_algebraic = vec(is_algebraic_M),  J0c = J0c, verbose = verbose, kwargs... )
        y_M = reshape(y_M, size(yend_M)...)
        _setindex!(y, y_M)
        if a !== nothing
            _setindex!(a, apm, stategrid, Tsolution, y_M, bc_M)
            a = merge(y, a)
        end
        return EconPDEResult(y, residual_norm, a)
    end
end

catlast(iter) = cat(iter...; dims = ndims(first(iter)) + 1)

_Array_bc(::Nothing, yend, grid) = zeros(size(first(values(yend)))..., length(yend))
function _Array_bc(bc, yend, grid)
    keys_grid = collect(keys(grid))
    bc_M = _Array_bc(nothing, yend, grid)
    for (k, yname) in enumerate(keys(yend))
        if length(keys_grid) == 1
            bc_M[1, k], bc_M[end, k] = bc[Symbol(yname, keys_grid[1])]
        elseif length(keys_grid) == 2
            bc_M[1, :,  k],  bc_M[end, :, k] = bc[Symbol(yname, keys_grid[1])]
            bc_M[:, 1,  k], bc_M[:, end,  k] = bc[Symbol(yname, keys_grid[2])]
        end
    end
    return bc_M
end

function get_a(apm, stategrid::StateGrid, Tsolution, y_M::AbstractArray, bc_M::AbstractArray)
    i0 = first(eachindex(stategrid))
    derivatives = differentiate(Tsolution, stategrid, y_M, i0, bc_M)
    result = apm(stategrid[i0], derivatives)
    if length(result) == 1
        return nothing
    else
        return OrderedDict(a_key => Array{Float64}(undef, size(stategrid)) for a_key in keys(result[2]))
    end
end

function sparsity_jac(stategrid::StateGrid, @nospecialize(yend))
    s = size(stategrid)
    l = prod(s)
    t = (ndims(stategrid), length(yend) > 1)
    if t == (1, 0)
        J = Tridiagonal(ones(l - 1), ones(l), ones(l -1))
        return J, matrix_colors(J)
    elseif t == (2, 0)
        J = BandedBlockBandedMatrix(Ones(l, l), fill(s[1], s[2]), fill(s[1], s[2]), (1, 1), (1, 1))
        return sparse(J), matrix_colors(J)
    elseif t == (1, 1)
        J = BandedBlockBandedMatrix(Ones(l * length(yend), l * length(yend)), fill(l, length(yend)) ,fill(l, length(yend)), (length(yend) - 1, length(yend) - 1), (1, 1))
        return sparse(J), matrix_colors(J)
    elseif t == (2, 1)
        J = BandedBlockBandedMatrix(Ones(l * length(yend), l * length(yend)), repeat(fill(s[1], s[2]), outer = length(yend)), repeat(fill(s[1], s[2]), outer = length(yend)), (s[2] * length(yend) - 1, s[2] * length(yend) - 1), (1, 1))
        return sparse(J), matrix_colors(J)
    else
        return nothing, nothing
    end
end

function _setindex!(@nospecialize(y), y_M::AbstractArray)
    for (i, v) in enumerate(values(y))
        v[:] = selectdim(y_M, ndims(y_M), i)
    end
end

function _setindex!(@nospecialize(a), apm, stategrid::StateGrid, Tsolution, y_M::AbstractArray, bc_M::AbstractArray)
    for i in eachindex(stategrid)
         solution = differentiate(Tsolution, stategrid, y_M, i, bc_M)
         outi = apm(stategrid[i], solution)[2]
         for (k, v) in zip(values(a), values(outi))
            k[i] = v
         end
     end
 end

 # create hjb! that accepts and returns AbstractVector rather than AbstractArrays
 function hjb!(apm, stategrid::StateGrid, Tsolution, ydot::AbstractVector, y::AbstractVector, bc_M::AbstractArray, ysize::NTuple)
     y_M = reshape(y, ysize...)
     ydot_M = reshape(ydot, ysize...)
     vec(hjb!(apm, stategrid, Tsolution, ydot_M, y_M, bc_M))
 end

 function hjb!(apm, stategrid::StateGrid, Tsolution, ydot_M::AbstractArray, y_M::AbstractArray, bc_M::AbstractArray)
     Tt = [Symbol(v, :t) for v in Tsolution.parameters[1]]
     for i in eachindex(stategrid)
         solution = differentiate(Tsolution, stategrid, y_M, i, bc_M)
         outi = apm(stategrid[i], solution)
         if isa(outi[1], Number)
            _setindex!(ydot_M, Tsolution, outi, i)
        else
            _setindex!(ydot_M, Tsolution, outi[1], i)
        end
     end
     return ydot_M
 end

 @generated function _setindex!(ydot_M::AbstractArray, ::Type{Tsolution}, outi::NamedTuple, i::CartesianIndex) where {Tsolution}
     N = length(Tsolution.parameters[1])
     quote
          $(Expr(:meta, :inline))
          $(Expr(:block, [Expr(:call, :setindex!, :ydot_M, Expr(:call, :getproperty, :outi, Meta.quot(Symbol(Tsolution.parameters[1][k], :t))), :i, k) for k in 1:N]...))
     end
 end

ext: jl


text:
# A modified next reaction method Heaps version of mnrm_v1
struct MNRM <: AbstractAggregatorAlgorithm end
mutable struct MNRMJumpAggregation{T,S,F1,F2,RNG,DEPGR,PQ} <: AbstractSSAJumpAggregator
  next_jump::Int
  prev_jump::Int
  next_jump_time::T
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  ma_jumps::S
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
  rng::RNG
  dep_gr::DEPGR
  pq::PQ
end
function MNRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; num_specs, dep_graph=nothing, kwargs...) where {T,S,F1,F2,RNG}

    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the Next Reaction Method (NRM) algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph

        # make sure each jump depends on itself
        add_self_dependencies!(dg)
    end
    pq = MutableBinaryHeap{T, DataStructures.FasterForward}()
    MNRMJumpAggregation{T,S,F1,F2,RNG,typeof(dg), typeof(pq)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng, dg, pq)
end

function aggregate(aggregator::MNRM, u, p, t, end_time, constant_jumps, ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using tuples
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)
    build_jump_aggregation(MNRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, save_positions, rng; num_specs = length(u), kwargs...)
end

function initialize!(p::MNRMJumpAggregation, integrator, u, params, t)
    fill_rates_and_get_times!(p, u, params, t)
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

@inline function execute_jumps!(p::MNRMJumpAggregation, integrator, u, params, t)
    update_state!(p, integrator, u)
    update_internal_times!(p, u, params, t)
    nothing
end
# calculate the next jump / jump time
function generate_jumps!(p::MNRMJumpAggregation, integrator, u, params, t)
    ttnj, p.next_jump = top_with_handle(p.pq)
    @fastmath p.next_jump_time = t + ttnj
    # TODO how to update the whole thing
    for rx in eachindex(p.cur_rates)
        update!(p.pq, rx, p.pq[rx] - ttnj)
    end
    nothing
end

@inline function update_internal_times!(p::MNRMJumpAggregation, u, params, t)
    @unpack pq, cur_rates, rates, ma_jumps, rng = p
    dep_rxs = p.dep_gr[p.next_jump]
    num_majumps = get_num_majumps(ma_jumps)
    @inbounds for rx in dep_rxs
        oldrate = cur_rates[rx]
        @inbounds cur_rates[rx] = calculate_jump_rate(ma_jumps,num_majumps,rates,u,params,t,rx)
        if rx != p.next_jump && oldrate > zero(oldrate)
            if cur_rates[rx] > zero(eltype(cur_rates))
                update!(pq, rx, pq[rx]*oldrate/cur_rates[rx])
            else
                update!(pq, rx, typemax(t))
            end
        else
            if cur_rates[rx] > zero(eltype(cur_rates))
                update!(pq, rx, randexp(rng)/cur_rates[rx])
            else
                update!(pq, rx, typemax(t))
            end
        end
    end
    nothing
end

# fill the propensity rates
# reevaulate all rates, recalculate all jump times, and reinit the priority queue
function fill_rates_and_get_times!(p::MNRMJumpAggregation, u, params, t)

    # mass action jumps
    majumps   = p.ma_jumps
    cur_rates = p.cur_rates
    pqdata = Vector{typeof(t)}(undef,length(cur_rates))
    @inbounds for i in 1:get_num_majumps(majumps)
        cur_rates[i] = evalrxrate(u, i, majumps)
        pqdata[i] = randexp(p.rng) / cur_rates[i]
    end

    # constant rates
    rates = p.rates
    idx   = get_num_majumps(majumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        pqdata[idx] = randexp(p.rng) / cur_rates[idx]
        idx += 1
    end

    # setup a new indexed priority queue to storing rx times
    p.pq = MutableBinaryHeap{typeof(t),DataStructures.FasterForward}(pqdata)
    nothing
end

# Change below


text:
using Drawing; const D = Drawing

fm = D.get_font_map_default()
ff = D.list_families(fm)
for f in ff
    m = D.is_monospace(f) ? "monospace" : ""
    println("$(D.get_name(f)) $m")
    fc = D.list_faces(f)
    for c in fc
        d1 = D.describe(c)
        println("  $(D.get_name(c)): $(d1)")
        println("  $(D.get_style(d1)), $(D.get_variant(d1)), $(D.get_weight(d1)), $(D.get_stretch(d1)), $(D.get_gravity(d1))")
        d2 = convert(D.FontDescription, "$(d1)")
        @assert d1 == d2
        @assert hash(d1) == hash(d2)
    end
    gc()  # call finalizer
end


text:
# using Pkg
# Pkg.add("Documenter")
using Documenter, PolyChaos

include("pages.jl")

makedocs(
    sitename = "PolyChaos.jl",
    format = Documenter.HTML(analytics = "UA-90474609-3",
                         assets = ["assets/favicon.ico"],
                         canonical="https://polychaos.sciml.ai/stable/"),
    modules = [PolyChaos],
    authors = "tillmann.muehlpfordt@kit.edu",
    doctest = false,
    pages = pages
)

deploydocs(
    repo = "github.com/SciML/PolyChaos.jl.git",
    target = "build",
)

max_stars_repo_name: timueh/PolyChaos

max_issues_repo_name: timueh/PolyChaos

max_forks_repo_name: timueh/PolyChaos


text:
include("../../external/mbh2dfmm/julia/modbhrouts.jl")

using FINUFFT

## FREESPACE SOLVER WITH NUFFT
"
fs_modstokes(F1, F2, L, lambda, xnu, ynu; precomp=GR)

Solves modified Stokes equation on NxN grid of length L with rhs with
freespace boundary conditions and compactly supported rhs (F1, F2).

GR is precomputed truncated Green's function to modified biharmonic:
GR = fs_modbh_precomp(N, L, lambda)
"
function fs_modstokes(F1, F2, L, lambda, xnu=[], ynu=[]; precomp=[])
    # Check input
    M = size(F1)
    @assert length(M)>1 && M[1]==M[2] "Grid must be square"
    @assert size(F1)==size(F2)
    N = M[1]
    # Get truncated Green's function
    if isempty(precomp)
        GRpre = fs_modbh_precomp(N, L, lambda)
    else
        GRpre = precomp
    end
    # Pad input to double size
    sf = 2
    N2 = sf*N
    info("FFT size is $(N2)x$(N2)")
    Fpad1 = zeros(Complex128, N2, N2)
    Fpad2 = zeros(Complex128, N2, N2)
    Fpad1[1:N, 1:N] = F1
    Fpad2[1:N, 1:N] = F2
    # Do FFTs in-place
    fft!(Fpad1)
    fft!(Fpad2)
    Fhat1, Fhat2 = Fpad1, Fpad2
    # Convolve with truncated greens function
    Uhat1, Uhat2 = Fhat1, Fhat2 # Operate in-place
    k1, k2 = k_vectors([N2, N2], [L*sf, L*sf])
    k1 = ifftshift(k1)
    k2 = ifftshift(k2)
    k1sq, k2sq = k1.^2, k2.^2
    for j=1:N2
        for i=1:N2
            K1 = k1[i]
            K2 = k2[j]
            Ksq = k1sq[i] + k2sq[j]
            KdotFhat = K1*Fhat1[i,j] + K2*Fhat2[i,j]
            Uhat1[i,j] = GRpre[i,j]*(Ksq*Fhat1[i,j] - K1*KdotFhat)
            Uhat2[i,j] = GRpre[i,j]*(Ksq*Fhat2[i,j] - K2*KdotFhat)
        end
    end
    if length(xnu) > 0
        # NUFFT for off-grid points
        @assert length(xnu)==length(ynu)
        # TODO: Get this scaling right so we can use modeord=0
        origin = pi/sf
        scale = 2*pi/(sf*L)
        xn = origin + xnu*scale
        yn = origin + ynu*scale
        Nnu = length(xnu)
        unu1 = Array{Complex128}(Nnu)
        unu2 = Array{Complex128}(Nnu)
        opts = finufft_default_opts()
        opts.modeord = 1
        nufft2d2!(xn, yn, unu1, 1, 1e-15, Uhat1, opts)
        nufft2d2!(xn, yn, unu2, 1, 1e-15, Uhat2, opts)
        uscale = 1/(N*sf)^2
        unu1 = real(unu1)*uscale
        unu2 = real(unu2)*uscale
    else
        unu1, unu2 = [], []
    end
    # Transform back
    ifft!(Uhat1)
    ifft!(Uhat2)
    # Truncate results to original grid size
    U1 = zeros(N,N)
    U2 = zeros(N,N)
    for j=1:N
        for i=1:N
            U1[i,j] = real(Uhat1[i,j])
            U2[i,j] = real(Uhat2[i,j])
        end
    end
    return U1, U2, unu1, unu2
end

"Precompute truncated Green's function to modified biharmonic equation"
function fs_modbh_precomp(N, L, lambda)
    sf = 1+sqrt(2) # Minimum oversampling in 2D
    Nf = ceil(Int, sf*N)
    if mod(Nf,2) != 0
        Nf+=1 # Even grids are always faster (right?)
    end
    info("Precompute FFT size is $(Nf)x$(Nf)")
    sf = Nf/N # Effective oversampling
    R = sqrt(2)*L
    logR = log(R)
    rscale = 1.0
    n = 1
    qn = modbh_qn([R],lambda,rscale,n)
    Q0 = qn[1,1]
    Q1 = qn[2,1]
    lam2 = lambda*lambda
    # Compute
    GR = Array{Complex128}(Nf, Nf)
    k1, k2 = k_vectors([Nf, Nf], [L*sf, L*sf])
    k1 = ifftshift(k1)
    k2 = ifftshift(k2)
    k1sq, k2sq = k1.^2, k2.^2
    for j=1:Nf
        for i=1:Nf
            K1 = k1[i]
            K2 = k2[j]
            Ksq = k1sq[i] + k2sq[j]
            K = sqrt(Ksq)
            J1 = besselj(1, K*R)
            J0M1 = besselj0m1(K*R)
            J0 = J0M1+1
            GR[i,j] = lam2*(J0M1 + K*R*logR*J1)/(Ksq*(lam2+Ksq))
            GR[i,j] += (Q0*K*R*J1-lambda*R*Q1*J0)/(lam2+Ksq)
        end
    end
    # Correct at k=0
    mid = 1 # grid is even
    GR[mid, mid] = R*R*(-0.25+0.5*logR)- R*Q1/lambda
    GR .*= 1/lam2
    ifft!(GR)
    # GR now in real space and has rubbish in the center,
    # Truncate by picking out corner values
    GRtrunc = Array{Complex128}(2*N, 2*N)
    idx = 1:N
    offset = Nf - N
    for i=idx
        for j=idx
            GRtrunc[i,  j  ] = GR[i,       j       ]
            GRtrunc[i+N,j  ] = GR[i+offset,j       ]
            GRtrunc[i,  j+N] = GR[i,       j+offset]
            GRtrunc[i+N,j+N] = GR[i+offset,j+offset]
        end
    end
    # Return Fourier transform
    fft!(GRtrunc)
    return GRtrunc
end

## PJ FUNCTIONS
function GR_purejulia(K, R)
    J0 = besselj(0, K*R)
    J1 = besselj(1, K*R)
    K0 = besselk(0, lambda*R)
    K1 = besselk(1, lambda*R)
    GR = -( J0*(lambda*K^2*R*K1 - K^2 - lambda^2) -
            K*R*J1*((K^2+lambda^2)*log(R) + K^2*K0) + lambda^2
            ) / (K^2*(lambda^2+K^2))
    return GR
end

function GR0_purejulia(R)
    GR = -(lambda^2*R^2 - 2*lambda^2*R^2*log(R) +
           4*lambda*R*besselk(1, R*lambda) - 4)/(4*lambda^2)
    return GR
end

## HELPERS
function k_vectors(M,box)
    function k_vec(M,L)
        if mod(M,2)==0
            MM = M/2;
            k = (2*pi/L)*collect(-MM:(MM-1));
        elseif mod(M-1,2)==0
            MM = (M-1)/2;
            k = (2*pi/L)*collect(-MM:MM);
        else
            error("k-vectors not computed");
        end
        return k
    end
    k1 = k_vec(M[1], box[1])
    k2 = k_vec(M[2], box[2])
    return k1, k2
end


text:
module TestSamples

using ..Utils
using Test

const SAMPLES = joinpath(@__DIR__, "../../samples")

function test_undefined_positive()
    proc = Utils.exec("using UndefinedPositive"; append_load_path = [SAMPLES])
    @test !success(proc)
    @test occursin("Following API is declared but not defined", proc.stderr)
    @test occursin("WILL_NOT_BE_DEFINED", proc.stderr)
    @test !occursin("WILL_BE_DEFINED", proc.stderr)
end

function test_strict_negative()
    proc = Utils.exec("using StrictNegative"; append_load_path = [SAMPLES])
    @test success(proc)
end

function test_strict_positive()
    proc = Utils.exec("using StrictPositive"; append_load_path = [SAMPLES])
    @test !success(proc)
    @test occursin("Non-public API imported", proc.stderr)
end

end  # module


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["VanDerPol"] = BenchmarkGroup()

# ==============================================================================
# Jet-based approach using Taylor Models
# ==============================================================================
include("vanderpol.jl")

# ----------------------------------------
# Case 1: μ = 1
# ----------------------------------------

# benchmark settings
𝑃, 𝑂 = vanderpol(μ=1)

# algorithm-specific options
𝑂jets = Options(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2, :max_steps=>500)

# first run
sol_1 = solve(𝑃, 𝑂, op=TMJets(𝑂jets))

# verify that specification holds
@assert all([ρ([0.0, 1.0], sol_1.Xk[i].X) < 2.75 for i in eachindex(sol_1.Xk)])

# benchmark
SUITE["VanDerPol"]["μ = 1: x[2] <= 2.75"] = @benchmarkable solve($𝑃, $𝑂, op=TMJets($𝑂jets))

# ----------------------------------------
# Case 2: μ = 2
# ----------------------------------------

X0_μ2 = Hyperrectangle(low=[1.55, 2.35], high=[1.85, 2.45])

# algorithm-specific options
𝑂jets = Options(:abs_tol=>1e-10, :orderT=>8, :orderQ=>2, :max_steps=>500)

# the idea is to split the initial states along the "x" direction
# n controls the splitting
nsplits_x = 8

function compute_μ2(;n::Int=8, validate=true)
    sol_2 = []
    for X0i in split(X0_μ2, n, 1)
        𝑃, 𝑂 = vanderpol(μ=2.0, T=8.0, X0=X0i, property=(t,x) -> x[2] < 4.0)
        𝑂jets = Options(:abs_tol=>1e-10, :orderT=>8, :orderQ=>1, :max_steps=>500)
        sol_2_i = solve(𝑃, 𝑂, op=TMJets(𝑂jets))
        push!(sol_2, sol_2_i)
        if validate
            # verify that specification holds
            @assert all([ρ([0.0, 1.0], sol_2_i.Xk[i].X) < 4.0 for i in eachindex(sol_2_i.Xk)])
        end
    end
    return sol_2
end

sol_2 = compute_μ2(n=nsplits_x, validate=true)

# benchmark
SUITE["VanDerPol"]["μ = 2: x[2] <= 4.0"] = @benchmarkable compute_μ2(n=$nsplits_x, validate=false)

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# ==============================================================================
# Create plots
# ==============================================================================

# --------------------------
# Case 1
# --------------------------

plot(sol_1,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"x\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"y\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[-3., -2., -1., 0., 1., 2., 3.], ytick=[-3., -2., -1., 0., 1., 2., 3.],
     xlims=(-3., 3.), ylims=(-3., 3.),
     bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000), linecolor="red", color="red")

plot!(x->x, x->2.75, -3., 3., line=2, color="red", linestyle=:dash, legend=nothing)
savefig(@relpath "vanderpol_case_1.png")

# --------------------------
# Case 2
# --------------------------

plot_2 = plot(x->x, x->4.0, -2.5, 3., line=2, color="red", linestyle=:dash, legend=nothing)

for i in 1:nsplits_x
    plot!(plot_2, sol_2[i], tickfont=font(30, "Times"), guidefontsize=45,
                   xlab=L"x\raisebox{-0.5mm}{\textcolor{white}{.}}",
                   ylab=L"y\raisebox{2mm}{\textcolor{white}{.}}",
                   xtick=[-2., -1., 0., 1., 2., 3.], ytick=[-4., -3., -2., -1., 0., 1., 2., 3., 4.],
                   xlims=(-2.5, 3.), ylims=(-4.5, 4.),
                   bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
                   size=(1000, 1000), color="blue", linewidth=0.0, linecolor="blue", alpha=.5)
end

savefig(plot_2, @relpath "vanderpol_case_2.png")

# --------------------------
# Cases 1 and 2 overlapped
# --------------------------

plot_all = plot(x->x, x->4.0, -2.5, 2.5, line=2, color="red", linestyle=:dash, legend=nothing)

for i in 1:nsplits_x
    plot!(plot_all, sol_2[i], tickfont=font(30, "Times"), guidefontsize=45,
                   xlab=L"x\raisebox{-0.5mm}{\textcolor{white}{.}}",
                   ylab=L"y\raisebox{1mm}{\textcolor{white}{.}}",
                   xtick=[-2., -1., 0., 1., 2.], ytick=[-4., -3., -2., -1., 0., 1., 2., 3., 4.],
                   xlims=(-2.5, 2.5), ylims=(-4.5, 4.),
                   bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
                   size=(1000, 1000), color="blue", linewidth=0., linecolor="blue")
end

plot!(plot_all, sol_1, tickfont=font(30, "Times"), guidefontsize=45,
                       xlab=L"x\raisebox{-0.5mm}{\textcolor{white}{.}}",
                       ylab=L"y\raisebox{1mm}{\textcolor{white}{.}}",
                       xtick=[-2., -1., 0., 1., 2.], ytick=[-4., -3., -2., -1., 0., 1., 2., 3., 4.],
                       xlims=(-2.5, 2.5), ylims=(-4.5, 4.),
                       bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
                       size=(1000, 1000), color="red", linewidth=0., linecolor="red")

plot!(plot_all, x->x, x->2.75, -2.5, 2.5, line=2, color="red", linestyle=:dash, legend=nothing)

savefig(plot_all, @relpath "vanderpol_case_all.png")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
# TODO: remove this testset when we dicide to drop 1.0 compatibility
if VERSION < v"1.3"
    @testset "safe_rand" begin
        mutex = Augmentor.rand_mutex[]
        typeof(mutex) <: Threads.Mutex
        # check that its not a null pointer
        @test reinterpret(Int, mutex.handle) > 0

        num = @inferred Augmentor.safe_rand()
        @test 0 <= num <= 1
        @test typeof(num) <: Float64
        num = @inferred Augmentor.safe_rand(2)
        @test all(0 .<= num .<= 1)
        @test typeof(num) <: Vector{Float64}
    end
end


text:
module GridLayoutBase

using GeometryBasics
using Observables
using Match

const DEFAULT_COLGAP = Ref{Any}(20.0)
const DEFAULT_ROWGAP = Ref{Any}(20.0)
# These function refs can be mutated by other packages to override the default
# way of retrieving default column and row gap sizes
const DEFAULT_ROWGAP_GETTER = Ref{Function}(() -> DEFAULT_ROWGAP[])
const DEFAULT_COLGAP_GETTER = Ref{Function}(() -> DEFAULT_COLGAP[])

include("types.jl")
include("gridlayout.jl")
include("layout_engine.jl")
include("layoutobservables.jl")
include("gridapi.jl")
include("helpers.jl")
include("geometry_integration.jl")
include("gridlayoutspec.jl")

export GridLayout, GridPosition
export GridLayoutSpec
export BBox
export LayoutObservables
export Inside, Outside, Mixed
export Fixed, Auto, Relative, Aspect
export width, height, top, bottom, left, right
export with_updates_suspended
export appendcols!, appendrows!, prependcols!, prependrows!, deletecol!, deleterow!, trim!
export gridnest!
export AxisAspect, DataAspect
export colsize!, rowsize!, colgap!, rowgap!
export Left, Right, Top, Bottom, TopLeft, BottomLeft, TopRight, BottomRight
export grid!, hbox!, vbox!
export swap!
export protrusionsobservable, suggestedbboxobservable, reportedsizeobservable, autosizeobservable, computedbboxobservable, gridcontent
export ncols, nrows

end


text:
module MHDFlows

using
  CUDA,
  Statistics,
  SpecialFunctions,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
import Base: show, summary

"Abstract supertype for problem."
abstract type AbstractProblem end
abstract type MHDVars <: AbstractVars end

include("DyeModule.jl")
include("Problems.jl")
include("pgen.jl")
include("HDSolver.jl")
include("MHDSolver.jl")
include("DiagnosticWrapper.jl")
include("integrator.jl")
include("datastructure.jl")
include("utils/VectorCalculus.jl")
include("utils/MHDAnalysis.jl")
include("utils/GeometryFunction.jl")

export Problem,
       TimeIntegrator!,
       Curl,
       Div,
       LaplaceSolver,
       Crossproduct,
       Dotproduct,
       xy_to_polar,
       ScaleDecomposition,
       h_k,
       VectorPotential,
       LaplaceSolver,
       getL

end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
const XY{T} = Tuple{T, T}
const Polyline{T} = Vector{XY{T}}
const Polylines{T} = Vector{Polyline{T}}
struct Contour{T}
    polylines::Polylines{T}
    level::T
end

Contour(level::T) where T = Contour(Vector{XY{T}}[], level)

"""
    tricontour(x, y, z, t, levels)

Return contour lines of unstructured triangular data. `x`, `y`, and `z` are one-dimensional
input arrays of the same length (the number of points). `t` is an integer array of size
`(3,nt)`, where `nt` is the number of triangles. The coordinates of the `i`th vertex of the
`j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. Function values at vertices are
specified by the `z` array, and contours levels are defined by `levels.`

If `levels` is an integer, then equally spaced levels between the extrema of `z` will be
generated. Otherwise, `levels` is a vector of increasing values, each value corresponding
to one contour line.
"""
function tricontour(x, y, z, t, levels)
    m = TriMesh(x, y, t)
    if levels isa Integer
        levels = range(minimum(z), maximum(z), length=levels)
    end
    tricontour(m, z, levels)
end

function tricontour(m::TriMesh, z, levels)
    @assert issorted(levels)
    contours = Contour{eltype(levels)}[]
    for level=levels
        push!(contours, generate_unfilled_contours(m, z, level))
    end
    contours
end

function generate_unfilled_contours(m::TriMesh, z, level)
    reset_visited!(m)
    contour = Contour(level)
    add_bdr_lines!(contour.polylines, m, z, level)
    add_interior_lines!(contour.polylines, m, z, level; on_upper=false, filled=false)
    contour
end

function interp(m::TriMesh, p1, p2, z, level)
    t = (z[p2] - level)/(z[p2] - z[p1])
    x = t*m.x[p1] + (1.0-t)*m.x[p2]
    y = t*m.y[p1] + (1.0-t)*m.y[p2]
    x,y
end

function interp(m::TriMesh, edge::TriEdge, z, level)
    p1 = m.t[edge.ie,edge.it]
    p2 = m.t[mod1(edge.ie+1,3),edge.it]
    interp(m, p1, p2, z, level)
end

const exit_edges = [0,3,1,3,2,2,1,0] # Lookup table for exit edges using encoded config

function get_exit_edge(m::TriMesh, it, z, level; on_upper)
    config = (z[m.t[1,it]] >= level) |
             (z[m.t[2,it]] >= level) << 1 |
             (z[m.t[3,it]] >= level) << 2
    on_upper && (config = 7-config)
    exit_edges[config+1]
end

function follow_interior!(polyline::Polyline, m::TriMesh, edge::TriEdge, z, level; end_on_bdr, on_upper=false)
    push!(polyline, interp(m, edge, z, level))
    it = edge.it
    while true
        vi = on_upper ? it + size(m.t,2) : it
        !end_on_bdr && m.visited[vi] && return edge
        m.visited[vi] = true
        ie = get_exit_edge(m, it, z, level; on_upper)
        @assert ie > 0
        edge = TriEdge(ie,it)
        push!(polyline, interp(m, edge, z, level))
        # Move to next triangle
        next_edge = get_neighbor_edge(m,edge)
        it = next_edge.it
        end_on_bdr && it == 0 && return edge
        edge = next_edge
    end
end

function follow_interior(m::TriMesh, edge::TriEdge, z, level::T; end_on_bdr, on_upper=false) where T
    polyline = Polyline{T}()
    edge = follow_interior!(polyline, m, edge, z, level; end_on_bdr, on_upper)
    polyline, edge
end

function add_bdr_lines!(polylines::Polylines, m::TriMesh, z, level)
    for bdr=m.bdr_loops
        nedges = length(bdr)
        end_above = false
        for i=1:nedges
            edge = bdr[i]
            if i == 1
                start_above = z[m.t[edge.ie,edge.it]] >= level
            else
                start_above = end_above
            end
            end_above = z[m.t[mod1(edge.ie+1,3),edge.it]] >= level
            if start_above && !end_above
                polyline,_ = follow_interior(m, edge, z, level; end_on_bdr=true)
                push!(polylines, polyline)
            end
        end
    end
end

function add_interior_lines!(polylines::Polylines, m::TriMesh, z, level; on_upper, filled)
    nt = size(m.t,2)
    for it=1:nt
        vi = on_upper ? it + size(m.t,2) : it
        m.visited[vi] && continue
        m.visited[vi] = true
        ie = get_exit_edge(m, it, z, level; on_upper)
        ie == 0 && continue
        # Found a new outgoing contour line, start from neighbor
        edge = get_neighbor_edge(m, TriEdge(ie,it))
        polyline,_ = follow_interior(m, edge, z, level; on_upper, end_on_bdr=false)
        push!(polyline, first(polyline)) # Make closed loop
        push!(polylines, polyline)
    end
end


text:
using ExoplanetsSysSim
using HDF5, JLD

stellar_catalog_file_in = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv")
stellar_catalog_file_out = replace(stellar_catalog_file_in,r".csv$",".jld")

stellar_catalog = setup_star_table(ascii(stellar_catalog_file_in))
save(stellar_catalog_file_out,"stellar_catalog", stellar_catalog, "stellar_catalog_usable", ExoplanetsSysSim.StellarTable.usable)


text:
function add_constrained_variables(
    dual_problem,
    primal_model,
    variable_parameters,
)
    single_or_vector_variables_types =
        MOIU.sorted_variable_sets_by_cost(dual_problem.dual_model, primal_model)
    params = Set(variable_parameters)
    for S in single_or_vector_variables_types
        if S <: MOI.AbstractVectorSet
            _add_constrained_variables(
                dual_problem.primal_dual_map,
                primal_model,
                S,
                params,
            )
        elseif S <: MOI.AbstractScalarSet
            _add_constrained_variable(
                dual_problem.primal_dual_map,
                primal_model,
                S,
                params,
            )
        end
    end
    return
end
const NO_CONSTRAINT = CI{Nothing,Nothing}(0)
# Function barrier for the type unstability of `F` and `S`.
function _add_constrained_variables(
    m::PrimalDualMap,
    primal_model,
    ::Type{S},
    params,
) where {S<:MOI.AbstractVectorSet}
    cis = MOI.get(
        primal_model,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}(),
    )
    for ci in cis
        f = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
        if all(
            vi -> !haskey(m.constrained_var_idx, vi) && !(vi in params),
            f.variables,
        )
            for (i, vi) in enumerate(f.variables)
                m.constrained_var_idx[vi] = (ci, i)
            end
            # Placeholder to indicate this constraint is part of constrained variables,
            # it will be replaced later with a dual constraints
            m.constrained_var_dual[ci] = NO_CONSTRAINT
        end
    end
    return
end
function _add_constrained_variable(
    m::PrimalDualMap,
    primal_model,
    ::Type{S},
    params,
) where {S<:MOI.AbstractScalarSet}
    cis = MOI.get(
        primal_model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,S}(),
    )
    for ci in cis
        f = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
        if !haskey(m.constrained_var_idx, f) && !(f in params)
            set = MOI.get(primal_model, MOI.ConstraintSet(), ci)
            if !iszero(MOI.constant(set))
                continue
            end
            m.constrained_var_idx[f] = (ci, 1)
            # Placeholder to indicate this constraint is part of constrained variables,
            # it will be replaced later with a dual constraints
            m.constrained_var_dual[ci] = NO_CONSTRAINT
        end
    end
    return
end


text:
# Utility #
# ========================== #
function evaluate_objective(objective::JuMP.GenericAffExpr, x::AbstractVector)
    val = objective.constant
    for (var, coeff) in objective.terms
        val += coeff*x[index(var).value]
    end
    return val
end

function invalidate_cache!(stochasticprogram::StochasticProgram)
    cache = problemcache(stochasticprogram)
    delete!(cache, :evp)
    delete!(cache, :dep)
    delete!(cache, :stage_1)
    cache = solutioncache(stochasticprogram)
    delete!(cache, :solution)
    return nothing
end

function remove_scenarios!(stochasticprogram::StochasticProgram, s::Integer = 2)
    remove_scenarios!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function remove_decision_variables!(stochasticprogram::StochasticProgram, s::Integer)
    remove_decision_variables!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function remove_subproblems!(stochasticprogram::StochasticProgram, s::Integer = 2)
    remove_subproblems!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function transfer_model!(dest::StochasticProgram, src::StochasticProgram)
    empty!(dest.generator)
    merge!(dest.generator, src.generator)
    return dest
end

function copy_decision_objective!(src::JuMP.Model, dest::JuMP.Model, vars::Vector{DecisionRef})
    src_obj = objective_function(src)
    src_obj_sense = objective_sense(src)
    dest_obj_sense = objective_sense(dest)
    for var in vars
        src_var = decision_by_name(src, name(var))
        src_var === nothing && error("Cannot copy objective function. Variable $var not in src model.")
        coeff = JuMP._affine_coefficient(src_obj, src_var)
        if dest_obj_sense == src_obj_sense
            set_objective_coefficient(dest, var, coeff)
        else
            set_objective_coefficient(dest, var, -coeff)
        end
    end
    dest_obj = objective_function(dest)
    if dest_obj_sense == src_obj_sense
        set_objective_function(dest, dest_obj + constant(src_obj))
    else
        set_objective_function(dest, -dest_obj + constant(src_obj))
        set_objective_sense(dest, src_obj_sense)
    end
    return nothing
end

function supports_zero(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(zero, (Type{vartype}, ))
            !provided_def && @warn "Zero not defined for $vartype. Cannot generate zero function."
            return false
        end
    end
    return true
end

function supports_expected(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(+, (vartype, vartype))
            !provided_def && @warn "Addition not defined for $vartype. Cannot generate expectation function."
            return false
        end
        if !hasmethod(*, (Float64, vartype)) || Base.code_typed(*, (Float64, vartype))[1].second != vartype
            !provided_def && @warn "Scalar multiplication with Float64 not defined for $vartype. Cannot generate expectation function."
            return false
        end
    end
    return true
end

problemcache(stochasticprogram::StochasticProgram) = stochasticprogram.problemcache
solutioncache(stochasticprogram::StochasticProgram) = stochasticprogram.solutioncache
function get_problem(stochasticprogram::StochasticProgram, key::Symbol)
    haskey(stochasticprogram.problemcache, key)|| error("No $key in problem cache")
    return stochasticprogram.problemcache[key]
end

typename(dtype::UnionAll) = dtype.body.name.name
typename(dtype::DataType) = dtype.name.name

function add_subscript(src::AbstractString, subscript::Integer)
    return @sprintf("%s%s", src, unicode_subscript(subscript))
end
add_subscript(src::Symbol, subscript::Integer) = add_subscript(String(src), subscript)

function unicode_subscript(subscript::Integer)
    if subscript < 0
        error("$subscript is negative")
    end
    return join('₀'+d for d in reverse(digits(subscript)))
end

function extract_set(expr)
    set = NoSpecifiedConstraint()
    found = false
    quit = false
    new_expr = prewalk(expr) do x
        if @capture(x, var_Symbol in constrset_) && !found && !quit
            set = constrset
            found = true
            return :($var)
        elseif @capture(x, var_Symbol[ids__] in constrset_) && !found && !quit
            set = constrset
            found = true
            return :($(x.args[2]))
        elseif @capture(x, [ids__] in constrset_) && !found && !quit
            set = constrset
            found = true
            return :($(x.args[2]))
        elseif @capture(x, set = constrset_) && !found && !quit
            set = constrset
            found = true
            return :()
        elseif @capture(x, var_Symbol[ids__])
            # Break here to prevent indices from being filtered
            quit = true
            return x
        else
            return x
        end
    end
    return set, new_expr, found
end

function decision_variables_at_stage(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    return all_decisions(get_decisions(proxy(stochasticprogram, s), s))
end

function decision_constraints_at_stage(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    return mapreduce(vcat, MOI.get(proxy(stochasticprogram, s), MOI.ListOfConstraints())) do (F, S)
        if F <: SingleDecision || F <: AffineDecisionFunction || F <: QuadraticDecisionFunction
            return MOI.get(proxy(stochasticprogram, s), MOI.ListOfConstraintIndices{F,S}())
        end
        return CI[]
    end
end

function attach_mocks!(structure::DeterministicEquivalent)
    MOIU.attach_optimizer(structure.model)
    return nothing
end

function attach_mocks!(structure::VerticalStructure)
    MOIU.attach_optimizer(structure.first_stage)
    for s in subproblems(structure, 2)
        MOIU.attach_optimizer(s)
    end
    return nothing
end

function attach_mocks!(structure::HorizontalStructure)
    return nothing
end

function decision_index(model::MOI.ModelLike, index::MOI.VariableIndex)
    ci = CI{MOI.SingleVariable,SingleDecisionSet{Float64}}(index.value)
    if MOI.is_valid(model, ci)
        return MOI.get(model, DecisionIndex(), ci)
    end
    # Locate multiple decision set
    F = MOI.VectorOfVariables
    S = MultipleDecisionSet{Float64}
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        f = MOI.get(model, MOI.ConstraintFunction(), ci)
        i = something(findfirst(vi -> vi == index, f.variables), 0)
        if i != 0
            return MOI.get(model, DecisionIndex(), ci)
        end
    end
end

function decision_index(model::MOI.ModelLike, ci::MOI.ConstraintIndex)
    return MOI.get(model, DecisionIndex(), ci)
end

function MOI.get(src::MOIU.MockOptimizer, attr::ScenarioDependentModelAttribute)
    return MOI.get(src, attr.attr)
end

function MOI.get(src::MOIU.MockOptimizer, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(src, attr.attr, index)
end

function MOI.get(src::MOIU.MockOptimizer, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex)
    return MOI.get(src, attr.attr, ci)
end

function MOI.set(src::MOIU.MockOptimizer, attr::ScenarioDependentModelAttribute, value)
    return MOI.set(src, attr.attr, value)
end

function MOI.set(src::MOIU.MockOptimizer, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex, value)
    return MOI.set(src, attr.attr, index, value)
end

function MOI.set(src::MOIU.MockOptimizer, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex, value)
    return MOI.set(src, attr.attr, ci, value)
end
# ========================== #


text:
function readcpp(io::IO, ::Type{Vector{T}}; kwargs...) where {T}
  format = get(kwargs, :format, "v3")
  v = Vector{T}()
  if format == "v3"
    size = read(io, UInt64)
    resize!(v, size)
    for n in 1:size
      v[n] = readcpp(io, T; kwargs...)
    end
  else
    throw(ArgumentError("read Vector: format=$format not supported"))
  end
  return v
end


text:
struct DualTermEvaluator{T, TPD<:PrimalData{T}}
    pd::TPD
end
function (dte::DualTermEvaluator)(λ, ji::Tuple)
    j, i = ji
    @unpack pd = dte
    @unpack p, q, ρ, σ, x1, x = pd
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    xj = x[j]
    Ujxj = Uj - xj
    xjLj = xj - Lj
    Δ = ρ[i]*σj/4
    return λ[i]*((p[j,i] + Δ)/Ujxj + (q[j,i] + Δ)/xjLj)
end
function (dte::DualTermEvaluator)(j::Int)
    pd = dte.pd
    @unpack p0, q0, σ, x1, x = pd
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    return p0[j]/(Uj-x[j]) + q0[j]/(x[j]-Lj)
end
struct DualObjVal{T, TPD<:PrimalData}
    pd::TPD
    dte::DualTermEvaluator{T, TPD}
    x_updater::XUpdater{TPD}
end
DualObjVal(pd::TPD, x_updater::XUpdater{TPD}) where {TPD <: PrimalData} = DualObjVal(pd, DualTermEvaluator(pd), x_updater)
function (dobj::DualObjVal{T, TPD})(λ) where {T, TPD<:PrimalData{T}}
    @unpack pd, dte = dobj
    @unpack p, r, r0 = pd
    # Updates x to the Lagrangian minimizer for the input λ
    #println("Updating x")
    dobj.x_updater(λ)
    nv, nc = size(p)
    φ = r0[] + dot(λ, r)
    all(isfinite, λ) || (@show λ; error("Wait a minute, you screwed up. 0"))

    isfinite(φ) || error("Wait a minute, you screwed up. 1")
    #println("Adding objective terms")
    φ += mapreduce(dte, +, 1:nv, init=T(0))
    isfinite(φ) || error("Wait a minute, you screwed up. 2")
    #println("Adding constraint terms")
    φ += mapreduce((ji)->dte(λ, ji), +, Base.Iterators.product(1:nv, 1:nc), init=T(0))
    #@show φ, λ
    isfinite(φ) || error("Wait a minute, you screwed up. 3")
    return -φ
end

struct DualGradTermEvaluator{TPD<:PrimalData}
    pd::TPD
end
function (gte::DualGradTermEvaluator)(ji::Tuple)
    j, i = ji
    pd = gte.pd
    @unpack p, q, ρ, σ, x1, x = pd
    σj = σ[j]
    Lj, Uj = minus_plus(x1[j], σj)
    Ujxj = Uj - x[j]
    xjLj = x[j] - Lj
    Δ = ρ[i]*σj/4
    return (p[j,i] + Δ)/Ujxj + (q[j,i] + Δ)/xjLj
end
struct DualObjGrad{TPD<:PrimalData}
    pd::TPD
    gte::DualGradTermEvaluator{TPD}
    x_updater::XUpdater{TPD}
end
DualObjGrad(pd::PrimalData, x_updater::XUpdater) = DualObjGrad(pd, DualGradTermEvaluator(pd), x_updater)
function (dgrad::DualObjGrad{TPD})(∇φ::AbstractVector{T}, λ) where {T, TPD<:PrimalData{T}}
    @unpack pd, gte = dgrad
    @unpack p, r, r0 = pd
    # Updates x to the Lagrangian minimizer for the input λ
    dgrad.x_updater(λ)
    nv, nc = size(p)
    # Negate since we have a maximization problem
    map!((i)->(-r[i] - mapreduce(gte, +, Base.Iterators.product(1:nv, i:i), init=T(0))),
        ∇φ, 1:nc)
    return ∇φ
end


text:
using Random
using JuMP
using DecisionProgramming

include("min_expected_value.jl")

rng = MersenneTwister(22)

S = States([2, 3])
C = [ChanceNode(2, Node[1])]
D = [DecisionNode(1, Node[])]
V = [ValueNode(3, Node[1, 2])]
X = [Probabilities(rng, c, S) for c in C]
Y = [Consequences(rng, v, S) for v in V]

validate_influence_diagram(S, C, D, V)
sort!.((C, D, V, X, Y), by = x -> x.j)

U = DefaultPathUtility(V, Y)
U⁺ = PositivePathUtility(S, U)

k = 1
ϵ = 0.1
deviations = [
    s_I => Deviation(X[k].data[s_I..., :], ϵ) for s_I in paths(S[C[k].I_j])
]

model = Model()
z = DecisionVariables(model, S, D)
(min_ev, mpevs) = min_expected_value(deviations, model, S, C, X, U⁺, z, k)
@objective(model, Max, min_ev)

using SCIP
optimizer = optimizer_with_attributes(
    () -> SCIP.Optimizer()
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Construct minimizing path probability")
X′ = min_probabilities(mpevs, X, k)
P = DefaultPathProbability(C, X′)

@info("Reoptimize model with minimum path probability")
model2 = Model()
z2 = DecisionVariables(model2, S, D)
π_s = PathProbabilityVariables(model2, z2, S, P; hard_lower_bound=false)
ev = expected_value(model2, π_s, U⁺)
@objective(model2, Max, ev)

optimizer2 = optimizer_with_attributes(
    () -> SCIP.Optimizer()
)
set_optimizer(model2, optimizer2)
optimize!(model2)

@info("Extracting results.")
Z = DecisionStrategy(z2)

@info("Printing decision strategy:")
print_decision_strategy(S, Z)

@info("Computing utility distribution.")
udist = UtilityDistribution(S, P, U, Z)

@info("Printing utility distribution.")
print_utility_distribution(udist)

@info("Printing statistics")
print_statistics(udist)

max_stars_repo_name: jaantollander/RobustDecisionProgramming

max_issues_repo_name: jaantollander/RobustDecisionProgramming

max_forks_repo_name: jaantollander/RobustDecisionProgramming


text:
using DeltaSigma
using Test

include("calculateSNR.jl")
include("exampleHBF.jl")
include("find2dPIS.jl")
include("mapABCD.jl")
include("peakSNR.jl")
include("predictSNR.jl")
include("realizeNTF.jl")
include("scaleABCD.jl")
include("simulateDSM.jl")
include("simulateHBF.jl")
include("simulateMS.jl")
include("simulateSNR.jl")
include("stuffABCD.jl")
include("synthesizeNTF.jl")

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
module JuliaSub

#--------------------------------------------------
# Imports
#--------------------------------------------------

using MacroTools
using Multisets

#--------------------------------------------------
# Files
#--------------------------------------------------

include("utils/lib.jl")
include("lb-analysis/lib.jl")

#--------------------------------------------------
# Exports
#--------------------------------------------------

#--------------------------------------------------
# Code
#--------------------------------------------------

VERBOSE = true
DEBUG = true

# print analysis info every this number of packages
const PKGS_NUM_STEP = 100

end

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
using JLD2, StaticNets,  Utilities
using PyCall; pygui(:qt); using PyPlot; pygui(true)

 N = 10#
 dgpType = "SIN"
 Nsample = 100
 targetStatic = true
 T=250
 save_fold = "./data/sdergm_estimates/"
 @load(save_fold*"FilterTestDynamicFixedSize_N_$(N)_Nsample_$(Nsample)_scaling_DGP_" * dgpType  *
      "targetStatic_$(targetStatic)"* ".jld",
      dynFitDgp,dynDegsSam,estFitSS,indsTVnodes,storeGasPar,rmseSSandGas,filFitGas, degb, unifDeg)

T = size(dynDegsSam.s)[2]
  for s=1:Nsample
    for t=1:T
          estFitSS.s[:,t,s] =  StaticNets.identify(StaticNets.ErgmDirBin1(dynDegsSam.s[:,t,s]),
                                                estFitSS.s[:,t,s];idType = "firstZero" )
          filFitGas.s[:,t,s] =  StaticNets.identify(StaticNets.ErgmDirBin1(dynDegsSam.s[:,t,s]),
                                                filFitGas.s[:,t,s];idType = "firstZero" )
        end
      rmseSSandGas.s[:,1,s] = sqrt.(meanSq((dynFitDgp - estFitSS.s[:,:,s]).^2,2))
      rmseSSandGas.s[:,2,s] = sqrt.(meanSq((dynFitDgp - filFitGas.s[:,:,s]).^2,2))
    end
    sum(estFitSS.s[1,:,:])==0 ? () : error()
    # AvgRMSE_SS_SD = (mean(rmseSSandGas.s[:,1,:]) ,   mean(rmseSSandGas.s[:,2,:]))
    #
    #  @save(save_fold*"FilterTestDynamicFixedSize_N_$(N)_Nsample_$(Nsample)_scaling_DGP_" * dgpType  * ".jld",
    #              dynFitDgp,dynDegsSam.s,estFitSS.s,indsTVnodes,storeGasPar,rmseSSandGas.s,filFitGas.s, degb,unifDeg)
    #
    #

 close()
 figure(figsize=(9,7))
 labSize = 26
 tickSize = 20
 legSize = 20
 for n=1:Nsample
    parInd =8
    #subplot(2,1,1);
    plot(1:T,filFitGas.s[ parInd,:,n],"r")
                   plot(1:T,estFitSS.s[parInd,:,n],".b",markersize = 2)
                   plot(1:T,dynFitDgp[parInd,:],"k",linewidth=5)

    #parInd = 17
    # subplot(2,1,2);plot(1:T,filFitGas.s[ parInd,:,n],"r")
    #                plot(1:T,estFitSS.s[parInd,:,n],".b")
    #                plot(1:T,dynFitDgp[parInd,:],"k",linewidth=5)
 end
 legTex = [ "SD- Beta Model "; "Beta Model "; "DGP"]
 legend(legTex,fontsize = legSize,framealpha = 0.95)
 xlabel("Time",fontsize = labSize)
 ylabel("\$ \\overleftarrow{\\theta} \$",fontsize = labSize)

 xticks(fontsize = tickSize)
 yticks(fontsize = tickSize)
 ylim([-2.5,5.5])

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
macro odbc(func,args,vals...)
    if is_windows()
        esc(quote
            ret = ccall( ($func, odbc_dm), stdcall, SQLRETURN, $args, $(vals...))
        end)
    else
        esc(quote
            ret = ccall( ($func, odbc_dm),          SQLRETURN, $args, $(vals...))
        end)
    end
end

macro check(handle, func)
    str = string(func)
    esc(quote
        ret = $func
        if ret ∉ [SQL_SUCCESS, SQL_SUCCESS_WITH_INFO] && ODBCError($handle)
            throw(ODBCError("$str failed; return code: $ret => $(RETURN_VALUES[ret])"))
        end
        nothing
    end)
end


text:
@testset "Flux summary" begin
    model = load_model(model_paths["e_coli_core.json"])

    sol = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [change_optimizer_attribute("IPM_IterationsLimit", 200)],
    )

    fr = flux_summary(sol; keep_unbounded = true, large_flux_bound = 25)

    @test isapprox(
        fr.biomass_fluxes["BIOMASS_Ecoli_core_w_GAM"],
        0.8739215022690006;
        atol = TEST_TOLERANCE,
    )
    @test isapprox(fr.export_fluxes["EX_co2_e"], 22.80983339307183; atol = TEST_TOLERANCE)
    @test isapprox(fr.import_fluxes["EX_o2_e"], -21.799492758430517; atol = TEST_TOLERANCE)
    @test isapprox(
        fr.unbounded_fluxes["EX_h2o_e"],
        29.175827202663395;
        atol = TEST_TOLERANCE,
    )
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function boxmaptest(fname)
    Drawing(850, 850, fname)
    origin()
    background("ivory")

    # fibonacci
    tiles = boxmap([1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89], O - (400, 400), 800, 800)
    for (n, t) in enumerate(tiles)
    randomhue()

    #box(t, :fill) # box() would work, but BoundingBoxes are fun too

    bb = BoundingBox(t)
    sethue("black")
    box(bb - 2, :stroke)

    randomhue()
    box(bb - 4, :fill)

    # text labels
    sethue("white")
    # rescale text to fit better
    fontsize(boxwidth(bb) > boxheight(bb) ? boxheight(bb)/4 : boxwidth(bb)/4)
    text(string(sort([1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89], rev=true)[n]),
        midpoint(bb[1], bb[2]),
        halign=:center,
        valign=:middle)
    end

    @test finish() == true
end

fname = "boxmaptest.pdf"
boxmaptest(fname)
println("...finished test: output in $(fname)")


text:
# This file contains additional custom distibutions using the interface of Distibutions.jl

"""
    Sine()

A custom univariate distribution.
The pdf is a sine function defined between 0 and π.
This is a common prior distribution used when fitting orbits to astrometry.

The full Distributions.jl interface is not yet defined for this distribution,
but the following methods work:
pdf, logpdf, minimum, maximum, insupport, mean, var, cdf, quantile
"""
struct Sine <: ContinuousUnivariateDistribution end
export Sine

function Distributions.pdf(d::Sine, x::Real)
    if 0 < x < π
        return sin(x)/2
    else
        return 0
    end
end
function Distributions.logpdf(d::Sine, x::Real)
    if 0 < x < π
        return log(sin(x)/2)
    else
        return -Inf
    end
end
Distributions.minimum(d::Sine) = 0
Distributions.maximum(d::Sine) = π
Distributions.insupport(d::Sine, x::Real) = 0 < x < π
Distributions.mean(d::Sine) = π/2

# See https://stats.libretexts.org/Bookshelves/Probability_Theory/Probability_Mathematical_Statistics_and_Stochastic_Processes_(Siegrist)/05%3A_Special_Distributions/5.27%3A_The_Sine_Distribution
Distributions.var(d::Sine) = 1/4 - 2/pi^2
Distributions.cdf(d::Sine, x::Real)= 1/2*(1-cos(x))
Distributions.quantile(d::Sine, p::Real) = acos(1-2p)


text:
using Statistics
using Images
using FFTW
using Plots
using DSP
using ImageFiltering
using OffsetArrays

function clamp_at_boundary(M, i, j)
	return M[
		clamp(i, 1, size(M, 1)),
		clamp(j, 1, size(M, 2)),
	]
end

function convolve(M, kernel, M_index_function=clamp_at_boundary)
	height = size(kernel, 1)
	width = size(kernel, 2)

	half_height = height ÷ 2
	half_width = width ÷ 2

	new_image = similar(M)

	# (i, j) loop over the original image
	for i in 1:size(M, 1)
		for j in 1:size(M, 2)
			# (k, l) loop over the neighbouring pixels
			new_image[i, j] = sum([
						kernel[k, l] * M_index_function(M, i - k, j - l)
						for k in -half_height:-half_height + height - 1
						for l in -half_width:-half_width + width - 1
					])
		end
	end
	new_image
end

K = OffsetArray(gaussian((3,3), 0.25), -1:1, -1:1)
U = rand(1.0:100.0, 6, 6);
convolve(U, K)

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
using DataTools
using Documenter

makedocs(;
    modules=[DataTools],
    authors="Takafumi Arakaki <aka.tkf@gmail.com> and contributors",
    repo="https://github.com/JuliaFolds/DataTools.jl/blob/{commit}{path}#L{line}",
    sitename="DataTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaFolds.github.io/DataTools.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    strict = get(ENV, "CI", "false") == "true",
)

deploydocs(;
    repo="github.com/JuliaFolds/DataTools.jl",
    push_preview = true,
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libbluray_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libbluray")
JLLWrappers.@generate_main_file("libbluray", UUID("4a6b8444-0b90-5145-9aa0-b7d88b695d06"))
end  # module libbluray_jll


text:
# MAT_v5.jl
# Tools for reading MATLAB v5 files in Julia
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# MATLAB's file format documentation can be found at
# http://www.mathworks.com/help/pdf_doc/matlab/matfile_format.pdf

module MAT_v5
using CodecZlib, BufferedStreams, HDF5, SparseArrays
import Base: read, write, close

round_uint8(data) = round.(UInt8, data)
complex_array(a, b) = complex.(a, b)

mutable struct Matlabv5File <: HDF5.H5DataStore
    ios::IOStream
    swap_bytes::Bool
    varnames::Dict{String, Int64}

    Matlabv5File(ios, swap_bytes) = new(ios, swap_bytes)
end

const miINT8 = 1
const miUINT8 = 2
const miINT16 = 3
const miUINT16 = 4
const miINT32 = 5
const miUINT32 = 6
const miSINGLE = 7
const miDOUBLE = 9
const miINT64 = 12
const miUINT64 = 13
const miMATRIX = 14
const miCOMPRESSED = 15
const miUTF8 = 16
const miUTF16 = 17
const miUTF32 = 18
const mxCELL_CLASS = 1
const mxSTRUCT_CLASS = 2
const mxOBJECT_CLASS = 3
const mxCHAR_CLASS = 4
const mxSPARSE_CLASS = 5
const mxDOUBLE_CLASS = 6
const mxSINGLE_CLASS = 7
const mxINT8_CLASS = 8
const mxUINT8_CLASS = 9
const mxINT16_CLASS = 10
const mxUINT16_CLASS = 11
const mxINT32_CLASS = 12
const mxUINT32_CLASS = 13
const mxINT64_CLASS = 14
const mxUINT64_CLASS = 15
const mxFUNCTION_CLASS = 16
const mxOPAQUE_CLASS = 17
const READ_TYPES = Type[
    Int8, UInt8, Int16, UInt16, Int32, UInt32, Float32, Union{},
    Float64, Union{}, Union{}, Int64, UInt64]
const CONVERT_TYPES = Type[
    Union{}, Union{}, Union{}, Union{},
    Union{}, Float64, Float32, Int8, UInt8,
    Int16, UInt16, Int32, UInt32, Int64, UInt64]

read_bswap(f::IO, swap_bytes::Bool, ::Type{T}) where T =
    swap_bytes ? bswap(read(f, T)) : read(f, T)
function read_bswap(f::IO, swap_bytes::Bool, ::Type{T}, dim::Union{Int, Tuple{Vararg{Int}}}) where T
    d = read!(f, Array{T}(undef, dim))
    if swap_bytes
        for i = 1:length(d)
            @inbounds d[i] = bswap(d[i])
        end
    end
    d
end
function read_bswap(f::IO, swap_bytes::Bool, d::AbstractArray{T}) where T
    readbytes!(f, reinterpret(UInt8, d))
    if swap_bytes
        for i = 1:length(d)
            @inbounds d[i] = bswap(d[i])
        end
    end
    d
end

skip_padding(f::IO, nbytes::Int, hbytes::Int) = if nbytes % hbytes != 0
    skip(f, hbytes-(nbytes % hbytes))
end

# Read data type and number of bytes at the start of a data element
function read_header(f::IO, swap_bytes::Bool)
    dtype = read_bswap(f, swap_bytes, UInt32)

    if (dtype & 0xFFFF0000) != 0
        # Small Data Element Format
        (dtype & 0x0000FFFF, convert(Int, dtype >> 16), 4)
    else
        # Data Element Format
        (dtype, convert(Int, read_bswap(f, swap_bytes, UInt32)), 8)
    end
end

# Read data element as a vector of a given type
function read_element(f::IO, swap_bytes::Bool, ::Type{T}) where T
    (dtype, nbytes, hbytes) = read_header(f, swap_bytes)
    data = read_bswap(f, swap_bytes, T, Int(div(nbytes, sizeof(T))))
    skip_padding(f, nbytes, hbytes)
    data
end

# Read data element as encoded type
function read_data(f::IO, swap_bytes::Bool)
    (dtype, nbytes, hbytes) = read_header(f, swap_bytes)
    read_type = READ_TYPES[dtype]
    data = read_bswap(f, swap_bytes, read_type, Int(div(nbytes, sizeof(read_type))))
    skip_padding(f, nbytes, hbytes)
    data
end

# Read data element as encoded type with given dimensions, converting
# to another type if necessary and collapsing one-element matrices to
# scalars
function read_data(f::IO, swap_bytes::Bool, ::Type{T}, dimensions::Vector{Int32}) where T
    (dtype, nbytes, hbytes) = read_header(f, swap_bytes)
    read_type = READ_TYPES[dtype]
    if (read_type === UInt8) && (T === Bool)
        read_type = Bool
    end

    read_array = any(dimensions .!= 1)
    if sizeof(read_type)*prod(dimensions) != nbytes
        error("Invalid element length")
    end
    if read_array
        data = read_bswap(f, swap_bytes, read_type, tuple(convert(Vector{Int}, dimensions)...))
    else
        data = read_bswap(f, swap_bytes, read_type)
    end
    skip_padding(f, nbytes, hbytes)

    read_array ? convert(Array{T}, data) : convert(T, data)
end

function read_cell(f::IO, swap_bytes::Bool, dimensions::Vector{Int32})
    data = Array{Any}(undef, convert(Vector{Int}, dimensions)...)
    for i = 1:length(data)
        (ignored_name, data[i]) = read_matrix(f, swap_bytes)
    end
    data
end

function read_struct(f::IO, swap_bytes::Bool, dimensions::Vector{Int32}, is_object::Bool)
    if is_object
        class = String(read_element(f, swap_bytes, UInt8))
    end
    field_length = read_element(f, swap_bytes, Int32)[1]
    field_names = read_element(f, swap_bytes, UInt8)
    n_fields = div(length(field_names), field_length)

    # Get field names as strings
    field_name_strings = Vector{String}(undef, n_fields)
    n_el = prod(dimensions)
    for i = 1:n_fields
        sname = field_names[(i-1)*field_length+1:i*field_length]
        index = something(findfirst(iszero, sname), 0)
        field_name_strings[i] = String(index == 0 ? sname : sname[1:index-1])
    end

    data = Dict{String, Any}()
    sizehint!(data, n_fields+1)
    if is_object
        data["class"] = class
    end

    if n_el == 1
        # Read a single struct into a dict
        for field_name in field_name_strings
            data[field_name] = read_matrix(f, swap_bytes)[2]
        end
    else
        # Read multiple structs into a dict of arrays
        for field_name in field_name_strings
            data[field_name] = Array{Any}(undef, dimensions...)
        end
        for i = 1:n_el
            for field_name in field_name_strings
                data[field_name][i] = read_matrix(f, swap_bytes)[2]
            end
        end
    end

    data
end

function plusone!(A)
    for i = 1:length(A)
        @inbounds A[i] += 1
    end
    A
end

function read_sparse(f::IO, swap_bytes::Bool, dimensions::Vector{Int32}, flags::Vector{UInt32})
    local m::Int, n::Int
    if length(dimensions) == 2
        (m, n) = dimensions
    elseif length(dimensions) == 1
        m = dimensions[1]
        n = 1
    elseif length(dimensions) == 0
        m = 0
        n = 0
    else
        error("invalid dimensions encountered for sparse array")
    end

    m = isempty(dimensions) ? 0 : dimensions[1]
    n = length(dimensions) <= 1 ? 0 : dimensions[2]
    ir = plusone!(convert(Vector{Int}, read_element(f, swap_bytes, Int32)))
    jc = plusone!(convert(Vector{Int}, read_element(f, swap_bytes, Int32)))
    if (flags[1] & (1 << 9)) != 0 # logical
        # WTF. For some reason logical sparse matrices are tagged as doubles.
        pr = read_element(f, swap_bytes, Bool)
    else
        pr = read_data(f, swap_bytes)
        if (flags[1] & (1 << 11)) != 0 # complex
            pr = complex_array(pr, read_data(f, swap_bytes))
        end
    end

    SparseMatrixCSC(m, n, jc, ir, pr)
end

truncate_to_uint8(x) = x % UInt8

function read_string(f::IO, swap_bytes::Bool, dimensions::Vector{Int32})
    (dtype, nbytes, hbytes) = read_header(f, swap_bytes)
    if dtype <= 2 || dtype == 16
        # If dtype <= 2, this may give an error on non-ASCII characters, since the string
        # would be ISO-8859-1 and not UTF-8. However, MATLAB 2012b always saves strings with
        # a 2-byte encoding in v6 format, and saves UTF-8 in v7 format. Thus, this may never
        # happen in the wild.
        chars = read!(f, Vector{UInt8}(undef, nbytes))
        if dimensions[1] <= 1
            data = String(chars)
        else
            data = Vector{String}(undef, dimensions[1])
            for i = 1:dimensions[1]
                data[i] = rstrip(String(chars[i:dimensions[1]:end]))
            end
        end
    elseif dtype <= 4 || dtype == 17
        # Technically, if dtype == 3 or dtype == 4, this is ISO-8859-1 and not Unicode.
        # However, the first 256 Unicode code points are derived from ISO-8859-1, so UCS-2
        # is a superset of 2-byte ISO-8859-1.
        chars = read_bswap(f, swap_bytes, UInt16, convert(Int, div(nbytes, 2)))
        bufs = [IOBuffer() for i = 1:dimensions[1]]
        i = 1
        while i <= length(chars)
            for j = 1:dimensions[1]
                char = convert(Char, chars[i])
                if 255 < convert(UInt32, char)
                    # Newer versions of MATLAB seem to write some mongrel UTF-8...
                    char = String([truncate_to_uint8(chars[i] >> 8), truncate_to_uint8(chars[i])])[1]
                end
                write(bufs[j], char)
                i += 1
            end
        end

        if dimensions[1] == 0
            data = ""
        elseif dimensions[1] == 1
            data = String(take!(bufs[1]))
        else
            data = String[rstrip(String(take!(buf))) for buf in bufs]
        end
    else
        error("Unsupported string type")
    end
    skip_padding(f, nbytes, hbytes)
    data
end

# Read matrix data
function read_matrix(f::IO, swap_bytes::Bool)
    (dtype, nbytes) = read_header(f, swap_bytes)
    if dtype == miCOMPRESSED
        return read_matrix(ZlibDecompressorStream(IOBuffer(read!(f, Vector{UInt8}(undef, nbytes)))), swap_bytes)
    elseif dtype != miMATRIX
        error("Unexpected data type")
    elseif nbytes == 0
        # If one creates a cell array using
        #     y = cell(m, n)
        # then MATLAB will save the empty cells as zero-byte matrices. If one creates a
        # empty cells using
        #     a = {[], [], []}
        # then MATLAB does not save the empty cells as zero-byte matrices. To avoid
        # surprises, we produce an empty array in both cases.
        return ("", Matrix{Union{}}(undef, 0, 0))
    end

    flags = read_element(f, swap_bytes, UInt32)
    class = flags[1] & 0xFF

    if class == mxOPAQUE_CLASS
        s0 = read_data(f, swap_bytes)
        s1 = read_data(f, swap_bytes)
        s2 = read_data(f, swap_bytes)
        arr = read_matrix(f, swap_bytes)
        return ("__opaque__", Dict("s0"=>s0, "s1"=>s1, "s2"=>s2, "arr"=>arr))
    end

    dimensions = read_element(f, swap_bytes, Int32)
    name = String(read_element(f, swap_bytes, UInt8))

    local data
    if class == mxCELL_CLASS
        data = read_cell(f, swap_bytes, dimensions)
    elseif class == mxSTRUCT_CLASS || class == mxOBJECT_CLASS
        data = read_struct(f, swap_bytes, dimensions, class == mxOBJECT_CLASS)
    elseif class == mxSPARSE_CLASS
        data = read_sparse(f, swap_bytes, dimensions, flags)
    elseif class == mxCHAR_CLASS && length(dimensions) <= 2
        data = read_string(f, swap_bytes, dimensions)
    elseif class == mxFUNCTION_CLASS
        data = read_matrix(f, swap_bytes)
    else
        if (flags[1] & (1 << 9)) != 0 # logical
            data = read_data(f, swap_bytes, Bool, dimensions)
        else
            convert_type = CONVERT_TYPES[class]
            data = read_data(f, swap_bytes, convert_type, dimensions)
            if (flags[1] & (1 << 11)) != 0 # complex
                data = complex_array(data, read_data(f, swap_bytes, convert_type, dimensions))
            end
        end
    end

    return (name, data)
end

# Open MAT file for reading
matopen(ios::IOStream, endian_indicator::UInt16) =
    Matlabv5File(ios, endian_indicator == 0x494D)

# Read whole MAT file
function read(matfile::Matlabv5File)
    seek(matfile.ios, 128)
    vars = Dict{String, Any}()
    while !eof(matfile.ios)
        (name, data) = read_matrix(matfile.ios, matfile.swap_bytes)
        vars[name] = data
    end
    vars
end
# Read only variable names from an HDF5 file
function getvarnames(matfile::Matlabv5File)
    if !isdefined(matfile, :varnames)
        seek(matfile.ios, 128)
        matfile.varnames = varnames = Dict{String, Int64}()
        while !eof(matfile.ios)
            offset = position(matfile.ios)
            (dtype, nbytes, hbytes) = read_header(matfile.ios, matfile.swap_bytes)
            if dtype == miCOMPRESSED
                f = ZlibDecompressorStream(matfile.ios)
                read_header(f, matfile.swap_bytes)
            elseif dtype == miMATRIX
                f = matfile.ios
            else
                error("Unexpected data type")
            end

            read_element(f, matfile.swap_bytes, UInt32)
            read_element(f, matfile.swap_bytes, Int32)
            varnames[String(read_element(f, matfile.swap_bytes, UInt8))] = offset

            seek(matfile.ios, offset+nbytes+hbytes)
        end
    end
    matfile.varnames
end

Base.haskey(matfile::Matlabv5File, varname::String) =
    haskey(getvarnames(matfile), varname)
Base.keys(matfile::Matlabv5File) =
    keys(getvarnames(matfile))

# Read a variable from a MAT file
function read(matfile::Matlabv5File, varname::String)
    varnames = getvarnames(matfile)
    if !haskey(varnames, varname)
        error("no variable $varname in file")
    end
    seek(matfile.ios, varnames[varname])
    (name, data) = read_matrix(matfile.ios, matfile.swap_bytes)
    data
end

# Complain about writing to a MAT file
function write(parent::Matlabv5File, name::String, s)
    error("Writing to a MATLAB v5 file is not currently supported. Create a new file instead.")
end

# Close MAT file
close(matfile::Matlabv5File) = close(matfile.ios)
end


text:
import DelimitedFiles
import Statistics
import NPZ
import RestrictedBoltzmannMachines as RBMs

using Test: @test, @testset
using Statistics: mean, std

#=
Comparisons against https://github.com/jertubiana/PGM.

I computed the energy, free_energy, and so on, using the PGM
code. So here we compare if using the same RBM I obtain the same
values of energies, etc. with my code.

The only difference is in the definition of some θ's (for Gaussian, ReLU, and dReLU layers).
I wanted to be consistent with the definition of the external fields, which appear with
a minus sign in the energy for discrete units. Therefore in the following tests pay
attention to these minus signs.
=#

@testset "Binary / Binary RBM" begin
    g = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_g_v.txt"))
    θ = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_g_h.txt"))
    w = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_w.txt")'

    v = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_v.txt")'
    h = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_h.txt")'
    E = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_E.txt"))
    F = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_bin_F.txt"))

    rbm = RBMs.RBM(RBMs.Binary(g), RBMs.Binary(θ), w)
    @test RBMs.energy(rbm, v, h) ≈ E rtol=1e-5
    @test RBMs.free_energy(rbm, v) ≈ F rtol=1e-5
end

@testset "Binary / Gaussian RBM" begin
    #= Note the minus sign in front of θ =#
    θ = -vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_theta_h.txt"))
    γ =  vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_gamma_h.txt"))
    g = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_g_v.txt"))
    w = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_w.txt")'

    v = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_v.txt")'
    h = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_h.txt")'
    Ev = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_Ev.txt"))
    Eh = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_Eh.txt"))
    E = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_E.txt"))
    F = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_gauss_F.txt"))

    rbm = RBMs.RBM(RBMs.Binary(g), RBMs.Gaussian(θ, γ), w)
    @test RBMs.energy(rbm.visible, v) ≈ Ev rtol=1e-5
    @test RBMs.energy(rbm.hidden,  h) ≈ Eh rtol=1e-5
    @test RBMs.energy(rbm, v, h) ≈ E rtol=1e-5
    @test RBMs.free_energy(rbm, v) ≈ F rtol=1e-5
end

@testset "Binary / dReLU RBM" begin
    #= Note the minus sign in front of θp. For Gaussian, ReLU, and dReLU layers
    (as well as for the reparameterized pReLU), my θ's are minus those of Jerome.
    I did this because I wanted to be consistent with the external field of Binary,
    Spin and Potts layer, for which the field appears with this sign. =#
    θp = -vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_thetap_h.txt"))
    θn =  vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_thetam_h.txt"))
    γp = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_gammap_h.txt"))
    γn = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_gammam_h.txt"))
    g = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_g_v.txt"))
    w = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_w.txt")'

    v = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_v.txt")'
    h = DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_h.txt")'
    Ev = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_Ev.txt"))
    Eh = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_Eh.txt"))
    E = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_E.txt"))
    F = vec(DelimitedFiles.readdlm("compare_to_pgm/RBM_bin_dReLU_F.txt"))

    rbm = RBMs.RBM(RBMs.Binary(g), RBMs.dReLU(θp, θn, γp, γn), w)
    @test RBMs.energy(rbm.visible, v) ≈ Ev rtol=1e-5
    @test RBMs.energy(rbm.hidden,  h) ≈ Eh rtol=1e-5
    @test RBMs.energy(rbm, v, h) ≈ E rtol=1e-5
    @test RBMs.free_energy(rbm, v) ≈ F rtol=1e-3
end

@testset "Binary pseudolikelihood" begin
    w = Matrix(DelimitedFiles.readdlm("compare_to_pgm/PL/RBM_Bernoulli_weights.txt")')
    visible_g = vec(DelimitedFiles.readdlm("compare_to_pgm/PL/RBM_Bernoulli_visible_fields.txt"))
    hidden_g = vec(DelimitedFiles.readdlm("compare_to_pgm/PL/RBM_Bernoulli_hidden_fields.txt"))
    v = BitMatrix(DelimitedFiles.readdlm("compare_to_pgm/PL/RBM_Bernoulli_data.txt", Bool)')
    pl = vec(DelimitedFiles.readdlm("compare_to_pgm/PL/RBM_Bernoulli_PL.txt"))

    rbm = RBMs.RBM(RBMs.Binary(visible_g), RBMs.Binary(hidden_g), w)
    pl_rbm = RBMs.log_pseudolikelihood(rbm, v)

    @test mean(pl_rbm) ≈ mean(pl) rtol=0.05
    @test std(pl_rbm)  ≈ std(pl)  rtol=0.1
end

@testset "Potts pseudolikelihood" begin
    # https://numpy.org/doc/stable/reference/generated/numpy.savez_compressed.html
    d = NPZ.npzread("compare_to_pgm/PL/RBM_Potts_PL.npz")
    rbm = RBMs.RBM(
        RBMs.Potts(Matrix(d["g"]')),
        RBMs.dReLU(-d["theta_plus"], d["theta_minus"], d["gamma_plus"], d["gamma_minus"]),
        permutedims(d["weights"], (3,2,1))
    )
    data = RBMs.onehot_encode(d["data"]')
    pl = RBMs.log_pseudolikelihood(rbm, data)

    # TODO: #7 see if we can decrease rtol
    @test mean(pl) ≈ mean(d["PL"]) rtol=0.1
    @test std(pl) ≈ std(d["PL"]) rtol=0.1
end


text:
using MLDatasets, JLD

function sigmoid(x)
    return 1 / (1 + exp(-x))
end

function softmax(a)
    c = maximum(a)
    exp_a = exp.(a .- c)  # オーバーフロー対策
    sum_exp_a = sum(exp_a)
    y = exp_a ./ sum_exp_a

    return y
end

function get_data()
    # x_train, t_train = MNIST.traindata()
    x_test, t_test = MNIST.testdata()
    x_test = reshape(x_test, (784, :))
    return x_test, t_test
end

function init_network()
    network = load("sample_weight.jld")
    return network
end

function predict(network, x)
    W1, W2, W3 = network["W1"], network["W2"], network["W3"]
    b1, b2, b3 = network["b1"], network["b2"], network["b3"]

    a1 = W1 * x .+ b1
    z1 = sigmoid.(a1)
    a2 = W2 * z1 .+ b2
    z2 = sigmoid.(a2)
    a3 = W3 * z2 .+ b3
    y = softmax(a3)

    return y
end

function main()
    x, t = get_data()
    network = init_network()

    accuracy_cnt = 0
    for i in 1:length(t)
        y = predict(network, x[:,i])
        p = argmax(y) - 1
        if p == t[i]
            accuracy_cnt += 1
        end
    end
    accuracy = accuracy_cnt / length(t)
    println("Accuracy: $accuracy")  # 0.9352
end

main()