text:
# # Vlasov-Poisson
#
#md # [notebook](@__NBVIEWER_ROOT_URL__notebooks/vlasov-poisson.ipynb)
#
# We consider the dimensionless Vlasov-Poisson equation for one species
# with a neutralizing background.
#
# ```math
#  \frac{∂f}{∂t}+ v⋅∇_x f + E(t,x) ⋅ ∇_v f = 0, \\
#  - Δϕ = 1 - ρ, E = - ∇ ϕ \\
#  ρ(t,x)  =  ∫ f(t,x,v) dv
# ```
#
#  - [Vlasov Equation - Wikipedia](https://en.wikipedia.org/wiki/Vlasov_equation)
#  - [Landau damping - Wikipedia](https://en.wikipedia.org/wiki/Landau_damping)
#

using Plots, LinearAlgebra
import Splittings:UniformMesh, BSpline
import Splittings:@Strang
import Splittings

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

function push_t!(f, mesh1, v, n2, dt)
    Splittings.advection!(f, mesh1, v, n2, dt, BSpline(5))
end

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

function push_v!(f, fᵗ, mesh1, mesh2, nrj, dt)
    rho = Splittings.compute_rho(mesh2, f)
    e   = Splittings.compute_e(mesh1, rho)
    push!(nrj, 0.5*log(sum(e.*e)*mesh1.step))
    transpose!(fᵗ, f)
    Splittings.advection!(fᵗ, mesh2, e, mesh1.length, dt, BSpline(5))
    transpose!(f, fᵗ)
end

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

function landau(tf, nt)

  n1, n2 = 32, 64
  x1min, x1max = 0.0, 4π
  x2min, x2max = -6., 6.
  mesh1 = UniformMesh(x1min, x1max, n1; endpoint=false)
  mesh2 = UniformMesh(x2min, x2max, n2; endpoint=false)
  x = mesh1.points
  v = mesh2.points

  ϵ, kx = 0.001, 0.5
  f = zeros(Complex{Float64},(n1,n2))
  f .= (1.0.+ϵ*cos.(kx*x))/sqrt(2π) * transpose(exp.(-0.5*v.^2))
  fᵗ = zeros(Complex{Float64},(n2,n1))

  dt = tf / nt

  nrj = Float64[]

  for it in 1:nt
      @Strang( push_t!(f, mesh1, v, n2, dt),
               push_v!(f, fᵗ, mesh1, mesh2, nrj, dt))
  end

  nrj

end

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

nt = 500
tf = 50.0
t  = range(0.0, stop=tf, length=nt)
@time nrj = landau(tf, nt)
plot( t, nrj)
plot!(t, -0.1533*t.-5.50)
savefig("landau-plot.png"); nothing # hide

using Test                       #src
@testset "Vlasov-Poisson" begin  #src
@test length(nrj) > 0            #src
end                              #src

#md # ![](landau-plot.png)


text:
abstract type AbstractSecondaryConnectivityModel end

@inline (mdl::AbstractSecondaryConnectivityModel)(n_vertices, n_pins, n_nets, n_local_nets, k) = mdl(n_vertices, n_pins, n_nets, n_local_nets)

struct AffineSecondaryConnectivityModel{Tv} <: AbstractSecondaryConnectivityModel
    α::Tv
    β_vertex::Tv
    β_pin::Tv
    β_local_net::Tv
    β_remote_net::Tv
end

function AffineSecondaryConnectivityModel(; α = false, β_vertex = false, β_pin = false, β_local_net = false, β_remote_net = false)
    AffineSecondaryConnectivityModel(promote(α, β_vertex, β_pin, β_local_net, β_remote_net)...)
end

@inline cost_type(::Type{AffineSecondaryConnectivityModel{Tv}}) where {Tv} = Tv

(mdl::AffineSecondaryConnectivityModel)(n_vertices, n_pins, n_local_nets, n_remote_nets, k) = mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin + n_local_nets * mdl.β_local_net + n_remote_nets * mdl.β_remote_net

function bound_stripe(A::SparseMatrixCSC, K, Π, mdl::AffineSecondaryConnectivityModel)
    @assert mdl.β_vertex >= 0
    @assert mdl.β_pin >= 0
    @assert mdl.β_local_net >= 0
    @assert mdl.β_remote_net >= 0
    adj_A = adjointpattern(A)
    c_hi = bottleneck_value(adj_A, Π, AffineConnectivityModel(mdl.α, mdl.β_vertex, mdl.β_pin, mdl.β_remote_net))
    c_lo = bottleneck_value(adj_A, Π, AffineWorkModel(mdl.α, mdl.β_vertex, mdl.β_pin))
    return (c_lo, c_hi)
end

struct SecondaryConnectivityOracle{Ti, Lcc, Mdl} <: AbstractOracleCost{Mdl}
    Π::SplitPartition{Ti}
    pos::Vector{Ti}
    πos::Vector{Ti}
    lcc::Lcc
    mdl::Mdl
end

oracle_model(ocl::SecondaryConnectivityOracle) = ocl.mdl

function bound_stripe(A::SparseMatrixCSC, K, Π::SplitPartition, ocl::SecondaryConnectivityOracle{<:Any, <:Any, <:AffineSecondaryConnectivityModel})
    @inbounds begin
        c_lo = 0
        c_hi = 0
        (m, n) = size(A)
        mdl = oracle_model(ocl)
        @assert mdl.β_vertex >= 0
        @assert mdl.β_pin >= 0
        @assert mdl.β_local_net >= 0
        @assert mdl.β_remote_net >= 0
        for k = 1:K
            n_vertices = Π.spl[k + 1] - Π.spl[k]
            n_pins = ocl.pos[Π.spl[k + 1]] - ocl.pos[Π.spl[k]]
            n_remote_nets = ocl.πos[k + 1] - ocl.πos[k]
            c_lo = max(c_lo, mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin)
            c_hi = max(c_hi, mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin + n_remote_nets * mdl.β_remote_net)
        end
        return (c_lo, c_hi)
    end
end

function oracle_stripe(hint::AbstractHint, mdl::AbstractSecondaryConnectivityModel, A::SparseMatrixCSC, Π::SplitPartition; adj_pos=nothing, adj_A=nothing, net=nothing, kwargs...)
    @inbounds begin
        m, n = size(A)

        if adj_A === nothing
            adj_A = adjointpattern(A)
        end

        if adj_pos === nothing
            adj_pos = adj_A.colptr
        end

        (n, K, n′, πos, prm), _ = partwise(A, convert(MapPartition, Π))
        lcc = partwisecount!(hint, n, K, n′, πos, prm, VertexCount(); kwargs...)

        return SecondaryConnectivityOracle(Π, adj_pos, πos, lcc, mdl)
    end
end

@inline function (cst::SecondaryConnectivityOracle{Ti, Mdl})(i::Ti, i′::Ti, k) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[cst.Π.spl[k + 1]] - cst.pos[cst.Π.spl[k]]
        d = cst.πos[k + 1] - cst.πos[k]
        l = cst.lcc(i, i′, k)
        return cst.mdl(cst.Π.spl[k + 1] - cst.Π.spl[k], w, l, d - l, k)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k) where {Ti, Mdl, Ocl <: SecondaryConnectivityOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = destep(_k)
        w = cst.pos[cst.Π.spl[k + 1]] - cst.pos[cst.Π.spl[k]]
        d = cst.πos[k + 1] - cst.πos[k]
        l = Step(cst.lcc)(_j, _j′, _k)
        return cst.mdl(cst.Π.spl[k + 1] - cst.Π.spl[k], w, l, d - l, k)
    end
end

function compute_objective(g::G, A, Π, Φ::SplitPartition, mdl::AffineSecondaryConnectivityModel) where {G}
    return compute_objective(g, adjointpattern(A), Φ, Π, AffinePrimaryConnectivityModel(mdl.α, mdl.β_vertex, mdl.β_pin, mdl.β_local_net, mdl.β_remote_net))
end
function compute_objective(g::G, A, Π, Φ, mdl::AffineSecondaryConnectivityModel) where {G}
    return compute_objective(g, adjointpattern(A), Φ, Π, AffinePrimaryConnectivityModel(mdl.α, mdl.β_vertex, mdl.β_pin, mdl.β_local_net, mdl.β_remote_net))
end

ext: jl


text:
using TriangleMesh
using Plots
using Triplot

poly = polygon_regular(5)
poly = polygon_Lshape()
mesh = create_mesh(poly, set_area_max=true, quality_meshing=true)

function plotTriMesh(mesh)
    x,y = mesh.point[1,:], mesh.point[2,:]

    plot() # init plot
    for vertex_ids in eachcol(mesh.cell)
        ids = vcat(vertex_ids, vertex_ids[1])
        plot!(x[ids],y[ids],linecolor=:black)
    end
    display(plot!(legend=false,ratio=1))
end
plotTriMesh(mesh)

function to_rgba(x::UInt32)
    a = ((x & 0xff000000)>>24)/255
    b = ((x & 0x00ff0000)>>16)/255
    g = ((x & 0x0000ff00)>>8)/255
    r = (x & 0x000000ff)/255
    RGBA(r, g, b, a)
end

# plotting
x,y = mesh.point[1,:], mesh.point[2,:]
z,t = (@. sin(pi*x)*sin(pi*y)), mesh.cell
img = to_rgba.(Triplot.rasterize(x,y,z,t)')
plot([extrema(x)...], [extrema(y)...], img, yflip=false)


text:
import PathWeightSampling
using PathWeightSampling: SMCEstimate, name, propagate, DirectMCEstimate, marginal_configuration, ConditionalEnsemble, MarginalEnsemble, energy_difference, trajectory_energy, generate_configuration, log_marginal, logpdf, simulate
using Test
using Statistics

system = PathWeightSampling.chemotaxis_system(dtimes = 0:0.05:2.0)

@test system.dist.aggregator.update_map == [0, 0, 0, 0, 1, 2]

conf = PathWeightSampling.generate_full_configuration(system)
mconf = marginal_configuration(conf)
cens = ConditionalEnsemble(system)
mens = MarginalEnsemble(system)

@test energy_difference(conf, cens) == energy_difference(mconf, mens)

merged = PathWeightSampling.merge_trajectories(conf.s_traj, conf.r_traj, conf.x_traj)
for (i, T) in enumerate(system.dtimes)
    @test energy_difference(conf, cens)[i] ≈ -trajectory_energy(cens.dist, merged, tspan = (0.0, T))
    @test energy_difference(conf, cens)[i] ≈ -trajectory_energy(mens.dist, merged, tspan = (0.0, T))

    if i >= 2
        E(i) = energy_difference(conf, cens)[i]
        T(i) = system.dtimes[i]
        dE = E(i - 1) - E(i)
        @test dE ≈ trajectory_energy(mens.dist, merged, tspan = (T(i - 1), T(i)))
        @test dE ≈ trajectory_energy(cens.dist, merged, tspan = (T(i - 1), T(i)))
    end
end

algs = [SMCEstimate(256), DirectMCEstimate(256)]

result = Dict(map(algs) do alg
    mi = map(1:20) do i
        @time begin
            cresult = simulate(alg, conf, cens)
            mresult = simulate(alg, mconf, mens)
            log_marginal(cresult) - log_marginal(mresult)
        end
    end
    name(alg) => hcat(mi...)
end)

# Test that all results lie within 2 standard deviations of each other.
# This should find coarse implementation errors for specific estimates.
for (name, mi) in result
    lower_bound = mean(mi, dims = 2) - 2 * std(mi, dims = 2)
    for (name2, mi2) in result
        upper_bound = mean(mi2, dims = 2) + 2 * std(mi2, dims = 2)
        @test all(lower_bound .<= upper_bound)
    end
end

# # Some plotting routines are commented out below to see what goes wrong.
# # Use them to see what's wrong when the tests above fail.
#--
# using Plots
# plot()
# colormap = Dict("SMC" => 1, "Direct MC" => 2)
# for (name, mi) in result
#     plot!(system.dtimes, mean(mi, dims=2), ribbon=3*std(mi, dims=2), color=colormap[name], label=name)
#     plot!(system.dtimes, mi, color=colormap[name], linewidth=0.3, label="")
# end
# plot!()
# --


text:
# isdefined(Base, :__precompile__) && __precompile__()

module DynamoDB

using AWS
using Requests
import JSON

include("crypto.jl")

# TODO -- kill this file when AWS.jl adds support for specified payloads
include("dynamo_requests.jl")

include("dynamo_json.jl")
include("dynamo_dsl.jl")

include("dynamo.jl")
include("dynamo_row_ops.jl")

# package code goes here

# core types for modeling dynamo tables/indexes
export DynamoTable, DynamoLocalIndex, DynamoGlobalIndex,
       dynamo_table, dynamo_local_index, dynamo_global_index,
       random_key

# core operations
export get_item, put_item, update_item, delete_item

# high-read operations (easy to go over your read budget using them)
export query, scan

# batched opperations (these have gotchas... maybe read the docs before using)
export batch_get_item, batch_write_item, batch_put_item, batch_delete_item

export attribute, attr, getindex
export no_conditions, size, begins_with, contains, is_string, is_string_set, is_real, is_real_set,
       is_binary, is_binary_set, is_bool, is_null, is_list, is_map, is_document, exists, not_exists,
       <, >, <=, >=, eq, !=, and, or, !, not, between
export get_or_else, append_to_list, assign, add_to_set, remove_from_set, delete

end # module


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415

# %%
using ForwardDiff
using PreallocationTools

function f(x, tmp)
    tmp = get_tmp(tmp, x)
    @. tmp = x + 1
    sum(tmp)
end

x = [1.0, 2.0, 3.0]
tmp = dualcache(similar(x))
@show f(x, tmp)
ForwardDiff.gradient(x -> f(x, tmp), x)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
struct Label{Name}; end

@pure Label(name::Symbol) = Label{name}()

macro l_str(ex)
    return Label(Symbol(ex::String))
end

@pure (==)(::Label{Name}, ::Label{Name}) where {Name} = true
@pure (==)(::Label{Name1}, ::Label{Name2}) where {Name1, Name2} = false

show(io::IO, ::Label{Name}) where {Name} = print(io, "l\"$(Name)\"")
print(io::IO, ::Label{Name}) where {Name} = print(io, Name)

@pure function label_to_index(names::Tuple{Vararg{Symbol}}, ::Label{name}) where name
    for i = 1:length(names)
        if names[i] == name
            return i
        end
    end
    error("Could not find $name in $names")
end

@inline getnames(x::Tuple{Vararg{Pair{<:Label}}}) = _getnames(x...)
@inline _getnames(::Pair{Label{name}}, x...) where {name} = (name, _getnames(x...)...)
@inline _getnames() = ()

@inline getvalues(x::Tuple{Vararg{Pair}}) = _getvalues(x...)
@inline _getvalues(x::Pair, y...)  = (last(x), _getvalues(y...)...)
@inline _getvalues() = ()

@pure function Base.issubset(x::Tuple{Vararg{Label}}, y::Tuple{Vararg{Label}})
    for l ∈ x
        if l ∉ y
            return false
        end
    end
    return true
end

@pure function Base.union(x::Tuple{Vararg{Label}}, y::Tuple{Vararg{Label}})
    out = x
    for l ∈ y
        if l ∉ x
            out = (out..., l)
        end
    end
    return out
end

@pure function Base.intersect(x::Tuple{Vararg{Label}}, y::Tuple{Vararg{Label}})
    out = ()
    for l ∈ x
        if l ∈ y
            out = (out..., l)
        end
    end
    return out
end


text:
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using ParametricOptInterface
using MathOptInterface

const POI = ParametricOptInterface
const MOI = MathOptInterface

using ParameterJuMP
using JuMP

using GLPK
const OPTIMIZER = GLPK.Optimizer;
using TimerOutputs

using LinearAlgebra
using Random

const N_Candidates = 200
const N_Observations = 2000
const N_Nodes = 200

const Observations = 1:N_Observations
const Candidates = 1:N_Candidates
const Nodes = 1:N_Nodes;

#' Initialize a random number generator to keep results deterministic

rng = Random.MersenneTwister(123);

#' Building regressors (explanatory) sinusoids

const X = zeros(N_Candidates, N_Observations)
const time = [obs / N_Observations * 1 for obs in Observations]
for obs in Observations, cand in Candidates
    t = time[obs]
    f = cand
    X[cand, obs] = sin(2 * pi * f * t)
end

#' Define coefficients

β = zeros(N_Candidates)
for i in Candidates
    if rand(rng) <= (1 - i / N_Candidates)^2 && i <= 100
        β[i] = 4 * rand(rng) / i
    end
end
# println("First coefs: $(β[1:min(10, N_Candidates)])")

const y = X' * β .+ 0.1 * randn(rng, N_Observations)

function full_model_regression()
    time_build = @elapsed begin # measure time to create a model

        # initialize a optimization model
        full_model = Model(OPTIMIZER)

        # create optimization variables of the problem
        @variables(full_model, begin
            ɛ_up[Observations] >= 0
            ɛ_dw[Observations] >= 0
            β[1:N_Candidates]
            # 0 <= β[Candidates] <= 8
        end)

        # define constraints of the model
        @constraints(
            full_model,
            begin
                ɛ_up_ctr[i in Observations],
                ɛ_up[i] >= +sum(X[j, i] * β[j] for j in Candidates) - y[i]
                ɛ_dw_ctr[i in Observations],
                ɛ_dw[i] >= -sum(X[j, i] * β[j] for j in Candidates) + y[i]
            end
        )

        # construct the objective function to be minimized
        @objective(
            full_model,
            Min,
            sum(ɛ_up[i] + ɛ_dw[i] for i in Observations)
        )
    end

    # solve the problem
    time_solve = @elapsed optimize!(full_model)

    println(
        "First coefficients in solution: $(value.(β)[1:min(10, N_Candidates)])",
    )
    println("Objective value: $(objective_value(full_model))")
    println("Time in solve: $time_solve")
    println("Time in build: $time_build")

    return nothing
end

function ObsSet(K)
    obs_per_block = div(N_Observations, N_Nodes)
    return (1+(K-1)*obs_per_block):(K*obs_per_block)
end

function slave_model(PARAM, K)

    # initialize the JuMP model
    slave = if PARAM == 0
        # special constructor exported by ParameterJuMP
        # to add the functionality to the model
        Model(OPTIMIZER)
    elseif PARAM == 1
        # POI constructor
        direct_model(POI.Optimizer(OPTIMIZER()))
        # Model(() -> POI.ParametricOptimizer(OPTIMIZER()))
    else
        # regular JuMP constructor
        Model(OPTIMIZER)
    end

    # Define local optimization variables for norm-1 error
    @variables(slave, begin
        ɛ_up[ObsSet(K)] >= 0
        ɛ_dw[ObsSet(K)] >= 0
    end)

    # create the regression coefficient representation
    if PARAM == 0
        # here is the main constructor of the Parameter JuMP packages
        # It will create model *parameters* instead of variables
        # Variables are added to the optimization model, while parameters
        # are not. Parameters are merged with LP problem constants and do not
        # increase the model dimensions.
        @variable(slave, β[i = 1:N_Candidates] == 0, Param())
    elseif PARAM == 1
        # Create parameters
        β = [@variable(slave, set = POI.Parameter(0)) for i in 1:N_Candidates]
        for (i, βi) in enumerate(β)
            set_name(βi, "β[$i]")
        end
    else
        # Create fixed variables
        @variables(slave, begin
            β[Candidates]
            β_fixed[1:N_Candidates] == 0
        end)
        @constraint(slave, β_fix[i in Candidates], β[i] == β_fixed[i])
    end

    # create local constraints
    # Note that *parameter* algebra is implemented just like variables
    # algebra. We can multiply parameters by constants, add parameters,
    # sum parameters and variables and so on.
    @constraints(
        slave,
        begin
            ɛ_up_ctr[i in ObsSet(K)],
            ɛ_up[i] >= +sum(X[j, i] * β[j] for j in Candidates) - y[i]
            ɛ_dw_ctr[i in ObsSet(K)],
            ɛ_dw[i] >= -sum(X[j, i] * β[j] for j in Candidates) + y[i]
        end
    )

    # create local objective function
    @objective(slave, Min, sum(ɛ_up[i] + ɛ_dw[i] for i in ObsSet(K)))

    # return the correct group of parameters
    return (slave, β)
end

function master_model(PARAM)
    master = Model(OPTIMIZER)
    @variables(master, begin
        ɛ[Nodes] >= 0
        β[1:N_Candidates]
    end)
    @objective(master, Min, sum(ɛ[i] for i in Nodes))
    sol = zeros(N_Candidates)
    return (master, ɛ, β, sol)
end

function master_solve(PARAM, master_model)
    model = master_model[1]
    β = master_model[3]
    optimize!(model)
    return (value.(β), objective_value(model))
end

#' ### Query dual of parameters in POI
function dual_of_parameter(vref::VariableRef)
    if mode(vref.model) == JuMP.AUTOMATIC
        backend_optimizer_model = JuMP.backend(vref.model).optimizer.model
        poi_optimizer = backend_optimizer_model.optimizer
        optimizer_map = backend_optimizer_model.model_to_optimizer_map
        c_index = optimizer_map[vref.index].value
        cp = MOI.ConstraintIndex{MOI.VariableIndex,POI.Parameter}(c_index)
        return MOI.get(poi_optimizer, MOI.ConstraintDual(), cp)
    elseif mode(vref.model) == JuMP.DIRECT
        poi_optimizer = JuMP.backend(vref.model)
        cp = MOI.ConstraintIndex{MOI.VariableIndex,POI.Parameter}(
            vref.index.value,
        )
        return MOI.get(poi_optimizer, MOI.ConstraintDual(), cp)
    end
    return throw(ErrorException("Only works for JuMP.Model"))
end

function slave_solve(PARAM, model, master_solution)
    β0 = master_solution[1]
    slave = model[1]

    # The first step is to fix the values given by the master problem
    @timeit "fix" if PARAM == 0
        # ParameterJuMP: *parameters* can be set to new values and the optimization
        # model will be automatically updated
        β = model[2]
        set_value.(β, β0)
    elseif PARAM == 1
        # POI: assigning new values to *parameters* and the optimization
        # model will be automatically updated
        β = model[2]
        MOI.set.(slave, POI.ParameterValue(), β, β0)
    else
        # JuMP: it is also possible to fix variables to new values
        β_fixed = slave[:β_fixed]
        fix.(β_fixed, β0)
    end

    # here the slave problem is solved
    @timeit "opt" optimize!(slave)

    # query dual variables, which are sensitivities
    # They represent the subgradient (almost a derivative)
    # of the objective function for infinitesimal variations
    # of the constants in the linear constraints
    @timeit "dual" if PARAM == 0
        # ParameterJuMP: we can query dual values of *parameters*
        π = dual.(β)
    elseif PARAM == 1
        # POI: we can query dual values of *parameters*
        π = -dual_of_parameter.(β)
    else
        # or, in pure JuMP, we query the duals form
        # constraints that fix the values of our regression
        # coefficients
        π = dual.(slave[:β_fix])
    end

    # π2 = shadow_price.(β_fix)
    # @show sum(π .- π2)
    obj = objective_value(slave)
    rhs = obj - dot(π, β0)
    return (rhs, π, obj)
end

function master_add_cut(PARAM, master_model, cut_info, node)
    master = master_model[1]
    ɛ = master_model[2]
    β = master_model[3]

    rhs = cut_info[1]
    π = cut_info[2]

    @constraint(master, ɛ[node] >= sum(π[j] * β[j] for j in Candidates) + rhs)
end
function decomposed_model(PARAM; print_timer_outputs::Bool = true)
    reset_timer!() # reset timer fo comparision
    time_init = @elapsed @timeit "Init" begin
        # println("Initialize decomposed model")

        # Create the master problem with no cuts
        # println("Build master problem")
        @timeit "Master" master = master_model(PARAM)

        # initialize solution for the regression coefficients in zero
        # println("Build initial solution")
        @timeit "Sol" solution = (zeros(N_Candidates), Inf)
        best_sol = deepcopy(solution)

        # Create the slave problems
        # println("Build slave problems")
        @timeit "Slaves" slaves =
            [slave_model(PARAM, i) for i in Candidates]

        # Save initial version of the slave problems and create
        # the first set of cuts
        # println("Build initial cuts")
        @timeit "Cuts" cuts =
            [slave_solve(PARAM, slaves[i], solution) for i in Candidates]
    end

    UB = +Inf
    LB = -Inf

    # println("Initialize Iterative step")
    time_loop = @elapsed @timeit "Loop" for k in 1:80

        # Add cuts generated from each slave problem to the master problem
        @timeit "add cuts" for i in Candidates
            master_add_cut(PARAM, master, cuts[i], i)
        end

        # Solve the master problem with the new set of cuts
        # Obtain new solution candidate for the regression coefficients
        @timeit "solve master" solution = master_solve(PARAM, master)

        # Pass the new candidate solution to each of the slave problems
        # Solve the slave problems and obtain cutting planes
        # @show solution[2]
        @timeit "solve nodes" for i in Candidates
            cuts[i] = slave_solve(PARAM, slaves[i], solution)
        end

        LB = solution[2]
        new_UB = sum(cuts[i][3] for i in Candidates)
        if new_UB <= UB
            best_sol = deepcopy(solution)
        end
        UB = min(UB, new_UB)
        # println("Iter = $k, LB = $LB, UB = $UB")

        if abs(UB - LB) / (abs(UB) + abs(LB)) < 0.05
            # println("Converged!")
            break
        end
    end
    # println(
    #     "First coefficients in solution: $(solution[1][1:min(10, N_Candidates)])",
    # )
    # println("Objective value: $(solution[2])")
    # println("Time in loop: $time_loop")
    # println("Time in init: $time_init")

    print_timer_outputs && print_timer()

    return best_sol[1]
end

# println("ParameterJuMP")
# GC.gc()
# β1 = decomposed_model(0; print_timer_outputs = false);
# GC.gc()
# β1 = decomposed_model(0);

println("POI")
GC.gc()
β2 = decomposed_model(1; print_timer_outputs = false);
GC.gc()
β2 = decomposed_model(1);

# println("Pure JuMP")
# GC.gc()
# β3 = decomposed_model(2; print_timer_outputs = false);
# GC.gc()
# β3 = decomposed_model(2);


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

include("BaseOptimizationTest.jl")
include("SingularOptTest.jl")
include("CircOptTest.jl")


text:
Topic = Tuple{String, QOS}
Base.show(io::IO, x::Topic) = print(io, "(", join(x, ", "), ")")

struct Subscribe <: HasId
    header::UInt8
    id::UInt16
    topics::Array{Topic}
end
Subscribe(topics::Array{Topic}) = Subscribe(convert(UInt8, SUBSCRIBE) | 0x02, 0x0000, topics)
Subscribe(packet::Subscribe, id::UInt16) = Subscribe(packet.header, id, packet.topics)

function write(s::IO, packet::Subscribe)
    mqtt_write(s, packet.id)
    for topic in packet.topics
        for i in topic
            mqtt_write(s, i)
        end
    end
end

Base.show(io::IO, x::Subscribe) = print(io, "subscribe[id: ", x.id, ", topics: ", join(x.topics, ", "), "]")

struct Suback <: Ack
    header::UInt8
    id::UInt16
end
Suback(id) = Suback(SUBACK, id)
Base.show(io::IO, x::Suback) = print(io, "suback[id: ", x.id, "]")

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
using BibInternal
using Test

# TODO: make error testing for new version
#@test true

@testset "isless() for BibInternal.Date" begin

    @test BibInternal.Date("","","2000") < BibInternal.Date("","","2001");
    @test BibInternal.Date("","10","2000") < BibInternal.Date("","11","2000");
    @test BibInternal.Date("5","01","2000") < BibInternal.Date("15","01","2000");

    @test BibInternal.Date("","","1970") > BibInternal.Date("","","1969");
    @test BibInternal.Date("","7","1970") > BibInternal.Date("","6","1970");
    @test BibInternal.Date("7","7","1970") > BibInternal.Date("6","7","1970");

    @test BibInternal.Date("","","1805") == BibInternal.Date("","","1805");
    @test BibInternal.Date("","4","1805") == BibInternal.Date("","4","1805");
    @test BibInternal.Date("3","4","1805") == BibInternal.Date("3","4","1805");

    @test BibInternal.Date("","","1805") !== BibInternal.Date("","","1905");
    @test BibInternal.Date("","4","1805") !== BibInternal.Date("","5","1805");
    @test BibInternal.Date("3","4","1805") !== BibInternal.Date("4","4","1805");

    @test_skip BibInternal.Date("","May","1805") == BibInternal.Date("","5","1805");
    @test_skip BibInternal.Date("1th","5","1805") !== BibInternal.Date("1","5","1805");

    @test_throws ArgumentError BibInternal.Date("","","1") < BibInternal.Date("","","2000a");
    @test_throws ArgumentError BibInternal.Date("","","1") < BibInternal.Date("","","");
    @test_throws ArgumentError BibInternal.Date("","","1") < BibInternal.Date("","","90ies");
    @test_throws ArgumentError BibInternal.Date("","","40k") < BibInternal.Date("","","40000");

end

@testset "isless() for BibInternal.Name" begin
    # TODO: consider testing for alphabetizing rules

    @test BibInternal.Name("","Cow","","John","") < BibInternal.Name("","Doe","","John","");
    @test BibInternal.Name("","Doe","","John","A.") < BibInternal.Name("","Doe","","John","B.");
    @test BibInternal.Name("","Doe","","Bronn","") < BibInternal.Name("","Doe","","John","");
    @test BibInternal.Name("","Bronn","","Would","") < BibInternal.Name("","Bronn","","Would","Not");

    @test BibInternal.Name("","Doe","","John","") == BibInternal.Name("","Doe","","John","");
    @test BibInternal.Name("","Doe","jun.","John","") == BibInternal.Name("","Doe","jun.","John","");
    @test BibInternal.Name("","Doe","","John","E.") == BibInternal.Name("","Doe","","John","E.");

    @test BibInternal.Name("","Doe","","Bronn","") !== BibInternal.Name("","Doe","","John","");
    @test BibInternal.Name("","Doe","jun.","John","") !== BibInternal.Name("","Doe","sen.","John","");
    @test BibInternal.Name("","Doe","","John","E.") !== BibInternal.Name("","Doe","","John","A.");

end


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "Coveralls" begin
    @testset "Plugin creation" begin
        p = Coveralls()
        @test p.gitignore == ["*.jl.cov", "*.jl.*.cov", "*.jl.mem"]
        @test p.src === nothing
        @test p.dest == ".coveralls.yml"
        @test p.badges == [
            Badge(
                "Coveralls",
                "https://coveralls.io/repos/github/{{USER}}/{{PKGNAME}}.jl/badge.svg?branch=master",
                "https://coveralls.io/github/{{USER}}/{{PKGNAME}}.jl?branch=master",
            )
        ]
        @test isempty(p.view)
        p = Coveralls(; config_file=nothing)
        @test p.src === nothing
        p = Coveralls(; config_file=test_file)
        @test p.src == test_file
        @test_throws ArgumentError Coveralls(; config_file=fake_path)
    end

    @testset "Badge generation" begin
        p = Coveralls()
        @test badges(p, me, test_pkg) == ["[![Coveralls](https://coveralls.io/repos/github/$me/$test_pkg.jl/badge.svg?branch=master)](https://coveralls.io/github/$me/$test_pkg.jl?branch=master)"]
    end

    @testset "File generation" begin
        p = Coveralls()
        @test isempty(gen_plugin(p, t, test_pkg))
        @test !isfile(joinpath(pkg_dir, ".coveralls.yml"))
        p = Coveralls(; config_file=test_file)
        @test gen_plugin(p, t, test_pkg) == [".coveralls.yml"]
        @test isfile(joinpath(pkg_dir, ".coveralls.yml"))
    end
end

rm(pkg_dir; recursive=true)


text:
#######################
# selection functions #
#######################

struct SelectionFunction
    julia_function::Function
end

macro sel(ast)
    if ast.head != :function
        error("syntax error at $(ast) in $(ast.head)")
    end
    if length(ast.args) != 2
        error("syntax error at $(ast) in $(ast.args)")
    end
    signature = ast.args[1]
    body = ast.args[2]
    if signature.head != :call
        error("syntax error at $(ast) in $(signature)")
    end
    function_name = signature.args[1]
    args = signature.args[2:end]
    escaped_args = map(esc, args)
    fn_args = [esc(state), escaped_args...]
    Expr(:block,
        Expr(:(=),
            esc(function_name),
            Expr(:call, :SelectionFunction,
                Expr(:function, Expr(:tuple, fn_args...), esc(body)))))
end

macro select(addr)
    Expr(:call, :select_addr, esc(state), esc(addr))
end

struct SelectionState{T}
    assignment::T
    selection::DynamicAddressSet
    visitor::AddressVisitor
end

function SelectionState(assignment::T) where {T}
    SelectionState(assignment, DynamicAddressSet(), AddressVisitor())
end

function exec(fn::SelectionFunction, state::SelectionState, args::Tuple)
    fn.julia_function(state, args...)
end

# NOTE: unlike for @bijective and @gen functions, invoking @addr
# changes both the namespace of @read (and @select)

function addr(state::SelectionState, fn::SelectionFunction, args, addr)
    visit!(state.visitor, addr)
    if has_internal_node(state.assignment, addr)
        node = get_internal_node(state.assignment, addr)
    else
        node = EmptyAssignment()
    end
    if has_internal_node(state.selection, addr) || has_leaf_node(state.selection, addr)
        error("Address $addr, or an address under $addr, was already selected")
    end
    (selection, value) = select(fn, args, node)
    set_internal_node!(state.selection, addr, selection)
    value
end

function read(state::SelectionState, addr)
    get_leaf_node(state.assignment, addr)
end

function select_addr(state::SelectionState, addr)
    visit!(state.visitor, addr)
    push_leaf_node!(state.selection, addr)
end

function splice(state::SelectionState, fn::SelectionFunction, args::Tuple)
    exec(fn, state, args)
end

function select(fn::SelectionFunction, args::Tuple, assignment::T) where {T}
    state = SelectionState(assignment)
    value = exec(fn, state, args)
    (state.selection, value)
end

# alternative syntax for constructing selection from a flat list of addresses

function select(addrs...)
    set = DynamicAddressSet()
    for addr in addrs
        push_leaf_node!(set, addr)
    end
    return set
end

export @sel
export @select
export select

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
module test_pkgs_pycall

using Test
using PyCall

numpy = pyimport("numpy")
@test numpy.math.__name__ == "math"
@test repr(pytypeof(numpy.math)) == "PyObject <class 'module'>"

sys = pyimport("sys")
pushfirst!(sys."path", @__DIR__)

ex = pyimport("ex")
@test ex.f() == 42

end # module test_pkgs_pycall


text:
"""
**Template Likelihood**

Template file for likelihood creation

```julia
TemplateLikelihood()
```
See all functions you need to implement
---

"""
struct TemplateLikelihood <: AbstractLikelihood
    ## Additional parameters can be added
end

function implemented(
    ::TemplateLikelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function pdf(l::TemplateLikelihood, y::Real, f::Real) end

function Base.show(io::IO, model::TemplateLikelihood)
    return print(io, "Template Likelihood")
end

function compute_proba(
    l::TemplateLikelihood, μ::AbstractVector{T}, σ²::AbstractVector{T}
) where {T}
    N = length(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    return pred, sig_pred
end

### Local Updates Section ###

function local_updates!(
    local_vars,
    l::TemplateLikelihood,
    y::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
)
    # Update the local variables here (in place if you want)
    return local_vars
end

function sample_local!(
    local_vars, l::TemplateLikelihood, y::AbstractVector, f::AbstractVector
)
    # Sample the local variables here (in place if you want)
    return local_vars
end

### Natural Gradient Section ###

∇E_μ(l::TemplateLikelihood, ::AOptimizer, y, state) = (nothing,)
∇E_Σ(l::TemplateLikelihood, ::AOptimizer, y, state) = (nothing,)


text:
opr_pm = [+|-]
opr_dm = [/|*]

oprlr_eq_dif = [==|!=]
oprlr_lgt = [<|>]
oprlr_lgt_eq = [<=|>=]

comma = ","
o_bracket = ["["]
c_bracket = ["]"]

oprl_algebc = {oprlr_eq_dif} | {oprlr_lgt} | {oprlr_lgt_eq}
opr_algebc = opr_pm | opr_dm

oprl_and_or = [:blank:]+"and"[:blank:]+ | [:blank:]+"or"[:blank:]+
opr_not = "not"[:blank:]*

oprlr = {oprlr_eq_dif} | {oprlr_lgt} | {oprlr_lgt_eq} | {orplr_and_or}

opr_concat = "++"+
ct_str = "\"".*"\""?[{opr_concat}{ct_string}?{ct_str}]

ct_int = ^{opr_pm}?[0-9]+$
ct_float = {ct_int}[.][0-9]*$
ctn = {ct_int}|{ct_float}
id = ^[:alpha:]+[:alnum:]*

loop_while = "while" [:blank:]+"("{expr_bool}")"[:blank:]+ "{"{cmd}"}"
loop_for = "for"[:blank:]+"("{int_'dclr'}{comma}{ct_int}{comma}{ct_int}")"[:blank:]*"{"{cmd}"}"
ctb = "1" | "0" | "true" | "false"
expr_bool =  {expr_arit} {opr_logic} {expr_arit} | {ctb}

term = {ctn}|{id}

expr_arit =  {term} {opr_algebc} {term} | "(" {expr_arit} [{opr_algebc}{expr_arit}]*")"[{opr_algebc}{expr_arit}]*
data = {expr_arit} | {term} | {expr_bool}
int_dclr = "int" {id} "=" [0-9]*

float_type = "float"
bool_type = "bool"
char_type = "char"
string_type = "string"
void_type = "void"

data_type = {float_type} | {bool_type} | {char_type} | {string_type} | {void}

var_dclr = {int_dclr} | [{bool_type} | {float_type} | {char_type} | {string_type}] {atrib}]
atrib = {id} "=" {data}
cmd = {loop} | {atrib} | {expr_arit} | {expr_bool} | {term}

params = {id}|{id}{comma}{id}?[{comma}{params}]
fn_dclr = {data_type} {id}"("{params}")"[:blank:]*"{"{cmd}*?["return" {data}]"}"
fn_call = {id}"("?{params}")"

comment_block = "/*"[:alnum:]*"*/" | "/*"{comment_block}"*/"
coment_line = "//"[:alnum:]*
comment = comment_block | comment_line

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
using DemoActions
using BenchmarkTools

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

test_algs = [eulers_inplace_f!, midpoint_inplace_f!, rk4_inplace_f!, dp5_inplace!]
test_models = [lotka!]
u0 = [1., 1.]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0., 10.)
dt = 0.25
# prob = MyProbInplace(test_algs[1], lotka!, u0, tspan, dt, p)

for alg in test_algs
    for f in test_models
        prob = MyProbInplace(alg, f, u0, tspan, dt, p)
        SUITE["inplace"][string(alg) * "_" * string(f), u0] = @benchmarkable mysolve_inplace($prob)
    end
end

run(SUITE)


text:
Unitful.register(@__MODULE__)
Base.eltype(::Type{Length{T}}) where T = T
@unit Ang "Ang" Angstrom              1e-1u"nm"               false
@unit e₀  "eₒ"  ElementaryCharge      1.602176620898e-19*u"C" false
@unit kₑ  "kₑ"  CoulombForceConstant  1/(4π)u"ϵ0"             false
@unit a₀  "a₀"  BohrRadius            1u"ħ^2/(1kₑ*me*e₀^2)"   false
@unit Eₕ  "Eₕ"  HartreeEnergy         1u"me*e₀^4*kₑ^2/(1ħ^2)" true
@unit Ry  "Ry"  RydbergEnergy         0.5Eₕ                   true

const localunits = Unitful.basefactors
const LengthType{T, A} = Quantity{T,𝐋,FreeUnits{A,𝐋,nothing}}
const ReciprocalType{T, A} = Quantity{T,𝐋^-1,FreeUnits{A,𝐋^-1,nothing}}

Base.show(io::IO, ::Type{Quantity{T,𝐋,FreeUnits{A,𝐋,nothing}}}) where {T,A} =
    dfprint(io, "LengthType{$T, $A}")

@inline function StaticArrays._inv(::StaticArrays.Size{(3,3)}, A::SMatrix{3,3, LT}) where {LT<:Union{LengthType, ReciprocalType}}

    @inbounds x0 = SVector{3}(A[1], A[2], A[3])
    @inbounds x1 = SVector{3}(A[4], A[5], A[6])
    @inbounds x2 = SVector{3}(A[7], A[8], A[9])

    y0 = cross(x1,x2)
    d  = StaticArrays.bilinear_vecdot(x0, y0)
    x0 = x0 / d
    y0 = y0 / d
    y1 = cross(x2,x0)
    y2 = cross(x0,x1)

    @inbounds return SMatrix{3, 3}((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3]))
end


text:
#=
Copyright 2019, E.A.

Contains code derived from the Python module "Newick" Copyright 2003-2008 by Thomas Mailund, released under the GPL v.2
=#

#==
    READING NEWICK STRINGS
==#

# TODO: Implement support for Newick comments in the format `[&...]`

const nullnode = Node()

#= Define token names & ther RegEx patterns.
Note that it is important the order in which the token definitions are given. The regexes for :number and :string will have matches in common, but the lexer will always process a :number token because it accepts the first match found in the token definition vector.
=#
const NEWICK_TOKENDEFS = [
    (:number,       r"^[+-]?\d+(\.\d+)?([eE][+-]?\d+)?"),
    (:colon,        r"^[\s\n]*:[\s\n]*"),
    (:semicolon,    r"^[\s\n]*;[\s\n]*"),
    (:comma,        r"^[\s\n]*,[\s\n]*"),
    (:rparen,       r"^[\s\n]*\)[\s\n]*"),
    (:lparen,       r"^[\s\n]*\([\s\n]*"),
    (:comment,      r"\[\[^\]]*\]"),
    (:string,       r"^((\"[^\"]+\")|('[^']+')|([^,:(); \r\t\n\[\]]+))"),
    (:whitespace,   r"^[\s\n\r]+")
    ]

function parsefork!(lxr::Lexer, p::AbstractNode)::Nothing
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    if next_token.name == :lparen
        parseinternal!(lxr, p)
    else
        parsetip!(lxr, p)
    end
    next_token = peektoken(lxr)
    skiptoken!(lxr, :comment)
    while next_token.name == :comma
        readtoken!(lxr, :comma)
        skiptoken!(lxr, :comment)
        next_token = peektoken(lxr)
        if next_token.name == :lparen
            parseinternal!(lxr, p)
        else
            parsetip!(lxr, p)
        end
        skiptoken!(lxr, :comment)
        next_token = peektoken(lxr)
    end

    return nothing
end

function parseinternal!(lxr::Lexer, p::AbstractNode)::AbstractNode
    readtoken!(lxr, :lparen)
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    q = Node()
    p ≠ nullnode && link!(p, q)
    parsefork!(lxr, q)
    readtoken!(lxr, :rparen)
    skiptoken!(lxr, :comment)
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    if next_token.name ∈ (:string, :number)
        q.label = next_token.value
        readtoken!(lxr, next_token.name)
        next_token = peektoken(lxr)
    end
    skiptoken!(lxr, :comment)
    if next_token.name == :colon
        p == nullnode && @error "Branch lengths at the root node are not supported."
        readtoken!(lxr, :colon)
        skiptoken!(lxr, :comment)
        next_token = readtoken!(lxr, :number)
        brlength!(p, q, parse(Float64, next_token.value))
    end

    return q
end

function parsetip!(lxr::Lexer, p::AbstractNode)::AbstractNode
    q = Node()
    next_token = peektoken(lxr)
    p ≠ nullnode && link!(p, q)

    if next_token.name ∈ (:string, :number)
        q.label = next_token.value
        next_token = readtoken!(lxr, next_token.name)
    end
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    if next_token.name == :colon
        p == nullnode && @error "Branch lengths at the root node are not supported."
        readtoken!(lxr, :colon)
        skiptoken!(lxr, :comment)
        next_token = readtoken!(lxr, :number)
        brlength!(p, q, parse(Float64, next_token.value))
    end
    skiptoken!(lxr, :comment)

    return q
end

function parse_newick(str::String)
    lxr = Lexer(str, nothing, NEWICK_TOKENDEFS)
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    if next_token.name == :lparen
        p = parseinternal!(lxr, nullnode)
    end
    readtoken!(lxr, :semicolon)

    return p
end

function parse_newick!(lxr::Lexer)
    skiptoken!(lxr, :comment)
    next_token = peektoken(lxr)
    if next_token.name == :lparen
        p = parseinternal!(lxr, nullnode)
    end
    readtoken!(lxr, :semicolon)

    return p
end

"""
    read_newick(filename::AbstractString; nhint::Int=1000)

Read trees in Newick format from a file.

Performance with large files can be improved by adjusting the `nhint` value to the approximate number of trees in the file (1000 by default).
"""
function read_newick(filename::AbstractString; nhint::Int=1000)::TreeVector
    file = open(filename)
    intrees = Vector{AbstractTree}(undef, nhint)
    treecount = 0
    notree = false
    while ! eof(file)
        instring = ""
        i = ' '
        while i ≠ ';'
            if eof(file)
                notree = occursin(r"^[\s\n\r]*$", instring)
                notree || @error "End of file reached before end of Newick string."
                break
            end
            i = read(file, Char)
            instring *= i
        end
        notree && break
        lxr = Lexer(instring, nothing, NEWICK_TOKENDEFS)
        tree_i = Tree(parse_newick!(lxr))
        treecount += 1
        if treecount <= nhint
            intrees[treecount] = tree_i
        else
            push!(intrees, tree_i)
        end
    end
    close(file)

    treecount < 1 && @error "No trees were found in the file."
    resize!(intrees, treecount)

    return intrees
end

"""
    read_newick(text=<NEWICK>)

Read trees in Newick format from a text string.
"""
function read_newick(;text::AbstractString)
    textlines = split(text, ';'; keepempty=false)
    ntrees = length(textlines)
    trees = Vector{Tree}(undef, ntrees)

    for (i, line) in enumerate(textlines)
        tree = Tree(parse_newick(line * ';'))
        trees[i] = tree
    end

    return trees
end

#==
    WRITING NEWICK STRINGS
==#

"""
Generate the Newick substring (not capped by \";\") of a non-origin `Node`. This function is intended to be always called from `_newick_string_origin`.
"""
function _newick_substring(p::AbstractNode, parent_p::AbstractNode)::String
    outstring = ""
    child_counter = 0
    for q in neighbours(p)
        q == parent_p && continue
        child_counter += 1
        outstring *= child_counter > 1 ? "," : ""
        outstring *= _newick_substring(q, p)
    end

    outstring = child_counter > 0 ? "(" * outstring * ")" : ""
    label_p = label(p)
    # TODO: Check that the label string does not break Newick and add quote marks or underscores if necessary.
    outstring *= label_p
    br_p = brlength(parent_p, p)
    outstring *= isnothing(br_p) ? "" : string(":", br_p)

    return outstring
end

"""
Generate the Newick substring (not capped by \";\") of the origin `Node`.
"""
function _newick_substring_origin(p::AbstractNode)::String
    outstring = ""
    child_counter = 0
    for q in neighbours(p)
        child_counter += 1
        outstring *= child_counter > 1 ? "," : ""
        outstring *= _newick_substring(q, p)
    end

    outstring = child_counter > 0 ? "(" * outstring * ")" : ""
    label_p = label(p)
    # TODO: Check that the label string does not break Newick and add quote marks or underscores if necessary.
    outstring *= label_p

    return outstring
end

"""
    newick_string(p::Node)::String

Return the Newick representation string of the tree subtended by node `p`.

If `p` belongs to a rooted tree, the function will return only the Newick string corresponding to the subtree subtended by `p`. With unrooted trees, the full Newick string is given, with the tree "pseudo-rooted" in `p`.
"""
newick_string(p::AbstractNode)::String = _newick_substring_origin(p) * ";"

"""
    newick_string(p::Node)::String

Return the Newick representation of tree `t`.
"""
newick_string(t::AbstractTree)::String =
    _newick_substring_origin(t.origin) * ";"

"""
    write_newick(filename::AbstractString, t::Union{Tree,TreeVector}; append::Bool=false)

Write the Newick representation of a tree or a vector of trees to a file.

Set `append` to `true` to add the tree(s) at the end instead of overwriting the contents of the file.
"""
function write_newick(filename::AbstractString, t::Union{AbstractTree,TreeVector}; append::Bool=false)
    open(filename, append ? "a" : "w") do f
        write_newick(f, t)
    end

    return nothing
end

function write_newick(io::IO, t::Tree)
    print(io, newick_string(t), "\n")

    return nothing
end

function write_newick(io::IO, t::TreeVector)
    for tt in t
        print(io, newick_string(tt), "\n")
    end

    return nothing
end

write_newick(t::Tree) = newick_string(t)

write_newick(t::Vector{Tree}) = newick_string.(t)


text:
My guess: 9517. How many bulls and cows?
0 1
My guess: 3285. How many bulls and cows?
1 1
My guess: 4298. How many bulls and cows?
0 2
My guess: 2465. How many bulls and cows?
0 2
My guess: 3926. How many bulls and cows?
2 0
My guess: 3724. How many bulls and cows?
3 0
My guess: 3124. How many bulls and cows?
4 0
The winning pick: 3124

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    a = [parse(Int, x) for x in split(readline())]
    b = [parse(Int, x) for x in split(readline())]

    dp = zeros(Int, n+1, m+1)
    for i=1:n+1
        dp[i, 1] = i-1
    end
    for j=1:m+1
        dp[1, j] = j-1
    end

    for i=2:n+1, j=2:m+1
        if a[i-1] == b[j-1]
            x = dp[i-1, j-1]
        else
            x = dp[i-1, j-1] + 1
        end
        dp[i, j] = min(dp[i-1, j]+1, dp[i, j-1]+1, x)
    end
    dp[n+1, m+1]
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using TSPLIB
using Test

@testset "TSPLIB.jl" begin
    exclude_list = [
        :rl11849
        :usa13509
        :brd14051
        :d15112
        :d18512
        :pla33810
        :pla85900
        ]
    for instance in keys(Optimals)
        if !(instance in exclude_list)
            @test readTSPLIB(instance).optimal == Optimals[instance]
        end
    end
end


text:
permute_error() = throw(ArgumentError("invalid permutation of dimensions"))
@generated function check_permute_parameters(::Val{N}, ::Val{perm}) where {N, perm}
    perm isa Tuple{Vararg{Int}} || return :(permute_error())
    N == length(perm) || return :(permute_error())
    for p in perm
        checkbounds(Bool, 1:N, p) || return :(permute_error())
    end
end

# Space

@generated function _permutedims(::Space{S}, ::Val{perm}) where {S, perm}
    # numbering each dimension
    num = Int[]
    for (i, s) in enumerate(S)
        append!(num, fill(i, length(s)))
    end
    num = num[collect(perm)] # allow to use invalid permutation in `_permutedims`
    # permute!(num, collect(perm))

    # collect the same number of group
    groups = [[num[1]]]
    for i in 2:length(num)
        if groups[end][end] == num[i] # check number
            push!(groups[end], num[i]) # add to the same group
        else
            push!(groups, [num[i]]) # create new group
        end
    end

    # apply `Symmetry` if needed
    newspace = map(groups) do group
        space_num = group[end]
        s = S[space_num]
        if length(group) == 1
            # if Symmetry, extract first dimension
            # if Real, just return it
            s[1]
        else
            if s isa Symmetry
                Symmetry(fill(s[1], length(group))...)
            else
                error() # unreachable
            end
        end
    end

    quote
        @_inline_meta
        Space($(newspace...))
    end
end

@inline function Base.permutedims(x::Space, perm::Val)
    check_permute_parameters(Val(tensororder(x)), perm)
    _permutedims(x, perm)
end

# AbstractTensor

@generated function ispermuted(::Val{N}, ::Val{perm}) where {N, perm}
    res = permute!(collect(1:N), collect(perm))
    quote
        @_inline_meta
        check_permute_parameters(Val(N), Val(perm))
        $(res == 1:N)
    end
end

@generated _revperm(::Val{perm}) where {perm} = :(@_inline_meta; tuple($(sortperm(collect(perm))...)))

@inline function Base.permutedims(x::AbstractTensor, perm::Val)
    ispermuted(Val(ndims(x)), perm) && return x
    S = permutedims(Space(x), perm)
    tensortype(S)() do ij...
        @_inline_meta
        @inbounds x[getindex.(Ref(ij), _revperm(perm))...]
    end
end


text:
export norminvAS241, norminvBSMS

const SPLIT1 = 0.425
const SPLIT2 = 5.0
const CONST1 = 0.180625
const CONST2 = 1.6

# Coefficients for P close to 0.5
const A0 = 3.3871328727963666080E0
const A1 = 1.3314166789178437745E+2
const A2 = 1.9715909503065514427E+3
const A3 = 1.3731693765509461125E+4
const A4 = 4.5921953931549871457E+4
const A5 = 6.7265770927008700853E+4
const A6 = 3.3430575583588128105E+4
const A7 = 2.5090809287301226727E+3
const B1 = 4.2313330701600911252E+1
const B2 = 6.8718700749205790830E+2
const B3 = 5.3941960214247511077E+3
const B4 = 2.1213794301586595867E+4
const B5 = 3.9307895800092710610E+4
const B6 = 2.8729085735721942674E+4
const B7 = 5.2264952788528545610E+3

#  Coefficients for P not close to 0, 0.5 or 1.
const C0 = 1.42343711074968357734E0
const C1 = 4.63033784615654529590E0
const C2 = 5.76949722146069140550E0
const C3 = 3.64784832476320460504E0
const C4 = 1.27045825245236838258E0
const C5 = 2.41780725177450611770E-1
const C6 = 2.27238449892691845833E-2
const C7 = 7.74545014278341407640E-4
const D1 = 2.05319162663775882187E0
const D2 = 1.67638483018380384940E0
const D3 = 6.89767334985100004550E-1
const D4 = 1.48103976427480074590E-1
const D5 = 1.51986665636164571966E-2
const D6 = 5.47593808499534494600E-4
const D7 = 1.05075007164441684324E-9

#  Coefficients for P near 0 or 1.
const E0 = 6.65790464350110377720E0
const E1 = 5.46378491116411436990E0
const E2 = 1.78482653991729133580E0
const E3 = 2.96560571828504891230E-1
const E4 = 2.65321895265761230930E-2
const E5 = 1.24266094738807843860E-3
const E6 = 2.71155556874348757815E-5
const E7 = 2.01033439929228813265E-7
const F1 = 5.99832206555887937690E-1
const F2 = 1.36929880922735805310E-1
const F3 = 1.48753612908506148525E-2
const F4 = 7.86869131145613259100E-4
const F5 = 1.84631831751005468180E-5
const F6 = 1.42151175831644588870E-7
const F7 = 2.04426310338993978564E-15

#InverseCdf computes the inverse cumulative normal distribution value using AS241 algorithm
function norminvAS241(p::Float64)::Float64
    #    DOUBLE PRECISION FUNCTION PPND16 (P, IFAULT)
    # ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
    # Produces the normal deviate Z corresponding to a given lower
    #   tail area of P Z is accurate to about 1 part in 10**16.
    #
    #   The hash sums below are the sums of the mantissas of the
    #   coefficients.   They are included for use in checking
    #  transcription.

    q = 0.0
    r = 0.0
    ppnd16 = 0.0

    q = p - 0.5
    if abs(q) <= SPLIT1
        r = CONST1 - q * q
        ppnd16 =
            q * (
                ((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r +
                A0
            ) /
            (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0)
        return ppnd16
    else
        if q < 0
            r = p
        else
            r = 1.0 - p
        end

        if r <= 0
            if p == 0.0
                ppnd16 = typemin(Float64)
            elseif p == 1.0
                ppnd16 = typemax(Float64)
            else
                ppnd16 = NaN
            end
            return ppnd16
        end

        r = sqrt(-log(r))
        if r <= SPLIT2
            r = r - CONST2
            ppnd16 =
                (
                    ((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) *
                    r + C0
                ) / (
                    ((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) *
                    r + 1.0
                )
        else
            r = r - SPLIT2
            ppnd16 =
                (
                    ((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) *
                    r + E0
                ) / (
                    ((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) *
                    r + 1.0
                )
        end
        return q < 0 ? -ppnd16 : ppnd16
    end
end

const BSMS_A0 = 2.5066282186077179315
const BSMS_A1 = -18.700226003845536893
const BSMS_A2 = 41.841057826285191701
const BSMS_A3 = -25.965959614983341521

const BSMS_B0 = -8.5075119518749713281
const BSMS_B1 = 23.298480931332967694
const BSMS_B2 = -21.419275595427556433
const BSMS_B3 = 3.2269921657498348465

const BSMS_C0 = 0.3374754822726147
const BSMS_C1 = 0.9761690190917186
const BSMS_C2 = 0.1607979714918209
const BSMS_C3 = 0.0276438810333863
const BSMS_C4 = 0.0038405729373609
const BSMS_C5 = 0.0003951896511919
const BSMS_C6 = 0.0000321767881768
const BSMS_C7 = 0.0000002888167364
const BSMS_C8 = 0.0000003960315187

function norminvBSMS(x::Float64)::Float64
    temp = x - 0.5
    local result
    if abs(temp) < 0.42
        result = temp*temp
        result =
            temp * (((BSMS_A3 * result + BSMS_A2) * result + BSMS_A1) * result + BSMS_A0) /
            (
                (((BSMS_B3 * result + BSMS_B2) * result + BSMS_B1) * result + BSMS_B0) *
                result + 1.0
            )
    else
        if (x < 0.5)
            result = x
        else
            result = 1.0 - x
        end
            result = log(-log(result))
            result =
                BSMS_C0 +
                result * (
                    BSMS_C1 +
                    result * (
                        BSMS_C2 +
                        result * (
                            BSMS_C3 +
                            result * (
                                BSMS_C4 +
                                result * (
                                    BSMS_C5 +
                                    result *
                                    (BSMS_C6 + result * (BSMS_C7 + result * BSMS_C8))
                                )
                            )
                        )
                    )
                )
        if (x < 0.5)
            result = -result
        end
    end

    return result
end


text:
# short named functions, not flexible, bad assumptions or improper namings might be commited

function convolution( image::Array{T1,N}, kernel::Array{T2,N}; typ=Float32 ) where {T1<:Real,T2<:Real,N}
	return FFTConvolution( image, kernel, typ=typ );
end

function convolution_b( image::Array{T1,2}, kernel::Array{T2,2}; typ=Float32 ) where {T1<:Real,T2<:Real}
	return FFTConvolution_b( image, kernel, typ=typ );
end

function smooth( image::Array{T1,2}, size ) where {T1<:Real}
	if size == 0
		return image
	end
	meanKern = meanKernel( size );
	return FFTConvolution( image, meanKern )
end

function gradient( image::Array{T1,2}; kernel="Prewitt" ) where {T1<:Real}

	kh = ( kernel=="Sobel" ) ? Sobel_h() : Prewitt_h();
    kv = ( kernel=="Sobel" ) ? Sobel_v() : Prewitt_v();

	Gy = convolution( image, kh )
	Gx = convolution( image, kv )

	return Gy[2:end-1,2:end-1], Gx[2:end-1,2:end-1]
end

# smoothing kernel

function meanKernel( size::Integer; typ=Float32 )
	return ones( typ, size, size )./( size*size )
end

# edge kernels

# Prewitt, separable
function Prewitt_h( ; typ=Float32 )
	p1, z0, n1 = typ.( ( 1, 0, -1 ) )
	return [ p1 z0 n1; p1 z0 n1; p1 z0 n1 ]
end

function Prewitt_v( ; typ=Float32 )
	p1, z0, n1 = typ.( ( 1, 0, -1 ) )
	return [ p1 p1 p1; z0 z0 z0; n1 n1 n1 ]
end

# Sobel, separable
function Sobel_h( ; typ=Float32 )
	p1, p2, z0, n1, n2 = typ.( ( 1, 2, 0, -1, -2 ) )
	return [ p1 z0 n1; p2 z0 n2; p1 z0 n1 ]
end

function Sobel_v( ; typ=Float32 )
	p1, p2, z0, n1, n2 = typ.( ( 1, 2, 0, -1, -2 ) )
	return [ p1 p2 p1; z0 z0 z0; n1 n2 n1 ]
end

# Cross-correlation and convolution

function paddedConvolution( image::Array{T1,2}, kernel::Array{T2,2}; typ=Float32 ) where {T1<:Real,T2<:Real}

	N1, N2 = size( image  );
	M1, M2 = size( kernel );

	paddedImage = similar(image, N1+2*M1-2, N2+2*M2-2);
	paddedImage[ M1:(M1+N1-1), M2:(M2+N2-1) ] .= image;

	corrmatrix = zeros( typ, N1 + M1 - 1, N2 + M2 - 1 )

	for col in 1:( N2 + M2 - 1)
		for row in 1:( N1 + M1 - 1 )
			corrmatrix[ row, col ] = sum( paddedImage[ row:row+M1-1, col:col+M2-1 ] .* kernel )
		end
	end

	return corrmatrix
end

function FFTConvolution!( image::Array{<:Real,N}, kernel::Array{<:Real,N}, conv::Array{T,N} ) where {T,N}

	@assert all( size(conv) .== size(image) .+ size(kernel) .- 1 );

	paddedkernel = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2), 1:size(kernel,3) ] .= convert.( Complex{T}, kernel )

	paddedimage  = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2), 1:size(image,3) ] .= convert.( Complex{T}, image )

	FFTW.fft!( paddedimage  );
	FFTW.fft!( paddedkernel );

	@inbounds @simd for idx in 1:length(paddedimage)
		paddedimage[idx] = paddedimage[idx] * paddedkernel[idx]
	end

	FFTW.ifft!( paddedimage );

	@inbounds @simd for idx in 1:length(paddedimage)
		conv[idx] = real( paddedimage[idx] )
	end

	return conv
end

function FFTConvolution( image::Array{<:Real,N}, kernel::Array{<:Real,N}; typ=Float32 ) where {N}

	return FFTConvolution!( image, kernel, zeros( typ, size(image).+size(kernel).-1) );
end

function FFTConvolution_crop!( image::Array{<:Real,N}, kernel::Array{<:Real,N}, conv::Array{T,N} ) where {T,N}

	@assert all( size(conv) .== size(image) );

	paddedkernel = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2), 1:size(kernel,3) ] .= convert.( Complex{T}, kernel )

	paddedimage  = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2), 1:size(image,3) ] .= convert.( Complex{T}, image )

	FFTW.fft!( paddedimage  );
	FFTW.fft!( paddedkernel );

	@inbounds @simd for idx in 1:length(paddedimage)
		paddedimage[idx] = paddedimage[idx] * paddedkernel[idx]
	end

	FFTW.ifft!( paddedimage );

	offs = div.( size( kernel ), 2 );
	if N == 2
		offs[3] = 0
	end

	for zet in 1+offs[3]:size(image,3)+offs[3]
	for col in 1+offs[2]:size(image,2)+offs[2]
		@simd for row in 1+offs[1]:size(image,1)+offs[1]
			conv[row-offs[1],col-offs[2],zet-offs[3]] = real( paddedimage[row,col,zet] )
	end end end

	return conv
end

function FFTConvolution_crop( image::Array{<:Real,N}, kernel::Array{<:Real,N}; typ=Float32 ) where {N}

	return FFTConvolution_crop!( image, kernel, zeros( typ, size(image) ) );
end

function FFTConvolution_b( image::Array{T1,2}, kernel::Array{T2,2}; typ=Float32 ) where {T1<:Real,T2<:Real}

	paddedkernel = zeros( typ, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2) ] .= convert.( typ, kernel )

	paddedimage  = zeros( typ, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2) ] .= convert.( typ, image )

	Ff = FFTW.fft( paddedimage  );
	Fg = FFTW.fft( paddedkernel );

	Fcorr =  Ff .* Fg;

	return real.( FFTW.bfft( Fcorr ) )
end


text:
struct Heisenberg{D} <: AbstractHeisenberg{D}
    size::NTuple{D, Int}
    periodic::Bool
    Heisenberg(size::Int...; periodic::Bool) = new{length(size)}(size, periodic)
end

Base.size(model::Heisenberg) = model.size

heisenberg_ij(nbit::Int, i::Int, j::Int=i+1) = put(nbit, i=>X)*put(nbit, j=>X) + put(nbit, i=>Y)*put(nbit, j=>Y) + put(nbit, i=>Z)*put(nbit, j=>Z)
const heisenberg_term = repeat(2, X, 1:2) + repeat(2, Y, 1:2) + repeat(2, Z, 1:2)

function hamiltonian(model::Heisenberg)
    nbit = nspin(model)
    sum(x->heisenberg_ij(nbit, x[1], x[2]), get_bonds(model))*0.25
end

function get_bonds(model::Heisenberg{2})
    m, n = model.size
    cis = LinearIndices(model.size)
    bonds = Tuple{Int, Int, Float64}[]
    for i=1:m, j=1:n
        (i!=m || model.periodic) && push!(bonds, (cis[i,j], cis[i%m+1,j], 1.0))
        (j!=n || model.periodic) && push!(bonds, (cis[i,j], cis[i,j%n+1], 1.0))
    end
    bonds
end

function get_bonds(model::Heisenberg{1})
    nbit, = model.size
    [(i, i%nbit+1, 1.0) for i in 1:(model.periodic ? nbit : nbit-1)]
end

max_stars_repo_name: Roger-luo/QuantumMPS

max_issues_repo_name: Roger-luo/QuantumMPS

max_forks_repo_name: Roger-luo/QuantumMPS


text:
# Parity operators based on MPB code:
# https://github.com/NanoComp/mpb/blob/master/src/maxwell/maxwell_constraints.c
using LoopVectorization

export 𝓟x!, 𝓟x̄!, 𝓟x, 𝓟x̄

function _𝓟x!(Ho::AbstractArray{Complex{T},3},Hi::AbstractArray{Complex{T},3},grid::Grid{2}) where T<:Real
    Nx,Ny = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iy ∈ 1:Ny
		for ix ∈ 1:(Nx÷2), l in 0:0 #1:((Nx÷2)-1), l in 0:0
			ix1 = ix
			ix2 = (ix > 1 ? Nx - (ix-2) : 1)
			temp1 		= 0.5*( Hi[1+l,ix1,iy] - Hi[1+l,ix2,iy] )
			temp2 		= 0.5*( Hi[2+l,ix1,iy] + Hi[2+l,ix2,iy] )
			Ho[1+l,ix1,iy] 	= temp1
			Ho[2+l,ix1,iy] 	= temp2
			Ho[1+l,ix2,iy] = -temp1
			Ho[2+l,ix2,iy] = temp2
		end
	end
	return Ho
end

function _𝓟x̄!(Ho::AbstractArray{Complex{T},3},Hi::AbstractArray{Complex{T},3},grid::Grid{2}) where T<:Real
    Nx,Ny = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iy ∈ 1:Ny
		for ix ∈ 1:(Nx÷2), l in 0:0 # 1:((Nx÷2)-1), l in 0:0
			ix1 = ix
			ix2 = (ix > 1 ? Nx - (ix-2) : 1)
			temp1 		= 0.5*( Hi[1+l,ix1,iy] + Hi[1+l,ix2,iy] )
			temp2 		= 0.5*( Hi[2+l,ix1,iy] - Hi[2+l,ix2,iy] )
			Ho[1+l,ix1,iy] 	= temp1
			Ho[2+l,ix1,iy] 	= temp2
			Ho[1+l,ix2,iy] = temp1
			Ho[2+l,ix2,iy] = -temp2
		end
	end
	return Ho
end

function _𝓟x!(Ho::AbstractArray{Complex{T},4},Hi::AbstractArray{Complex{T},4},grid::Grid{3}) where T<:Real
    Nx,Ny,Nz = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iz ∈ 1:Nz, iy ∈ 1:Ny
		for ix ∈ 1:(Nx÷2), l in 0:0 # 1:((Nx÷2)-1), l in 0:0
			ix1 = ix
			ix2 = (ix > 1 ? Nx - (ix-2) : 1)
			temp1 		= 0.5*( Hi[1+l,ix1,iy,iz] - Hi[1+l,ix2,iy,iz] )
			temp2 		= 0.5*( Hi[2+l,ix1,iy,iz] + Hi[2+l,ix2,iy,iz] )
			Ho[1+l,ix1,iy,iz] 	= temp1
			Ho[2+l,ix1,iy,iz] 	= temp2
			Ho[1+l,ix2,iy,iz] = -temp1
			Ho[2+l,ix2,iy,iz] = temp2
		end
	end
	return Ho
end

function _𝓟x̄!(Ho::AbstractArray{Complex{T},4},Hi::AbstractArray{Complex{T},4},grid::Grid{3}) where T<:Real
    Nx,Ny,Nz = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iz ∈ 1:Nz, iy ∈ 1:Ny
		for ix ∈ 1:(Nx÷2), l in 0:0 # 1:((Nx÷2)-1), l in 0:0
			ix1 = ix
			ix2 = (ix > 1 ? Nx - (ix-2) : 1)
			temp1 		= 0.5*( Hi[1+l,ix1,iy,iz] + Hi[1+l,ix2,iy,iz] )
			temp2 		= 0.5*( Hi[2+l,ix1,iy,iz] - Hi[2+l,ix2,iy,iz] )
			Ho[1+l,ix1,iy,iz] 	= temp1
			Ho[2+l,ix1,iy,iz] 	= temp2
			Ho[1+l,ix2,iy,iz] = temp1
			Ho[2+l,ix2,iy,iz] = -temp2
		end
	end
	return Ho
end

𝓟x!(grid::Grid{2}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x!(reshape(Ho,(2,grid.Nx,grid.Ny)),reshape(Hi,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=true)
𝓟x!(grid::Grid{3}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x!(reshape(Ho,(2,grid.Nx,grid.Ny,grid.Nz)),reshape(Hi,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=true)
𝓟x̄!(grid::Grid{2}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x̄!(reshape(Ho,(2,grid.Nx,grid.Ny)),reshape(Hi,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=true)
𝓟x̄!(grid::Grid{3}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x̄!(reshape(Ho,(2,grid.Nx,grid.Ny,grid.Nz)),reshape(Hi,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=true)

𝓟x(grid::Grid{2}) = LinearMap{ComplexF64}(H -> vec(_𝓟x!(reshape(copy(H),(2,grid.Nx,grid.Ny)),reshape(H,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=false)
𝓟x(grid::Grid{3}) = LinearMap{ComplexF64}(H -> vec(_𝓟x!(reshape(copy(H),(2,grid.Nx,grid.Ny,grid.Nz)),reshape(H,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=false)
𝓟x̄(grid::Grid{2}) = LinearMap{ComplexF64}(H -> vec(_𝓟x̄!(reshape(copy(H),(2,grid.Nx,grid.Ny)),reshape(H,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=false)
𝓟x̄(grid::Grid{3}) = LinearMap{ComplexF64}(H -> vec(_𝓟x̄!(reshape(copy(H),(2,grid.Nx,grid.Ny,grid.Nz)),reshape(H,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=false)

# import Base: * #, transpose, adjoint
# function Base.:(*)(M::LinearMaps.FunctionMap,X::Matrix)
# 	#if isequal(size(M),size(X)) # size check?
# 	ncolsX = size(X)[2]
# 	# @assert ncolsX == size(M)[1]
# 	Y = similar(X)
# 	for i in 1:ncolsX
# 		@views Y[:,i] = M * X[:,i]
# 	end
# 	return Y
# end

#
# function LinearAlgebra.mul!(y::AbstractVecOrMat, M̂::HelmholtzMap, x::AbstractVector)
#     LinearMaps.check_dim_mul(y, M̂, x)
# 	M̂(y, x)
# end
#
# function LinearAlgebra.mul!(y::AbstractVecOrMat, P̂::HelmholtzPreconditioner, x::AbstractVector)
#     LinearMaps.check_dim_mul(y, P̂, x)
# 	P̂(y, x)
# end
#
# Base.adjoint(A::HelmholtzMap) = A
# Base.transpose(P̂::HelmholtzPreconditioner) = P̂
# LinearAlgebra.ldiv!(c,P̂::HelmholtzPreconditioner,b) = mul!(c,P̂,b) # P̂(c, b) #
# LinearAlgebra.ldiv!(P̂::HelmholtzPreconditioner,b) = mul!(b,P̂,b)


text:
function FiniteDifferences.to_vec(x::RegularSpacing)
    function from_vec_RegularSpacing(x_vec)
        return RegularSpacing(x_vec[1], x_vec[2], x.N)
    end
    return [x.t0, x.Δt], from_vec_RegularSpacing
end

@testset "regular_data" begin
    t0 = randn()
    Δt = randn()
    N = 5
    x = RegularSpacing(t0, Δt, N)
    x_range = range(t0; step=Δt, length=N)

    @test size(x) == size(x_range)
    @test getindex(x, 3) ≈ getindex(x_range, 3)
    @test collect(x) ≈ collect(x_range)
    @test step(x) == step(x_range)
    @test length(x) == length(x_range)

    let
        x, back = Zygote.pullback(RegularSpacing, t0, Δt, N)

        Δ_t0 = randn()
        Δ_Δt = randn()
        @test back((t0 = Δ_t0, Δt = Δ_Δt, N=nothing)) == (Δ_t0, Δ_Δt, nothing)

        adjoint_test((t0, Δt) -> RegularSpacing(t0, Δt, 10), (randn(), randn()))
    end
end


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

mutable struct GrantAccessData <: SwaggerModel
    access::Any # spec type: Union{ Nothing, String } # spec name: access
    durationInSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: durationInSeconds

    function GrantAccessData(;access=nothing, durationInSeconds=nothing)
        o = new()
        validate_property(GrantAccessData, Symbol("access"), access)
        setfield!(o, Symbol("access"), access)
        validate_property(GrantAccessData, Symbol("durationInSeconds"), durationInSeconds)
        setfield!(o, Symbol("durationInSeconds"), durationInSeconds)
        o
    end
end # type GrantAccessData

const _property_map_GrantAccessData = Dict{Symbol,Symbol}(Symbol("access")=>Symbol("access"), Symbol("durationInSeconds")=>Symbol("durationInSeconds"))
const _property_types_GrantAccessData = Dict{Symbol,String}(Symbol("access")=>"String", Symbol("durationInSeconds")=>"Int32")
Base.propertynames(::Type{ GrantAccessData }) = collect(keys(_property_map_GrantAccessData))
Swagger.property_type(::Type{ GrantAccessData }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_GrantAccessData[name]))}
Swagger.field_name(::Type{ GrantAccessData }, property_name::Symbol) =  _property_map_GrantAccessData[property_name]

const _allowed_GrantAccessData_access = ["None", "Read", "Write"]

function check_required(o::GrantAccessData)
    (getproperty(o, Symbol("durationInSeconds")) === nothing) && (return false)
    true
end

function validate_property(::Type{ GrantAccessData }, name::Symbol, val)
    if name === Symbol("access")
        Swagger.validate_param(name, "GrantAccessData", :enum, val, _allowed_GrantAccessData_access)
    end
end


text:
cfld(x::NTuple{2,Int}, n) = (cld(x[1], n), fld(x[2], n))

function sub_extdegree(X::AbstractVector{<:MP.AbstractMonomial}, vars)
    if isempty(X)
        return (0, 0)
    else
        return extrema(map(mono -> sum(var -> MP.degree(mono, var), vars), X))
    end
end

# Cheap approximation of the convex hull as the approximation of:
#
# z such that mindegree < sum(z) < maxdegree
# |\
# |#\ <-------- sum(z) = maxdegree
# |##\
# |\<-\-------- sum(z) = mindegree
# | \##\
# +---------
#
# and:
#
# z such that minmultideg < z < maxmultideg
# | +----+ <--- maxmultidegree
# | |####|
# | |####|
# | +----+
# | ^---------- minmultidegree
# +---------

function _filter(X::AbstractVector{<:MP.AbstractMonomial}, extdeg, exp, n)
    mindeg, maxdeg = cfld(extdeg, 2)
    minmultideg, maxmultideg = Vector{Int}(undef, n), Vector{Int}(undef, n)
    for i in 1:n
        exponent_i(mono) = exp(mono, i)
        minmultideg[i] = cld(mapreduce(exponent_i, min, X), 2)
        maxmultideg[i] = fld(mapreduce(exponent_i, max, X), 2)
    end
    return mindeg, maxdeg, mono -> begin
        all(i -> minmultideg[i] <= exp(mono, i) <= maxmultideg[i], 1:n)
    end
end
function _full_filter(X::AbstractVector{<:MP.AbstractMonomial}, extdeg, exp, n)
    mindeg, maxdeg, filter = _filter(X, extdeg, exp, n)
    return mono -> mindeg <= MP.degree(mono) <= maxdeg && filter(mono)
end

function _sub_half_newton_polytope(X::AbstractVector{<:MP.AbstractMonomial},
                                   extdeg, exp, vars)
    mindeg, maxdeg, filter = _filter(X, extdeg, exp, length(vars))
    return MP.monomials(vars, mindeg:maxdeg, filter)
end

function sub_half_newton_polytope(X::AbstractVector{<:MP.AbstractMonomial},
                                  vars)
    _sub_half_newton_polytope(X, sub_extdegree(X, vars),
                              (mono, i) -> MP.degree(mono, vars[i]), vars)
end

is_commutative(vars) = length(vars) < 2 || prod(vars[1:2]) == prod(reverse(vars[1:2]))

# Multipartite
# TODO we might do this recursively : do 2 parts, merge them, merge with next
#      one and so on so that the filter at the end prunes more.
function half_newton_polytope(X::AbstractVector, parts::Tuple; apply_post_filter = true)
    if !is_commutative(MP.variables(X))
        throw(ArgumentError("Multipartite Newton polytope not supported with noncommutative variables."))
    end
    if !all(i -> all(j -> i == j || isempty(parts[i] ∩ parts[j]),
                     1:length(parts)),
            1:length(parts))
        throw(ArgumentError("Parts are not disjoint in multipartite Newton polytope estimation: $parts."))
    end
    # Some variables might be in no part...
    missing = setdiff(MP.variables(X), reduce(union, parts))
    if isempty(missing)
        all_parts = parts
    else
        # in that case, create a part with the missing ones
        all_parts = (parts..., missing)
    end
    if length(all_parts) == 1
        # all variables on same part, fallback to shortcut
        return half_newton_polytope(X, tuple(); apply_post_filter = apply_post_filter)
    end
    monovecs = map(vars -> sub_half_newton_polytope(X, vars), all_parts)
    # Cartesian product of the newton polytopes of the different parts
    product = [prod(monos) for monos in Iterators.product(monovecs...)]
    mindeg, maxdeg = cfld(MP.extdegree(X), 2)
    # We know that the degree inequalities are satisfied variable-wise and
    # part-wise but for all variables together so we filter with that
    monos = MP.monovec(filter(mono -> mindeg <= MP.degree(mono) <= maxdeg, product))
    if apply_post_filter
        return post_filter(monos, X)
    else
        return monos
    end
end

function __chip(cur, i, vars, exps, n, op)
    if n > 0
        exp = min(exps[i], n)
        next = iszero(exp) ? cur : op(cur, vars[i]^exp)
        return __chip(next, i + 1, vars, exps, n - exp, op)
    else
        return cur
    end
end
function _chip(mono, n)
    vars = MP.variables(mono)
    exps = MP.exponents(mono)
    if n < 0
        vars = reverse(vars)
        exps = reverse(exps)
        op(a, b) = b * a
    else
        op = *
    end
    return __chip(MP.constantmonomial(mono), 1, vars, exps, abs(n), op)
end

function _is_hermitian_square(mono)
    d = MP.degree(mono)
    if isodd(d)
        return false
    end
    n = div(d, 2)
    return _chip(mono, n) == _chip(mono, -n)
end

# Shortcut for more efficient `extdeg` and `exp` function in case all the
# variables are in the same part
function half_newton_polytope(X::AbstractVector, parts::Tuple{}; apply_post_filter = true)
    vars = MP.variables(X)
    if is_commutative(vars)
        # Commutative variables
        exp(mono, i) = MP.exponents(mono)[i]
        monos = _sub_half_newton_polytope(X, MP.extdegree(X), exp, vars)
    else
        # Non-commutative variables
        # We use Newton chip method of [Section 2.3, BKP16].
        #
        # [BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.
        # *Optimization of polynomials in non-commuting variables*.
        # Berlin: Springer, 2016.
        vars = unique!(sort(vars))
        function ncexp(mono, i)
            mvars = MP.variables(mono)
            return mapreduce(j -> mvars[j] == vars[i] ? MP.exponents(mono)[j] : 0,
                             +, eachindex(mvars), init=0)
        end
        filter = _full_filter(X, MP.extdegree(X), ncexp, length(vars))
        _monos = eltype(X)[]
        for mono in X
            if _is_hermitian_square(mono)
                for i in 1:div(MP.degree(mono), 2)
                    w = _chip(mono, -i)
                    if filter(w)
                        push!(_monos, w)
                    end
                end
            end
        end
        monos = MP.monovec(_monos)
    end
    if apply_post_filter
        return post_filter(monos, X)
    else
        return monos
    end
end

# If `mono` is such that there is no other way to have `mono^2` by multiplying
# two different monomials of `monos` and `mono` is not in `X` then, the corresponding
# diagonal entry of the Gram matrix will be zero hence the whole column and row
# will be zero hence we can remove this monomial.
# See [Theorem 2, L09] or [Section 2.4, BKP16].
#
# [L09] Lofberg, Johan.
# *Pre-and post-processing sum-of-squares programs in practice*.
# IEEE transactions on automatic control 54.5 (2009): 1007-1011.
#
# [BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.
# *Optimization of polynomials in non-commuting variables*.
# Berlin: Springer, 2016.
function post_filter(monos, X)
    num = Dict(mono => 1 for mono in X)
    function _increase(mono)
        num[mono] = get(num, mono, 0) + 1
    end
    function _decrease(mono)
        value = num[mono] - 1
        if iszero(value)
            delete!(num, mono)
        else
            num[mono] = value
        end
        return iszero(value)
    end
    for a in monos
        for b in monos
            if a != b
                _increase(a * b)
            end
        end
    end
    back = Dict{eltype(monos), Int}()
    keep = ones(Bool, length(monos))
    function _delete(i)
        keep[i] = false
        a = monos[i]
        for (j, b) in enumerate(monos)
            if keep[j] && a != b
                for w in [a * b, b * a]
                    if _decrease(w) && haskey(back, w)
                        _delete(back[w])
                    end
                end
            end
        end
    end
    for (i, mono) in enumerate(monos)
        w = mono^2
        if haskey(num, w)
            back[w] = i
        else
            _delete(i)
        end
    end
    return monos[findall(keep)]
end

function monomials_half_newton_polytope(X::AbstractVector, parts; apply_post_filter = true)
    half_newton_polytope(MP.monovec(X), parts; apply_post_filter = apply_post_filter)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
include("stacknet.jl")
include("logger.jl")

function testing()

    function train(X, Y, net::StackNet)
        const MAXITER = 20000
        const BATCHSIZE = 5

        for epoch = 1:MAXITER
            loss = 0
            for i = 1:size(X)[2]
                # every single sample
                outputs = forward(net, X[:, i])
                loss += net.loss_func(outputs, Y[i])
                backward(net, outputs, Y[i])

                if i % BATCHSIZE == 0
                    batch_update(net, eta=0.03)
                    #dbglog(shownet(net))
                end
            end
            batch_update(net, eta=0.03)
            #dbglog("loss $(loss[1])")
        end
    end

    function test(X, Y, net::StackNet)
        for i = 1:size(X)[2]
            # every single sample
            outputs = forward(net, X[:, i])
            dbglog("loss $(net.loss_func(outputs, Y[i]))")

            infolog(outputs[1], Y[i])
        end
    end

    X = [
    1 2 0 0 1 0 0 -1 -2 -1;
    0 0 -1 -2 -1 1 2 0 0 1
    ]

    testX = [
    2 1 -1 -2;
    -1 -2 2 1
    ]

    testY = [ 1 1 -1 -1 ]
    testY = [ 1 1 0 0 ]

    Y = [ 1 1 1 1 1 -1 -1 -1 -1 -1 ]
    Y = [ 1 1 1 1 1 0 0 0 0 0 ]

    net = StackNet(square_error)
    add_layer(net, PerceptronLayer(2, 3, tanh))
    add_layer(net, PerceptronLayer(3, 1, sigmoid))

    dbglog(shownet(net))
    test(testX, testY, net)

    train(X, Y, net)
    dbglog(shownet(net))

    test(testX, testY, net)

end

# homework main execution
function main()
    X = Array[
    [ 1.58,  2.32,  -5.8],   [ 0.67,  1.58,  -4.78],  [ 1.04,  1.01,  -3.63],
    [-1.49,  2.18,  -3.39],  [-0.41,  1.21,  -4.73],  [1.39,  3.16,  2.87],

    [ 0.21,  0.03,  -2.21],   [ 0.37,  0.28,  -1.8],  [ 0.18,  1.22,  0.16],
    [-0.24,  0.93,  -1.01],  [-1.18,  0.39,  -0.39],  [0.74,  0.96,  -1.16],

    [-1.54,  1.17,  0.64],   [5.41,  3.45,  -1.33],  [ 1.55,  0.99,  2.69],
    [1.86,  3.19,  1.51],    [1.68,  1.79,  -0.87],  [3.51,  -0.22,  -1.39],

    ]

    testX = Array[
    [ 1.20,  1.40,  -1.89],  [-0.92,  1.44,  -3.22],  [ 0.45,  1.33,  -4.38],
    [-0.76,  0.84,  -1.96],
    [-0.38,  1.94,  -0.48],  [0.02,  0.72,  -0.17],  [ 0.44,  1.31,  -0.14],
    [ 0.46,  1.49,  0.68],
    [1.40,  -0.44,  -0.92],  [0.44,  0.83,  1.97],  [ 0.25,  0.68,  -0.99],
    [ 0.66,  -0.45,  0.08]
    ]

    Y = Array[
    [1, 0, 0], [1, 0, 0], [1, 0, 0],
    [1, 0, 0], [1, 0, 0], [1, 0, 0],

    [0, 1, 0], [0, 1, 0], [0, 1, 0],
    [0, 1, 0], [0, 1, 0], [0, 1, 0],

    [0, 0, 1], [0, 0, 1], [0, 0, 1],
    [0, 0, 1], [0, 0, 1], [0, 0, 1]
    ]

    testY = Array[
    [1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0],
    [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0],
    [0, 0, 1], [0, 0, 1], [0, 0, 1], [0, 0, 1]
    ]

    X = [X; testX]
    Y = [Y; testY]

    const MAXITER = 30000
    const BATCHSIZE = 2
    const ETA = 0.1
    const HIDDEN = 7

    net = StackNet(square_error, batch_size=BATCHSIZE)
    add_layer(net, PerceptronLayer(3, HIDDEN, tanh))
    add_layer(net, PerceptronLayer(HIDDEN, 3, sigmoid))

    net.layers[1].weights = [0.34220853408533164 0.3245861926046375 0.3332052733100308
    0.22859717196907367 0.46779743149877756 0.3036053965321488
    0.10892493071988898 0.5410961612858103 0.3499789079943008
    0.4423683719806599 0.2032223830720614 0.3544092449472787
    0.270107374555619 0.16764827123181456 0.5622443542125665
    0.02227635058706275 0.8910201267019922 0.08670352271094506
    0.4865634514215009 0.4626169052427591 0.05081964333573993]
    net.layers[1].bias = [0.9952117507065281,0.388087128339186,0.8725765350929733,0.41961076086301596,0.1167031740076856,0.9087224278524852,0.05980293821390692]
    net.layers[2].weights = [0.0526149662305822 0.2031494320476992 0.1346865896996537 0.11584528187906622 0.17054130500350598 0.15808138735561705 0.16508103778387567
    0.24831808672833983 0.2403426545152851 0.11561300113671563 0.10207716354505683 0.14601009172035231 0.0947610652356002 0.052877937118650215
    0.2356175051594146 0.008793297652765042 0.07416296981297327 0.08912410166055519 0.15911479330166178 0.19253821277987473 0.2406491196327555]
    net.layers[2].bias = [0.7230274076027488,0.34008837831362837,0.8944012775586168]

    function train()
        for epoch = 1:MAXITER
            loss = 0
            for i = 1:length(X)

                outputs = forward(net, X[i])
                loss += net.loss_func(outputs, Y[i])
                backward(net, outputs, Y[i])

                if i % BATCHSIZE == 0
                    batch_update(net, eta=ETA)
                end

            end
            batch_update(net, eta=ETA)

            dbglog("train loss: $HIDDEN, $BATCHSIZE, $ETA, $epoch, $loss")
        end
    end

    function test()
        for i = 1:length(testX)
            dbglog("------------\n test sample number: $i $(testX[i])")
            outputs = forward(net, testX[i])
            dbglog("loss $(net.loss_func(outputs, testY[i]))")
            infolog(outputs, testY[i])
        end
    end

    #dbglog(shownet(net))
    dbglog(repr(net))
    #test()

    train()

    #dbglog(shownet(net))
    #test()
end

if !isinteractive()
    main()
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -17.39516631142554 <= q <= 12.050938789121952)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.1892165831624366 + -0.6797303891528506*1/(1 + exp(-(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])))) + -0.5563706493225857*1/(1 + exp(-(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])))) + -0.5881166064998418*1/(1 + exp(-(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])))) + -0.5324749908347566*1/(1 + exp(-(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))))) + 1/(1 + exp(-(-0.7087227741366906 + -0.3339546766337098*1/(1 + exp(-(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])))) + -0.8987724666577948*1/(1 + exp(-(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])))) + 0.3907505378456695*1/(1 + exp(-(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])))) + -0.7983638213649877*1/(1 + exp(-(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))))) + 1/(1 + exp(-(0.16171979073463483 + 0.31171990946083294*1/(1 + exp(-(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])))) + 0.2757921127232579*1/(1 + exp(-(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])))) + -0.732639957443523*1/(1 + exp(-(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])))) + -0.947011017876227*1/(1 + exp(-(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))))) + 1/(1 + exp(-(-0.17479526548317725 + -0.7121643152938684*1/(1 + exp(-(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])))) + -0.20573796119512977*1/(1 + exp(-(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])))) + 0.920351043144791*1/(1 + exp(-(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])))) + -0.8770524767395065*1/(1 + exp(-(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
abstract BamIO <: IO

type BamIOS <: BamIO
    bgzf::Ptr{BGZF}
    h::Ptr{Bam_hdr_t}
    b::Ptr{Bam1_t}
    kstr::Ptr{KStr}
    end_file::Bool
end

function bam_open{T<:DirectIndexString}(path::T, mode::Union{ASCIIString,Char})
    @assert isfile(path)
    bai_path = string(path,".bai")
    @assert isfile(bai_path)

    ppath = get_pointer(path)
    pmode = get_pointer(mode)
    bgzf_is_bgzf(ppath)
    bgzf = bgzf_open(ppath, pmode)
    h = bam_hdr_read(bgzf)
    b = bam_init1()

    kstr = get_ptr_of_null_kstr(100000)

    BamIOS(bgzf, h, b, kstr, false)
end

function eof(bios::BamIOS)
    bios.end_file
end

function close(bios::BamIOS)
    bam_destroy1(bios.b)
    bam_hdr_destroy(bios.h)
    bgzf_close(bios.bgzf)
end

function readline(bios::BamIOS)
    flag = bam_read1!(bios.bgzf, bios.b)
    if flag == -1
        bios.end_file = true
    end
    sam_format1!(bios.h, bios.b, bios.kstr)
    str = kstrToASCII(bios.kstr)
    str = deepcopy(str)
    split(str,"\t")
end

function readlines(bios::BamIOS)
    i = 0
    while !eof(bios)
        i += 1
        i % 1000_000 == 0 && info(i)
        line = readline(bios)
        i % 1000_000 ==0 && println(line)
    end
end

function test_file(file)
    bios = bam_open(file,"rb")
    readlines(bios)
    close(bios)
end

function tests()
    test_file("/home/guo/haplox/private/DarkVC/data/test_Bams/test.bam")
    test_file("/home/guo/haplox/private/DarkVC/data/test_Bams/test2.bam")
end


text:
export ScheduleOnPipelineStage, schedule_updates

import Rocket: release!

"""
    ScheduleOnPipelineStage{S} <: AbstractPipelineStage

Applies the `schedule_on()` operator from `Rocket.jl` library to the given pipeline with a provided `scheduler`

# Arguments
- `scheduler`: scheduler to schedule updates on. Must be compatible with `Rocket.jl` library and `schedule_on()` operator.

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
struct ScheduleOnPipelineStage{S} <: AbstractPipelineStage
    scheduler :: S
end

apply_pipeline_stage(stage::ScheduleOnPipelineStage, factornode, tag, stream) = stream |> schedule_on(stage.scheduler)

Rocket.release!(stage::ScheduleOnPipelineStage)                        = release!(stage.scheduler)
Rocket.release!(stages::NTuple{N, <: ScheduleOnPipelineStage}) where N = foreach(release!, stages)
Rocket.release!(stages::AbstractArray{ <: ScheduleOnPipelineStage })   = foreach(release!, stages)

update!(stage::ScheduleOnPipelineStage)                        = release!(stage.scheduler)
update!(stages::NTuple{N, <: ScheduleOnPipelineStage}) where N = foreach(update!, stages)
update!(stages::AbstractArray{ <: ScheduleOnPipelineStage })   = foreach(update!, stages)

function _schedule_updates end

__schedule_updates(var::AbstractVariable)                         = __schedule_updates((var, ))
__schedule_updates(vars::NTuple{N, <: AbstractVariable }) where N = __schedule_updates(ScheduleOnPipelineStage(PendingScheduler()), vars)
__schedule_updates(vars::AbstractArray{ <: AbstractVariable })    = __schedule_updates(ScheduleOnPipelineStage(PendingScheduler()), vars)

__schedule_updates(pipeline_stage::ScheduleOnPipelineStage, var::AbstractVariable) = __schedule_updates(pipeline_stage, (var, ))

function __schedule_updates(pipeline_stage::ScheduleOnPipelineStage, vars::NTuple{N, <: AbstractVariable }) where N
    foreach((v) -> add_pipeline_stage!(v, pipeline_stage), vars)
    return pipeline_stage
end

function __schedule_updates(pipeline_stage::ScheduleOnPipelineStage, vars::AbstractArray{ <: AbstractVariable })
    foreach((v) -> add_pipeline_stage!(v, pipeline_stage), vars)
    return pipeline_stage
end

"""
    schedule_updates(variables...; pipeline_stage = ScheduleOnPipelineStage(PendingScheduler()))

Schedules posterior marginal updates for given variables using `stage`. By default creates `ScheduleOnPipelineStage` with `PendingScheduler()` from `Rocket.jl` library.
Returns a scheduler with `release!` method available to release all scheduled updates.

See also: [`ScheduleOnPipelineStage`](@ref), [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref)
"""
function schedule_updates(args...; pipeline_stage = ScheduleOnPipelineStage(PendingScheduler()))
    return map((arg) -> __schedule_updates(pipeline_stage, arg), args)
end


text:
#
# This example assumes below packages are available in your environment
#

using Pkg, DrWatson

#@quickactivate "StatisticalRethinkingStan"
using StatsModelComparisons, StanSample
using StatsFuns, CSV, Random
using ParetoSmooth
using StatisticalRethinking

ProjDir = @__DIR__

df = CSV.read(joinpath(ProjDir, "cars.csv"), DataFrame)

cars_stan = "
data {
    int N;
    vector[N] speed;
    vector[N] dist;
}
parameters {
    real a;
    real b;
    real sigma;
}
transformed parameters{
    vector[N] mu;
    mu = a + b * speed;
}
model {
    a ~ normal(0, 100);
    b ~ normal(0, 10);
    sigma ~ exponential(1);
    dist ~ normal(mu, sigma)    ;
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
}
"

Random.seed!(1)
cars_stan_model = SampleModel("cars.model", cars_stan)
data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
rc = stan_sample(cars_stan_model; data)
println()

if success(rc)
    #stan_summary(cars_stan_model, true)
    nt_cars = read_samples(cars_stan_model);
end

log_lik = nt_cars.log_lik';
n_sam, n_obs = size(log_lik)
lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs);

pwaic = [var(log_lik[:, i]) for i in 1:n_obs];
@show -2(sum(lppd) - sum(pwaic))
println()

@show waic(log_lik)
println()

#Random.seed!(123)
loo, loos, pk = psisloo(log_lik);
@show -2loo

#Random.seed!(123)
ll = reshape(nt_cars.log_lik, 50, 1000, 4);
psis_ll = psis(ll);

lwp = deepcopy(ll);
lwp += psis_ll.weights;
lwpt = Matrix(reshape(lwp, 50, 4000)');
loos = reshape(logsumexp(lwpt; dims=1), size(lwpt, 2));

@show loo = sum(loos)
@show 2loo

pk_plot(pk)
savefig(joinpath(ProjDir, "pk.png"))

pk_plot(psis_ll.pareto_k)
savefig(joinpath(ProjDir, "pareto_k.png"))


text:
```
MIT License

Copyright (c) 2021 Qi Bao

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.
```

include("Graph_class.jl")

using DataStructures # contain function: counter
using SpecialFunctions # contains function: logfactorial()
using Random
using SparseArrays
using Distributed

function powerLaw(n; alp = 2.5, xmin = 1)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    for i = 1 : n
        x[i] = xmin * ((1 - x[i])^(-1.0/(alp - 1.0)))
    end
    xm = x[argmax(x)]
    x ./= xm
    return x
end

function Uniform(n)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    return x
end

function Exponential(n; lmd = 1, xmin = 1)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    for i = 1 : n
        x[i] = xmin - (1.0/lmd)*log(1-x[i])
    end
    xm = x[argmax(x)]
    x ./= xm
    return x
end

# (I + L)⁻¹
function getW(G)
	L = zeros(G.n, G.n)
	for (ID, u, v, w) in G.E
		L[u, u] += w
		L[v, v] += w
		L[u, v] -= w
		L[v, u] -= w
	end
	for i = 1 : G.n
		L[i, i] += 1.0
	end
	return inv(L)
end

function getL(G)
	L = zeros(G.n, G.n)
	for (ID, u, v, w) in G.E
		L[u, u] += w
		L[v, v] += w
		L[u, v] -= w
		L[v, u] -= w
	end
	return L
end

function getSparseIpL(G)
    d = ones(G.n)
    for (ID, u, v, w) in G.E
        d[u] += w
        d[v] += w
    end
    Is = zeros(Int32, G.m*2+G.n)
    Js = zeros(Int32, G.m*2+G.n)
    Vs = zeros(G.m*2+G.n)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = -w
        Is[ID + G.m] = v
        Js[ID + G.m] = u
        Vs[ID + G.m] = -w
    end
    for i = 1 : G.n
        Is[G.m + G.m + i] = i
        Js[G.m + G.m + i] = i
        Vs[G.m + G.m + i] = d[i]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function getSparseL(G)
    d = zeros(G.n)
    for (ID, u, v, w) in G.E
        d[u] += w
        d[v] += w
    end
    Is = zeros(Int32, G.m*2+G.n)
    Js = zeros(Int32, G.m*2+G.n)
    Vs = zeros(G.m*2+G.n)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = -w
        Is[ID + G.m] = v
        Js[ID + G.m] = u
        Vs[ID + G.m] = -w
    end
    for i = 1 : G.n
        Is[G.m + G.m + i] = i
        Js[G.m + G.m + i] = i
        Vs[G.m + G.m + i] = d[i]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

## NEW

# (1) Trivalency Model
function Trivalency(w, ratio=5)
    Random.seed!(round(Int, time() * 10000))
    x = rand([0.1, 0.01, 0.001]*ratio, w)
    return x
end

# (2) transform set to vector
@everywhere toarray(s::Union{Set, Vector}) = [toarray.(s)...]
@everywhere toarray(v::Number) = v

"""
    (3) Find the most common element in a vector
"""
function most_common(a)
    c = collect(counter(a))
    pair = c[findmax(map(x->x[2], c))[2]]
    return pair[1], pair[2]
end

# (4) logfactorial
logfact(x, y) = logfact(x) / logfact(y)
logfact(x) = logfactorial(x)

function logchoose(N, k)
    return logfact(N) - logfact(k) - logfact(N-k)
end

# (5) topk, find the top k values in a vector
function topk(a, k)
    sort_a = sortperm(a, rev=true)
    if k <= length(a)
        return sort_a[1:k]
    else
        return sort_a[1:length(a)]
    end
end


text:
using Term.box
import Term: Segment

@testset "\e34mBOX" begin
    testbox = Term.box.Box(
        "ASCII",
        """
        +--+
        | ||
        |-+|
        | ||
        |-+|
        |-+|
        | ||
        +--+
        """,
    )

    @test string(testbox) == "Box\e[2m(ASCII)\e[0m"

    @suppress_out begin
        @test_nothrow println(stdout, testbox)
    end

    _bstring = Term.box.fit(Term.box.ASCII, [1, 1, 1, 1, 1, 1, 1, 1])
    @test _bstring == "+---------------+\n| | | | | | | | |\n|-+-+-+-+-+-+-+-|\n| | | | | | | | |\n|-+-+-+-+-+-+-+-|\n|-+-+-+-+-+-+-+-|\n| | | | | | | | |\n+---------------+"

    # test get row
    @test get_row(ROUNDED, [3], :top) ==  "╭───╮"
    @test get_row(ROUNDED, 3, :top) == "╭─╮"
    @test get_row(ROUNDED, [3, 5], :top) == "╭───┬─────╮"

    # test get title row with no justification
    tr = get_title_row(:top, Term.box.ROUNDED, nothing, width=22, style="red")
    @test typeof(tr) == Segment
    @test tr.measure.w == 22

    # test title row with justification
    left = get_title_row(:top, Term.box.ROUNDED, "test", width=22,
                            justify=:left, title_style="blue", style="red")
    @test left.measure.w == 22
    @test left.text == "\e[31m╭──── \e[34mtest\e[39m\e[31m\e[31m ──────────╮\e[39m\e[0m\e[39m\e[31m"

    right = get_title_row(:top, Term.box.ROUNDED, "test", width=22,
                        justify=:right, title_style="blue", style="red")
    @test right.measure.w == 22
    @test right.text == "\e[31m╭─────────── \e[34mtest\e[39m\e[31m\e[31m ───╮\e[39m\e[0m\e[39m\e[31m"

    center = get_title_row(:top, Term.box.ROUNDED, "test", width=22,
                        justify=:center, title_style="blue", style="red")
    @test center.measure.w == 22
    @test center.text == "\e[31m╭───────\e[39m \e[34mtest\e[39m\e[31m ───────╮\e[39m\e[0m"

    for width in (15, 21, 33, 58), justify in (:left, :center, :right)
        line = get_title_row(:top, Term.box.DOUBLE, "test", width=width,
                        justify=justify, title_style="on_green", style="blue bold")

        @test line.measure.w == width
    end
end