text:
# This does simple Poisson estimation using the meps1996.data

# The MEPS data

# Define dep and expl vbls
#
# 	The dep. vbls, in corresponding column
# 	Office based doctor visits	1
# 	Outpatient doctor visits	2
# 	Emergency room visits		3
# 	Inpatient visits		4
# 	Dental visits			5
# 	Prescriptions			6
#
using Econometrics, DelimitedFiles
which_dep = 1
if (which_dep == 1) dep = "OBDV"   end
if (which_dep == 2) dep = "OPV"    end
if (which_dep == 3) dep = "IPV"    end
if (which_dep == 4) dep = "ERV"    end
if (which_dep == 5) dep = "DV"     end
if (which_dep == 6) dep = "PRESCR" end
data = readdlm("../Data/meps1996.data")
y = data[:,which_dep]
x = data[:,7:12]
n = size(x,1)
x = [ones(n,1) x]
# to check effects of poor scaling, try commenting the next line
x[:,end] = x[:,end]/1000.0 # put income in thousands
names = ["constant", "pub. ins.","priv. ins.", "sex", "age","edu","inc"]
title = "Poisson model, "*dep* ",  MEPS 1996 full data set"
model = theta -> poisson(theta, y, x)
theta = zeros(size(x,2)) # start values for estimation
# try adding the option vc=2 for "Hessian" or vc=3 for OPG
thetahat, objvalue, V, converged = mleresults(model, theta, title, names, vc=1)
nothing

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using Distributed

function typeof_future(r::Future)
    if isready(r)
        pid = r.where
        f   = @spawnat pid typeof(fetch(r))
        T   = fetch(f)
        return T
    else
        error("future is not ready")
        return nothing
    end
end

"""
    struct RemoteHMatrix{S,T}

A light wrapper for a `Future` storing an `HMatrix`.
"""
struct RemoteHMatrix{S,T}
    future::Future
end

Base.fetch(r::RemoteHMatrix) = fetch(r.future)

function Base.getindex(r::RemoteHMatrix,i::Int,j::Int)
    pid = r.future.where
    @fetchfrom pid getindex(fetch(r),i,j)
end

"""
    mutable struct DHMatrix{R,T} <: AbstractHMatrix{T}

Concrete type representing a hierarchical matrix with data distributed amongst
various workers. Its structure is very similar to `HMatrix`, except that the
leaves store a [`RemoteHMatrix`](@ref) object.

The `data` on the leaves of a `DHMatrix` may live on a different worker, so
calling `fetch` on them should be avoided whenever possible.
"""
mutable struct DHMatrix{R,T} <: AbstractHMatrix{T}
    rowtree::R
    coltree::R
    # admissible:: Bool --> false
    data::Union{RemoteHMatrix{R,T},Nothing}
    children::Matrix{DHMatrix{R,T}}
    parent::DHMatrix{R,T}
    # inner constructor which handles `nothing` fields.
    function DHMatrix{R,T}(rowtree,coltree,data,children,parent) where {R,T}
        dhmat = new{R,T}(rowtree,coltree,data)
        dhmat.children = isnothing(children) ? Matrix{DHMatrix{R,T}}(undef,0,0) : children
        dhmat.parent   = isnothing(parent) ? dhmat : parent
        return dhmat
    end
end

"""
    DHMatrix{T}(rowtree,coltree;partition_strategy=:distribute_columns)

Construct the block structure of a distributed hierarchical matrix covering
`rowtree` and `coltree`. Returns a `DHMatrix` with leaves that are empty.

The `partition_strategy` keyword argument determines how to partition the blocks
for distributed computing. Currently, the only available options is
`distribute_columns`, which will partition the columns of the underlying matrix
into `floor(log2(nw))` parts, where `nw` is the number of workers available.
"""
function DHMatrix{T}(rowtree::R, coltree::R; partition_strategy=:distribute_columns) where {R,T}
    #build root
    root  = DHMatrix{R,T}(rowtree,coltree,nothing,nothing,nothing)
    # depending on the partition strategy, dispatch to appropriate (recursive)
    # method
    if partition_strategy == :distribute_columns
        nw        = nworkers()
        dmax      = floor(Int64,log2(nw))
        _build_block_structure_distribute_cols!(root,dmax)
    else
        error("unrecognized partition strategy")
    end
    return root
end

function _build_block_structure_distribute_cols!(current_node::DHMatrix{R,T},dmax) where {R,T}
    if Trees.depth(current_node) == dmax
        return current_node
    else
        X = rowtree(current_node)
        Y = coltree(current_node)
        (isleaf(X) || isleaf(Y)) && error("you should split the tree further")
        # do not recurse on row, only on columns
        row_children = [X]
        col_children = Y.children
        children     = [DHMatrix{R,T}(r,c,nothing,nothing,current_node) for r in row_children, c in col_children]
        current_node.children = children
        for child in children
            _build_block_structure_distribute_cols!(child,dmax)
        end
        return current_node
    end
end

Base.size(H::DHMatrix) = length(rowrange(H)), length(colrange(H))

function Base.show(io::IO,::MIME"text/plain",hmat::DHMatrix)
    # isclean(hmat) || return print(io,"Dirty DHMatrix")
    println(io,"Distributed HMatrix of $(eltype(hmat)) with range $(rowrange(hmat)) × $(colrange(hmat))")
    nodes = collect(PreOrderDFS(hmat))
    println(io,"\t number of nodes in tree: $(length(nodes))")
    leaves = collect(Leaves(hmat))
    @printf(io,"\t number of leaves: %i\n",length(leaves))
    for (i,leaf) in enumerate(leaves)
        r   = leaf.data.future
        pid = r.where
        irange,jrange = @fetchfrom pid rowrange(fetch(r)), colrange(fetch(r))
        println("\t\t leaf $i on process $pid spanning $irange × $jrange")
    end
end

"""
    _assemble_hmat_distributed(K,rtree,ctree;adm=StrongAdmissibilityStd(),comp=PartialACA();global_index=true,threads=false)

Internal methods called **after** the `DHMatrix` structure has been initialized
in order to construct the `HMatrix` on each of the leaves of the `DHMatrix`.
"""
function _assemble_hmat_distributed(K,rtree,ctree;adm=StrongAdmissibilityStd(),comp=PartialACA(),global_index=true,threads=false)
    #
    R = typeof(rtree)
    T = eltype(K)
    wids   = workers()
    root   = DHMatrix{T}(rtree,ctree;partition_strategy=:distribute_columns)
    leaves = Leaves(root) |> collect
    @info "Assembling distributed HMatrix on $(length(leaves)) processes"
    @sync for (k, leaf) in enumerate(leaves)
        pid = wids[k] # id of k-th worker
        r   = @spawnat pid assemble_hmat(K, rowtree(leaf), coltree(leaf); adm, comp,global_index,threads,distributed=false)
        leaf.data = RemoteHMatrix{R,T}(r)
    end
    return root
end

function mul!(y::AbstractVector,A::DHMatrix,x::AbstractVector,a::Number,b::Number;
                            global_index=true,threads=true)
    # since the HMatrix represents A = Pr*H*Pc, where Pr and Pc are row and column
    # permutations, we need first to rewrite C <-- b*C + a*(Pc*H*Pb)*B as
    # C <-- Pr*(b*inv(Pr)*C + a*H*(Pc*B)). Following this rewrite, the
    # multiplication is performed by first defining B <-- Pc*B, and C <--
    # inv(Pr)*C, doing the multiplication with the permuted entries, and then
    # permuting the result  C <-- Pr*C at the end.
    ctree     = A.coltree
    rtree     = A.rowtree
    # permute input
    if global_index
        x         = x[ctree.loc2glob]
        y         = permute!(y,rtree.loc2glob)
        rmul!(x,a) # multiply in place since this is a new copy, so does not mutate exterior x
    else
        x = a*x # new copy of x
    end
    rmul!(y,b)
    leaves = filter_tree(x->Trees.isleaf(x),A)
    nb  = length(leaves)
    acc = Vector{Future}(undef,nb)
    @sync for i in 1:nb
        r        = leaves[i].data.future
        pid      = r.where
        jrange   = @fetchfrom pid colrange(fetch(r))
        T,n      = eltype(y), length(y)
        acc[i] = @spawnat pid mul!(zeros(T,n),fetch(r),view(x,jrange),1,0;global_index=false,threads)
    end
    # reduction stage: fetch everybody to worker 1 and sum them up
    for yi in acc
        axpy!(1,fetch(yi),y)
    end
    # permute output
    global_index && invpermute!(y,loc2glob(rtree))
    return y
end

function isclean(H::DHMatrix)
    for node in PreOrderDFS(H)
        if isleaf(node)
            if !hasdata(node)
                @warn "leaf node without data found"
                return false
            end
            if !isclean(data(node))
                @warn "dirty `HMatrix` leaf found"
                return false
            end
        else
            if hasdata(node)
                @warn "data found on non-leaf node"
                return false
            end
        end
    end
    return true
end

ext: jl

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
module ServiceInvocationClients
using gRPCClient
using ProtoBuf:google

include("dapr.jl")
using .dapr.proto.internals.v1

import Base: show

# begin service: dapr.proto.internals.v1.ServiceInvocation

export ServiceInvocationBlockingClient, ServiceInvocationClient

struct ServiceInvocationBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::ServiceInvocationBlockingStub

    function ServiceInvocationBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = ServiceInvocationBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct ServiceInvocationClient
    controller::gRPCController
    channel::gRPCChannel
    stub::ServiceInvocationStub

    function ServiceInvocationClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = ServiceInvocationStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::ServiceInvocationBlockingClient) = print(io, "ServiceInvocationBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::ServiceInvocationClient) = print(io, "ServiceInvocationClient(", client.channel.baseurl, ")")

import .dapr.proto.internals.v1: CallActor
"""
    CallActor

- input: dapr.proto.internals.v1.InternalInvokeRequest
- output: dapr.proto.internals.v1.InternalInvokeResponse
"""
CallActor(client::ServiceInvocationBlockingClient, inp::dapr.proto.internals.v1.InternalInvokeRequest) = CallActor(client.stub, client.controller, inp)
CallActor(client::ServiceInvocationClient, inp::dapr.proto.internals.v1.InternalInvokeRequest, done::Function) = CallActor(client.stub, client.controller, inp, done)

import .dapr.proto.internals.v1: CallLocal
"""
    CallLocal

- input: dapr.proto.internals.v1.InternalInvokeRequest
- output: dapr.proto.internals.v1.InternalInvokeResponse
"""
CallLocal(client::ServiceInvocationBlockingClient, inp::dapr.proto.internals.v1.InternalInvokeRequest) = CallLocal(client.stub, client.controller, inp)
CallLocal(client::ServiceInvocationClient, inp::dapr.proto.internals.v1.InternalInvokeRequest, done::Function) = CallLocal(client.stub, client.controller, inp, done)

# end service: dapr.proto.internals.v1.ServiceInvocation

end # module DaprClients


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

####  Create data for testing the case k=1, m>1
# `t` is the number of samples,
# `m` is the number of datasets,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise,
#  the more data are correlated.
function getData(t, m, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    X=randn(n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(n, t)) for i=1:m]
    return 𝐗
end

function getData(::Type{Complex{T}}, t, m, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    X=randn(ComplexF64, n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(ComplexF64, n, t)) for i=1:m]
    return 𝐗
end

# REAL data: check that for the case m=2 gCCA gives the same result as CCA
t, m, n, noise = 20, 2, 6, 0.1
Xset=getData(t, m, n, noise)
Cx=(Xset[1]*Xset[1]')/t
Cy=(Xset[2]*Xset[2]')/t
Cxy=(Xset[1]*Xset[2]')/t

gc=gcca(Xset; simple=true)

c=cca(Hermitian(Cx), Hermitian(Cy), Cxy; simple=true)

@test (c.F[1]'*Cxy*c.F[2]) ≈ (gc.F[1]'*Cxy*gc.F[2])
@test gc.F[1]'*Cx*gc.F[1]≈I
@test gc.F[2]'*Cy*gc.F[2]≈I
D=gc.F[1]'*Cxy*gc.F[2]
@test norm(D-Diagonal(D))+1≈1.

# COMPLEX data: check that for the case m=2 gCCA gives the same result as CCA
t, m, n, noise = 20, 2, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)
Ccx=(Xcset[1]*Xcset[1]')/t
Ccy=(Xcset[2]*Xcset[2]')/t
Ccxy=(Xcset[1]*Xcset[2]')/t

gcc=gcca(Xcset; simple=true, maxiter=10000)
cc=cca(Hermitian(Ccx), Hermitian(Ccy), Ccxy; simple=true)

# for complex data just do a sanity check as the order of vectors
# is arbitrary. The following two tests currently fail
# @test spForm(cc.F[1]'gcc.F[1])<0.001
# @test spForm(cc.F[2]'gcc.F[2])<0.001

@test gcc.F[1]'*Ccx*gcc.F[1]≈I
@test gcc.F[2]'*Ccy*gcc.F[2]≈I
# sanity check only as there is noise in the complex case
D=gcc.F[1]'*Ccxy*gcc.F[2]
@test norm(D-Diagonal(D))/(n^2-n)<0.001

# REAL data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xset=getData(t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gc=gcca(Xset, eVar=0.9)

# name of the filter
gc.name

𝒞=Array{Matrix}(undef, 1, m, m)
for i=1:m, j=1:m 𝒞[1, i, j]=(Xset[i]*Xset[j]')/t end

using Plots
# plot regularized accumulated eigenvalues
plot(gc.arev)

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all cross-covariances in a single matrix of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

 C=𝒞2Mat(𝒞, m, 1)
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-covariances")
 𝒮=_rotate_crossCov(gc.F, 𝒞, m, 1)
 S=𝒞2Mat(𝒮, m, 1)
 h2 = heatmap(S, clim=(0, 1), yflip=true, c=:amp, title="all rotated cross-covariances")
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FiggCCA.png")

# COMPLEX data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gcc=gcca(Xcset, eVar=0.9)


text:
using Photometry.Aperture: edges,
                           bbox

@testset "Apertures" begin
    ap_rect = RectangularAperture(50, 40, 10, 10, 0)
    @test edges(ap_rect) == (-5.5, 5.5, -5.5, 5.5)
    @test bbox(ap_rect) == (45, 55, 35, 45)
    @test size(ap_rect) == (11, 11)
    @test RectangularAperture([50, 40], 10, 10, 0) == ap_rect

    ap_ann = RectangularAnnulus(50, 40, 5, 10, 10, 0)
    @test edges(ap_ann) == (-5.5, 5.5, -5.5, 5.5)
    @test bbox(ap_ann) == (45, 55, 35, 45)
    @test size(ap_ann) == (11, 11)
    @test RectangularAnnulus([50, 40], 5, 10, 10, 0) == ap_ann
end

@testset "Rectangle Aperture" begin
    ap0 = RectangularAperture(0, 0, 0, 0, 0)

    @test sprint(show, ap0) == "RectangularAperture(0, 0, w=0, h=0, θ=0°)"

    ap1 = RectangularAperture(0, 0, 1, 1, 0)

    @test sprint(show, ap1) == "RectangularAperture(0, 0, w=1, h=1, θ=0°)"

    @test_throws ErrorException RectangularAperture(0, 0, -12, 2, 0)
    @test_throws ErrorException RectangularAperture(0, 0, 5, -2, 0)
end

@testset "Rectangle Annulus" begin

    ap1 = RectangularAnnulus(0, 0, 1, 1, 1, 0)

    @test sprint(show, ap1) == "RectangularAnnulus(0.0, 0.0, w_in=1.0, w_out=1.0, h_in=1.0, h_out=1.0, θ=0.0°)"

    @test_throws ErrorException RectangularAnnulus(0, 0, 0, 0, 0, 0)
    @test_throws ErrorException RectangularAnnulus(0, 0, -12, 2, 3, 0)
    @test_throws ErrorException RectangularAnnulus(0, 0, 4, 2, 3, 0)
    @test_throws ErrorException RectangularAnnulus(0, 0, 5, -2, 4, 0)
    @test_throws ErrorException RectangularAnnulus(0, 0, 5, 8, -4, 0)
end


text:
__precompile__()

module MDP

include("Preprocessing/Preprocessing.jl");
include("Gencode/Gencode.jl");
include("Postprocessing/Postprocessing.jl");

export initializeapp, initializeconfig, initializemdp, preprocessing, readconfig, readEXTXYZ, readJSON, readJSONfromFolder
export setlattice, setregion, setdomain

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
"""
Tools for the statistical analysis of Monte Carlo data.

### Exports:
- [`blocking_analysis`](@ref)
- [`ratio_of_means`](@ref)
- [`growth_witness`](@ref)
- [`smoothen`](@ref)
- [`shift_estimator`](@ref)
- [`projected_energy`](@ref)
- [`growth_estimator`](@ref)
- [`mixed_estimator`](@ref)
- [`val_and_errs`](@ref)
- [`val`](@ref)
- [`mean_and_se`](@ref)
"""
module StatsTools

using Statistics, MonteCarloMeasurements, Distributions, DataFrames
using StrLiterals, StrFormat # for Base.show() methods
using Random, SpecialFunctions, LinearAlgebra
import MacroTools
import Measurements

import Statistics: cov
import Measurements: measurement
import MonteCarloMeasurements: Particles
import Base: show, NamedTuple

export growth_witness, smoothen
export blocking_analysis, mean_and_se
export ratio_of_means
export replica_fidelity, med_and_errs, ratio_with_errs, to_measurement
export growth_estimator, mixed_estimator, w_lin, w_exp, projected_energy, shift_estimator
export pmedian, pquantile, pmiddle, piterate, pextrema, pminimum, pmaximum, pmean, pcov
export val_and_errs, errs, val

include("growth_witness.jl")
include("variances.jl")
include("blocking.jl")
include("ratio_of_means.jl")
include("convenience.jl")
include("fidelity.jl")
include("reweighting.jl")

end  # module StatsTools


text:
ContractAttributes() = ContractAttributes(Dict{String, Any}())

Base.setindex!(a::ContractAttributes, x, y) = setindex!(a.attributes, x, y)
Base.getindex(a::ContractAttributes, k) = getindex(a.attributes, k)
Base.haskey(a::ContractAttributes, k) = haskey(a.attributes, k)
Base.keys(a::ContractAttributes) = keys(a.attributes)

function ContractAttributes(x::Pair)
    result = ContractAttributes()
    result[string(x.first)] = x.second
    return result
end

function ContractAttributes(x...)
    result = ContractAttributes()

    for i in x
        @assert isa(i, Pair)
        result[string(i.first)] = i.second
    end

    return result
end

function Base.:(==)(a1::ContractAttributes, a2::ContractAttributes)
    if isempty(a1.attributes) && isempty(a2.attributes)
        return true
    end

    if (isempty(a1.attributes) && !isempty(a2.attributes)) || (!isempty(a1.attributes) && isempty(a2.attributes))
        return false
    end

    for k in keys(a1.attributes)
        if !haskey(a2.attributes, k)
            return false
        end
    end

    for k in keys(a2.attributes)
        if !haskey(a1.attributes, k)
            return false
        end
    end

    for k in keys(a1.attributes)
        if a1.attributes[k] != a2.attributes[k]
            return false
        end
    end

    return true
end

Base.hash(a::ContractAttributes) = 101 + hash(a.attributes)

# attributes with dedicated API

get_currency(a::ContractAttributes) :: Currencies.Currency = Currencies.Currency(a["currency"])
get_maturity(a::ContractAttributes) :: Date = Date(a["maturity"])
get_issue(a::ContractAttributes) :: Date = Date(a["issue"])

"Returns `true` if `sc` has a carry curve defined in the contract attributes."
function has_carry_curve(attr::ContractAttributes, sc::SpotCurrency) :: Bool
    @assert haskey(attr, "carry_curves") "Contract has no `carry_curves` attribute."
    currency_name = string(Currencies.currency_symbol(sc))
    carry_curves = attr["carry_curves"]
    @assert isa(carry_curves, Dict) # maps Currency -> Curve

    # returns true if there's a carry curve for the given currency, or to any currency
    return haskey(carry_curves, currency_name) || haskey(carry_curves, "Any")
end

function get_carry_curve_symbol(attr::ContractAttributes, sc::SpotCurrency) :: Symbol
    @assert haskey(attr, "carry_curves") "Contract has no `carry_curves` attribute."
    carry_curves = attr["carry_curves"]
    cur_sym = Currencies.currency_symbol(sc)
    cur_str = String(cur_sym)
    @assert haskey(carry_curves, cur_str) || haskey(carry_curves, "Any") "Contract has no carry curve defined for the currency $cur_sym."

    if haskey(carry_curves, cur_str)
        return Symbol(carry_curves[cur_str])
    else
        return Symbol(carry_curves["Any"])
    end
end


text:
@testset "Weighting" begin
  @testset "UniformWeighting" begin
    d = RegularGrid(100,100)
    s = georef((z=rand(100,100),))
    for o in [d, s]
      w = weight(o, UniformWeighting())
      @test length(w) == nelms(o)
      @test all(w .== 1)
    end
  end

  @testset "BlockWeighting" begin
    d = RegularGrid(100,100)
    s = georef((z=rand(100,100),))
    for o in [d, s]
      w = weight(o, BlockWeighting(10,10))
      @test length(unique(w)) == 1
      @test w[1] == 1 / 100
    end
  end

  @testset "DensityRatioWeighting" begin
    n = 1000; Random.seed!(123)

    r1 = Normal(0, 2)
    r2 = MixtureModel([Normal(-2,1), Normal(2,2)], [0.2, 0.8])

    z1 = sort(rand(r1, n))
    z2 = sort(rand(r2, n))

    d1 = georef(DataFrame(z=z1), PointSet(reshape(1:n,1,:)))
    d2 = georef(DataFrame(z=z2), PointSet(reshape(1:n,1,:)))

    w = weight(d1, DensityRatioWeighting(d2))

    if visualtests
      @plottest begin
        plot( z1, pdf.(r1, z1), size=(800,400), label="source")
        plot!(z1, pdf.(r2, z1), label="target")
        plot!(z1, w .* pdf.(r1, z1), label="approx")
      end joinpath(datadir,"density-ratio.png") !istravis
    end
  end
end


text:
using BinaryTreeTassellation
using Test

@testset "BinaryTreeTassellation.jl" begin
    # Write your tests here.
end


text:
let
    epc   = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe    = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    R     = rECItoECEF(epc)

    a_grav_point = accel_point_mass(x[1:3])

    a_grav_00 = accel_gravity(x, R, 0, 0)

    tol=1e-9
    @test isapprox(a_grav_00[1], a_grav_point[1], atol=tol)
    @test isapprox(a_grav_00[2], a_grav_point[2], atol=tol)
    @test isapprox(a_grav_00[3], a_grav_point[3], atol=tol)

    a_grav_20 = accel_gravity(x, R, 20, 20)
    @test norm(a_grav_20) > norm(a_grav_00)
end

let
    # Define Initial State
    epc   = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe    = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    r_sun = sun_position(epc)

    # Call function
    a_sun = accel_thirdbody_sun(epc, x)

    @test a_sun[1] < 1e-5
    @test a_sun[2] < 1e-5
    @test a_sun[3] < 1e-5

    a_sun2 = accel_thirdbody_sun(x, r_sun)

    @test a_sun[1] == a_sun2[1]
    @test a_sun[2] == a_sun2[2]
    @test a_sun[3] == a_sun2[3]
end

let
    # Define Initial State
    epc    = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe     = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x      = sOSCtoCART(oe, use_degrees = true)
    r_moon = moon_position(epc)

    # Call function
    a_moon = accel_thirdbody_moon(epc, x)

    @test a_moon[1] < 1e-5
    @test a_moon[2] < 1e-5
    @test a_moon[3] < 1e-5

    a_moon2 = accel_thirdbody_moon(x, r_moon)

    @test a_moon[1] == a_moon2[1]
    @test a_moon[2] == a_moon2[2]
    @test a_moon[3] == a_moon2[3]
end

let
    # Define Initial State
    epc   = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe    = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    T     = rECItoECEF(epc)

    # Call function
    rho    = density_harris_priester(epc, x)
    a_drag = accel_drag(x, rho, 1, 1, 2.3, T)

    @test a_drag[1] < 1e-5
    @test a_drag[2] < 1e-5
    @test a_drag[3] < 1e-5
end

let
    # Define Initial State
    epc   = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe    = [R_EARTH + 500e3, 0, 0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    r_sun = sun_position(epc)

    # Call function
    nu = eclipse_cylindrical(x, r_sun)
    @test nu == 1.0
    @test nu == eclipse_cylindrical(epc, x)

    oe    = [R_EARTH + 500e3, 0, 0, 180.0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    nu = eclipse_cylindrical(x, r_sun)
    @test nu == 0.0
    @test nu == eclipse_cylindrical(epc, x)
end

let
    # Define Initial State
    epc   = Epoch(2018, 3, 20, 16, 15, 0) # Test on Vernal equinox
    oe    = [R_EARTH + 500e3, 0, 0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    r_sun = sun_position(epc)

    # Call function
    nu = eclipse_conical(x, r_sun)
    @test nu == 0.0
    @test nu == eclipse_conical(epc, x)

    oe    = [R_EARTH + 500e3, 0, 0, 180.0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    nu = eclipse_conical(x, r_sun)
    @test nu == 1.0
    @test nu == eclipse_conical(epc, x)
end

let
    # Define Initial State
    epc   = Epoch(2018, 1, 1, 0, 0, 0)
    oe    = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x     = sOSCtoCART(oe, use_degrees = true)
    r_sun = sun_position(epc)

    # Call function
    a_srp = accel_srp(x, r_sun, 1, 1)

    @test a_srp[1] < 1e-5
    @test a_srp[2] < 1e-5
    @test a_srp[3] < 1e-5
end

let
    # Define Initial State
    oe = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    x  = sOSCtoCART(oe, use_degrees=true)

    # Call function
    a_rel = accel_relativity(x)

    @test a_rel[1] < 1e-7
    @test a_rel[2] < 1e-7
    @test a_rel[3] < 1e-7
end


text:
using Pkg
Pkg.build()
Pkg.test(; coverage = true)


text:
using SpecialFunctions, BesselK, ForwardDiff, StaticArrays
import ForwardDiff: derivative

# Convenience tool that will remove allocs:
@inline pv(scale, range, v) = @SVector [scale, range, v]

# Raw Bessel functions. Note also that there is the `adbesselkxv` method in
# BesselK.jl that gives you (x^v)*besselk(v,x) directly, sometimes with at least
# a slight gain to accuracy and speed. That isn't bound here, but obviously you
# could just slightly tweak this version.
R_besselk(v, x)    = BesselK.adbesselk(v, x)
R_besselk_dx(v, x) = derivative(_x->R_besselk(v, _x), x)
R_besselk_dv(v, x) = derivative(_v->R_besselk(_v, x), v)
R_besselk_dx_dx(v, x) = derivative(_x->R_besselk_dx(v, _x), x)
R_besselk_dx_dv(v, x) = derivative(_x->R_besselk_dv(v, _x), x)
R_besselk_dv_dv(v, x) = derivative(_v->R_besselk_dv(_v, x), v)

# Unlike the Julia-native version, this takes a pre-computed distance instead of
# two coordinates x and y, because if we have to take them as straight
# Vector{Float64} items then that will make allocations in the kernel calls and
# totally kill performance. In Julia this is avoided by using StaticArrays, but
# I don't think I can ask R users to deal with that interface. If you care
# enough about performance and need the additional flexibility that that
# provides, you could extend this code...or just switch to Julia more fully.
#
# Note here that the parameter order is:
# (sigma (scale), rho (range), nu (smoothness))
#
# You can off course change all of this up, but then be careful to also update
# the book-keeping stuff in the derivative functions.
#
# TODO (cg 2022/01/02 17:54): what is the most common signature for Julia stuff
# here? As it stands, the parameters can be passed in as a tuple, but should the
# R version take all parameters as scalars to again avoid using straight
# Vector{Float64}?
function matern(dist, params)
  (sg, rho, nu) = params
  iszero(dist) && return sg*sg
  arg = sqrt(2*nu)*dist/rho
  (sg*sg*(2^(1-nu))/gamma(nu))*BesselK.adbesselkxv(nu, arg)
end

# First derivatives:
matern_d1(dist, p) = matern(dist, (sqrt(2*p[1]), p[2], p[3]))
matern_d2(dist, p) = derivative(_p->matern(dist,pv(p[1],   _p, p[3])), p[2])
matern_d3(dist, p) = derivative(_p->matern(dist,pv(p[1], p[2],   _p)), p[3])

# Second derivatives:
matern_d1_d1(dist, p) = derivative(_p->matern_d1(dist, pv(_p,   p[2], p[3])), p[1])
matern_d1_d2(dist, p) = derivative(_p->matern_d1(dist, pv(p[1],   _p, p[3])), p[2])
matern_d1_d3(dist, p) = derivative(_p->matern_d1(dist, pv(p[1], p[2],   _p)), p[3])
matern_d2_d2(dist, p) = derivative(_p->matern_d2(dist, pv(p[1],   _p, p[3])), p[2])
matern_d2_d3(dist, p) = derivative(_p->matern_d2(dist, pv(p[1], p[2],   _p)), p[3])
matern_d3_d3(dist, p) = derivative(_p->matern_d3(dist, pv(p[1], p[2],   _p)), p[3])


text:
import Random # overload Random.rand!

abstract type AbstractSystem{T} <: AbstractArray{T, 4} end
Base.IndexStyle(::Type{<:AbstractSystem}) = IndexLinear()
Base.size(sys::S) where {S <: AbstractSystem} = size(sys.sites)
Base.getindex(sys::S, i::Int) where {S <: AbstractSystem} = sys.sites[i]
Base.setindex!(sys::S, v, i::Int) where {S <: AbstractSystem} = setindex!(sys.sites, v, i)

@inline function eachcellindex(sys::S) where {S <: AbstractSystem}
    return eachcellindex(sys.lattice)
end

"""
    nbasis(sys::SpinSystem)
"""
@inline function nbasis(sys::S) where {S <: AbstractSystem}
    return nbasis(sys.lattice)
end

"""
Defines a collection of charges. Currently primarily used to test ewald
 summation calculations.
"""
mutable struct ChargeSystem <: AbstractSystem{Float64}
    lattice :: Lattice             # Definition of underlying lattice
    sites   :: Array{Float64, 4}   # Holds charges at each site
end

"""
Defines a collection of spins, as well as the Hamiltonian they interact under.
 This is the main type to interface with most of the package.
"""
mutable struct SpinSystem <: AbstractSystem{Vec3}
    lattice     :: Lattice            # Definition of underlying lattice
    hamiltonian :: HamiltonianCPU     # Contains all interactions present
    sites       :: Array{Vec3, 4}     # Holds actual spin variables: Axes are [Basis, CellA, CellB, CellC]
    sites_info  :: Vector{SiteInfo}   # Characterization of each basis site
end

"""
    ChargeSystem(lat::Lattice)

Construct a `ChargeSystem` on the given lattice, initialized to all zero charges.
"""
function ChargeSystem(lat::Lattice)
    sites_size = (nbasis(lat), lat.size...)
    sites = zeros(sites_size)

    return ChargeSystem(lat, sites)
end

function ChargeSystem(cryst::Crystal, latsize)
    lattice = Lattice(crystal, latsize)
    return ChargeSystem(lattice)
end

"""
    rand!(sys::ChargeSystem)

Sets charges to random values uniformly drawn from ``[-1, 1]``,
then shifted to charge-neutrality.
"""
function Random.rand!(sys::ChargeSystem)
    sys.sites .= 2 .* rand(Float64, size(sys.sites)) .- 1.
    sys.sites .-= sum(sys.sites) / length(sys.sites)
    return
end

"""
    propagate_site_info(cryst::Crystal, sites_info::Vector{SiteInfo})

Given an incomplete list of site information, propagates spin magnitudes and
symmetry-transformed g-tensors to all symmetry-equivalent sites. If SiteInfo is
not provided for a site, sets S=1/2 and g=2 for that site. Throws an error if
two symmetry-equivalent sites are provided in `sites_info`.
"""
function propagate_site_info(crystal::Crystal, sites_info::Vector{SiteInfo})
    # All sites not explicitly provided are by default S = 1/2, g = 2
    all_sites_info = [SiteInfo(i, 1/2, 2) for i in 1:nbasis(crystal)]

    specified_atoms = Int[]
    for siteinfo in sites_info
        @unpack site, S, g = siteinfo
        (sym_bs, sym_gs) = all_symmetry_related_couplings(crystal, Bond(site, site, [0,0,0]), g)
        for (sym_b, sym_g) in zip(sym_bs, sym_gs)
            sym_atom = sym_b.i
            if sym_atom in specified_atoms
                # Perhaps this should only throw if two _conflicting_ SiteInfo are passed?
                # Then propagate_site_info can be the identity on an already-filled list.
                error("Provided two `SiteInfo` which describe symmetry-equivalent sites!")
            else
                push!(specified_atoms, sym_atom)
            end

            all_sites_info[sym_atom] = SiteInfo(sym_atom, S, sym_g)
        end
    end

    return all_sites_info
end

"""
    SpinSystem(crystal::Crystal, ints::Vector{<:AbstractInteraction}, latsize, sites_info::Vector{SiteInfo}=[];
               μB, μ0)

Construct a `SpinSystem` with spins of magnitude `S` residing on the lattice sites
 of a given `crystal`, interactions given by `ints`, and the number of unit cells along
 each lattice vector specified by `latsize`. Initialized to all spins pointing along
 the ``+𝐳̂`` direction. μB and μ0 set the Bohr magneton and vacuum permeability. By
 default, these are set so that the unit system is (meV, T, Å).
"""
function SpinSystem(crystal::Crystal, ints::Vector{<:AbstractInteraction}, latsize,
                    sites_info::Vector{SiteInfo}=SiteInfo[]; μB=BOHR_MAGNETON, μ0=VACUUM_PERM)
    latsize = collect(Int64.(latsize))
    lattice = Lattice(crystal, latsize)

    all_sites_info = propagate_site_info(crystal, sites_info)
    ℋ_CPU = HamiltonianCPU(ints, crystal, latsize, all_sites_info; μB, μ0)

    # Initialize sites to all spins along +z
    sites_size = (nbasis(lattice), lattice.size...)
    sites = fill(SA[0.0, 0.0, 1.0], sites_size)

    # Default unit system is (meV, K, Å, T)
    SpinSystem(lattice, ℋ_CPU, sites, all_sites_info)
end

function Base.show(io::IO, ::MIME"text/plain", sys::SpinSystem)
    printstyled(io, "Spin System\n"; bold=true, color=:underline)
    sz = size(sys.sites)
    println(io, "Basis $(sz[1]), Lattice dimensions $(sz[2:end])")
end

"""
    rand!(sys::SpinSystem)

Sets spins randomly sampled on the unit sphere.
"""
function Random.rand!(sys::SpinSystem)
    sys.sites .= randn(Vec3, size(sys.sites))
    @. sys.sites /= norm(sys.sites)
    return
end

"""
    randflips!(sys::SpinSystem)

Sets spins randomly either aligned or anti-aligned
with their original direction.
"""
function randflips!(sys::SpinSystem)
    sys.sites .*= rand((-1, 1), size(sys))
end

"""
    energy(sys::SpinSystem)

Computes the energy of the system under `sys.hamiltonian`.
"""
energy(sys::SpinSystem) = energy(sys.sites, sys.hamiltonian)

"""
    field!(B::Array{Vec3}, sys::SpinSystem)

Updates B in-place to contain the local field at each site in the
system under `sys.hamiltonian`. The "local field" is defined as

``𝐁_i = -∇_{𝐬_i} ℋ / S_i``

with ``𝐬_i`` the unit-vector variable at site i, and ``S_i`` is
the magnitude of the associated spin.
"""
field!(B::Array{Vec3}, sys::SpinSystem) = field!(B, sys.sites, sys.hamiltonian)

"""
    field(sys::SpinSystem)

Compute the local field B at each site of the system under
`sys.hamiltonian`.
"""
@inline function field(sys::SpinSystem)
    B = zero(sys)
    field!(B, sys)
    B
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
model = get_model("hopper_2D")
nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.b
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

# get trajectory
ref_traj = get_trajectory("hopper_2D", "gait_in_place_high", load_type=:joint_traj)

H = ref_traj.H
h = ref_traj.h
κ = 1.0e-8

# Cost function
obj = TrackingObjective(H, model.dim,
    q = [[Diagonal(1.0e-5 * [0.01, 1,1,1]) for t = 1:H-1]; [Diagonal(1.0e-4 * ones(nq))]],
    u = [Diagonal(1.0e-5 * [0.1, 1.0]) for t = 1:H],
    γ = [Diagonal(1.0e-100 * ones(nc)) for t = 1:H],
    b = [Diagonal(1.0e-100 * ones(nb)) for t = 1:H])
n_opts = NewtonOptions(r_tol=3e-8, max_iter=100)
im_traj = ImplicitTraj(ref_traj, model, κ=κ)
core = Newton(H, h, model, ref_traj, im_traj, obj = obj, opts = n_opts)

ϕ_fast(model, ref_traj.q[end])

q0_dist = deepcopy(ref_traj.q[1] + [-0.25,0.0,0,0])
q1_dist = deepcopy(ref_traj.q[2] + [-0.25,0.0,0,0])
newton_solve!(core, model, im_traj, ref_traj, q0=q0_dist, q1=q1_dist, verbose=true)

plot_surface!(vis, model.env, n=200)
visualize_robot!(vis, model, core.traj.q)
visualize_robot!(vis, model, ref_traj.q)

plot(hcat(Vector.(ref_traj.u)...)')
plot(hcat(Vector.(core.traj.u)...)')
plot(hcat(Vector.(ref_traj.q)...)')
plot(hcat(Vector.(core.traj.q)...)')
# Check ~perfect loop
@show round.(core.traj.q[end-2] - ref_traj.q[1], digits=3)
@show round.(core.traj.q[end-1] - ref_traj.q[2], digits=3)

# Save trajectory
traj = deepcopy(core.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_2D", "gaits", "gait_forward_high.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = get_trajectory("hopper_2D", "gait_forward_high", load_type=:joint_traj)
plot(hcat(Vector.(traj.q)...)')


text:
```
ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок и маркеров.
РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.

Рассмотреть отдельно еще случай, когда изначально в некоторых клетках поля могут находиться маркеры.
```

```
Если на поле будут стоять маркеры, то они не повлияют на ход программы. Иначе можно было бы добавить проверку
на наличие маркера в клетке, стоящей на пути Робота.
```

using HorizonSideRobots
originRobot = Robot(animate = true, "1example.sit")

# get_line! создаёт линию из маркеров в заданном направлении и возвращает Робота в исходное положение
function get_line!(robot, side)::Nothing
    steps = 0
    while ((isborder(robot,side))==false)
        move!(robot,side)
        putmarker!(robot)
        steps += 1
    end
    for j in 1:steps
        move!(robot, HorizonSide((Int(side)+2)%4))
    end
end

# Point1_master! главная функция программы
function Point1_master!(robot::Robot)::Nothing
    for side in [Nord, Ost, Sud, West]
        get_line!(robot, side)
    end
    putmarker!(robot)
end

Point1_master!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
function bott(eigvec::AbstractMatrix{<:Number}, tb::TBModel,
              outputRange::AbstractRange=Base.OneTo(size(eigvec,2)))

    # computes exp(iθx) & exp(iθy) with θx ≡ x/2π & θy ≡ y/2π and [x y] = tb.pos[:,1:2],
    # and then block-repeats it for cases with non-unity basis sizes (assumed split over a
    # similar block structure). This implicitly assumes x and y to be ∈[0,1]
    expx = repeat(cis.(tb.pos[:,1].*(2π)), tb.basisdim)
    expy = repeat(cis.(tb.pos[:,2].*(2π)), tb.basisdim)

    if outputRange != 1:size(eigvec, 2)
        eigvec = @view eigvec[:,outputRange]
    end

    X = eigvec'*(expx.*eigvec)
    Y = eigvec'*(expy.*eigvec)

    β = zeros(Float64, length(outputRange))
    @showprogress 2 "Computing Bott indices ..." for idx in eachindex(outputRange)
        # Not worth it to do views here, I think; BLAS is more important than GC...
        X′ = X[Base.OneTo(idx), Base.OneTo(idx)]
        Y′ = Y[Base.OneTo(idx), Base.OneTo(idx)]
        β[idx] = imag(sum(log.(eigvals(Y′*X′*Y′'*X′'))))/(2*pi)
    end

    return β
end

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
nwk = readstring("src/hd_consensus_tree.nwk")

function find_start(str::String, taxon::String, levels_up::Int64)
    start = findfirst(taxon, str)[1]
    start_ = findprev("(", str, start)[1]
    counter = 0
    while counter < levels_up
        if str[start_-1] == ')'
            counter -= 1
            start_ -= 1
        elseif str[start_-1] == '('
            counter += 1
            start_ -= 1
        else
            start_ -= 1
        end
    end
    return start_
end

function find_end(str::String, taxon::String, levels_up::Int64)
    start = findfirst(taxon, str)[1]
    start_search_1 = findnext(")", str, start)[1]
    start_search_2 = findnext("(", str, start)[1]
    if start_search_1 < start_search_2
        end_ = start_search_1
        counter = 0
        while counter < levels_up
            if str[end_+1] == '('
                counter -= 1
                end_ += 1
            elseif str[end_+1] == ')'
                counter += 1
                end_ += 1
            else
                end_ += 1
            end
        end
    else
        end_ = start_search_2 - 2
        counter = 0
        while counter < levels_up + 1
            if str[end_+1] == '('
                counter -= 1
                end_ += 1
            elseif str[end_+1] == ')'
                counter += 1
                end_ += 1
            else
                end_ += 1
            end
        end
    end

    return end_

end

function get_subset_string(str::String, taxon::String, levels_up::String)
    levels_up_int = parse(Int64, levels_up)
    start_ = find_start(str, taxon, levels_up_int)
    end_ = find_end(str, taxon, levels_up_int)
    return str[start_:end_] * ';'
end

function replace_symbols!(str::String)
    str = replace(str, "(", "")
    str = replace(str, ")", "")
    str = replace(str, ",", "")
    str = replace(str, "\'", "")
end

drugFile = open("src/drugid_names.txt")

dname_dict  = Dict{Any, Any}()
for ln in eachline(drugFile)
    ln   = chomp(ln)
    txn  = split(ln, "\$")[1]
    name = split(ln, "\$")[2]
    name = replace(name, r"\,", "")
    dname_dict[txn] = replace_symbols!(ucfirst(name))
end

function replace_name!(s::String, map::Dict{Any, Any})

    begin_ = findnext("_", s, 1)[1]

    while begin_ < length(s)
        end_ = findnext(":", s, begin_)
        global str_ = s[begin_-2:end_[1]-1]
        name = map[str_]

        s = s[1:begin_[1]-3] * name * s[end_[1]:length(s)]
        try
            begin_ = findnext("_", s, end_[1])[1]
        catch
            break
        end
    end

    return s

end

path_file = open("src/hd_tree_taxon_path.txt")

path_dict  = Dict{Any, Any}()
for ln in eachline(path_file)
    ln   = chomp(ln)
    txn  = split(ln, "\$")[1]
    path_ = length(split(split(ln, "\$")[2], "\."))
    max_level = path_ - 8

    path_dict[txn] = path_
end

function check_level(str::String, levels_up::String)
    max_level = path_dict[str] - 9 / 2
    if levels_up > max_level
        levels_up = max_level
    end
    return level_up
end

function main(str::String, taxon::String, levels_up::String, map::Dict{Any, Any})
    final_string = get_subset_string(str, taxon, levels_up)
    return final_string#replace_name!(final_string, map)
end

max_stars_repo_name: bcbi/PhytoKB-API

max_issues_repo_name: bcbi/PhytoKB-API

max_forks_repo_name: bcbi/PhytoKB-API


text:
# using PyCall
# using Pkg; Pkg.add("PyPlot")

push!(LOAD_PATH, "../src/")
using Documenter, AdFem
makedocs(sitename="AdFem", modules=[AdFem],
    pages = Any[
        "index.md",
        "Tutorial"=>Any["gallery.md"],
        "Forward Computation"=>Any["coupled.md", "staticelasticity.md", "plasticity.md",
            "viscoelasticity.md", "viscoelasticity_earth.md", "earthquake.md",
            "heatequation.md","elastodynamics.md", "twophaseflow.md", "NavierStokes2.md",
            "SteadyStateNavierStokes.md","SteadyStateNavierStokes3D.md", "CHTCoupled.md",
            "fwd_mixed_poisson.md",
            "fwd_linear_elasticity.md",
            "fwd_stress_based_viscoelasticity.md", "pml.md"],
        "Inverse Modeling"=>Any["inverse.md", "inv_viscoelasticity.md", "coupled_viscoelasticity.md",
            "inv_twophaseflow.md", "inv_viscoelasticity_nonparametric.md", "inv_viscoelasticity_earth.md"],
        "Advanced Topics"=>Any["mfem_tutorial.md", "mfem_mesh.md", "dev_unstructured.md", "BDMElement.md", "mfem3d.md"],
        "api.md",
        "docker_install_guide.md"
    ],
    authors = "Kailai Xu"
)

deploydocs(
    repo = "github.com/kailaix/AdFem.jl.git",
)


text:
function molecule(::Molecule"H₂")
  return """
  H 0.0 0.0 0.0
  H 0.76 0.0 0.0
  """
end
diis(::Molecule"H₂") = false
oda(::Molecule"H₂") = false


text:
include("../src/RandomField.jl")
using .RandomField

# choose parameters:
## remember, we need that 4|n_i for i=1,...,d
N = [10,20,30]		# number pixels per dimension
L = [80.0,80.0,80.0]	# real space box is [-l,l]^d
f(V::Vector{T} where T<:Real) = rect_rf(V)

# derived parameters
d = dim(N) # dimension
H = (pi/2) .* N ./L

# allocate memory
V = Cloud(d,N)
X = Cloud(d,N)
F = Cscalar(V)
G = Cscalar(X)

init!(X,L)
init!(V,H)

tovtk(X, "../data/x.cl.vtk")

Xx = slice(X,1)
tovtk(Xx, "../data/xx.cl.vtk")

Xy = slice(X,2)
tovtk(Xy, "../data/xy.cl.vtk")

Xz = slice(X,3)
tovtk(Xz, "../data/xz.cl.vtk")

sample!(F,f)

tovtk(F, "../data/v.F.vtk")

Fx = slice(F,1)
tovtk(Fx, "../data/v.f1.vtk")

Fy = slice(F,2)
tovtk(Fy, "../data/v.f2.vtk")

Fz = slice(F,3)
tovtk(Fz, "../data/v.f3.vtk")
tocsv(Fz, "../data/v.f3.csv")


text:
# # Daisyworld
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../daisyworld.mp4" type="video/mp4">
# </video>
# ```
#
# Study this example to learn about
# - Simple agent properties with complex model interactions
# - Diffusion of a quantity in a `GridSpace`
# - Including a "surface property" in the model
# - counting time in the model and having time-dependent dynamics
# - performing interactive scientific research
#
# ## Overview of Daisyworld
#
# This model explores the [Gaia hypothesis](https://en.wikipedia.org/wiki/Gaia_hypothesis),
# which considers the Earth as a single, self-regulating system including both living and
# non-living parts.
#
# Daisyworld is filled with black and white daisies.
# Their albedo's differ, with black daisies absorbing light and heat,
# warming the area around them; white daisies doing the opposite.
# Daisies can only reproduce within a certain temperature range, meaning too much
# (or too little) heat coming from the sun and/or surrounds will ultimately halt daisy
# propagation.
#
# When the climate is too cold it is necessary for the black daisies to propagate in order
# to raise the temperature, and vice versa -- when the climate is too warm, it is
# necessary for more white daisies to be produced in order to cool the temperature.
# The interplay of the living and non living aspects of this world manages to find an
# equilibrium over a wide range of parameter settings, although with enough external
# forcing, the daisies will not be able to regulate the temperature of the planet
# and eventually go extinct.

# ## Defining the agent types

# `Daisy` has three values (other than the required
# `id` and `pos` for an agent that lives on a [`GridSpace`](@ref). Each daisy has an `age`,
# confined later by a maximum age set by the user, a `breed` (either `:black` or `:white`)
# and an associated `albedo` value, again set by the user.
# `Land` represents the surface. We could make `Land` also have an albedo field, but
# in this world, the entire surface has the same albedo and thus we make it a model parameter.

# Notice that the surface is not an agent, as per design of Agents.jl, but rather a
# standard Julia array.

using Agents
using Statistics: mean
using Random # hide

mutable struct Daisy <: AbstractAgent
    id::Int
    pos::Dims{2}
    breed::Symbol
    age::Int
    albedo::Float64 # 0-1 fraction
end

const DaisyWorld = ABM{<:GridSpace, Daisy};

# ## World heating

# The surface temperature of the world is heated by its sun, but daisies growing upon it
# absorb or reflect the starlight -- altering the local temperature.

function update_surface_temperature!(pos, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    absorbed_luminosity = if isempty(ids) # no daisy
        ## Set luminosity via surface albedo
        (1 - model.surface_albedo) * model.solar_luminosity
    else
        ## Set luminosity via daisy albedo
        (1 - model[ids[1]].albedo) * model.solar_luminosity
    end
    ## We expect local heating to be 80 ᵒC for an absorbed luminosity of 1,
    ## approximately 30 for 0.5 and approximately -273 for 0.01.
    local_heating = absorbed_luminosity > 0 ? 72 * log(absorbed_luminosity) + 80 : 80
    ## Surface temperature is the average of the current temperature and local heating.
    model.temperature[pos...] = (model.temperature[pos...] + local_heating) / 2
end

# In addition, temperature diffuses over time
function diffuse_temperature!(pos, model::DaisyWorld)
    ratio = get(model.properties, :ratio, 0.5) # diffusion ratio
    npos = nearby_positions(pos, model)
    model.temperature[pos...] =
        (1 - ratio) * model.temperature[pos...] +
        ## Each neighbor is giving up 1/8 of the diffused
        ## amount to each of *its* neighbors
        sum(model.temperature[p...] for p in npos) * 0.125 * ratio
end

# ## Daisy dynamics

# The final piece of the puzzle is the life-cycle of each daisy. This method defines an
# optimal temperature for growth. If the temperature gets too hot or too cold, daisies
# will not wish to propagate. So long as the temperature is favorable,
# daisies compete for land and attempt to spawn a new plant of their `breed` in locations
# close to them.

function propagate!(pos, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    if !isempty(ids)
        daisy = model[ids[1]]
        temperature = model.temperature[pos...]
        ## Set optimum growth rate to 22.5 ᵒC, with bounds of [5, 40]
        seed_threshold = (0.1457 * temperature - 0.0032 * temperature^2) - 0.6443
        if rand(model.rng) < seed_threshold
            ## Collect all adjacent position that have no daisies
            empty_neighbors = Tuple{Int,Int}[]
            neighbors = nearby_positions(pos, model)
            for n in neighbors
                if isempty(ids_in_position(n, model))
                    push!(empty_neighbors, n)
                end
            end
            if !isempty(empty_neighbors)
                ## Seed a new daisy in one of those position
                seeding_place = rand(model.rng, empty_neighbors)
                add_agent!(seeding_place, model, daisy.breed, 0, daisy.albedo)
            end
        end
    end
end

# And if the daisies cross an age threshold, they die out.
# Death is controlled by the `agent_step` function
function agent_step!(agent::Daisy, model::DaisyWorld)
    agent.age += 1
    agent.age >= model.max_age && kill_agent!(agent, model)
end

# The model step function advances Daisyworld's dynamics:

function model_step!(model)
    for p in positions(model)
        update_surface_temperature!(p, model)
        diffuse_temperature!(p, model)
        propagate!(p, model)
    end
    model.tick += 1
    solar_activity!(model)
end

# Notice that `solar_activity!` changes the incoming solar radiation over time,
# if the given "scenario" (a model parameter) is `:ramp`.
# The parameter `tick` of the model keeps track of time.

function solar_activity!(model::DaisyWorld)
    if model.scenario == :ramp
        if model.tick > 200 && model.tick <= 400
            model.solar_luminosity += model.solar_change
        end
        if model.tick > 500 && model.tick <= 750
            model.solar_luminosity -= model.solar_change / 2
        end
    elseif model.scenario == :change
        model.solar_luminosity += model.solar_change
    end
end

# ## Initialising Daisyworld

# Here, we construct a function to initialize a Daisyworld. We use [`fill_space!`](@ref)
# to fill the space with `Land` instances. Then, we need to know how many
# daisies of each type to seed the planet with and what their albedo's are.
# We also want a value for surface albedo, as well as solar intensity
# (and we also choose between constant or time-dependent intensity with `scenario`).

import StatsBase
import DrWatson: @dict
using Random

function daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 0.2, # % cover of the world surface of white breed
    init_black = 0.2, # % cover of the world surface of black breed
    albedo_white = 0.75,
    albedo_black = 0.25,
    surface_albedo = 0.4,
    solar_change = 0.005,
    solar_luminosity = 1.0, # initial luminosity
    scenario = :default,
    seed = 165,
)

    rng = MersenneTwister(seed)
    space = GridSpace(griddims)
    properties = @dict max_age surface_albedo solar_luminosity solar_change scenario
    properties[:tick] = 0
    properties[:temperature] = zeros(griddims)

    model = ABM(Daisy, space; properties, rng)

    ## Populate with daisies: each position has only one daisy (black or white)
    grid = collect(positions(model))
    num_positions = prod(griddims)
    white_positions =
        StatsBase.sample(grid, Int(init_white * num_positions); replace = false)
    for wp in white_positions
        wd = Daisy(nextid(model), wp, :white, rand(model.rng, 0:max_age), albedo_white)
        add_agent_pos!(wd, model)
    end
    allowed = setdiff(grid, white_positions)
    black_positions =
        StatsBase.sample(allowed, Int(init_black * num_positions); replace = false)
    for bp in black_positions
        wd = Daisy(nextid(model), bp, :black, rand(model.rng, 0:max_age), albedo_black)
        add_agent_pos!(wd, model)
    end

    ## Adjust temperature to initial daisy distribution
    for p in positions(model)
        update_surface_temperature!(p, model)
    end

    return model
end

# ## Visualizing & animating
# %% #src
# Lets run the model with constant solar isolation and visualize the result

using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

model = daisyworld()

# To visualize we need to define the necessary functions for [`abmplot`](@ref).
# We will also utilize its ability to plot an underlying heatmap,
# which will be the model surface temperature,
# while daisies will be plotted in black and white as per their breed.
# Notice that we will explicitly provide a `colorrange` to the heatmap keywords,
# otherwise the colormap will be continuously and automatically updated to match
# the underlying temperature values while we are animating the time evolution.

daisycolor(a::Daisy) = a.breed

plotkwargs = (
    ac=daisycolor, as = 12, am = '♠',
    heatarray = :temperature,
    heatkwargs = (colorrange = (-20, 60),),
)
fig, _ = abmplot(model; plotkwargs...)
fig

# And after a couple of steps
Agents.step!(model, agent_step!, model_step!, 5)
fig, _ = abmplot(model; heatarray = model.temperature, plotkwargs...)
fig

# Let's do some animation now
model = daisyworld()
abmvideo(
    "daisyworld.mp4",
    model,
    agent_step!,
    model_step!;
    title = "Daisy World",
    plotkwargs...,
)

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../daisyworld.mp4" type="video/mp4">
# </video>
# ```
# Running this animation for longer hints that this world achieves quasi-equilibrium
# for some input parameters, where one `breed` does not totally dominate the other.
# Of course we can check this easily through data collection.
# Notice that here we have to define a function `breed` that returns the daisy's `breed`
# field. We cannot use just `:breed` to automatically find it, because in this mixed
# agent model, the `Land` doesn't have any `breed`.
# %% #src
black(a) = a.breed == :black
white(a) = a.breed == :white
adata = [(black, count), (white, count)]

model = daisyworld(; solar_luminosity = 1.0)

agent_df, model_df = run!(model, agent_step!, model_step!, 1000; adata)
figure = Figure(resolution = (600, 400))
ax = figure[1, 1] = Axis(figure, xlabel = "tick", ylabel = "daisy count")
blackl = lines!(ax, agent_df[!, :step], agent_df[!, :count_black], color = :red)
whitel = lines!(ax, agent_df[!, :step], agent_df[!, :count_white], color = :blue)
figure[1, 2] = Legend(figure, [blackl, whitel], ["black", "white"], textsize = 12)
figure

# ## Time dependent dynamics
# %% #src

# To use the time-dependent dynamics we simply use the keyword `scenario = :ramp` during
# model creation. However, we also want to see how the planet surface temperature changes
# and would be nice to plot solar luminosity as well.
# Thus, we define in addition
temperature(model) = mean(model.temperature)
mdata = [temperature, :solar_luminosity]

# And we run (and plot) everything
model = daisyworld(solar_luminosity = 1.0, scenario = :ramp)
agent_df, model_df =
    run!(model, agent_step!, model_step!, 1000; adata = adata, mdata = mdata)

figure = CairoMakie.Figure(resolution = (600, 600))
ax1 = figure[1, 1] = Axis(figure, ylabel = "daisy count", textsize = 12)
blackl = lines!(ax1, agent_df[!, :step], agent_df[!, :count_black], color = :red)
whitel = lines!(ax1, agent_df[!, :step], agent_df[!, :count_white], color = :blue)
figure[1, 2] = Legend(figure, [blackl, whitel], ["black", "white"], textsize = 12)

ax2 = figure[2, 1] = Axis(figure, ylabel = "temperature", textsize = 12)
ax3 = figure[3, 1] = Axis(figure, xlabel = "tick", ylabel = "L", textsize = 12)
lines!(ax2, model_df[!, :step], model_df[!, :temperature], color = :red)
lines!(ax3, model_df[!, :step], model_df[!, :solar_luminosity], color = :red)
for ax in (ax1, ax2); ax.xticklabelsvisible = false; end
figure


text:
using LocalStore

# Define an item type. This should contain enough information
# to be able to reconstruct the desired object.

struct A
    x::Int
end

# Since store location is based on hash, it is recommended to
# define a hash function for custom types.

Base.hash(a::A, h::UInt) = hash(:A, hash(a.x, h))

# Define how the item is stored on disk. The `path` is a dir, handled
# automatically by the LocalStore package.

LocalStore.save(a::A, path::String) = open("$path/file.txt", "w") do f
    ## some difficult computation you only want to do once
    x2 = a.x^2
    ## maybe download some remote files

    ## now save to disk the results
    write(f, string(x2))
end

# Define how the item is read from disk.

LocalStore.load(a::A, path::String) = open("$path/file.txt", "r") do f
    parse(Int, readline(f))
end

# Automatically checks if the requested object is stored. If not,
# it creates a local directory and saves it there and loads it.
# Next time `load` is called, the stored object is returned.

a = LocalStore.load(A(5))


text:
@testset "Screening functions" begin
    m = test_toyModel()

    # nothing to analyze
    @test_throws DomainError screen(m; analysis = identity)

    # array dimensionalities must match (this is actually caught by typechecker)
    @test_throws MethodError screen(
        m;
        analysis = identity,
        variants = [[], []],
        args = [() ()],
    )

    # sizes must match
    @test_throws DomainError screen(
        m;
        analysis = identity,
        variants = [[], []],
        args = [()],
    )

    # argument handling
    @test screen(m, analysis = identity, variants = [[]]) == [m]
    @test screen(m, analysis = identity, args = [()]) == [m]
    @test screen(m, analysis = (a, b) -> b, args = [(1,), (2,)]) == [1, 2]

    # test modifying some reactions
    quad_rxn(i) = (m::CoreModel) -> begin
        mm = copy(m)
        mm.S = copy(m.S)
        mm.S[:, i] .^= 2
        return mm
    end

    @test screen_variants(
        m,
        [[quad_rxn(i)] for i = 1:3],
        m -> flux_balance_analysis_vec(m, Tulip.Optimizer);
        workers = W,
    ) == [
        [250.0, -250.0, -1000.0, 250.0, 1000.0, 250.0, 250.0],
        [500.0, 500.0, 1000.0, 500.0, -1000.0, 500.0, 500.0],
        [500.0, 500.0, 1000.0, -500.0, 1000.0, 500.0, 500.0],
    ]

    # test solver modifications
    @test screen(
        m;
        analysis = (m, sense) -> flux_balance_analysis_vec(
            m,
            Tulip.Optimizer;
            modifications = [change_sense(sense)],
        ),
        args = [(MIN_SENSE,), (MAX_SENSE,)],
    ) == [
        [-500.0, -500.0, -1000.0, 500.0, 1000.0, -500.0, -500.0],
        [500.0, 500.0, 1000.0, -500.0, -1000.0, 500.0, 500.0],
    ]
end


text:
using Documenter
using NPFinancial

makedocs(
    modules = [NPFinancial],
    sitename="NPFinancial"
)

deploydocs(
    repo = "github.com/tk3369/NPFinancial.jl.git",
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    make   = () -> run(`mkdocs build`),
    target = "site"
)


text:
module API

import LLVM.API: LLVMValueRef, LLVMModuleRef, LLVMTypeRef, LLVMContextRef
using Enzyme_jll
using Libdl
using LLVM
using CEnum

const EnzymeLogicRef = Ptr{Cvoid}
const EnzymeTypeAnalysisRef = Ptr{Cvoid}
const EnzymeAugmentedReturnPtr = Ptr{Cvoid}

struct IntList
    data::Ptr{Int64}
    size::Csize_t
end
IntList() = IntList(Ptr{Int64}(0),0)

@cenum(CConcreteType,
  DT_Anything = 0,
  DT_Integer = 1,
  DT_Pointer = 2,
  DT_Half = 3,
  DT_Float = 4,
  DT_Double = 5,
  DT_Unknown = 6
)

struct EnzymeTypeTree end
const CTypeTreeRef = Ptr{EnzymeTypeTree}

EnzymeNewTypeTree() = ccall((:EnzymeNewTypeTree, libEnzyme), CTypeTreeRef, ())
EnzymeNewTypeTreeCT(T, ctx) = ccall((:EnzymeNewTypeTreeCT, libEnzyme), CTypeTreeRef, (CConcreteType, LLVMContextRef), T, ctx)
EnzymeNewTypeTreeTR(tt) = ccall((:EnzymeNewTypeTreeTR, libEnzyme), CTypeTreeRef, (CTypeTreeRef,), tt)

EnzymeFreeTypeTree(tt) = ccall((:EnzymeFreeTypeTree, libEnzyme), Cvoid, (CTypeTreeRef,), tt)
EnzymeSetTypeTree(dst, src) = ccall((:EnzymeSetTypeTree, libEnzyme), Cvoid, (CTypeTreeRef, CTypeTreeRef), dst, src)
EnzymeMergeTypeTree(dst, src) = ccall((:EnzymeMergeTypeTree, libEnzyme), Cvoid, (CTypeTreeRef, CTypeTreeRef), dst, src)
EnzymeTypeTreeOnlyEq(dst, x) = ccall((:EnzymeTypeTreeOnlyEq, libEnzyme), Cvoid, (CTypeTreeRef, Int64), dst, x)
EnzymeTypeTreeShiftIndiciesEq(dst, dl, offset, maxSize, addOffset) =
    ccall((:EnzymeTypeTreeShiftIndiciesEq, libEnzyme), Cvoid, (CTypeTreeRef, Cstring, Int64, Int64, UInt64),
        dst, dl, offset, maxSize, addOffset)

EnzymeTypeTreeToString(tt) = ccall((:EnzymeTypeTreeToString, libEnzyme), Cstring, (CTypeTreeRef,), tt)
EnzymeTypeTreeToStringFree(str) = ccall((:EnzymeTypeTreeToStringFree, libEnzyme), Cvoid, (Cstring,), str)

struct CFnTypeInfo
    arguments::Ptr{CTypeTreeRef}
    ret::CTypeTreeRef

    known_values::Ptr{IntList}
end

@cenum(CDIFFE_TYPE,
  DFT_OUT_DIFF = 0,  # add differential to an output struct
  DFT_DUP_ARG = 1,   # duplicate the argument and store differential inside
  DFT_CONSTANT = 2,  # no differential
  DFT_DUP_NONEED = 3 # duplicate this argument and store differential inside,
                     # but don't need the forward
)

# Create the derivative function itself.
#  \p todiff is the function to differentiate
#  \p retType is the activity info of the return
#  \p constant_args is the activity info of the arguments
#  \p returnValue is whether the primal's return should also be returned
#  \p dretUsed is whether the shadow return value should also be returned
#  \p additionalArg is the type (or null) of an additional type in the signature
#  to hold the tape.
#  \p typeInfo is the type info information about the calling context
#  \p _uncacheable_args marks whether an argument may be rewritten before loads in
#  the generated function (and thus cannot be cached).
#  \p augmented is the data structure created by prior call to an augmented forward
#  pass
#  \p AtomicAdd is whether to perform all adjoint updates to memory in an atomic way
#  \p PostOpt is whether to perform basic optimization of the function after synthesis
function EnzymeCreatePrimalAndGradient(logic, todiff, retType, constant_args, TA,
                                       returnValue, dretUsed, topLevel, additionalArg, typeInfo,
                                       uncacheable_args, augmented, atomicAdd, postOpt)
    ccall((:EnzymeCreatePrimalAndGradient, libEnzyme), LLVMValueRef,
        (EnzymeLogicRef, LLVMValueRef, CDIFFE_TYPE, Ptr{CDIFFE_TYPE}, Csize_t,
         EnzymeTypeAnalysisRef, UInt8, UInt8, UInt8, LLVMTypeRef, CFnTypeInfo,
         Ptr{UInt8}, Csize_t, EnzymeAugmentedReturnPtr, UInt8, UInt8),
        logic, todiff, retType, constant_args, length(constant_args), TA, returnValue,
        dretUsed, topLevel, additionalArg, typeInfo, uncacheable_args, length(uncacheable_args),
        augmented, atomicAdd, postOpt)
end

# Create an augmented forward pass.
#  \p todiff is the function to differentiate
#  \p retType is the activity info of the return
#  \p constant_args is the activity info of the arguments
#  \p returnUsed is whether the primal's return should also be returned
#  \p typeInfo is the type info information about the calling context
#  \p _uncacheable_args marks whether an argument may be rewritten before loads in
#  the generated function (and thus cannot be cached).
#  \p forceAnonymousTape forces the tape to be an i8* rather than the true tape structure
#  \p AtomicAdd is whether to perform all adjoint updates to memory in an atomic way
#  \p PostOpt is whether to perform basic optimization of the function after synthesis
function EnzymeCreateAugmentedPrimal(logic, todiff, retType, constant_args, TA,  returnUsed,
                                     typeInfo, uncacheable_args, forceAnonymousTape, atomicAdd, postOpt)
    ccall((:EnzymeCreateAugmentedPrimal, libEnzyme), EnzymeAugmentedReturnPtr,
        (EnzymeLogicRef, LLVMValueRef, CDIFFE_TYPE, Ptr{CDIFFE_TYPE}, Csize_t,
         EnzymeTypeAnalysisRef, UInt8,
         CFnTypeInfo, Ptr{UInt8}, Csize_t, UInt8, UInt8, UInt8),
        logic, todiff, retType, constant_args, length(constant_args), TA,  returnUsed,
        typeInfo, uncacheable_args, length(uncacheable_args), forceAnonymousTape, atomicAdd, postOpt)
end

# typedef bool (*CustomRuleType)(int /*direction*/, CTypeTree * /*return*/,
#                                CTypeTree * /*args*/, size_t /*numArgs*/,
#                                LLVMValueRef);
const CustomRuleType = Ptr{Cvoid}

function CreateTypeAnalysis(triple, rulenames, rules)
    @assert length(rulenames) == length(rules)
    ccall((:CreateTypeAnalysis, libEnzyme), EnzymeTypeAnalysisRef, (Cstring, Ptr{Cstring}, Ptr{CustomRuleType}, Csize_t), triple, rulenames, rules, length(rules))
end

function ClearTypeAnalysis(ta)
    ccall((:ClearTypeAnalysis, libEnzyme), Cvoid, (EnzymeTypeAnalysisRef,), ta)
end

function FreeTypeAnalysis(ta)
    ccall((:FreeTypeAnalysis, libEnzyme), Cvoid, (EnzymeTypeAnalysisRef,), ta)
end

function CreateLogic()
    ccall((:CreateEnzymeLogic, libEnzyme), EnzymeLogicRef, ())
end

function ClearLogic(logic)
    ccall((:ClearEnzymeLogic, libEnzyme), Cvoid, (EnzymeLogicRef,), logic)
end

function FreeLogic(logic)
    ccall((:FreeEnzymeLogic, libEnzyme), Cvoid, (EnzymeLogicRef,), logic)
end

function EnzymeExtractReturnInfo(ret, data, existed)
    @assert length(data) == length(existed)
    ccall((:EnzymeExtractReturnInfo, libEnzyme),
           Cvoid, (EnzymeAugmentedReturnPtr, Ptr{Int64}, Ptr{UInt8}, Csize_t),
           ret, data, existed, length(data))
end

function EnzymeExtractFunctionFromAugmentation(ret)
    ccall((:EnzymeExtractFunctionFromAugmentation, libEnzyme), LLVMValueRef, (EnzymeAugmentedReturnPtr,), ret)
end

function EnzymeExtractTapeTypeFromAugmentation(ret)
    ccall((:EnzymeExtractTapeTypeFromAugmentation, libEnzyme), LLVMTypeRef, (EnzymeAugmentedReturnPtr,), ret)
end

import Libdl
function EnzymeSetCLBool(name, val)
    handle = Libdl.dlopen(libEnzyme)
    ptr = Libdl.dlsym(handle, name)
    ccall((:EnzymeSetCLBool, libEnzyme), Cvoid, (Ptr{Cvoid}, UInt8), ptr, val)
end
function EnzymeGetCLBool(name)
    handle = Libdl.dlopen(libEnzyme)
    ptr = Libdl.dlsym(handle, name)
    ccall((:EnzymeGetCLBool, libEnzyme), UInt8, (Ptr{Cvoid},), ptr)
end
# void EnzymeSetCLInteger(void *, int64_t);

function printperf!(val)
    ptr = cglobal((:EnzymePrintPerf, libEnzyme))
    ccall((:EnzymeSetCLBool, libEnzyme), Cvoid, (Ptr{Cvoid}, UInt8), ptr, val)
end

function printtype!(val)
    ptr = cglobal((:EnzymePrintType, libEnzyme))
    ccall((:EnzymeSetCLBool, libEnzyme), Cvoid, (Ptr{Cvoid}, UInt8), ptr, val)
end

function printall!(val)
    ptr = cglobal((:EnzymePrint, libEnzyme))
    ccall((:EnzymeSetCLBool, libEnzyme), Cvoid, (Ptr{Cvoid}, UInt8), ptr, val)
end

function maxtypeoffset!(val)
    ptr = cglobal((:MaxTypeOffset, libEnzyme))
    ccall((:EnzymeSetCLInteger, libEnzyme), Cvoid, (Ptr{Cvoid}, Int64), ptr, val)
end

function looseTypeAnalysis!(val)
    ptr = cglobal((:looseTypeAnalysis, libEnzyme))
    ccall((:EnzymeSetCLInteger, libEnzyme), Cvoid, (Ptr{Cvoid}, UInt8), ptr, val)
end

function EnzymeRemoveTrivialAtomicIncrements(func)
    ccall((:EnzymeRemoveTrivialAtomicIncrements, libEnzyme), Cvoid, (LLVMValueRef,), func)
end

function EnzymeAddAttributorLegacyPass(PM)
    ccall((:EnzymeAddAttributorLegacyPass, libEnzyme),Cvoid,(LLVM.API.LLVMPassManagerRef,), PM)
end

end


text:
using JSON, RelocatableFolders

include("./bundlepaths.jl")
download_js_bundles()

const CONFIG_BEGIN_MARKER = "###JULIA-WEBIO-CONFIG-BEGIN"
const CONFIG_END_MARKER = "###JULIA-WEBIO-CONFIG-END"

macro meh(expr)
    quote
        try
            $(esc(expr))
        catch exc
            @info "Surpressed exception (this is probably harmless)" expr=$(string(expr)) exception=exc
        end
    end
end

try
    @info "Attempting to uninstall old Jupyter integrations (it's probably okay if error messages appear below this point)"
    jupyter_config_dir = get(ENV, "JUPYTER_CONFIG_DIR", joinpath(homedir(), ".jupyter"))

    config_file_py = joinpath(jupyter_config_dir, "jupyter_notebook_config.py")
    if isfile(config_file_py)
        config_regex = Regex(
            "\n?" * CONFIG_BEGIN_MARKER * ".*" * CONFIG_END_MARKER * "\n?",
            "s",
        )
        config_py = read(config_file_py, String)
        config_py = replace(config_py, config_regex => "")
        write(config_file_py, config_py)
    end

    config_file_json = joinpath(jupyter_config_dir, "jupyter_notebook_config.json")
    if isfile(config_file_json)
        config_json = JSON.parse(read(config_file_json, String))
        @meh delete!(config_json["NotebookApp"]["nbserver_extensions"], "jlstaticserve")
        write(config_file_json, JSON.json(config_json, 4))
    end

    notebook_config_json = joinpath(jupyter_config_dir, "nbconfig", "notebook.json")
    if isfile(notebook_config_json)
        config_json = JSON.parse(read(notebook_config_json, String))
        @meh delete!(config_json["load_extensions"], "webio-jupyter-notebook")
        write(notebook_config_json, JSON.json(config_json, 4))
    end

    # try to uninstall labextension using system Jupyter installation
    @meh run(`jupyter labextension uninstall @webio/jupyter-lab-provider`)

    # try to uninstall labextension using Conda Jupyter installation
    conda_root = joinpath(first(Base.DEPOT_PATH), "conda", "3")
    conda_py = normpath(joinpath(conda_root, @static(
        Sys.iswindows()
        ? joinpath("Scripts", "python.exe")
        : joinpath("bin", "python")
    )))
    @meh run(`$conda_py -m jupyterlab.labextensions uninstall @webio/jupyter-lab-provider`)
catch exc
    # pass
    @warn "Error while updating Jupyter config" exception=exc
end


text:
# TODO: copy some Transducers.jl internals here, instead of importing them.

abstract type AbstractScheduler end
# TODO: Find a better name. It's a bit strange to call it a "scheduler" when it
# can't handle concurrency(?).

struct ThreadsScheduler <: AbstractScheduler end

spawn(@nospecialize(f), ::ThreadsScheduler) = @spawn f()

struct DynamicFunction <: Function
    f::Any
end

(f::DynamicFunction)(x) = f.f(x)

function transduce_dac(sch::AbstractScheduler, ctx, rf, init, xs)
    output = Promise()
    transduce_cps_dac(sch, ctx, rf, init, xs) do acc
        tryput!(output, acc)
    end
    return fetch(output)[]
end

# Using CPS to avoid deadlock (= don't require full concurrency handling for
# the scheduler) and minimize spawn.
function transduce_cps_dac(_return_, sch::AbstractScheduler, ctx, rf, init, xs)
    _return_ = DynamicFunction(_return_)  # avoid type explosion
    if should_abort(ctx)
        _return_(Ok(init))
        return
    end
    if issmall(xs)
        spawn(sch) do
            result = try
                if should_abort(ctx)
                    _return_(Ok(init))
                    return
                end
                acc = _reduce_basecase(rf, init, xs)
                if acc isa Reduced
                    cancel!(ctx)
                end
                Ok(acc)
            catch err
                cancel!(ctx)
                Err(err)
            end
            _return_(result)
        end
    else
        left, right = _halve(xs)
        fg, bg = splitcontext(ctx)
        bridge = Promise()
        transduce_cps_dac(sch, fg, rf, init, left) do a
            b0 = tryput!(bridge, a)
            b0 === nothing && return
            b = something(b0)
            _return_(_combine(ctx, rf, a, b))
        end
        transduce_cps_dac(sch, bg, rf, init, right) do b
            a0 = tryput!(bridge, b)
            a0 === nothing && return
            a = something(a0)
            _return_(_combine(ctx, rf, a, b))
        end
    end
end

function _combine(ctx, rf, a0, b0)
    result = try
        # TODO: merge errors
        if a0 isa Err
            a0
        elseif b0 isa Err
            b0
        else
            Ok(__combine(rf, a0[], b0[]))
        end
    catch err
        Err(err)
    end
    if result isa Err
        cancel!(ctx)
    end
    return result
end

function __combine(rf, a, b)
    a isa Reduced && return a
    b isa Reduced && return combine_right_reduced(rf, a, b)
    return combine(rf, a, b)
end


text:
using SemanticWords
using JSON
using KernelMethods.KMap
using KernelMethods.Scores
using KernelMethods.Kernels
using KernelMethods.Supervised
using SimilaritySearch
using NearNeighborGraph
using TextModel

function read_dense(filename; header=true, delim=' ')
    words = String[]
    vectors = DenseCosine{Float32}[]

    for (lineno, line) in enumerate(readlines(filename))
        if header && lineno == 1
            continue
        end
        arr = split(line, delim)
        m = length(arr)
        vec = Vector{Float32}(m - 1)
        @inbounds for i in 2:m
            num = parse(Float32, arr[i])
            isnan(num) && error("parsed NaN object in $filename, lineno: $lineno -- $line")
            isinf(num) && error("parsed Inf object in $filename, lineno: $lineno -- $line")
            vec[i-1] = num
        end

        try
            v = DenseCosine(vec)
            push!(vectors, v)
        catch
            info("WARNING ignoring vector with zero norm; filename: $filename, lineno: $lineno")
            continue
        end
        push!(words, arr[1])

        if (length(vectors) % 50000) == 1
            info("advance $lineno -- #vectors: $(length(vectors))")
        end
    end

    info("a vocabulary of $(length(words)) has loaded")
    words, vectors
end

function read_sparse(filename; wdelim='\t', delim=' ')
    words = String[]
    vectors = VBOW[]

    f = open(filename)
    for (lineno, line) in enumerate(readlines(filename))
        if lineno == 1
            config = JSON.parse(line)
            continue
        end

        w, line = split(line, wdelim, limit=2)
        push!(words, w)
        arr = split(line, delim)
        m = length(arr)
        vec = Vector{WeightedToken}(m)
        @inbounds for i in 1:m
            a, b = split(arr[i], ':')
            num = parse(Float32, b)
            isnan(num) && error("parsed NaN object in $filename, lineno: $lineno -- $line")
            isinf(num) && error("parsed Inf object in $filename, lineno: $lineno -- $line")
            vec[i] = WeightedToken(parse(Int, a), num)
        end

        push!(vectors, VBOW(vec))

        if (length(vectors) % 50000) == 1
            info("advance $lineno -- #vectors: $(length(vectors))")
        end
    end
    close(f)

    info("a vocabulary of $(length(words)) has loaded")
    words, vectors
end

function create_config(nlist, qlist, skiplist=[])
    config = TextConfig()
    config.del_usr = false
    config.del_punc = false
    config.del_num = false
    config.del_url = true
    config.del_diac = false
    config.nlist = collect(Int, nlist)
    config.qlist = collect(Int, qlist)
    config.skiplist = collect(Int, skiplist)
    config
end

function create_encoder(config, words, centroids, codes)
    H = Dict{String,String}()
    for i in 1:length(words)
        # w = tokenize(words[i], config)
        w = words[i]

        # the encoding is designed to preserve these special chars
        # '!' = Char(33) => missing / unknown token
        # ' ' = Char(32) => space, token separator
        if length(w) > 0
            #sort!(codes[i])
            # H[w[1]] = [Char(c+512) for c in codes[i]] |> join
            H[w] = [Char(c+512) for c in codes[i]] |> join
        end
    end

    function encode(data)::Vector{String}
        if typeof(data) <: Vector
            L = String[]
            for text in data
                for w in tokenize(text::String, config)
                    ww = get(H, w, "")
                    if length(ww) > 0
                        push!(L, w)
                    end
                end
            end

            return L
        else
            return [get(H, w, w) for w in tokenize(data::String, config)]
        end
    end

    encode
end

function encode_corpus(encode::Function, corpus, key_text)
    for tweet in corpus
        tweet["original-text"] = tweet[key_text]
        tweet[key_text] = join(encode(tweet[key_text]), ' ')
    end

    corpus
end

function vectorize_collection(model, col, key_klass, key_text, allow_zero_length)
	X = Vector{VBOW}(length(col))
	y = Vector{String}(length(col))

	i = 0
	for tweet in col
		vec = vectorize(tweet[key_text], model)
		if allow_zero_length || length(vec) > 0
			i += 1
			X[i] = vec
			y[i] = tweet[key_klass]
		end
	end

	if !allow_zero_length
		resize!(X, i)
		resize!(y, i)
	end
	X, y
end

function rocchio_model(model, train, test; key_klass="klass", key_text="text")
    info("vectorizing train")
	X, y = vectorize_collection(model, train, key_klass, key_text, false)

    info("computing centroids")
    ŷ = unique(y)
    X̂ = [centroid!(X[y .== label]) for label in ŷ]
    index = Sequential(X̂, angle_distance)
    info("vectorizing test")
	Xtest, ytest = vectorize_collection(model, test, key_klass, key_text, true)

    info("predicting test")
    ypred = String[]
    for (i, x) in enumerate(Xtest)
        res = search(index, x, KnnResult(1))
        push!(ypred, ŷ[first(res).objID])
    end

    _scores = scores(ytest, ypred)
    _scores
end

function knnclassifier(model, train, test; key_klass="klass", key_text="text")
    info("*** starting kclassifier procedure ***")
    X, y = vectorize_collection(model, train, key_klass, key_text, false)
    info("*** computing kclassifier ***")
    nnc = NearNeighborClassifier(X, y, angle_distance)
    KernelMethods.Supervised.optimize!(nnc, accuracy, folds=3)
    info("vectorizing test")
	Xtest, ytest = vectorize_collection(model, test, key_klass, key_text, true)

    ypred = predict(nnc, Xtest)
    _scores = scores(ytest, ypred)
    _scores
end

function kclassifier(model, train, test; key_klass="klass", key_text="text")
    info("*** starting kclassifier procedure ***")
    X, y = vectorize_collection(model, train, key_klass, key_text, false)
    _X = VBOW[]
    _y = Vector{typeof(y[1])}()
    ŷ = unique(y)
    for label in ŷ
        info("*** computing centroids for label '$label' ***")
        cc = ApproxKDCentroids(31, 3, maxiter=10, tol=0.001, recall=0.9)
        centroids, codes, rlist = codebook(cc, X[y .== label])
        append!(_X, centroids)
        for i in 1:length(centroids)
            push!(_y, label)
        end
    end

    info("*** computing kclassifier ***")
    nnc = NearNeighborClassifier(_X, _y, angle_distance)
    KernelMethods.Supervised.optimize!(nnc, accuracy, folds=3)
    info("vectorizing test")
	Xtest, ytest = vectorize_collection(model, test, key_klass, key_text, true)

    ypred = predict(nnc, Xtest)
    _scores = scores(ytest, ypred)
    _scores
end

function create_codebook()
    numcenters = parse(Int, get(ENV, "numcenters", "10000"))
    k = parse(Int, get(ENV, "k", "7"))
    # kcut = parse(Int, get(ENV, "kcut", "")
    kind = get(ENV, "kind", "random")
    maxiter = parse(Int, get(ENV, "maxiter", "10"))
    tol = parse(Float64, get(ENV, "tol", "0.001"))
    @assert (kind in ("random", "fft", "dnet")) "kind=random|fft|dnet; specifies the algorithm to compute the codebook"
    vecfile = get(ENV, "vectors", "")
    @assert (length(vecfile) > 0) "vectors=file; the embedding file is mandatory"
    modelname = get(ENV, "model", "")
    @assert (length(modelname) > 0) "model=outfile"
    recall = parse(Float64, get(ENV, "recall", "0.9"))

    @show numcenters, k, kind, recall, maxiter, tol, vecfile
    cc = ApproxKDCentroids(numcenters, k, maxiter=maxiter, tol=tol, recall=recall)
    #cc = FFTraversal(numcenters, k)
    if get(ENV, "vformat", "dense") == "dense"
        begin
            words, X = read_dense(vecfile, header=false)
            centroids, codes, rlist = codebook(cc, X)
            # centroids = []
            # codes = [rand(1:numcenters, k) for i in 1:length(X)]
            @save modelname words centroids codes
        end
    else
        begin
            words, X = read_sparse(vecfile)
            centroids, codes, rlist = codebook(cc, X)
            @save modelname words centroids codes
        end
    end
end

function main()
    action = get(ENV, "action", "")
    key_klass = get(ENV, "klass", "klass")
    key_text = get(ENV, "text", "text")
    weighting = get(ENV, "weighting", "tfidf")
    @assert (weighting in ["tfidf", "tf", "idf", "freq", "entropy"]) "unknown weighting scheme $weighting"
    smoothing = parse(Int, get(ENV, "smoothing", "3"))
    nlist = [parse(Int, i) for i in split(get(ENV, "nlist", "1"), ',') if length(i) > 0]
    qlist = [parse(Int, i) for i in split(get(ENV, "qlist", "1,3,5"), ',') if length(i) > 0]
    filter_low = parse(Int, get(ENV, "filter_low", "1"))
    filter_high = parse(Float64, get(ENV, "filter_high", "1.0"))

    if action == "model"
        create_codebook()
    elseif action == "translate"
        modelname=get(ENV, "model", "")
        @assert (length(modelname) > 0) "argument model=modelfile is mandatory"
        @load modelname words centroids codes
        config = create_config([1], [])
        encode = create_encoder(config, words, centroids, codes)
        info("encoding $ARGS")
        key_text = get(ENV, "text", "text")
        @show modelname, key_text, ARGS
        for arg in ARGS
            outname = replace(modelname, ".jld2", "") * "." * basename(arg)
            f = open(outname, "w")
            for tweet in encode_corpus(encode, [JSON.parse(line) for line in readlines(arg)], key_text)
                println(f, JSON.json(tweet))
            end
            close(f)
        end
    elseif action == "evaluate"
        train = [JSON.parse(line) for line in readlines(ARGS[1])]
        test = [JSON.parse(line) for line in readlines(ARGS[2])]

		config = create_config(nlist, qlist)
        X = [item[key_text] for item in train]
        if weighting == "entropy"
            y = [item[key_klass] for item in train]
            le = LabelEncoder(y)
            model = DistModel(config, X, transform.(le, y))
            model = EntModel(model, smoothing)
        else
            vmodel = VectorModel(config)
            vmodel.filter_low = filter_low
            vmodel.filter_high = filter_high
            TextModel.fit!(vmodel, X)
            info("vocabulary size: ", length(vmodel.W))
            if weighting == "tfidf"
                model = TfidfModel(vmodel)
            elseif weighting == "tf"
                model = TfModel(vmodel)
            elseif weighting == "idf"
                model = IdfModel(vmodel)
            elseif weighting == "freq"
                model = FreqModel(vmodel)
            end
        end

        _scores = rocchio_model(model, train, test, key_klass=key_klass, key_text=key_text)
        #_scores = knnclassifier(model, train, test, key_klass=key_klass, key_text=key_text)
		println(JSON.json(_scores))
    elseif action == "semantic-vocabulary"
        train = [JSON.parse(line) for line in readlines(ARGS[1])]

		config = create_config(nlist, qlist)
        X = []
        for item in train
            for text in item[key_text]
                push!(X, text)
            end
        end
        # X = [item[key_text] for item in train]
        if weighting == "entropy"
            y = [item[key_klass] for item in train]
            le = LabelEncoder(y)
            model = DistModel(config, X, transform.(le, y))
            model = EntModel(model, 3)
        else
            vmodel = VectorModel(config)
            vmodel.filter_low = filter_low
            vmodel.filter_high = filter_high
            TextModel.fit!(vmodel, X)
            info("vocabulary size: ", length(vmodel.W))
            if weighting == "tfidf"
                model = TfidfModel(vmodel)
            elseif weighting == "tf"
                model = TfModel(vmodel)
            elseif weighting == "idf"
                model = IdfModel(vmodel)
            elseif weighting == "freq"
                model = FreqModel(vmodel)
            end
        end

        XX = [vectorize(text, model) for text in X]
        # clear X
        tokenmap = id2token(vmodel)
        tX = dtranspose(XX)
        @assert length(model.vmodel.W) == length(tX)
        header = Dict(k => getfield(config, k) for k in fieldnames(config) if k != :normalize)
        header[:normalize] = config.normalize != identity
        println(JSON.json(header))
        for (keyid, tokens) in tX
            #info("word $keyid - $(tokenmap[keyid]): ", [(a.id, a.weight) for a in tokens])
            if length(tokens) >= filter_low
                println(tokenmap[keyid], "\t", join([string(a.id |> Int, ":",a.weight) for a in tokens], ' '))
            end
        end

	elseif action == "index"
		filename = ARGS[1]
		indexname = replace(filename, ".json", "") * ".index.jld2"
		if isfile(indexname)
			info("loading index from $indexname")
			@load indexname train model nns
		else
			train = [JSON.parse(line) for line in readlines(filename)]
			config = create_config(nlist, qlist)
			model = VectorModel(config)
			model.filter_low = filter_low
			model.filter_high = filter_high
			TextModel.fit!(model, train, get_text=x -> x[key_text])
			X, y = vectorize_collection(model, train, key_klass, key_text, false)
			N = LogSatNeighborhood()
			nns = LocalSearchIndex(VBOW, cosine_distance, search=BeamSearch(), neighborhood=N)
			NearNeighborGraph.fit!(nns, X)
			NearNeighborGraph.optimize!(nns, 0.9)
			@save indexname train model nns
		end

		for p in search(nns, vectorize("a veces me dan ganas pero luego se me pasa", model), KnnResult(5))
			@show p, train[p.objID][key_klass]
		end
    else
        exit(1)
    end
end

if !isinteractive()
	main()
end


text:
function solve_master(node, data, tour_pool, branch_pool, relax)

    ### Constants ###
    tour_number = length(node["tours"])
    tour_costs = [tour_pool[i]["cost"] for i in node["tours"]]
    vehicle_number = data[1,:nb_vehicles]

    ### Model definiton ###
    model = Model(with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    ### Variables ###
    if relax
        @variable(model, alpha[1:tour_number] >= 0)
    else
        @variable(model, alpha[1:tour_number], Bin)
    end

    ### Objective ###
    @objective(model, Min, sum(alpha.*tour_costs))

    ### Constraints ###
    # Add the branch constraints depending on node type
    sum_branchs = []
    for i in node["branch"]
        sum_branch = 0
        for (j,k) in enumerate(node["tours"])
            if branch_pool[i][1] == tour_pool[k]["depot"]
                sum_branch += alpha[j]*(tour_pool[k]["X"][branch_pool[i][2]][branch_pool[i][3]])
            end
        end
        push!(sum_branchs, sum_branch)
    end
    if node["type"]=="odd"
        @constraint(model, branch_lower[i=1:length(node["branch"])],
            -sum_branchs[i]<=-0)
        @constraint(model, branch_upper[i=1:length(node["branch"])],
            sum_branchs[i]<=0)
    elseif node["type"]=="even"
        @constraint(model, branch_lower[i=1:length(node["branch"])],
            -sum_branchs[i]<=-1)
        @constraint(model, branch_upper[i=1:length(node["branch"])],
            sum_branchs[i]<=vehicle_number)
    end

    # Respect the number of vehicles
    @constraint(model, respect_number_vehicle,
        sum(alpha)<=vehicle_number)

    # Service on required arcs (at least one tour service each arc)
    @constraint(model, service_arcs[arc=1:data[1,:req_arcs]],
        sum(alpha.*[tour_pool[i]["X"]["req_arcs"][arc] for i in node["tours"]])>=1
        )

    # Service on required edges (at least one tour service each edge)
    @constraint(model, service_edges[edge=1:data[1,:req_edges]],
        sum(alpha.*[tour_pool[i]["X"]["req_edges"][edge] + tour_pool[i]["X"]["req_edges_swap"][edge] for i in node["tours"]])>=1
        )

    ### Solver master problem ###
    optimize!(model)

    ### Query dual values ###
    dual_values = Dict{String, Any}()
    # Dual value associated to vehicule number constraint
    try
        dual_values["vehicle"] = dual(respect_number_vehicle)
    catch
        dual_values["vehicle"] = 0
    end
    # Dual values associate to the arc servicing constraint
    dual_values_req_arcs = []
    for arc in 1:data[1,:req_arcs]
        try
            push!(dual_values_req_arcs, dual(service_arcs[arc]))
        catch
            push!(dual_values_req_arcs, 0)
        end
    end
    dual_values["req_arcs"] = dual_values_req_arcs
    # Dual values associate to the edge servicing constraint
    dual_values_req_edges = []
    for edge in 1:data[1,:req_edges]
        try
            push!(dual_values_req_edges, dual(service_edges[edge]))
        catch
            push!(dual_values_req_edges, 0)
        end
    end
    dual_values["req_edges"] = dual_values_req_edges
    # Dual values of branch upper
    dual_values_branch_upper = []
    for i in 1:length(node["branch"])
        try
            push!(dual_values_branch_upper, dual(branch_upper[i]))
        catch
            push!(dual_values_branch_upper, 0)
        end
    end
    dual_values["branch_upper"] = dual_values_branch_upper
    # Dual values of branch lower
    dual_values_branch_lower = []
    for i in 1:length(node["branch"])
        try
            push!(dual_values_branch_lower, dual(branch_lower[i]))
        catch
            push!(dual_values_branch_lower, 0)
        end
    end
    dual_values["branch_lower"] = dual_values_branch_lower

    node["alpha"] = value.(alpha)

    return(dual_values)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
module FeynmanKleinert

# Bring out the major leagues... https://www.youtube.com/watch?v=_E6DDktoPhg
using Optim

using QuadGK # For numerical integration

export K,Va2,Wtilde,verboseWtilde,W # export all the things.

"Integrand of (4) in Feynman and Kleinert"
K(xp,x,a2,V)=1/sqrt(2*π*a2)*exp(-(x-xp)^2/(2*a2))*V(xp)
"(4) in Feynman and Kleinert"
Va2(x,a2,V)=quadgk(xp->K(xp,x,a2,V), -Inf, +Inf) #, reltol=0, abstol=1e-5)

# (5) naively as written, in Feynman and Kleinert
#Wtilde(x0,a2,Ω,g,β)=(1/β) * log( sinh(β*Ω/2)/(β*Ω/2) ) - (Ω^2/2) * a2 + Va2(x0,a2,g)[1]
# Large β limit; by taking sinh(x)=e^x-e^-x , and then dropping e^-x
#Wtilde(x0,a2,Ω,g,β)=(1/β) * (β*Ω/2 - log(β*Ω/2)) - (Ω^2/2) * a2 + Va2(x0,a2,g)[1]

# Branch depending on β; for accuracy + lack of infinities
"""
    Wtilde(x0,a2,Ω,V,β)

(5) in Feynman and Kleinert. Switches depending on D=β*Ω/2 to use either an
exponential approximation for sinh(x) (D>7.5) or an explicit evaluation and
then take logarithms. Without doing this you get collapse to infinity for large
β.
"""
function Wtilde(x0,a2,Ω,V,β)
    D=β*Ω/2
    if D>7.5
        F=(1/β) * (D - log(2*D)) # Large limit;  by taking sinh(x)=e^x-e^-x , and then dropping e^-x
    else
        F=(1/β)*( log(sinh(D)/D) ) # explicit evaluation, naively as written in (5)
    end
    V0=Va2(x0,a2,V)[1]
    mid=(Ω^2/2) * a2
    return F-mid+V0
end

"""
   verboseWtilde(x0,a2,Ω,g,β)

Expanded version of Wtilde, reporting the decomposition - this was used to explain what was going wrong
"""
function verboseWtilde(x0,a2,Ω,g,β)
    D=β*Ω/2
    println("D=β*Ω/2 = $D")
    F= (1/β)*( log(sinh(D)) - log(D) )
    #F= (1/β)*( (log(2D)+D^2/6+D^4/180) - log(D) ) # generalised Puiseux series
    # Large D limit, through sinh(D) --> e^D [through away e^-D part]
    F = (1/β) * (D - log(2*D))
    V=Va2(x0,a2,g)[1]
    mid=(Ω^2/2) * a2
    println("Wtilde(x0=$x0,a2=$a2,Ω=$Ω,g=$g,β=$β) \n\t= F - mid + V")
    println("\t= $F - $mid + $V = ",F-mid+V)
    return F-mid+V
end

"""
    W(x0,g,β)

    Implementing optimisation / minimisation described in Feynman-Kleinert (6)
    Calculates optimal W for x0, with given g and Beta
    Use's (7) (a2 in terms of omega) to fix numerical noise from trying to do
    it directly...
"""
function W(x0,g,β)
    #println("W(x0=$x0,g=$g,β=$β)")

    a2(Ω)=1/(β*Ω^2)*( (β*Ω)/2 * coth(β*Ω/2)-1 ) #(7)
    myf(x)=Wtilde(x0,a2(x[1]),x[1],g,β)

    # Define bounds on solution
    lower=[0.0]
    upper=[5.0]
    initial=[1.0]

    res=optimize(OnceDifferentiable(myf, initial; autodiff = :forward), initial, lower, upper, Fminbox();
    optimizer=GradientDescent)

    minimum=Optim.minimum(res)
    Ω=Optim.minimizer(res)[1]
    mya2=a2(Ω)[1]
    #println("Optimised for W(x0=$x0,g=$g,β=$β). \n\tΩ=$Ω, a2=$mya2, minimum=$minimum.")
    #show(res)
    return minimum
end

end # module


text:
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "types.jl" begin
    include("types.jl")
end

@testset "query.jl" begin
    include("query.jl")
end

@testset "testing.jl" begin
    include("testing.jl")
end

@testset "jupyter.jl" begin
    include("jupyter.jl")
end

@testset "mockserver.jl" begin
    include("mockserver.jl")
end


text:
module RocketDiscontinueOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: discontinue()" begin

    println("Testing: operator discontinue()")

    run_proxyshowcheck("Discontinue", discontinue())

    @testset begin
        values    = [ ]
        errors    = [ ]
        completes = [ ]

        subject = BehaviorSubject(0)

        subscription1 = subscribe!(subject, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (d) -> push!(errors, d),
            on_complete = ()  -> push!(completes, 0)
        ))

        @test values    == [ 0 ]
        @test errors    == [ ]
        @test completes == [ ]

        subscription2 = subscribe!(subject |> map(Int, d -> d + 1) |> discontinue(), lambda(
            on_next     = (d) -> next!(subject, d),
            on_error    = (e) -> error!(subject, e),
            on_complete = ()  -> complete!(subject)
        ))

        @test values    == [ 0, 1 ]
        @test errors    == [ ]
        @test completes == [ ]

        error!(subject, "err")

        @test values    == [ 0, 1 ]
        @test errors    == [ "err" ]
        @test completes == [ ]
    end

    @testset begin
        values    = [ ]
        errors    = [ ]
        completes = [ ]

        subject = BehaviorSubject(0)

        subscription1 = subscribe!(subject, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (d) -> push!(errors, d),
            on_complete = ()  -> push!(completes, 0)
        ))

        @test values    == [ 0 ]
        @test errors    == [ ]
        @test completes == [ ]

        subscription2 = subscribe!(subject |> map(Int, d -> d + 1) |> discontinue(), lambda(
            on_next     = (d) -> next!(subject, d),
            on_error    = (e) -> error!(subject, e),
            on_complete = ()  -> complete!(subject)
        ))

        @test values    == [ 0, 1 ]
        @test errors    == [ ]
        @test completes == [ ]

        complete!(subject)

        @test values    == [ 0, 1 ]
        @test errors    == [ ]
        @test completes == [ 0 ]
    end

end

end


text:
using ApproxFun, SingularIntegralEquations, Base.Test

# Scalar case

G=Fun(z->2+cos(z+1/z),Circle()) # the symbol of the Toeplitz operator
T=ToeplitzOperator(G)

C  = Cauchy(-1)
V  = (I+(1-G)*C)\(G-1)

Φmi = 1+C*V
Φp = V+Φmi

L  = ToeplitzOperator(1/Φmi)
U  = ToeplitzOperator(Φp)

@test norm((T-U*L)[1:10,1:10]) < 10eps()  # check the accuracy

# Matrix case

G=Fun(z->[-1 -3; -3 -1]/z +
         [ 2  2;  1 -3] +
         [ 2 -1;  1  2]*z,Circle())

C  = Cauchy(-1)

A=(I+(I-G)*C)

F̃ = (G-I)[:,1]
F=Fun((G-I)[:,1])

@test Fun(F̃,space(F)) == F

V1  = A\F
Ṽ1 = A\F̃

A1=ApproxFun.choosespaces(A,(G-I)[:,1])
A2=ApproxFun.choosespaces(A,Fun((G-I)[:,1]))

@test A1\Fun((G-I)[:,1])  == V1
@test A1\(G-I)[:,1]  == V1

QR=qrfact(A1)

@test QR\Fun((G-I)[:,1]) == V1
@test QR\(G-I)[:,1] == V1

@test norm((V1-Ṽ1).coefficients) == 0
@test norm((A*V1-F).coefficients) < 100eps()

@test norm((F-Fun((G-I)[:,1])).coefficients) == 0

V  = (I+(I-G)*C)\(G-I)

@test norm((V1-Fun(V[:,1])).coefficients) == 0

Φmi = I+C*V
Φp = V+Φmi

T=ToeplitzOperator(G)

L  = ToeplitzOperator(inv(Φmi))
U  = ToeplitzOperator(Φp)

@test norm((T-U*L)[1:10,1:10]) < 100eps()  # check the accuracy


text:
##################################################################
#
#   Maximal cliques of undirected graph
#   Derived from Graphs.jl: https://github.com/julialang/Graphs.jl
#
##################################################################

@testset "Cliques" begin
    function setofsets(array_of_arrays)
        Set(map(Set, array_of_arrays))
    end

    function test_cliques(graph, expected)
        # Make test results insensitive to ordering
        setofsets(@inferred(maximal_cliques(graph))) == setofsets(expected)
    end

    @testset "simple maximal clique cases" begin
        gx = SimpleGraph(3)
        add_edge!(gx, 1, 2)
        @testset "$g simple single maximal clique" for g in testgraphs(gx)
            @test test_cliques(g, Array[[1, 2], [3]])
        end
        add_edge!(gx, 2, 3)
        @testset "$g" for g in testgraphs(gx)
            @test test_cliques(g, Array[[1, 2], [2, 3]])
        end
    end
    @testset "pivotdonenbrs defined" begin
        h = SimpleGraph(6)
        add_edge!(h, 1, 2)
        add_edge!(h, 1, 3)
        add_edge!(h, 1, 4)
        add_edge!(h, 2, 5)
        add_edge!(h, 2, 6)
        add_edge!(h, 3, 4)
        add_edge!(h, 3, 6)
        add_edge!(h, 5, 6)

        @testset "$g" for g in testgraphs(h)
            @test !isempty(@inferred(maximal_cliques(g)))
        end
    end

    # test for extra cliques bug
    @testset "does not find extra cliques" begin
        h = SimpleGraph(7)
        add_edge!(h, 1, 3)
        add_edge!(h, 2, 6)
        add_edge!(h, 3, 5)
        add_edge!(h, 3, 6)
        add_edge!(h, 4, 5)
        add_edge!(h, 4, 7)
        add_edge!(h, 5, 7)
        @testset "$g" for g in testgraphs(h)
            @test test_cliques(h, Array[[7, 4, 5], [2, 6], [3, 5], [3, 6], [3, 1]])
        end
    end
end