text:
input = joinpath(@__DIR__, "input")
lines = readlines(input)

function cubes(lines)
    acc = []
    for line in lines
        mode = startswith(line, "on") ? true : false
        l, r = parse.(Int, match(r"x=([-]?\d+)..([-]?\d+)", line).captures)
        xi = min(l, r)
        xa = max(l, r)
        l, r = parse.(Int, match(r"y=([-]?\d+)..([-]?\d+)", line).captures)
        yi = min(l, r)
        ya = max(l, r)
        l, r = parse.(Int, match(r"z=([-]?\d+)..([-]?\d+)", line).captures)
        zi = min(l, r)
        za = max(l, r)

        cube = (mode, xi:xa, yi:ya, zi:za)

        push!(acc, cube)
    end
    return acc
end

function boundcheck(cube)
    (_, x, y, z) = cube
    return issubset(x, -50:50) && issubset(y, -50:50) && issubset(z, -50:50)
end

function flip(switches)
    acc = Set()
    for cube in switches
        (mode, x, y, z) = cube
        for switch in collect(Iterators.product(x, y, z))
            if mode
                push!(acc, switch)
            else
                delete!(acc, switch)
            end
        end
    end
    return acc
end

function expand(cubes)
    acc = []
    for cube in cubes
        (emode, ex, ey, ez) = cube
        nacc = []
        for (nmode, nx, ny, nz) in acc
            x = intersect(ex, nx)
            y = intersect(ey, ny)
            z = intersect(ez, nz)
            if length(x) > 0 && length(y) > 0 && length(z) > 0
                push!(nacc, (!nmode, x, y, z))
            end
        end
        if emode
            push!(nacc, cube)
        end
        append!(acc, nacc)
    end
    return acc
end

function value(cube)
    (mode, x, y, z) = cube
    n = length(x) * length(y) * length(z)
    return mode ? n : -1 * n
end

cs = cubes(lines)
cs = filter(boundcheck, cs)
p1 = length(flip(cs))

cs = cubes(lines)
cs = expand(cs)
p2 = sum(value, cs)

println("-----------------------------------------------------------------------")
println("reactor reboot -- part one :: $p1")
println("reactor reboot -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 567496)
@assert(p2 == 1355961721298916)


text:
export update_y_l

"""
  update l and compute y . This is a subfunction for PARSDMM.jl
"""
function update_y_l(
  x                     ::Vector{TF},
  p                     ::Integer,
  i                     ::Integer,
  Blas_active           ::Bool,
  y                     ::Vector{Vector{TF}},
  y_old                 ::Vector{Vector{TF}},
  l                     ::Vector{Vector{TF}},
  l_old                 ::Vector{Vector{TF}},
  rho                   ::Vector{TF},
  gamma                 ::Vector{TF},
  prox                  ::Vector{Any},
  TD_OP                 ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
  log_PARSDMM,
  P_sub,
  counter               ::Integer,
  x_hat                 ::Vector{Vector{TF}},
  r_pri                 ::Vector{Vector{TF}},
  s                     ::Vector{Vector{TF}},
  feasibility_only=false::Bool
  ) where {TF<:Real,TI<:Integer}

#x_hat   = Vector{Vector{Float64}}(p);
#r_pri   = Vector{Vector{Float64}}(p);
#s       = Vector{Vector{Float64}}(p)
lx   = length(x);
rho1 = Vector{TF}(undef,p)
rho1 = TF(1.0) ./ rho;

for ii=1:p
  #y_old[ii]=deepcopy(y[ii])
  #l_old[ii]=deepcopy(l[ii])
  copy!(y_old[ii],y[ii]);
  copy!(l_old[ii],l[ii]);

  s[ii] = TD_OP[ii]*x;
  if Blas_active  #&& i>5
    if gamma[ii]==1 #without relaxation
      copyto!(y[ii],s[ii])
      #y[ii]       = prox[ii]( BLAS.axpy!(-rho1[ii],l[ii],y[ii]) )
      BLAS.axpy!(-rho1[ii],l[ii],y[ii])
      y[ii] = prox[ii](y[ii])
      r_pri[ii]   .= -s[ii] .+ y[ii];
      BLAS.axpy!(rho[ii],r_pri[ii],l[ii]);
    else #relaxed iterations
      x_hat[ii]   .= gamma[ii].*s[ii]
      BLAS.axpy!(TF(1.0) .- gamma[ii],y[ii],x_hat[ii]);
      #y[ii]       = copy(x_hat[ii]);
      copyto!(y[ii],x_hat[ii])
      #y[ii]       = prox[ii]( BLAS.axpy!(-rho1[ii],l[ii],y[ii]))
      BLAS.axpy!(-rho1[ii],l[ii],y[ii])
      y[ii] = prox[ii](y[ii])
      r_pri[ii]   .= -s[ii] .+ y[ii];
      BLAS.axpy!(rho[ii],y[ii] .- x_hat[ii],l[ii]);
    end
  else
    if gamma[ii]==1 #without relaxation
      #y[ii]       = prox[ii]( s[ii].-l[ii].*rho1[ii] );
      y[ii]      .= s[ii] .- l[ii].*rho1[ii]
      y[ii] = prox[ii](y[ii]);
      r_pri[ii]   .= -s[ii] .+ y[ii];
      l[ii]       .= l[ii] .+ rho[ii] .* r_pri[ii];
    else #relaxed iterations
      x_hat[ii]   .= gamma[ii].*s[ii] .+ ( TF(1.0) .- gamma[ii] ) .* y[ii];
      #y[ii]       = prox[ii]( x_hat[ii].-l[ii].*rho1[ii] );
      y[ii]       .= x_hat[ii] .- l[ii] .* rho1[ii]
      y[ii] = prox[ii]( y[ii] );
      r_pri[ii]   .= -s[ii] .+ y[ii];
      l[ii]       .= l[ii] .+ rho[ii] .* ( -x_hat[ii] .+ y[ii] );
    end
  end #end blas

log_PARSDMM.r_pri[i,ii]  = norm(r_pri[ii])#./(norm(y[ii])+(100*eps(TF))); #add 1-14, because y may be 0 for certain initializations of x,y,l
if Blas_active
  log_PARSDMM.r_dual[i,ii] = norm( BLAS.scal!(lx,rho[ii],TD_OP[ii]'*(y[ii] .- y_old[ii]),1 ))
else
 log_PARSDMM.r_dual[i,ii] = norm( rho[ii] .* (TD_OP[ii]'*(y[ii] .- y_old[ii])) )
end
#log feasibility
if feasibility_only==false
  if mod(i,10)==0 && ii<p#log every 10 it, or whatever number is suitable
    temp = P_sub[ii](deepcopy(s[ii]))
    #copy!(x_hat[ii],s[ii])
    log_PARSDMM.set_feasibility[counter,ii] = norm( temp .- s[ii]) ./ (norm(s[ii])+(100*eps(TF))); #use x_hat[ii] as output, because it is already allocated and not used afterwards. The next iteration it will be overwritten anyway before use
  end
else
  if mod(i,10)==0#log every 10 it, or whatever number is suitable
    #copy!(x_hat[ii],s[ii])
    temp = P_sub[ii](deepcopy(s[ii]))
    log_PARSDMM.set_feasibility[counter,ii] = norm( temp .- s[ii]) ./ (norm(s[ii])+(100*eps(TF))); #use x_hat[ii] as output, because it is already allocated and not used afterwards. The next iteration it will be overwritten anyway before use
  end
end

end #end y and l update loop
if mod(i,10)==0
  counter += 1;
end

return y,l,r_pri,s,log_PARSDMM,counter,y_old,l_old

end #end of function

ext: jl


text:
"""
    module ActorInterfacesTests

Test suite for ActorInterfaces.

To test if an actor library is compatible with ActorInterfaces.jl, you have to provide a proxy type,
subtype of ActorLib, that implements functions defined in this suite (prefixed with ex_) to allow
external manipulating of the actor library, namely spawning, sending messages, counting actors and
running until no messages remain to process.

The proxy type must have a default constructor that creates an empty instance of the actor system.
It is allowed to reuse the same instance after cleaning if there is only a single global instance.
No more than one test will be running at a time.

To run the suite, call `run_suite()`` with the proxy type.
"""
module ActorInterfacesTests

export ActorLib, run_suite

abstract type ActorLib end

"""
    function ex_spawn!(lib, bhv, args...) end

Spawn a new actor with the given behavior and args, and return its address.
"""
function ex_spawn!(lib, bhv, args...) end

"""
    function ex_send!(lib, addr, msg...) end

Send a message to the given addr.
"""
function ex_send!(lib, addr, msg...) end

"""
    function ex_actorcount(lib) end

Return the number of spawned actors in the actor system `lib`, not counting errored ones.
"""
function ex_actorcount(lib) end

"""
    function ex_runtofinish(lib) end

Block and run the actor system until all the messages get processed.
"""
function ex_runtofinish(lib) end

include("suite/suite.jl")

end # module


text:
@testset "Transition flux plots" begin

transFluxPlot( mpSim, 12, ("", "A junior"), ("", "B junior"),
    ("", "Reserve junior") )
transFluxPlot( mpSim, 12, ("A junior", "A senior"), ("B junior", "B senior"),
    ("Reserve junior", "A senior"), ("Reserve junior", "B senior") )
transFluxPlot( mpSim, 12, ("A junior", ""), ("B junior", ""),
    ("Reserve junior", "") )
transFluxPlot( mpSim, 12, ("A senior", "Master"), ("B senior", "Master") )
transFluxPlot( mpSim, 12, ("A senior", ""), ("B senior", ""), ("Master", "") )
#=
transFluxPlot( mpSim, 12, ("", "A junior"), ("", "B junior"),
    ("", "Reserve junior"), showPlot=false, filename="plot/transition test 1" )
transFluxPlot( mpSim, 12, ("A junior", "A senior"), ("B junior", "B senior"),
    ("Reserve junior", "A senior"), ("Reserve junior", "B senior"),
    showPlot=false, filename="plot/transition test 2" )
transFluxPlot( mpSim, 12, ("A junior", ""), ("B junior", ""),
    ("Reserve junior", ""), showPlot=false, filename="plot/transition test 3" )
transFluxPlot( mpSim, 12, ("A senior", "Master"), ("B senior", "Master"),
    showPlot=false, filename="plot/transition test 4" )
transFluxPlot( mpSim, 12, ("A senior", ""), ("B senior", ""), ("Master", ""),
    showPlot=false, filename="plot/transition test 5" )
=#
end  # @testset "Transition flux plots"


text:
module DiscriminativeCircuits

using Reexport

using LogicCircuits
# only reexport selectively from LogicCircuits
export pos_literals, neg_literals
# circuit queries
export issmooth, isdecomposable, isstruct_decomposable,
       isdeterministic, iscanonical
# circuit status
export num_edges, num_parameters
# datasets
export twenty_datasets

using ProbabilisticCircuits

include("discriminative/logistic_nodes.jl")
include("discriminative/parameters.jl")

# TODO to fix later after adding regression circuits, right now only have logistic circuit
const DiscriminativeCircuit = LogisticCircuit

include("bitcircuits/param_bit_circuit.jl")
include("bitcircuits/param_bit_circuit_pair.jl")

include("queries/expectation_rec.jl")
include("queries/expectation_graph.jl")
include("queries/expectation_bit.jl")
include("queries/queries.jl")

# include("discriminative/parameters.jl")
include("io/io.jl")

end # module

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
using ACE, Test, Printf, LinearAlgebra, StaticArrays, BenchmarkTools

##
@testset "ACE.jl" begin
    # ------------------------------------------
    #   basic polynomial basis building blocks
    include("polynomials/test_ylm.jl")
    include("testing/test_wigner.jl")
    # include("polynomials/test_rylm.jl")  # not used right now
    include("polynomials/test_transforms.jl")
    include("polynomials/test_orthpolys.jl")

    # --------------------------------------------
    # core permutation-invariant functionality
    include("test_1pbasis.jl")
    include("test_pibasis.jl")
    # include("test_pipot.jl")   # TODO...

    # ------------------------
    #   rotation_invariance
    include("test_cg.jl")
    include("test_symmbasis.jl")
    include("test_euclvec.jl")

    # # ----------------------
    # #   miscallaneous ...
    # include("compat/test_compat_v05.jl")
    # include("compat/test_compat.jl")
    # include("test_any.jl")

    # ----------------------------------
    #    old tests to be re-introduced
    # include("test_real.jl")
    # include("test_orth.jl")
    # include("test_descriptor.jl")
    # include("bonds/test_cylindrical.jl")
    # include("bonds/test_fourier.jl")
    # include("bonds/test_envpairbasis.jl")
end


text:
# Renyi entropies are not defined for nearest neighbor estimators, so need an extra method
# that doesn't feed the q keyword to genentropy
function _transferentropy(joint, ST, T𝒯, T, est::NearestNeighborEntropyEstimator;
        base = 2, q = 1)

    te = genentropy(T𝒯, est, base = base) +
        genentropy(ST, est, base = base) -
        genentropy(T, est, base = base) -
        genentropy(joint, est, base = base)
end


text:
pam_implied_pop_m, pam_implied_pop_f = eqm_consistency(pam.surplus, pam.msingle, pam.fsingle)
nam_implied_pop_m, nam_implied_pop_f = eqm_consistency(nam.surplus, nam.msingle, nam.fsingle)

# ensure convergence to valid equilibrium
@testset "validity of equilibrium population measures" begin
	match_tol = 1e-12
	@test pam.matches ≈ match_matrix(pam.surplus, pam.msingle, pam.fsingle) atol = match_tol
	@test pam_implied_pop_m ≈ pam.mdist atol = match_tol
	@test pam_implied_pop_f ≈ pam.fdist atol = match_tol

	@test nam.matches ≈ match_matrix(nam.surplus, nam.msingle, nam.fsingle) atol = match_tol
	@test nam_implied_pop_m ≈ nam.mdist atol = match_tol
	@test nam_implied_pop_f ≈ nam.fdist atol = match_tol
end

# verify that super/sub-modularity of production implies positive/negative assortative matching.
# sum of diag corners greater than sum of anti-diag corners
@testset "positive and negative assortativity" begin
	@test pam.matches[1,1] + pam.matches[end,end] > pam.matches[1,end] + pam.matches[end,1]
	@test nam.matches[1,1] + nam.matches[end,end] < nam.matches[1,end] + nam.matches[end,1]
end

# symmetry
@testset "symmetry of matches and singles" begin
	@test pam.matches ≈ pam.matches'
	@test pam.msingle ≈ pam.fsingle
end


text:
# constant loss function for comparability with other packages

############################################################################################
### Types
############################################################################################
"""
Constant loss term. Can be used for comparability to other packages.

# Constructor

    SemConstant(;constant_loss, kwargs...)

# Arguments
- `constant_loss::Number`: constant to add to the objective

# Examples
```julia
    my_constant = SemConstant(constant_loss = 42.0)
```

# Interfaces
Analytic gradients and hessians are available.

# Extended help
## Implementation
Subtype of `SemLossFunction`.
"""
struct SemConstant{C} <: SemLossFunction
    c::C
end

############################################################################################
### Constructors
############################################################################################

function SemConstant(;constant_loss, kwargs...)
    return SemConstant(constant_loss)
end

############################################################################################
### methods
############################################################################################

objective!(constant::SemConstant, par, model) = constant.c
gradient!(constant::SemConstant, par, model) = zero(par)
hessian!(constant::SemConstant, par, model) = zeros(eltype(par), length(par), length(par))

############################################################################################
### Recommended methods
############################################################################################

update_observed(loss_function::SemConstant, observed::SemObserved; kwargs...) = loss_function


text:
using PhomericsTest
using Test

@testset "PhomericsTest.jl" begin
    @test f2(1,2) == 6
    @test f2(10,10) == 40
    @test f_diff2(1,2) == 2
    @test f1(1,1) == 2
end


text:
# struct MosaicPlot{T, S} <: OnlineStat{VectorOb}
#     stat::GroupBy{T,CountMap{S,OrderedDict{S,Int}}}
# end
# MosaicPlot(T, S) = MosaicPlot(GroupBy{T}(CountMap(S)))
# nobs(o::MosaicPlot) = nobs(o.stat)
# _fit!(o::MosaicPlot, xy) = _fit!(o.stat, xy)

"""
    Mosaic(T::Type, S::Type)

Data structure for generating a mosaic plot, a comparison between two categorical variables.

# Example

    using OnlineStats, Plots
    x = [rand() > .8 for i in 1:10^5]
    y = rand([1,2,2,3,3,3], 10^5)
    o = fit!(Mosaic(Bool, Int), zip(x, y))
    plot(o)
"""
mutable struct Mosaic{T, C<:CountMap} <: OnlineStat{VectorOb}
    value::OrderedDict{T, C}
    n::Int
end
Mosaic(T::Type, S::Type) = Mosaic(OrderedDict{T, CountMap{S, OrderedDict{S,Int}}}(), 0)
Base.show(io::IO, o::Mosaic{T,S}) where {T, S} = print(io, "Mosaic: $T × $S")
function _fit!(o::Mosaic{T, C}, xy) where {T, S, C<:CountMap{S, OrderedDict{S,Int}}}
    o.n += 1
    if haskey(o.value, first(xy))
        _fit!(o.value[first(xy)], last(xy))
    else
        stat = CountMap(S)
        _fit!(stat, last(xy))
        o.value[first(xy)] = stat
    end
end
value(o::Mosaic) = sort!(o.value)
Base.keys(o::Mosaic) = sort!(collect(keys(o.value)))
subkeys(o::Mosaic) = sort!(mapreduce(x->collect(keys(x)), union, values(o.value)))
_merge!(o::Mosaic, o2::Mosaic) = (o.n += o2.n; merge!(merge!, o.value, o2.value))

@recipe function f(o::Mosaic{T,S}) where {T,S}
    kys = sort!(collect(keys(o.value)))
    n = nobs(o)
    xwidths = [nobs(o.value[ky]) / n for ky in kys]
    xedges = vcat(0.0, cumsum(xwidths))

    subkys = subkeys(o)
    y = zeros(length(subkys), length(kys))
    for (j, ky) in enumerate(kys)
        y[:, j] = probs(o.value[ky], subkys)
        y[:, j] = 1.0 .- vcat(0.0, cumsum(y[1:(end-1), j]))
    end

    seriestype := :bar
    bar_widths := xwidths
    labels := subkys
    xticks := (midpoints(xedges), kys)
    xlim := (0, 1)
    ylim := (0, 1)

    xedges, y'
end


text:
## One-electron integrals
### Overlap matrix elements

function overlap( a::PGBF, b::PGBF )
    s = overlap( a.expn, a.center, a.power, b.expn, b.center, b.power)
    return a.NORM*b.NORM*s
end

# overlap with itself
function overlap(aexpn::Float64, acenter::Tuple3F64, apower::Tuple3I64)
    return overlap(aexpn,acenter,apower,aexpn,acenter,apower)
end

# Low level method
function overlap(aexpn::Float64, acenter::Tuple3F64, apower::Tuple3I64,
    bexpn::Float64, bcenter::Tuple3F64, bpower::Tuple3I64)

    ax, ay, az = acenter[1],acenter[2],acenter[3]
    bx, by, bz = bcenter[1],bcenter[2],bcenter[3]

    aI,aJ,aK = apower[1],apower[2],apower[3]
    bI,bJ,bK = bpower[1],bpower[2],bpower[3]

    gamma = aexpn + bexpn
    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    rab2 = dist2(ax-bx,ay-by,az-bz)
    pre = (pi/gamma)^1.5*exp(-aexpn*bexpn*rab2/gamma)
    wx = overlap1d(aI,bI,px-ax,px-bx,gamma)
    wy = overlap1d(aJ,bJ,py-ay,py-by,gamma)
    wz = overlap1d(aK,bK,pz-az,pz-bz,gamma)
    return pre*wx*wy*wz
end

function gaussian_product_center(a::PGBF,b::PGBF)
    return (a.expn*[a.x,a.y,a.z]+b.expn*[b.x,b.y,b.z])/(a.expn+b.expn)
end

function gaussian_product_center(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,
                  bexpn::Float64,bx::Float64,by::Float64,bz::Float64)
    return (aexpn*[ax,ay,az]+bexpn*[bx,by,bz])/(aexpn+bexpn)
end

function overlap1d(la::Int64,lb::Int64,ax::Float64,bx::Float64,gamma::Float64)
    total = 0
    for i in 0:div(la+lb,2)
        total += binomial_prefactor(2i,la,lb,ax,bx)*factorial2(2i-1)/(2gamma)^i
    end
    return total
end

function binomial_prefactor( s::Int64, ia::Int64, ib::Int64, xpa::Float64, xpb::Float64 )
    total = 0
    for t in 0:s
        if (s-ia) <= t <= ib
          total += binomial(ia,s-t)*binomial(ib,t)*xpa^(ia-s+t)*xpb^(ib-t)
        end
    end
    return total
end

# where should we put this ??
overlap( a::CGBF, b::CGBF ) = contract( overlap, a, b )


text:
struct MadNLPExecutionStats{T} <: AbstractExecutionStats
    status::Status
    solution::StrideOneVector{T}
    objective::T
    constraints::StrideOneVector{T}
    dual_feas::T
    primal_feas::T
    multipliers::StrideOneVector{T}
    multipliers_L::StrideOneVector{T}
    multipliers_U::StrideOneVector{T}
    iter::Int
    counters::NLPModelsCounters
    elapsed_time::Real
end

struct InvalidNumberException <: Exception end
struct NotEnoughDegreesOfFreedomException <: Exception end

MadNLPExecutionStats(ips::InteriorPointSolver) =MadNLPExecutionStats(
    ips.status,view(ips.x,1:get_nvar(ips.nlp)),ips.obj_val,ips.c,
    ips.inf_du, ips.inf_pr,
    ips.l,view(ips.zl,1:get_nvar(ips.nlp)),view(ips.zu,1:get_nvar(ips.nlp)),
    ips.cnt.k, ips.nlp.counters,ips.cnt.total_time)
getStatus(result::MadNLPExecutionStats) = STATUS_OUTPUT_DICT[result.status]

# Utilities
has_constraints(ips) = ips.m != 0

# Print functions -----------------------------------------------------------
function print_init(ips::AbstractInteriorPointSolver)
    @notice(ips.logger,@sprintf("Number of nonzeros in constraint Jacobian............: %8i",get_nnzj(ips.nlp)))
    @notice(ips.logger,@sprintf("Number of nonzeros in Lagrangian Hessian.............: %8i\n",get_nnzh(ips.nlp)))

    num_fixed = length(ips.ind_fixed)
    num_var = get_nvar(ips.nlp) - num_fixed
    num_llb_vars = length(ips.ind_llb)
    num_lu_vars = sum((get_lvar(ips.nlp).!=-Inf).*(get_uvar(ips.nlp).!=Inf)) - num_fixed
    num_uub_vars = length(ips.ind_uub)
    num_eq_cons = sum(get_lcon(ips.nlp).==get_ucon(ips.nlp))
    num_ineq_cons = sum(get_lcon(ips.nlp).!=get_ucon(ips.nlp))
    num_ue_cons = sum((get_lcon(ips.nlp).!=get_ucon(ips.nlp)).*(get_lcon(ips.nlp).==-Inf).*(get_ucon(ips.nlp).!=Inf))
    num_le_cons = sum((get_lcon(ips.nlp).!=get_ucon(ips.nlp)).*(get_lcon(ips.nlp).!=-Inf).*(get_ucon(ips.nlp).==Inf))
    num_lu_cons = sum((get_lcon(ips.nlp).!=get_ucon(ips.nlp)).*(get_lcon(ips.nlp).!=-Inf).*(get_ucon(ips.nlp).!=Inf))
    get_nvar(ips.nlp) < num_eq_cons && throw(NotEnoughDegreesOfFreedomException())

    @notice(ips.logger,@sprintf("Total number of variables............................: %8i",num_var))
    @notice(ips.logger,@sprintf("                     variables with only lower bounds: %8i",num_llb_vars))
    @notice(ips.logger,@sprintf("                variables with lower and upper bounds: %8i",num_lu_vars))
    @notice(ips.logger,@sprintf("                     variables with only upper bounds: %8i",num_uub_vars))
    @notice(ips.logger,@sprintf("Total number of equality constraints.................: %8i",num_eq_cons))
    @notice(ips.logger,@sprintf("Total number of inequality constraints...............: %8i",num_ineq_cons))
    @notice(ips.logger,@sprintf("        inequality constraints with only lower bounds: %8i",num_le_cons))
    @notice(ips.logger,@sprintf("   inequality constraints with lower and upper bounds: %8i",num_lu_cons))
    @notice(ips.logger,@sprintf("        inequality constraints with only upper bounds: %8i\n",num_ue_cons))
    return
end

function print_iter(ips::AbstractInteriorPointSolver;is_resto=false)
    mod(ips.cnt.k,10)==0&& @info(ips.logger,@sprintf(
        "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls"))
    @info(ips.logger,@sprintf(
        "%4i%s% 10.7e %6.2e %6.2e %5.1f %6.2e %s %6.2e %6.2e%s  %i",
        ips.cnt.k,is_resto ? "r" : " ",ips.obj_val/ips.obj_scale[],
        is_resto ? ips.RR.inf_pr_R : ips.inf_pr,
        is_resto ? ips.RR.inf_du_R : ips.inf_du,
        is_resto ? log(10,ips.RR.mu_R) : log(10,ips.mu),
        ips.cnt.k == 0 ? 0. : norm(primal(ips.d),Inf),
        ips.del_w == 0 ? "   - " : @sprintf("%5.1f",log(10,ips.del_w)),
        ips.alpha_z,ips.alpha,ips.ftype,ips.cnt.l))
    return
end

function print_summary_1(ips::AbstractInteriorPointSolver)
    @notice(ips.logger,"")
    @notice(ips.logger,"Number of Iterations....: $(ips.cnt.k)\n")
    @notice(ips.logger,"                                   (scaled)                 (unscaled)")
    @notice(ips.logger,@sprintf("Objective...............:  % 1.16e   % 1.16e",ips.obj_val,ips.obj_val/ips.obj_scale[]))
    @notice(ips.logger,@sprintf("Dual infeasibility......:   %1.16e    %1.16e",ips.inf_du,ips.inf_du/ips.obj_scale[]))
    @notice(ips.logger,@sprintf("Constraint violation....:   %1.16e    %1.16e",norm(ips.c,Inf),ips.inf_pr))
    @notice(ips.logger,@sprintf("Complementarity.........:   %1.16e    %1.16e",
                                ips.inf_compl*ips.obj_scale[],ips.inf_compl))
    @notice(ips.logger,@sprintf("Overall NLP error.......:   %1.16e    %1.16e\n",
                                max(ips.inf_du*ips.obj_scale[],norm(ips.c,Inf),ips.inf_compl),
                                max(ips.inf_du,ips.inf_pr,ips.inf_compl)))
    return
end

function print_summary_2(ips::AbstractInteriorPointSolver)
    ips.cnt.solver_time = ips.cnt.total_time-ips.cnt.linear_solver_time-ips.cnt.eval_function_time
    @notice(ips.logger,"Number of objective function evaluations             = $(ips.cnt.obj_cnt)")
    @notice(ips.logger,"Number of objective gradient evaluations             = $(ips.cnt.obj_grad_cnt)")
    @notice(ips.logger,"Number of constraint evaluations                     = $(ips.cnt.con_cnt)")
    @notice(ips.logger,"Number of constraint Jacobian evaluations            = $(ips.cnt.con_jac_cnt)")
    @notice(ips.logger,"Number of Lagrangian Hessian evaluations             = $(ips.cnt.lag_hess_cnt)")
    @notice(ips.logger,@sprintf("Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  = %6.3f",
                                ips.cnt.solver_time))
    @notice(ips.logger,@sprintf("Total wall-clock secs in linear solver                      = %6.3f",
                                ips.cnt.linear_solver_time))
    @notice(ips.logger,@sprintf("Total wall-clock secs in NLP function evaluations           = %6.3f",
                                ips.cnt.eval_function_time))
    @notice(ips.logger,@sprintf("Total wall-clock secs                                       = %6.3f\n",
                                ips.cnt.total_time))
end

function print_ignored_options(logger,option_dict)
    @warn(logger,"The following options are ignored: ")
    for (key,val) in option_dict
        @warn(logger," - "*string(key))
    end
end
function string(ips::AbstractInteriorPointSolver)
    """
                Interior point solver

                number of variables......................: $(get_nvar(ips.nlp))
                number of constraints....................: $(get_ncon(ips.nlp))
                number of nonzeros in lagrangian hessian.: $(get_nnzh(ips.nlp))
                number of nonzeros in constraint jacobian: $(get_nnzj(ips.nlp))
                status...................................: $(ips.status)
                """
end
print(io::IO,ips::AbstractInteriorPointSolver) = print(io, string(ips))
show(io::IO,ips::AbstractInteriorPointSolver) = print(io,ips)


text:
"""
    AbstractDistributionFitOptimizer
    OptimOptimizer <: AbstractDistributionFitOptimizer
    set_optimizer(::AbstractDistributionFitOptimizer)

`DistributionFits.jl` uses the following interface to opimize an univariate
function `f` on bounded interval `[lower,upper]`:

    optimize(f, ::AbstractDistributionFitOptimizer, lower, upper)

Returning an object with fields `minimizer` and `converged`.

`DistributionFits.set_optimizer(::AbstractDistributionFitOptimizer)`
sets the optimizer used in calling the optimize function.
Specializing this function with a concrete type, allows using different
optimization packages.

By default, when the [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/)
package is in scope, the
`OptimOptimizer` is set, which implements the interface by using Optim's
[optimize](https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/#minimizing-a-univariate-function-on-a-bounded-interval) function.

Hence, the module using DistribuitonFits.jl has to either

- explicitly invoke `using Optim`, or
- call `set_optimizer` with the concrete subtype of `AbstractDistributionFitOptimizer`
  for which the corresponding optimize method is implemented.
"""
abstract type AbstractDistributionFitOptimizer end
struct OptimOptimizer <:  AbstractDistributionFitOptimizer; end
struct NotSetOptimizer <: AbstractDistributionFitOptimizer; end

optimize(f, ::NotSetOptimizer, lower, upper) = error(
    "Optimizer not set yet. Either invoke 'using Optim' or 'DistributionFits.set_optimizer(...)'.")

optimizer = NotSetOptimizer();
optimize(f, lower, upper) = optimize(f, optimizer, lower, upper)
set_optimizer(opt::AbstractDistributionFitOptimizer) = (global optimizer = opt)


text:
#dof_satter.jl

function gradc(lmm::LMM{T}, theta) where T
    if !lmm.result.fit error("Model not fitted!") end
    if !isnothing(lmm.result.grc) return lmm.result.grc end
    vloptf(x) = sweep_β_cov(lmm, x, lmm.result.beta)
    chunk  = ForwardDiff.Chunk{1}()
    jcfg   = ForwardDiff.JacobianConfig(vloptf, theta, chunk)
    jic    = ForwardDiff.jacobian(vloptf, theta, jcfg, Val{false}())
    grad   = Vector{Matrix{T}}(undef, thetalength(lmm))
    for i in 1:thetalength(lmm)
        gic     = reshape(view(jic, :, i), rankx(lmm), rankx(lmm)) #<Opt
        grad[i] = - lmm.result.c * gic * lmm.result.c
    end
    lmm.result.grc = grad
    grad
end

function getinvhes(lmm::LMM{T}) where T
    local A
    if isnothing(lmm.result.h)
        lmm.result.h = hessian(lmm)
        H = copy(lmm.result.h)
    else
        H = copy(lmm.result.h)
    end
    theta = copy(lmm.result.theta)
    qrd   = qr(H, Val(true))
    vals  = falses(thetalength(lmm))
    for i = 1:thetalength(lmm)
        if lmm.covstr.ct[qrd.jpvt[i]] == :var
            if abs(qrd.R[i, i]) > 1E-8
                vals[qrd.jpvt[i]] = true
            else
                theta[qrd.jpvt[i]] = zero(T)
                H[:,qrd.jpvt[i]]  .= zero(T)
                H[qrd.jpvt[i],:]  .= zero(T)
            end
        elseif lmm.covstr.ct[qrd.jpvt[i]] == :rho
            if 1.0 - abs(lmm.result.theta[qrd.jpvt[i]])  > 1E-6
                vals[qrd.jpvt[i]] = true
            else
                if lmm.result.theta[qrd.jpvt[i]] > 0 lmm.result.theta[qrd.jpvt[i]] = 1.0 else lmm.result.theta[qrd.jpvt[i]] = -1.0 end
                H[:,qrd.jpvt[i]] .= zero(T)
                H[qrd.jpvt[i],:] .= zero(T)
            end
        end
    end
    try
        vh  = view(H, vals, vals)
        vh .= inv(Matrix(vh))
        A = H * 2
    catch
        A = pinv(H) * 2
    end
    A, theta
end
"""
    dof_satter(lmm::LMM{T}, l) where T

Return Satterthwaite approximation for the denominator degrees of freedom, where `l` is a contrast vector (estimable linear combination
of fixed effect coefficients vector (`β`).

```math
df = \\frac{2(LCL')^{2}}{g'Ag}
```

Where: ``A = 2H^{-1}``, ``g = \\triangledown_{\\theta}(LC^{-1}_{\\theta}L')``

"""
function dof_satter(lmm::LMM{T}, l::AbstractVector) where T
    A, theta = getinvhes(lmm)
    grad  = gradc(lmm, theta)
    g  = Vector{T}(undef, length(grad))
    for i = 1:length(grad)
        #g[i] = (l' * grad[i] * l)[1]
        g[i] = mulαtβα(l, grad[i])
    end
    #d = g' * A * g
    d = mulαtβα(g, A)
    df = 2*(mulαtβα(l, lmm.result.c))^2 / d
    if df < 1.0 return 1.0 elseif df > dof_residual(lmm) return dof_residual(lmm) else return df end
end
"""
    dof_satter(lmm::LMM{T}, n::Int) where T

Return Satterthwaite approximation for the denominator degrees of freedom, where `n` - coefficient number.
"""
function dof_satter(lmm::LMM{T}, n::Int) where T
    l = zeros(Int, length(lmm.result.beta))
    l[n] = 1
    dof_satter(lmm, l)
end
"""
    dof_satter(lmm::LMM{T}) where T

Return Satterthwaite approximation for the denominator degrees of freedom for all coefficients.

"""
function dof_satter(lmm::LMM{T}) where T
    isfitted(lmm) || error("Model not fitted")
    lb       = length(lmm.result.beta)
    A, theta = getinvhes(lmm)
    grad     = gradc(lmm, theta)
    dof      = Vector{Float64}(undef, lb)
    for gi = 1:lb
        l     = zeros(Int, lb)
        l[gi] = 1
        g     = Vector{T}(undef, length(grad))
        for i = 1:length(grad)
            #g[i] = (l' * grad[i] * l)[1]
            g[i] = mulαtβα(l, grad[i])
        end
        #d = g' * A * g
        d = mulαtβα(g, A)
        df = 2*(mulαtβα(l, lmm.result.c))^2 / d
        if df < 1.0 dof[gi] = 1.0 elseif df > dof_residual(lmm) dof[gi] = dof_residual(lmm) else dof[gi] = df end
    end
    dof
end

"""
    dof_satter(lmm::LMM{T}, l::Matrix) where T

Return Satterthwaite approximation for the denominator degrees of freedom for conrast matrix `l`.

For `size(l, 1)` > 1:

```math
df = \\frac{2E}{E - rank(LCL')}
```

where:

* let ``LCL' = QΛQ^{-1}``, where ``QΛQ^{-1}`` - spectral decomposition of ``LCL'``
* ``Lq_i`` is the i-th row of ``Q^{-1}L``
* ``A = 2H^{-1}``, ``g = \\triangledown_{\\theta}(Lq_i C^{-1}_{\\theta} Lq_i')``
* ``v_i = \\frac{2*Λ_{i,i}^2}{g' * A * g}``
* ``E = \\sum_{i=1}^n {\\frac{v_i}(v_i - 2)}`` for ``v_i > 2``
"""
function dof_satter(lmm::LMM{T}, l::AbstractMatrix) where T
    if lmm.rankx != size(l, 2) error("size(l, 2) not equal rank X!") end
    A, theta = getinvhes(lmm)
    grad  = gradc(lmm, theta)
    g     = Vector{T}(undef, length(grad))
    lcl   = l*lmm.result.c*l'
    lclr  = rank(lcl)
    #if lclr != size(l, 1) error() end
    lcle  = eigen(lcl)
    pl    = lcle.vectors'*l
    vm    = Vector{T}(undef, lclr)
    em    = 0
    for i = 1:lclr
        plm = pl[i,:]
        for i2 = 1:length(grad)
            #g[i2] = (plm' * grad[i2] * plm)[1]
            g[i2] = mulαtβα(plm, grad[i2])
        end
        #d = g' * A * g
        d = mulαtβα(g, A)
        vm[i] = 2*lcle.values[i]^2 / d
        if vm[i] > 2.0 em += vm[i] / (vm[i] - 2.0) end
    end
    df = 2em/(em - lclr)
    if df < 1.0 return 1.0 elseif df > dof_residual(lmm) return dof_residual(lmm) else return df end
end


text:
module Politics

using MySQL
using ProgressMeter
using Requests

include("consts.jl")
include("politicians.jl")

if length(ARGS) > 0
  api_key = ARGS[1]
  conn = mysql_connect("localhost", "root", "cityhub", "cityhub")

  progress = Progress(length(NYC_ZIPCODES), 1, "Retrieving politicians...", 32)

  for zipcode in NYC_ZIPCODES
      get_politicians(api_key, conn, zipcode)
      next!(progress)
  end

  mysql_disconnect(conn)
else
  println("A key for the Google Civic Information API must be passed as an argument in order to retrieve politician data.")
end

end


text:
__precompile__()
module Mosek
  if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
  else
    error("Mosek not properly installed. Please run Pkg.build(\"Mosek\")")
  end

  import MathProgBase

  export
    makeenv, maketask, maketask_ptr,
    MosekError

  # A macro to make calling C API a little cleaner
  macro msk_ccall(func, args...)
    f = Base.Meta.quot(Symbol("MSK_$(func)"))
    args = [esc(a) for a in args]
    quote
      ccall(($f,libmosek), $(args...))
    end
  end

  # -----
  # Types
  # -----
  type MosekError <: Exception
    rcode :: Int32
    msg   :: String
  end

  # Environment: typedef void * MSKenv_t;
  type MSKenv
    env::Ptr{Void}
    streamcallbackfunc::Any
  end

  # Task: typedef void * MSKtask_t;
  type MSKtask
    env::MSKenv
    task::Ptr{Void}
    borrowed::Bool
    # need to keep a reference to callback funcs for GC
    streamcallbackfunc:: Any
    userstreamcallbackfunc:: Any
    callbackfunc:: Any
    usercallbackfunc:: Any
    nlinfo:: Any

    function MSKtask(env::MSKenv)
      temp = Array(Ptr{Void}, 1)
      res = @msk_ccall(maketask, Int32, (Ptr{Void}, Int32, Int32, Ptr{Void}), env.env, 0, 0, temp)

      if res != MSK_RES_OK
        throw(MosekError(res,""))
      end

      task = new(env,temp[1],false,nothing,nothing,nothing,nothing,nothing)
      task.callbackfunc = cfunction(msk_info_callback_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
      task.usercallbackfunc = nothing
      finalizer(task,deletetask)

      r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), task.task, task.callbackfunc, task)
      if r != MSK_RES_OK
          throw(MosekError(r,getlasterror(t)))
      end

      task
    end

    function MSKtask(t::MSKtask)
      temp = Array(Ptr{Void}, 1)
      res = @msk_ccall(clonetask, Int32, (Ptr{Void}, Ptr{Void}), t.task, temp)

      if res != MSK_RES_OK
        throw(MosekError(res,""))
      end

      task = new(env,temp[1],false,nothing,nothing,nothing,nothing,nothing)
      task.callbackfunc = cfunction(msk_info_callback_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
      task.usercallbackfunc = nothing

      finalizer(task,deletetask)

      r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), task.task, task.callbackfunc, task)
      if r != MSK_RES_OK
          throw(MosekError(r,getlasterror(t)))
      end

      task
    end

    function MSKtask(t::Ptr{Void},borrowed::Bool)
      task = new(msk_global_env,t,borrowed,nothing,nothing,nothing,nothing,nothing)
      task.callbackfunc = cfunction(msk_info_callback_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
      task.usercallbackfunc = nothing

      finalizer(task,deletetask)

      r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), task.task, task.callbackfunc, task)
      if r != MSK_RES_OK
          throw(MosekError(r,getlasterror(t)))
      end

      task
    end
  end

  # ------------
  # API wrappers
  # ------------
  # TODO: Support other argument
  function makeenv()
    temp = Array(Ptr{Void}, 1)
    res = @msk_ccall(makeenv, Int32, (Ptr{Ptr{Void}}, Ptr{UInt8}), temp, C_NULL)
    if res != 0
      # TODO: Actually use result code
      error("MOSEK: Error creating environment")
    end
    MSKenv(temp[1],nothing)
  end

  function makeenv(func::Function)
      temp = Array(Ptr{Void}, 1)
      res = @msk_ccall(makeenv, Int32, (Ptr{Ptr{Void}}, Ptr{UInt8}), temp, C_NULL)
      if res != 0
          # TODO: Actually use result code
          error("MOSEK: Error creating environment")
      end
      env = MSKenv(temp[1],nothing)

      try
          func(env)
      finally
          deleteenv(env)
      end
  end

  # Note on initialization of msk_global_env:
  #
  #  When loading Mosek from source this works fine, but when loading
  #  precompiled module, makeenv() is not called (and some garbage
  #  value is put in msk_global_env). It appears that static
  #  initializers must be called from __init__(). That is called a bad solution here:
  #    https://github.com/JuliaLang/julia/issues/12010
  msk_global_env = makeenv() :: MSKenv
  __init__() = (global msk_global_env = makeenv())

  function maketask(env::MSKenv)

      MSKtask(env)
  end

  function maketask(task::MSKtask)
    MSKtask(task)
  end

  function maketask()
    MSKtask(msk_global_env)
  end

  function maketask_ptr(t::Ptr{Void},borrowed::Bool)
    MSKtask(t,borrowed)
  end

  function maketask(func::Function, env::MSKenv)
      t = MSKtask(env)
      try
          func(t)
      finally
          deletetask(t)
      end
  end

  function maketask(func::Function)
      t = MSKtask(msk_global_env)
      try
          func(t)
      finally
          deletetask(t)
      end
  end

  function deletetask(t::MSKtask)
    if t.task != C_NULL
        if ! t.borrowed
            temp = Array(Ptr{Void},1)
            temp[1] = t.task
            @msk_ccall(deletetask,Int32,(Ptr{Ptr{Void}},), temp)
        end
        t.task = C_NULL
    end
  end

  function deleteenv(e::MSKenv)
    if e.env != C_NULL
      temp = Array(Ptr{Void},1)
      temp[1] = t.env
      @msk_ccall(deleteenv,Int32,(Ptr{Ptr{Void}},), temp)
      e.env = C_NULL
    end
  end

  function getlasterror(t::MSKtask)
    lasterrcode = Array(Cint,1)
    lastmsglen = Array(Cint,1)

    @msk_ccall(getlasterror,Cint,(Ptr{Void},Ptr{Cint},Cint,Ptr{Cint},Ptr{UInt8}),
               t.task, lasterrcode, 0, lastmsglen, C_NULL)
    lastmsg = Array(UInt8,lastmsglen[1])
    @msk_ccall(getlasterror,Cint,(Ptr{Void},Ptr{Cint},Cint,Ptr{Cint},Ptr{UInt8}),
               t.task, lasterrcode, lastmsglen[1], lastmsglen, lastmsg)
    convert(String,lastmsg[1:lastmsglen[1]-1])
  end

  include("msk_enums.jl")
  include("msk_functions.jl")
  include("msk_callback.jl")
  include("msk_geco.jl")

  include("MosekSolverInterface.jl")
  using Mosek.MosekMathProgSolverInterface
  export MosekSolver
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 VirtualMachineCaptureParameters <: SwaggerModel
    vhdPrefix::Any # spec type: Union{ Nothing, String } # spec name: vhdPrefix
    destinationContainerName::Any # spec type: Union{ Nothing, String } # spec name: destinationContainerName
    overwriteVhds::Any # spec type: Union{ Nothing, Bool } # spec name: overwriteVhds

    function VirtualMachineCaptureParameters(;vhdPrefix=nothing, destinationContainerName=nothing, overwriteVhds=nothing)
        o = new()
        validate_property(VirtualMachineCaptureParameters, Symbol("vhdPrefix"), vhdPrefix)
        setfield!(o, Symbol("vhdPrefix"), vhdPrefix)
        validate_property(VirtualMachineCaptureParameters, Symbol("destinationContainerName"), destinationContainerName)
        setfield!(o, Symbol("destinationContainerName"), destinationContainerName)
        validate_property(VirtualMachineCaptureParameters, Symbol("overwriteVhds"), overwriteVhds)
        setfield!(o, Symbol("overwriteVhds"), overwriteVhds)
        o
    end
end # type VirtualMachineCaptureParameters

const _property_map_VirtualMachineCaptureParameters = Dict{Symbol,Symbol}(Symbol("vhdPrefix")=>Symbol("vhdPrefix"), Symbol("destinationContainerName")=>Symbol("destinationContainerName"), Symbol("overwriteVhds")=>Symbol("overwriteVhds"))
const _property_types_VirtualMachineCaptureParameters = Dict{Symbol,String}(Symbol("vhdPrefix")=>"String", Symbol("destinationContainerName")=>"String", Symbol("overwriteVhds")=>"Bool")
Base.propertynames(::Type{ VirtualMachineCaptureParameters }) = collect(keys(_property_map_VirtualMachineCaptureParameters))
Swagger.property_type(::Type{ VirtualMachineCaptureParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineCaptureParameters[name]))}
Swagger.field_name(::Type{ VirtualMachineCaptureParameters }, property_name::Symbol) =  _property_map_VirtualMachineCaptureParameters[property_name]

function check_required(o::VirtualMachineCaptureParameters)
    (getproperty(o, Symbol("vhdPrefix")) === nothing) && (return false)
    (getproperty(o, Symbol("destinationContainerName")) === nothing) && (return false)
    (getproperty(o, Symbol("overwriteVhds")) === nothing) && (return false)
    true
end

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


text:
using ConstrainedDynamics

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1
box = Box(width, depth, length1, length1, color = RGBA(1., 1., 0.))

p2 = [0.0;0.0;length1 / 2] # joint connection point

# Initial orientation
ϕ1 = 0
q1 = Quaternion(RotX(ϕ1))

# Links
origin = Origin{Float64}()
link1 = Body(box)

# Constraints
joint_between_origin_and_link1 = EqualityConstraint(Revolute(origin, link1, zeros(3), p2, joint_axis))

links = [link1]
constraints = [joint_between_origin_and_link1]
shapes = [box]

mech = Mechanism(origin, links, constraints, shapes = shapes, tend = 10.)
setPosition!(mech,origin,link1,p2 = p2,Δq = q1)

pid = PID(2, pi / 2, P = 10., I = 10., D = 5.)

simulate!(mech,pid,save = true)
visualize!(mech)


text:
############################################################
# judiModeling ###############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiModeling, judiModelingException, judiSetupModeling, judiModelingAdjoint, judiModelingAdjointException, judiSetupModelingAdjoint, subsample

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections
struct judiModeling{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    options::Options
	fop::Function              # forward
	fop_T::Function  # transpose
end

mutable struct judiModelingException <: Exception
    msg :: String
end

struct judiModelingAdjoint{DDT,RDT} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    options::Options
	fop::Function              # forward
	fop_T::Function  # transpose
end

mutable struct judiModelingAdjointException <: Exception
    msg :: String
end

############################################################
## outer constructors
"""
    judiModeling(info, model; options=Options())
    judiModeling(info, model, src_geometry, rec_geometry; options=Options())

Create seismic modeling operator for a velocity model given as a `Model` structure. `info` is an `Info` structure\\
containing necessary dimensions to set up the operator. The function also takes the source and receiver geometries\\
as additional input arguments, which creates a combined operator `judiProjection*judiModeling*judiProjection'`.

Example
=======

`Pr` and `Ps` are projection operatos of type `judiProjection` and\\
`q` is a data vector of type `judiVector`:

    F = judiModeling(info, model)
    dobs = Pr*F*Ps'*q

    F = judiModeling(info, model, q.geometry, rec_geometry)
    dobs = F*q

"""
function judiModeling(info::Info, model::Model; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
	(DDT == Float32 && RDT == Float32) || throw(judiModelingException("Domain and range types not supported"))
	m = info.n * sum(info.nt)
	n = m
	srcnum = 1:info.nsrc

	return F = judiModeling{Float32,Float32}("forward wave equation", m, n, info, model, options,
							  args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options),
							  args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options)
							  )
end

function judiModelingAdjoint(info::Info, model::Model; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
	(DDT == Float32 && RDT == Float32) || throw(judiModelingAdjointException("Domain and range types not supported"))
	m = info.n * sum(info.nt)
	n = m
	srcnum = 1:info.nsrc

	return F = judiModelingAdjoint{Float32,Float32}("adjoint wave equation", m, n, info, model, options,
							  args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options),
							  args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options)
							  )
end

############################################################
## overloaded Base functions

# conj(jo)
conj(A::judiModeling{DDT,RDT}) where {DDT,RDT} =
    judiModeling{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.options,A.fop, A.fop_T)

# transpose(jo)
transpose(A::judiModeling{DDT,RDT}) where {DDT,RDT} =
    judiModelingAdjoint{DDT,RDT}("adjoint wave equation",A.n,A.m,A.info,A.model,A.options,A.fop_T, A.fop)

adjoint(A::judiModeling{DDT,RDT}) where {DDT,RDT} =
	judiModelingAdjoint{DDT,RDT}("adjoint wave equation",A.n,A.m,A.info,A.model,A.options,A.fop_T, A.fop)

# conj(jo)
conj(A::judiModelingAdjoint{DDT,RDT}) where {DDT,RDT} =
    judiModelingAdjoint{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.options,A.fop, A.fop_T)

# transpose(jo)
transpose(A::judiModelingAdjoint{DDT,RDT}) where {DDT,RDT} =
    judiModeling{DDT,RDT}("forward wave equation",A.n,A.m,A.info,A.model,A.options,A.fop_T, A.fop)

adjoint(A::judiModelingAdjoint{DDT,RDT}) where {DDT,RDT} =
	judiModeling{DDT,RDT}("forward wave equation",A.n,A.m,A.info,A.model,A.options,A.fop_T, A.fop)

############################################################
## Additional overloaded functions

# *(judiModeling,judiWavefield)
function *(A::judiModeling{ADDT,ARDT}, v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiModelingException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiModeling,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (nothing,v.data,nothing,nothing,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiModeling,judiWavefield):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(judiModelingAdjoint,judiWavefield)
function *(A::judiModelingAdjoint{ADDT,ARDT}, v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiModelingAdjointException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiModeling,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (nothing,nothing,nothing,v.data,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiModeling,judiWavefield):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(judiModeling,judiRHS)
function *(A::judiModeling{ADDT,ARDT}, v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiModelingException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiModeling,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (v.geometry,v.data,nothing,nothing,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiModeling,judiRHS):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(judiModelingAdjoint,judiRHS)
function *(A::judiModelingAdjoint{ADDT,ARDT}, v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiModelingException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiModeling,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (nothing, nothing, v.geometry,v.data,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiModeling,judiRHS):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

############################################################
## Additional overloaded functions

# Subsample Modeling function
function subsample(F::judiModeling, srcnum)
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiModeling(info, F.model, options=F.options)
end

function subsample(F::judiModelingAdjoint, srcnum)
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiModelingAdjoint(info, F.model, options=F.options)
end

getindex(F::judiModeling,a) = subsample(F,a)
getindex(F::judiModelingAdjoint,a) = subsample(F,a)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# Autogenerated wrapper script for CMake_jll for x86_64-linux-musl-cxx11
export cmake

using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("CMake")
JLLWrappers.@declare_executable_product(cmake)
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        cmake,
        "bin/cmake",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MomentClosure
using MomentClosure: define_μ
using ModelingToolkit: value
using Catalyst
using Test

@parameters t, σ_b, σ_u, ρ_b, ρ_u, σ_b_LMA
@variables g(t)

# simple feedback loop
rn_nonlinear = @reaction_network begin
      σ_b, g + p → 0
      σ_u*(1-g), 0 ⇒ g + p
      ρ_u, g → g + p
      ρ_b*(1-g), 0 ⇒ p
      1, p → 0
end σ_b σ_u ρ_b ρ_u

rn_linear = @reaction_network begin
      σ_b_LMA, g → 0
      σ_u*(1-g), 0 ⇒ g
      ρ_u, g → g+p
      (ρ_b*(1-g)), 0 ⇒ p
      1, p → 0
end σ_b_LMA σ_u ρ_b ρ_u

binary_vars = [speciesmap(rn_nonlinear)[g]]
LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false)

μ = define_μ(2,3)
expr1 = LMA_eqs.odes.eqs[1].rhs
expr2 = σ_u - σ_b*μ[1,1] - σ_u*μ[1,0]
@test isequal(expr1, value.(expr2))
expr1 = LMA_eqs.odes.eqs[2].rhs
expr2 = ρ_b + ρ_u*μ[1,0] - ρ_b*μ[1,0] - μ[0,1]
@test isequal(expr1, value.(expr2))
expr1 = LMA_eqs.odes.eqs[3].rhs
expr2 = ρ_u*μ[1,0] + σ_u*μ[0,1] - μ[1,1] - σ_u*μ[1,1] - σ_b*μ[1,1]^2/μ[1,0]
@test isequal(expr1, value.(expr2))

# cooperativity cp=2
rn_nonlinear = @reaction_network begin
      σ_b, g + 2*p → 0
      σ_u*(1-g), 0 ⇒ g + 2*p
      ρ_u, g → g + p
      ρ_b*(1-g), 0 ⇒ p
      1, p → 0
end σ_b σ_u ρ_b ρ_u

binary_vars = [speciesmap(rn_nonlinear)[g]]
_, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false)
expr1 = effective_params[value(σ_b_LMA)]
expr2 = (σ_b*μ[1,2] - σ_b*μ[1,1]) / μ[1,0]
@test isequal(expr1, value.(expr2))


text:
# ---
# title: 1318. Minimum Flips to Make a OR b Equal to c
# id: problem1318
# author: Indigo
# date: 2021-01-23
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/minimum-flips-to-make-a-or-b-equal-to-c/description/>
# hidden: true
# ---
#
# Given 3 positives numbers `a`, `b` and `c`. Return the minimum flips required
# in some bits of `a` and `b` to make ( `a` OR `b` == `c` ). (bitwise OR
# operation).
# Flip operation consists of change  **any**  single bit 1 to 0 or change the
# bit 0 to 1 in their binary representation.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/06/sample_3_1676.png)
#
#
#
#     Input: a = 2, b = 6, c = 5
#     Output: 3
#     Explanation: After flips a = 1 , b = 4 , c = 5 such that (a OR b == c)
#
# **Example 2:**
#
#
#
#     Input: a = 4, b = 2, c = 7
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: a = 1, b = 2, c = 3
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= a <= 10^9`
#   * `1 <= b <= 10^9`
#   * `1 <= c <= 10^9`
#
#
## @lc code=start
using LeetCode

min_flip(a, b, c) = min_flip(a|>Int32, b|>Int32, c|>Int32)
function min_flip(a::Int32, b::Int32, c::Int32)
    sa, sb, sc = bitstring.([a, b, c])
    res = 0
    for i in 1:32
        if sc[i] == '1'
            (sa[i] == '0' == sb[i]) && (res += 1)
        else
            res += (sa[i] == '1') + (sb[i] == '1')
        end
    end
    res
end
## @lc code=end


text:
#===================================================================================================
Infrastructure of multiobjective evolution.
===================================================================================================#

"""
    dominate(f1::AV{T}, f2::AV{T}, [senses::AV{<:Integer}]) where {T<:Real}

Test whether `f1` dominates `f2` w.r.t the given `senses`. If the argument `senses` is not provided,
then maximization is assumed.
"""
function dominate(f1::AV{T}, f2::AV{T}, senses::AV{<:Integer}) where {T<:Real}
    @assert length(f1) == length(f2) && length(f2) == length(senses)
    one_better = false
    m = length(f1)
    @inbounds for i = 1:m
        fw1 = f1[i] * senses[i]
        fw2 = f2[i] * senses[i]
        if fw1 < fw2
            return false
        elseif fw1 > fw2
            one_better = true
        end
    end
    return one_better
end

function dominate(f1::AbstractVector{T}, f2::AbstractVector{T}) where {T<:Real}
    @assert length(f1) == length(f2)
    one_better = false
    m = length(f1)
    @inbounds for i = 1:m
        fw1 = f1[i]
        fw2 = f2[i]
        if fw1 < fw2
            return false
        elseif fw1 > fw2
            one_better = true
        end
    end
    return one_better
end


text:
import FLoops
using Documenter: doctest

# Workaround `UndefVarError: FLoops not defined`
@eval Main import FLoops

doctest(FLoops, manual = true)


text:
#showutils.jl

function addspace(s::String, n::Int; first = false)::String
    if n > 0
        for i = 1:n
            if first s = Char(' ') * s else s = s * Char(' ') end
        end
    end
    return s
end

function printmatrix(io::IO, m::Matrix; header = false)
    sm = string.(m)
    lv = maximum(length.(sm), dims = 1)
    for r = 1:size(sm, 1)
        for c = 1:size(sm, 2)
            sm[r,c] = addspace(sm[r,c], lv[c] - length(sm[r,c]))*"   "
        end
    end
    i = 1
    if header
        line = "⎯" ^ sum(length.(sm[1,:]))
        println(io, line)
        for c = 1:size(sm, 2)
            print(io, sm[1,c])
        end
        print(io, "\n")
        println(io, line)
        i = 2
    end
    for r = i:size(sm, 1)
        for c = 1:size(sm, 2)
            print(io, sm[r,c])
        end
        print(io, "\n")
    end
end


text:
using ExplainableAI: MaxActivationSelector, IndexSelector
using Random

ns_max = @inferred MaxActivationSelector()
ns_idx = @inferred IndexSelector(4)

# Array
A = [-2.1694243, 2.4023275, 0.99464744, -0.1514646, 1.0307171]
R = similar(A)

@test_throws ArgumentError @inferred ns_max(A)
@test_throws ArgumentError @inferred ns_idx(A)

# 5x1 input
A = reshape(A, 5, 1)
R = similar(A)

I_max = @inferred ns_max(A)
I_idx = @inferred ns_idx(A)
@test I_max isa Vector{CartesianIndex{2}}
@test I_idx isa Vector{CartesianIndex{2}}
@test I_max == [CartesianIndex(2, 1)]
@test I_idx == [CartesianIndex(4, 1)]

# 5x3 input
A = [
    0.0564584 0.736398 0.370134
    0.112634 0.0193744 0.679662
    0.942529 0.335366 0.0132769
    0.170132 0.474743 0.590655
    0.218707 0.0440574 0.962128
]
R = similar(A)

I_max = @inferred ns_max(A)
I_idx = @inferred ns_idx(A)
@test I_max isa Vector{CartesianIndex{2}}
@test I_idx isa Vector{CartesianIndex{2}}
@test I_max == [
    CartesianIndex(3, 1)
    CartesianIndex(1, 2)
    CartesianIndex(5, 3)
]
@test I_idx == [
    CartesianIndex(4, 1)
    CartesianIndex(4, 2)
    CartesianIndex(4, 3)
]

# 4x3x2 input with Tuple IndexSelector
ns_idx = @inferred IndexSelector((4, 2))
A = rand(MersenneTwister(1234), Float32, 4, 3, 2)
R = similar(A)

I_max = @inferred ns_max(A)
I_idx = @inferred ns_idx(A)
@test I_max isa Vector{CartesianIndex{3}}
@test I_idx isa Vector{CartesianIndex{3}}
@test I_max == [
    CartesianIndex(2, 1, 1)
    CartesianIndex(4, 1, 2)
]
@test I_idx == [
    CartesianIndex(4, 2, 1)
    CartesianIndex(4, 2, 2)
]


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

struct MPCInputs <: AbstractInputs
    s::MPCScenario
    techs::Techs
    existing_sizes::DenseAxisArray{Float64, 1}  # (techs.all)
    time_steps::UnitRange
    time_steps_with_grid::Array{Int, 1}
    time_steps_without_grid::Array{Int, 1}
    hours_per_time_step::Float64
    months::UnitRange
    production_factor::DenseAxisArray{Float64, 2}  # (techs.all, time_steps)
    levelization_factor::Dict{String, Float64}  # (techs.all)
    value_of_lost_load_per_kwh::Array{R, 1} where R<:Real #default set to 1 US dollar per kwh
    pwf_e::Float64
    pwf_om::Float64
    pwf_fuel::Dict{String, Float64}
    third_party_factor::Float64
    ratchets::UnitRange
    techs_by_exportbin::DenseAxisArray{Array{String,1}}  # indexed on [:NEM, :WHL]
    export_bins_by_tech::Dict{String, Array{Symbol, 1}}
    cop::Dict{String, Float64}  # (techs.cooling)
    thermal_cop::Dict{String, Float64}  # (techs.absorption_chiller)
end

function MPCInputs(fp::String)
    s = MPCScenario(JSON.parsefile(fp))
    MPCInputs(s)
end

function MPCInputs(s::MPCScenario)

    time_steps = 1:length(s.electric_load.loads_kw)
    hours_per_time_step = 1 / s.settings.time_steps_per_hour
    techs, production_factor, existing_sizes = setup_tech_inputs(s)
    months = 1:length(s.electric_tariff.monthly_demand_rates)

    techs_by_exportbin = DenseAxisArray([ techs.all, techs.all, techs.all], s.electric_tariff.export_bins)
    # TODO account for which techs have access to export bins (when we add more techs than PV)

    levelization_factor = Dict(t => 1.0 for t in techs.all)
    pwf_e = 1.0
    pwf_om = 1.0
    pwf_fuel = Dict{String, Float64}()
    pwf_fuel["Generator"] = 1.0
    third_party_factor = 1.0

    time_steps_with_grid, time_steps_without_grid, = setup_electric_utility_inputs(s)

    export_bins_by_tech = Dict{String, Array{Symbol, 1}}()
    for t in techs.elec
        export_bins_by_tech[t] = s.electric_tariff.export_bins
    end
    # TODO implement export bins by tech (rather than assuming that all techs share the export_bins)

    #Placeholder COP because the REopt model expects it
    cop = Dict("ExistingChiller" => s.cooling_load.cop)
    thermal_cop = Dict{String, Float64}()

    MPCInputs(
        s,
        techs,
        existing_sizes,
        time_steps,
        time_steps_with_grid,
        time_steps_without_grid,
        hours_per_time_step,
        months,
        production_factor,
        levelization_factor,  # TODO need this?
        typeof(s.financial.value_of_lost_load_per_kwh) <: Array{<:Real, 1} ? s.financial.value_of_lost_load_per_kwh : fill(s.financial.value_of_lost_load_per_kwh, length(time_steps)),
        pwf_e,
        pwf_om,
        pwf_fuel,
        third_party_factor,
        # maxsize_pv_locations,
        1:length(s.electric_tariff.tou_demand_ratchet_time_steps),  # ratchets
        techs_by_exportbin,
        export_bins_by_tech,
        cop,
        thermal_cop
        # s.site.min_resil_time_steps,
        # s.site.mg_tech_sizes_equal_grid_sizes,
        # s.site.node
    )
end

function setup_tech_inputs(s::MPCScenario)

    techs = Techs(s)

    time_steps = 1:length(s.electric_load.loads_kw)

    # REoptInputs indexed on techs:
    existing_sizes = DenseAxisArray{Float64}(undef, techs.all)
    production_factor = DenseAxisArray{Float64}(undef, techs.all, time_steps)

    if !isempty(techs.pv)
        setup_pv_inputs(s, existing_sizes, production_factor)
    end

    if "Generator" in techs.all
        setup_gen_inputs(s, existing_sizes, production_factor)
    end

    return techs, production_factor, existing_sizes
end

function setup_pv_inputs(s::MPCScenario, existing_sizes, production_factor)
    for pv in s.pvs
        production_factor[pv.name, :] = pv.prod_factor_series
        existing_sizes[pv.name] = pv.size_kw
    end
    return nothing
end

function setup_gen_inputs(s::MPCScenario, existing_sizes, production_factor)
    existing_sizes["Generator"] = s.generator.size_kw
    production_factor["Generator", :] = ones(length(s.electric_load.loads_kw))
    return nothing
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
# To run this script, `cd` to the `./test/fixtures` directory and then, from the Julia terminal, `include("./runner.jl")`.

import JSON

# Positive normal values:
x = linspace( 0, 1001, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./positive_normal.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Negative normal values:
x = linspace( -1001, 0, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./negative_normal.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Positive small values:
x = linspace( 0, 1, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./positive_small.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Negative small values:
x = linspace( -1, 0, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./negative_small.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Positive tiny values:
x = linspace( 5e-39, 1e-30, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./positive_tiny.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Negative tiny values:
x = linspace( -5e-39, -1e-30, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./negative_tiny.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Positive subnormal values:
x = linspace( 1e-40, 1e-45, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./positive_subnormal.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Negative subnormal values:
x = linspace( -1e-40, -1e-45, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./negative_subnormal.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Large positive values:
x = linspace( 1e36, 1e38, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./positive_large.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )

# Large negative values:
x = linspace( -1e36, -1e38, 500 )

y = Array( Any, length(x) )
for i in eachindex(x)
	# Mimic implicit type promotion in JavaScript:
	y[i] = bits(convert(Float32, x[i]))
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("./negative_large.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )
close( outfile )


text:
function mm_spec_kyparisis3_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(2, 0, 0, 0, 3, 2, 0)
    return dimspec
end

function mm_spec_kyparisis3_defn(x, pr, ps)

    f = (1/2)*(x[1]-pr[1])^2 + (1/2)*(x[2]-pr[2])^2 + x[1] + x[2]

    ce = Vector{Num}()

    ci = Vector{Num}(undef, 3)
    ci[1] = x[1]
    ci[2] = x[2]
    ci[3] = x[1] + x[2] - pr[1] - pr[2]

    F = Matrix{Num}(undef, 0, 0)
    # F[1,1] = x[1]
    # F[2,1] = x[2]

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_kyparisis3_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_kyparisis3_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_kyparisis3_knownsols()

    function local_kyparisis3_knownsol(pr::Vector{R}, ps::Vector{Int64}) where {R <: Real}
        @assert ( 2 == length(pr) && 0 == length(ps) )
        e1 = pr[1]
        e2 = pr[2]

        if ( e1 >= 0 && e2 >= 0 )
            x_sol = (e1, e2)
        elseif ( e1 >= 0 && e2 < 0 )
            x_sol = (max(e1-1, e1 + e2, 0), 0)
        elseif ( e1 < 0 && e2 >= 0 )
            x_sol = (0, max(e2-1, e1 + e2, 0))
        else
            x_sol = (0, 0)
        end
        f_sol = (1/2)*(x_sol[1]-e1)^2 + (1/2)*(x_sol[2]-e2)^2 + x_sol[1] + x_sol[2]

        return (x_sol, f_sol)
    end

    return [ local_kyparisis3_knownsol ]
end

function mm_spec_kyparisis3_parameterisations(t)

    defn1 = MPCCParameterisationDefn(
            Vector{Num}([ 3.0, t ]),
            (-3.0, 3.0),
            "Vertical, rising"
        )

    return Vector{MPCCParameterisationDefn}( [ defn1 ] )
end


text:
# convert_a3d

# Function that allows federation by setting the `output_format`  in the Stanmodel().

function convert_a3d(a3d_array, cnames, ::Val{:dataframe})
  snames = [Symbol(cnames[i]) for i in 1: length(cnames)]
  [DataFrame(a3d_array[:,:,i], snames) for i in 1:size(a3d_array, 3)]
end


text:
"""
    KRLS(y, x, K, M, λ, s2n, forgetting)

Kernel Recursive Least Squares (Tracker) algorithm with a AFxed budget size via growing and pruning.

# Arguments
- `y::REALVEC`
- `x::REALMAT`
- `K::Kernel`
- `M::Integer`
- `λ::Real = 1`
- `s2n::Real = 1`
- `forgetting::String = "B2P"`

# Examples
```julia
2 + 3
```
"""
function KRLS(y::REALVEC,
              x::REALMAT,
              K::Kernel,
              M::Integer = size(y),
              λ::Real = 1,
              s2n::Real = 1,
              forgetting::String = "B2P")

    @assert s2n > 0 "Observation noise (regularization), s2n, must be strictly positive."
    @assert λ <= 1 && λ > 0 "Forgetting factor, λ, must be ]0;1]."
    @assert M <= length(y) && M > 0 "Budget size, M, must be [1;T]."
    @assert (forgetting == "B2P" || forgetting == "UI") "Input, forget, must be B2P or UI."

    T  = length(y)
    ŷ  = zeros(T)
    σ² = zeros(T)
    basis = [1]
    xb = x[basis, :]

    # Initialization
    ktt = kernel(x[1:1, :], K)[1] + eps()
    μ = [y[1] * ktt]  / (s2n + ktt) # predictive mean at time 1
    Σ = ktt - ktt^2 / (s2n + ktt)
    Q = reshape([1/ktt], 1, 1)
    s0n = y[1]^2 / (s2n + ktt)
    s0d = 1.0
    s02 = s0n / s0d
    σ²[1] = s02 * s2n              # predictive variance at time 1
    m = 1

    for t in 2:T
        # Forget
        forget!(μ, Σ, λ, kernel(xb, K), forgetting)

        ktt = kernel(x[t:t, :], K)[1] + eps()
        kbt = kernel(xb, x[t:t, :], K)

        q = vec(Q * kbt)
        h = vec(Σ * q)

        # Projection uncertainty
        γ² = ktt - dot(kbt, q)
        println("gamma: $γ²")
        #γ² < 0.0 && (γ² = eps()) # Ensure that gamma squared is not negative

        # Noiseless prediction variance
        s2f = γ² + dot(q, h)
        #s2f < 0.0 && (s2f = eps()) # Ensure that s2f is not negative

        # Predictive mean
        ȳ = dot(q, μ)

        # Predictive variance
        s2y = s2n + s2f

        ŷ[t]  = ȳ
        σ²[t] = s2y * s02

        # Observe yₜ

        # Estimate of s02 via Maximum Likelihood
        s0n = s0n + λ * (y[t] - ȳ) / s2y
        s0d = s0d + λ
        s02 = s0n / s0d

        # Get posteriors of N(f|μ, Σ)
        push!(μ, ȳ)
        p = vcat(h, s2f)
        μ = μ + (y[t] - ȳ)/s2y * p
        s = reshape(p, :, 1)
        Σ = hcat(vcat(Σ, h'), s)
        Σ = Σ - (p * p') / s2y

        # Ensure that gamma is not too small
        if γ² < (10*eps())
            μ  = μ[1:end-1]
            Σ  = Σ[1:end-1, 1:end-1]
        else
            # Dictionary update
            push!(basis, t)
            xb = x[basis, :]
            m = m + 1

            p = vcat(q, [-1.0])
            Q = hcat(vcat(Q, zeros(1, m - 1)), zeros(m, 1)) + p * p' / γ²

            if m > M
                # MSE pruning criterion
                errors = ((Q * μ) ./ diag(Q)).^2
                dd, r = findmin(errors)

                # Remove element r
                d  = setdiff(1:m, r)
                Qs = Q[d, :]
                qs = Q[r, r]
                Q  = Q[d, d]
                Q  = Q - (Qs * Qs')/qs
                μ  = μ[d]
                Σ  = Σ[d, d]
                m  = m - 1
                xb = xb[d, :]
                basis = basis[d]
            end
        end
    end
    return (predictions = ŷ, variances = σ², basis = basis, xb = xb, μ = μ, Σ = Σ, Q = Q)
end

function KRLS(y::REALVEC,
              x::REALVEC,
              K::Kernel,
              M::Integer = size(y),
              λ::Real = 1,
              s2n::Real = 1,
              forgetting::String = "B2P")
    KRLS(y, reshape(x, :, 1), K, M, λ, s2n, forgetting)
end

function forget!(μ, Σ, λ, kb, forgetting)
    if forgetting == "B2P"
        Σ = Σ * λ .+ kb * (1 - λ)
        μ = sqrt(λ) * μ
    elseif forgetting == "UI"
        Σ = Σ / λ
    end
end

predict(μ::Vector, Q::Matrix, x::REALMAT, xb::REALMAT, K::Kernel) = kernel(xb, x, K)' * Q * μ


text:
export Basis1d, Element1d

type Basis1d
    Q::Int
    z::Array{Float64,1}
    w::Array{Float64,1}
    D::Array{Float64,2}

    function Basis1d(q::Int)
        Q = q
        z = zeros_glj(Q, T=Float64)
        w = weights_glj(z)
        D = diffmat_glj(z)
        new(Q, z, w, D)
    end
end

type Element1d

    basis::Basis1d
    Q::Int
    dx::Float64
    x::Array{Float64,1}
    dzdx::Float64
    jac::Float64
    mass::Array{Float64,1}

    function Element1d(b::Basis1d, xnodes::Array{Float64,1})
        basis = b
        Q = b.Q
        dx = xnodes[2] - xnodes[1]
        x = (1-b.z)/2 * xnodes[1] + (1+b.z)/2 * xnodes[2]
        dzdx = 2 / dx
        jac = dx / 2
        mass = jac * basis.w
        new(basis, Q, dx, x, dzdx, jac, mass)
    end

end

using Base.LinAlg.BLAS

function deriv_z(elem::Element1d, f::Array{Float64,1})
    elem.basis.D * f
end

function deriv_x(elem::Element1d, f::Array{Float64,1})
    elem.dzdx * deriv_z(elem, f)
end

function convective(elem::Element1d, f::Array{Float64,1}, u::Array{Float64,1})

    udf = deriv_x(elem, f)
    for i::Int = 1:elem.Q
        udf[i] *= u[i]
    end
    return udf
end

function integrate(elem::Element1d, f::Array{Float64,1})

    dot(elem.mass, f)
end

function laplace_matrix(elem::Element1d)

    local Q::Int = elem.Q, i::Int, j::Int

end


text:
using Pkg
Pkg.activate("/home/caseykneale/.julia/dev/ObjectDetectionStats")
using ObjectDetectionStats

ods_machine = ObjectDetectionScore( 4 )

ods_machine.TP[:] .= [90,90,90,90]
ods_machine.FP[:] .= [10,10,10,10]
ods_machine.FN[:] .= [10,10,10,10]

classwise_precision( ods_machine )
classwise_recall( ods_machine )
macro_precision(ods_machine)
macro_recall(ods_machine)


text:
using Gadfly, DataFrames

plot(readtable(Pkg.dir("Gadfly", "test", "issue120.csv")), x=:x1, y=:x2)


text:
# ---
# title: 1354. Construct Target Array With Multiple Sums
# id: problem1354
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/construct-target-array-with-multiple-sums/description/>
# hidden: true
# ---
#
# Given an array of integers `target`. From a starting array, `A` consisting of
# all 1's, you may perform the following procedure :
#
#   * let `x` be the sum of all elements currently in your array.
#   * choose index `i`, such that `0 <= i < target.size` and set the value of `A` at index `i` to `x`.
#   * You may repeat this procedure as many times as needed.
#
# Return True if it is possible to construct the `target` array from `A`
# otherwise return False.
#
#
#
# **Example 1:**
#
#
#
#     Input: target = [9,3,5]
#     Output: true
#     Explanation: Start with [1, 1, 1]
#     [1, 1, 1], sum = 3 choose index 1
#     [1, 3, 1], sum = 5 choose index 2
#     [1, 3, 5], sum = 9 choose index 0
#     [9, 3, 5] Done
#
#
# **Example 2:**
#
#
#
#     Input: target = [1,1,1,2]
#     Output: false
#     Explanation: Impossible to create target array from [1,1,1,1].
#
#
# **Example 3:**
#
#
#
#     Input: target = [8,5]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `N == target.length`
#   * `1 <= target.length <= 5 * 10^4`
#   * `1 <= target[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

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


text:
##################################################
# Type promotion/conversion
##################################################

Base.convert(::Type{Sym}, x) = Sym(x)
Base.convert(::Type{Sym}, x::Sym) = x
Base.convert(::Type{Sym{T}}, x) where {T} = Sym{T}(x)
Base.convert(::Type{Sym{T}}, x::Sym) where {T} = Sym{T}(Val(:noargs), gethead(x), getargs(x)...)

Base.convert(::Type{Tuple{Vararg{Sym,N} where N}}, xs::Tuple) = Tuple(Sym(x) for x in xs)
Base.convert(::Type{Tuple{Vararg{Sym{T},N} where N}}, xs::Tuple) where {T} = Tuple(Sym{T}(x) for x in xs)

Base.promote_rule(::Type{Sym{Any}}, ::Type{Any}) = Sym
Base.promote_rule(::Type{Sym{T}}, ::Type{Any}) where {T} = Sym
Base.promote_rule(::Type{Sym{T}}, ::Type{S}) where {T,S} = Sym
Base.promote_rule(::Type{Sym{T}}, ::Type{T}) where {T} = Sym{T}
Base.promote_rule(::Type{Sym{T}}, ::Type{Sym{S}}) where {T,S} = Sym{promote_type(T, S)}

Base.promote_op(f, a::Type{<:Sym{T}}) where {T} = Sym{promote_tag(:call, f, T)}

Base.promote_op(f, a::Type{Sym{T}}, b::S) where {T,S} = Sym{promote_tag(:call, f, T, S)}
Base.promote_op(f, a::T, b::Type{S}) where {T,S} = Sym{promote_tag(:call, f, T, S)}
Base.promote_op(f, a::Type{Sym{T}}, b::Type{Sym{S}}) where {T,S} = Sym{promote_tag(:call, f, T, S)}

##################################################
# Tag promotion
##################################################

promote_tag(head::Symbol, args...) = promote_tag(Val(head), args...)
promote_tag(::Val{:protected}, argtag) = argtag
promote_tag(::Val{:fn}, argtag) = argtag
promote_tag(::Val{:call}, op::Sym, argtags...) = promote_tag(Val(:call), firstarg(op), argtags...)
promote_tag(::Val{:call}, op, argtags...) = Base.promote_op(op, argtags...)
promote_tag(::Val{:(.)}, x::Type, argtags...) = Base.promote_op(getproperty, x, argtags...)
promote_tag(::Val{:ref}, x::Type, argtags...) = Base.promote_op(getindex, x, argtags...)
promote_tag(::Val{:ref}, x::Type{<:Matrix}, argtags...) = eltype(x)
promote_tag(::Val{Symbol("'")}, argtags...) = Base.promote_op(adjoint, argtags...)

# Improve stability, assuming all numbers form a ring
promote_tag(::Val{:call}, ::typeof(Base.:+), argtags::Type{T}...) where {T<:Number} =
    Base.promote_type(T, typeof(zero(T) + zero(T)))
promote_tag(::Val{:call}, ::typeof(Base.:-), argtags::Type{T}...) where {T<:Number} =
    Base.promote_type(T, typeof(zero(T) - zero(T)))
promote_tag(::Val{:call}, ::typeof(Base.:*), argtags::Type{T}...) where {T<:Number} =
    T
promote_tag(::Val{:call}, ::typeof(Base.:/), argtags::Type{T}...) where {T<:Number} =
    Base.promote_type(T, typeof(zero(T)/oneunit(T)))


text:
using Gtk.ShortNames, DataFrames, DefaultApplication, Mustache, Luxor, Printf
using JLD

# Path to index logo
if Sys.iswindows()
    global indexImgPath =
        joinpath(dirname(Base.source_path()), "figures\\index_logo.png")
    global indexImgPathT =
        joinpath(dirname(Base.source_path()), "C:\\Windows\\Temp\\index_logo_temp.png")
end

if Sys.islinux()
    global indexImgPath =
        joinpath(dirname(Base.source_path()), "figures/index_logo.png")
    global indexImgPathT =
        joinpath(dirname(Base.source_path()), "/tmp/index_logo_temp.png")
end

function PIAssistantGUI()
    global wSc, hSc = screen_size()

    img = readpng(indexImgPath)
    wImg = img.width
    hImg = img.height

    # Ratio for proportion
    imgP = hImg/wImg

    x = round((hSc/2)/imgP)/wImg
    y = round(hSc/2)/hImg

    # Creates a temporary Image as index.
    Drawing(round((hSc/2)/imgP)+60, round(hSc/2)+60, indexImgPathT)

    origin()
    background("white")

    scale(x, y)
    placeimage(img, -wImg/2, -hImg/2, 1)
    finish()
    preview()

    indexImg = Image()
    set_gtk_property!(indexImg, :file, indexImgPathT)

    winIndex = Window(indexImg)
    set_gtk_property!(winIndex, :decorated, false)
    set_gtk_property!(winIndex, :window_position, 3)
    set_gtk_property!(winIndex, :accept_focus, true)
    set_gtk_property!(winIndex, :width_request, round((hSc/2)/imgP)+60)
    set_gtk_property!(winIndex, :height_request, round(hSc/2)+60)

    Gtk.showall(winIndex)

    sleep(3)
    destroy(winIndex)

    mainPI()
end


text:
julia_include_dir = joinpath(match(r"(.*)(bin)", Sys.BINDIR).captures[1],"include","julia")
if !isdir(julia_include_dir)  # then we're running directly from build
	julia_base_dir_aux = splitdir(splitdir(Sys.BINDIR)[1])[1]  # useful for running-from-build
	julia_include_dir = joinpath(julia_base_dir_aux, "usr", "include" )
	julia_include_dir *= ";" * joinpath(julia_base_dir_aux, "src", "support" )
	julia_include_dir *= ";" * joinpath(julia_base_dir_aux, "src" )
end
print(julia_include_dir)

max_stars_repo_name: dmikushin/libcxxwrap-julia

max_issues_repo_name: dmikushin/libcxxwrap-julia

max_forks_repo_name: dmikushin/libcxxwrap-julia


text:
using BinaryBuilder, BinaryBuilderBase, Downloads, Pkg

verbose = "--verbose" in ARGS

# Read in input `.json` file
json = String(read(ARGS[1]))
buff = IOBuffer(strip(json))
objs = []
while !eof(buff)
    push!(objs, BinaryBuilder.JSON.parse(buff))
end

# Merging modifies `obj`, so let's keep an unmerged version around
objs_unmerged = deepcopy(objs)

# Merge the multiple outputs into one
merged = BinaryBuilder.merge_json_objects(objs)
BinaryBuilder.cleanup_merged_object!(merged)
BinaryBuilder.cleanup_merged_object!.(objs_unmerged)

# Determine build version
name = merged["name"]
version = merged["version"]
# Filter out build-time dependencies that will not go into the dependencies of
# the JLL packages.
dependencies = Dependency[dep for dep in merged["dependencies"] if !isa(dep, BuildDependency)]
lazy_artifacts = merged["lazy_artifacts"]
build_version = BinaryBuilder.get_next_wrapper_version(name, version)
repo = "JuliaBinaryWrappers/$(name)_jll.jl"
code_dir = joinpath(Pkg.devdir(), "$(name)_jll")
julia_compat = merged["julia_compat"]

# Register JLL package using given metadata
BinaryBuilder.init_jll_package(
    name,
    code_dir,
    repo,
)

function download_cached_binaries(download_dir, platforms)
    # Grab things out of the aether for maximum consistency
    bb_hash = ENV["BB_HASH"]
    proj_hash = ENV["PROJ_HASH"]

    for platform in platforms
        url = "https://julia-bb-buildcache.s3.amazonaws.com/$(bb_hash)/$(proj_hash)/$(triplet(platform)).tar.gz"
        filename = "$(name).v$(version).$(triplet(platform)).tar.gz"
        Downloads.download(url, joinpath(download_dir, filename))
    end
end

function download_binaries_from_release(download_dir)
    function do_download(download_dir, info)
        url = info["url"]
        hash = info["sha256"]
        filename = basename(url)
        BinaryBuilderBase.download_verify(url, hash, joinpath(download_dir, filename))
    end

    # Doownload the tarballs reading the information in the current `Artifacts.toml`.
    artifacts = Pkg.Artifacts.load_artifacts_toml(joinpath(code_dir, "Artifacts.toml"))[name]
    if artifacts isa Dict
        # If it's a Dict, that means this is an AnyPlatform artifact, act accordingly.
        info = artifacts["download"][1]
        do_download(download_dir, info)
    else
        # Otherwise, it's a Vector, and we must iterate over all platforms.
        for artifact in artifacts
            info = artifact["download"][1]
            do_download(download_dir, info)
        end
    end
end

# Filter out build-time dependencies also here
for json_obj in [merged, objs_unmerged...]
    json_obj["dependencies"] = Dependency[dep for dep in json_obj["dependencies"] if !isa(dep, BuildDependency)]
end
skip_build = get(ENV, "SKIP_BUILD", "false") == "true"
mktempdir() do download_dir
    # Grab the binaries for our package
    if skip_build
        # We only want to update the wrappers, so download the tarballs from the
        # latest build.
        download_binaries_from_release(download_dir)
    else
        # We are going to publish the new binaries we've just baked, take them
        # out of the cache while they're hot.
        download_cached_binaries(download_dir, merged["platforms"])
    end

    # Push up the JLL package (pointing to as-of-yet missing tarballs)
    tag = "$(name)-v$(build_version)"
    upload_prefix = "https://github.com/$(repo)/releases/download/$(tag)"

    # If we didn't rebuild the tarballs, save the original Artifacts.toml
    artifacts_toml = skip_build ? read(joinpath(code_dir, "Artifacts.toml"), String) : ""
    # This loop over the unmerged objects necessary in the event that we have multiple packages being built by a single build_tarballs.jl
    for (i,json_obj) in enumerate(objs_unmerged)
        from_scratch = (i == 1)
        BinaryBuilder.rebuild_jll_package(json_obj; download_dir=download_dir, upload_prefix=upload_prefix, verbose=verbose, lazy_artifacts=json_obj["lazy_artifacts"], from_scratch=from_scratch)
    end

    if skip_build
        # Restore Artifacts.toml
        write(joinpath(code_dir, "Artifacts.toml"), artifacts_toml)
    else
        # Upload the tarballs to GitHub releases
        BinaryBuilder.upload_to_github_releases(repo, tag, download_dir; verbose=verbose)
    end
end
BinaryBuilder.push_jll_package(name, build_version)
BinaryBuilder.register_jll(name, build_version, dependencies, julia_compat)


text:
include("BaryTrees/BaryTrees.jl")
using Main.BaryTrees, Plots

# construct an empty tree from a Barycentric coordinate system defined
# by the triangle with vertices (0,sqrt(3)/2), (-0.5,0), (0.5,0)
T = BaryTree(
    BaryCentric(Point((0,sqrt(3)/2)),Point((-0.5,0)),Point((0.5,0.0)))
)

# random points in a square: (-1,-1), (-1,1), (1,1), (1,-1)
p = (rand(2^9,2) .- 0.5).*2;

pin = zeros(1,2)
pout = zeros(1,2)

for i in 1:size(p,1)
    b = Cartesian2BaryCentric(Point((p[i,:]...)),T.boundary.b)

    # test if in trig
    if b[1] >= 0 && b[2] >= 0 && b[1]+b[2] < 1
        pin = cat(pin,p[i:i,:],dims=1)
    else
        pout = cat(pout,p[i:i,:],dims=1)
    end
end

# insert only points inside trig
for i in 1:size(pin,1)
    insert!(Point((pin[i,:]...)),T)
end

scatter(pin[:,1],pin[:,2],label="",dpi=300,
    aspect_ratio=:equal,xaxis=nothing,yaxis=nothing,axis=:off
)

Draw!(T)
plot!()

savefig("barytree.png")

max_stars_repo_name: Jerboa-app/BaryTrees

max_issues_repo_name: Jerboa-app/BaryTrees

max_forks_repo_name: Jerboa-app/BaryTrees