text:
using SaddleSearch.TestSets: hessprecond
using Isaac
using SaddleSearch: ODE12r, odesolve, IterationLog

verbose = 1

# MullerPotential with good IC
V = MullerPotential()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
dimer = SuperlinearDimer(maximum_translation=0.2, max_num_rot=1, maxnumdE=100, verbose=verbose)
x1, v, res = run!(dimer, E, dE, x0, v0)
@show numE(res)[end], numdE(res)[end]

# # MullerPotential with good IC
# V = MullerPotential()
# x0, v0 = ic_dimer(V, :near)
# E, dE = objective(V)
# dimer = DimerMethod(verbose=verbose, a0_trans = 0.001, a0_rot = 0.001)
# x2, v, res = run!(dimer, E, dE, x0, v0)
# @show numE(res)[end], numdE(res)[end]

# @show x1, x2

# Lennard-Jones Cluster
V = LJcluster()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
dimer = SuperlinearDimer(maximum_translation=0.2, max_num_rot=1, maxnumdE=500, verbose=verbose)
x, v, res = run!(dimer, E, dE, x0, v0)
@show numE(res)[end], numdE(res)[end]

V = LJVacancy2D()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
dimer = SuperlinearDimer(maximum_translation=0.2, max_num_rot=1, maxnumdE=500, verbose=verbose)
x, v, res = run!(dimer, E, dE, x0, v0)
@show numE(res)[end], numdE(res)[end]

# # MullerPotential with good IC
# V = MullerPotential()
# x0, v0 = ic_dimer(V, :near)
# E, dE = objective(V)
#
# function Fdimer(t, xv)
#    n = length(xv) ÷ 2
#    x, v = xv[1:n], xv[n+1:end]
#    v /= norm(v)
#    fx = - (eye(n) - 2 * v * v') * dE(x)
#    fv = - (eye(n) - v * v') * ((dE(x + 1e-3 * v) - dE(x))/1e-3)
#    return [fx; fv]
# end
#
# ode = ODE12r(atol = 1e-3, rtol = 1e-3)
# odesolve(ode, Fdimer, [x0;v0], 2,
#          IterationLog(:numE => Int, :numdE => Int, :res => Float64),
#          verbose = 3)

# MullerPotential with good IC
V = MullerPotential()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
sol, it_hist, ierr, x_hist = nsoli(x0, dE)
@show length(x_hist)
@show it_hist[:, 2][end]

V = LJcluster()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
sol, it_hist, ierr, x_hist = nsoli(x0, dE)
@show length(x_hist)
@show it_hist[:, 2][end]

V = LJVacancy2D()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
sol, it_hist, ierr, x_hist = nsoli(x0, dE)
@show length(x_hist)
@show it_hist[:, 2][end]


text:
using BinaryBuilder

# Collection of sources required to build Gettext
name = "Gettext"
version = v"0.20.1"
sources = [
    ArchiveSource("https://ftp.gnu.org/pub/gnu/gettext/gettext-$(version).tar.xz",
                  "53f02fbbec9e798b0faaf7c73272f83608e835c6288dd58be6c9bb54624a3800"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gettext-*/
export CFLAGS="-O2"
./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgettextlib", "libgettextlib-0-20"], :libgettext)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Libiconv_jll"),
    Dependency("XML2_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl


text:
using Flux: param

struct WordVectors{A<:AbstractArray{<:AbstractFloat, 2}, B<:AbstractArray{<:AbstractString, 1}, C<:AbstractDict{<:AbstractString, <:Integer}}
    vectors::A # the vectors computed from word2vec
    vocab::B # id to vocab
    vocab_hash::C # vocabulary to id
end

function WordVectors(wordvector_dims::Integer, vocab::AbstractArray, vocab_hash)
    return WordVectors(param(randn(wordvector_dims, length(vocab))), vocab, vocab_hash)
end

(layer::WordVectors)(x::Integer, y::Integer) = layer.vectors[x, :], y

(layer::WordVectors)(x::AbstractArray{<:Real, 1}, y::Integer) = layer.vectors*x, y

(layer::WordVectors)(x::AbstractArray{<:Real, 2}, y::AbstractArray{<:Real, 1}) = layer.vectors*x, y

(layer::WordVectors)(data::Tuple) = layer(data...)

import Base: show, size

function Base.show(io::IO, wv::WordVectors)
    len_vecs, num_words = size(wv.vectors)
    print(io, "WordVectors $(num_words) words, $(len_vecs)-element vectors")
end

"""
    size(wv)

Return the word vector length and the number of words as a tuple.
"""
size(wv::WordVectors) = size(wv.vectors)

"""
    vocabulary(wv)

Return the vocabulary as a vector of words of the WordVectors `wv`.
"""
vocabulary(wv::WordVectors) = wv.vocab

"""
    in_vocabulary(wv, word)

Return `true` if `word` is part of the vocabulary of the WordVector `wv` and
`false` otherwise.
"""
in_vocabulary(wv::WordVectors, word) = word in wv.vocab

"""
    index(wv, word)

Return the index of `word` from the WordVectors `wv`.
"""
index(wv::WordVectors, word) = wv.vocab_hash[word]

"""
    get_vector(wv, word)

Return the vector representation of `word` from the WordVectors `wv`.
"""
get_vector(wv::WordVectors, word) =
      (idx = wv.vocab_hash[word]; wv.vectors[:,idx])

"""
    cosine(wv, word, n=10)

Return the position of `n` (by default `n = 10`) neighbors of `word` and their
cosine similarities.
"""
function cosine(wv::WordVectors, word, n=10)
    metrics = wv.vectors'*get_vector(wv, word)
    topn_positions = sortperm(metrics[:], rev = true)[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    similarity(wv, word1, word2)

Return the cosine similarity value between two words `word1` and `word2`.
"""
function similarity(wv::WordVectors, word1, word2)
    return get_vector(wv, word1)'*get_vector(wv, word2)
end

"""
    cosine_similar_words(wv, word, n=10)

Return the top `n` (by default `n = 10`) most similar words to `word`
from the WordVectors `wv`.
"""
function cosine_similar_words(wv::WordVectors, word, n=10)
    indx, metr = cosine(wv, word, n)
    return vocabulary(wv)[indx]
end

"""
    analogy(wv, pos, neg, n=5)

Compute the analogy similarity between two lists of words. The positions
and the similarity values of the top `n` similar words will be returned.
For example,
`king - man + woman = queen` will be
`pos=[\"king\", \"woman\"], neg=[\"man\"]`.
"""
function analogy(wv::WordVectors, pos::AbstractArray, neg::AbstractArray, n= 5)
    m, n_vocab = size(wv)
    n_pos = length(pos)
    n_neg = length(neg)
    anal_vecs = Array{AbstractFloat}(m, n_pos + n_neg)

    for (i, word) in enumerate(pos)
        anal_vecs[:,i] = get_vector(wv, word)
    end
    for (i, word) in enumerate(neg)
        anal_vecs[:,i+n_pos] = -get_vector(wv, word)
    end
    mean_vec = mean(anal_vecs, 2)
    metrics = wv.vectors'*mean_vec
    top_positions = sortperm(metrics[:], rev = true)[1:n+n_pos+n_neg]
    for word in [pos;neg]
        idx = index(wv, word)
        loc = findfirst(top_positions, idx)
        if loc != 0
            splice!(top_positions, loc)
        end
    end
    topn_positions = top_positions[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    analogy_words(wv, pos, neg, n=5)

Return the top `n` words computed by analogy similarity between
positive words `pos` and negaive words `neg`. from the WordVectors `wv`.
"""
function analogy_words(wv::WordVectors, pos, neg, n=5)
    indx, metr = analogy(wv, pos, neg, n)
    return vocabulary(wv)[indx]
end


text:
# the code was tested on v0.4.2
# update HetaSimulator
# with ] add https://github.com/hetalang/HetaSimulator.jl

using HetaSimulator, Plots
using Distributed # to use parallel simulations

platform = load_platform(".")
model = platform.models[:nameless]

### Single scenario simulations

mcscn1 = Scenario(
  model;
  tspan = (0., 200.),
  parameters = [:k1=>0.01],
  saveat = [50., 80., 150.]
)

mcscn2 = Scenario(
  model;
  tspan = (0., 200.),
  parameters = [:k1=>0.02],
  saveat = [50., 100., 200.]
)

mcsim1 = mc(
  mcscn1,
  [:k1=>Uniform(1e-3,1e-2), :k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)],
  100
)

plotd = plot(mcsim1, vars=[:b])
# savefig(plotd, "fig.png")

mc_df1 = DataFrame(mcsim1, vars=[:a, :b])

### Multiple scenarios simulations

mcsim2 = mc(
  [:mc1=>mcscn1,:mc2=>mcscn2],
  [:k1=>0.01, :k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)],
  100
)
plotd = plot(mcsim2)
# savefig(plotd,"file.png")

mc_df2 = DataFrame(mcsim2)

### Monte-Carlo for whole platform

scn_csv = read_scenarios("./scenarios.csv")
add_scenarios!(platform, scn_csv)

mcplat = mc(
  platform,
  [:k1=>0.01, :k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)],
  100
)
plotd = plot(mcplat)
# savefig(plotd,"file.png")

### Using pre-generated parameter set

#=
mcvecs = DataFrame(
  k1=0.01,
  k2=rand(Normal(1e-3,1e-4), 50),
  k3=rand(Uniform(1e-4,1e-2), 50)
)
=#

mcvecs = read_mcvecs("./params.csv")

mcv1 = mc(mcscn1, mcvecs)
plotd = plot(mcv1)
# savefig(plotd,"file.png")

# parallel simulations

addprocs(2)
@everywhere using HetaSimulator

mcv2 = mc(
  mcscn1,
  mcvecs;
  parallel_type=EnsembleDistributed()
)

### Monte-Carlo statistics
# see https://diffeq.sciml.ai/stable/features/ensemble/#Summary-Statistics

# mean
timestep_mean(mcv1,2)
timepoint_mean(mcv1,80)

# median
timestep_median(mcv1,2)
timepoint_median(mcv1,80)

# meanvar
timestep_meanvar(mcv1,2)
timepoint_meanvar(mcv1,80)

# meancov
timestep_meancov(mcv1,2,3)
timepoint_meancov(mcv1,80.,150.)

# meancor
timestep_meancor(mcv1,2,3)
timepoint_meancor(mcv1,80.,150.)

# !!!quantile
timestep_quantile(mcv1,0.95,2)
timepoint_quantile(mcv1,0.95,80.)

# full time steps statistics
timeseries_steps_mean(mcv1) # Computes the mean at each time step
timeseries_steps_median(mcv1) # Computes the median at each time step
timeseries_steps_quantile(mcv1,0.95) # Computes the quantile q at each time step
timeseries_steps_meanvar(mcv1) # Computes the mean and variance at each time step
timeseries_steps_meancov(mcv1) # Computes the covariance matrix and means at each time step
timeseries_steps_meancor(mcv1) # Computes the correlation matrix and means at each time step

# Ensemble Summary
ens = EnsembleSummary(mcsim1;quantiles=[0.05,0.95])
plot(ens)
# plotd = plot(ens)
# savefig(plotd,"file.png")


text:
# ---
# title: 795. Number of Subarrays with Bounded Maximum
# id: problem795
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/number-of-subarrays-with-bounded-maximum/description/>
# hidden: true
# ---
#
# We are given an array `A` of positive integers, and two positive integers `L`
# and `R` (`L <= R`).
#
# Return the number of (contiguous, non-empty) subarrays such that the value of
# the maximum array element in that subarray is at least `L` and at most `R`.
#
#
#
#     **Example :**
#     Input:
#     A = [2, 1, 4, 3]
#     L = 2
#     R = 3
#     Output: 3
#     Explanation: There are three subarrays that meet the requirements: [2], [2, 1], [3].
#
#
# **Note:**
#
#   * L, R  and `A[i]` will be an integer in the range `[0, 10^9]`.
#   * The length of `A` will be in the range of `[1, 50000]`.
#
#
## @lc code=start
using LeetCode

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


text:
module SpeedTests

using Test, Random
using DataFrames
using DataFramesMeta
using Devectorize

Random.seed!(1)
const n = 5_000_000
a = rand(n)
b = rand(n)
df = DataFrame(a = da, b = db)
df2 = DataFrame(Any[a, b])
names!(df2, [:a, :b])

function dot1(a::Vector, b::Vector)
    x = 0.0
    for i in 1:length(a)
        x += a[i] * b[i]
    end
    return x
end

function dot3(df::DataFrame)
    da, db = df[:a], df[:b]
    T = eltype(da)
    x = 0.0
    for i in 1:length(da)
        x += da[i]::T * db[i]::T
    end
    return x
end

function dot4(df::DataFrame)
    da, db = df[:a], df[:b]
    return dot2(da, db)
end

function dot8(a::Vector, b::Vector)
    x = 0.0
    for i in 1:length(a)
        if !(isnan(a[i]) || isnan(b[i]))
            x += a[i] * b[i]
        end
    end
    return x
end

function dot9(df::DataFrame)
    @with df begin
        x = 0.0
        for i in 1:length(:a)
            x += values(:a)[i] * values(:b)[i]
        end
        x
    end
end

t1 = @elapsed dot1(a, b)
t3 = @elapsed dot3(df)
t4 = @elapsed dot4(df)
t8 = @elapsed dot8(a, b)
t9 = @elapsed dot9(df)

end


text:
global const marmousi_folder=joinpath(@__DIR__, "marmousi2")

"""
```julia
mod=Medium(attrib)
```
As of now, only seismic models are predefined in this package. Choose `attrib::Symbol`

* `=:acou_homo1` : a test homogeneous acoustic model
* `=:acou_homo2` : a test homogeneous acoustic model, but with coarser spatial sampling (faster testing)
* `=:marmousi2` : marmousi2 model with lower resolution; useful for surface seismic experiments
* `=:pizza` : `:acous_homo2` with some perturbations
"""
function Medium(attrib::Symbol, δ::Real=0.0; verbose=false)
	bfrac=0.1;
	δ=Float64(δ)
	if((attrib == :acou_homo1))
		vp0 = [1500., 3500.] # bounds for vp
		rho0 = [1500., 3500.] # density bounds
		mgrid = repeat([range(-1000.0,stop=1000.0,length=201)],2)
		nz,nx=length.(mgrid)
		model=Medium(mgrid,[:vp,:rho])
		update!(model,[:vp,:rho],[vp0,rho0])
		fill!(model)
	elseif((attrib==:pizza))
		vp0 = [1700., 2300.] # bounds for vp
		rho0 = [1700., 2300.] # density bounds
		mgrid = repeat([range(-1000.0,stop=1000.0,length=51)],2)
		nz,nx=length.(mgrid)
		model=Medium(mgrid,[:vp,:rho])
		update!(model,[:vp,:rho],[vp0,rho0])
		fill!(model)

		# add some noise to starting model
		update!(model, [:vp,:rho], randn_perc=0.5)

		# add perturbations
		for ellip_loc in [[500.,0.], [0.,500.], [-500.,0.], [0.,-500.]]
			update!(model, [:vp,:rho], ellip_rad=50., ellip_loc=ellip_loc,
				ellip_pert=100.)
		end

	elseif((attrib == :acou_homo2))
		vp0 = [1700., 2300.] # bounds for vp
		rho0 = [1700., 2300.] # density bounds
		mgrid = repeat([range(-1000.0,stop=1000.0,length=51)],2)
		nz,nx=length.(mgrid)
		model=Medium(mgrid,[:vp,:rho])
		update!(model,[:vp,:rho],[vp0,rho0])
		fill!(model)

	elseif(attrib == :marmousi2)
		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii_0.1.su"))
		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii_0.1.su"))
		rho,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii_0.1.su"))
		vp .*= 1000.; vs .*= 1000.; #rho .*=1000
		mgrid=[range(0.,stop=3500.,length=size(vp,1)),range(0., stop=17000., length=size(vp,2))]
		model=Medium(mgrid,[:vp,:rho,:vs])
		copyto!(model[:vp],vp)
		copyto!(model[:rho],rho)
		copyto!(model[:vs],vs)
		update!(model,bfrac);
#		vs0=Models.bounds(vs,bfrac);
#		rho0=Models.bounds(rho, bfrac);
#		update!(model,[:vp,:rho,:vs],[vp0,rho0,vs0])
#=
	elseif(attrib == :seismic_marmousi2_high_res)
		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii.su"))
		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii.su"))
		rho,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii.su"))
		vp .*= 1000.; vs .*= 1000.; #rho .*=1000
		vp0=Models.bounds(vp,bfrac);
		vs0=Models.bounds(vs,bfrac);
		rho0=Models.bounds(rho, bfrac);
		mgrid=[range(0.,stop=3500.,length=size(vp,1)),range(0., stop=17000., length=size(vp,2))]
		model=Medium(mgrid,[:vp,:rho,:vs])
		update!(model,[:vp,:rho,:vs],[vp0,rho0,vs0])
		copyto!(model[:vp],vp)
		copyto!(model[:rho],rho)
		copyto!(model[:vs],vs)

	elseif(attrib == :seismic_marmousi2_xwell)
		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
				     zmin=1000., zmax=2000., xmin=8500., xmax=9500.,)
		update!(model, bfrac) # adjuts bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_surf)
		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=6000., xmax=12000.,)
		update!(model, bfrac) # adjust bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_downhole)
		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=9025., xmax=9125., zmin=1400., zmax=1600.,)
		update!(model, bfrac) # adjust bounds just inside the bounds
	elseif(attrib == :seismic_marmousi2_rvsp)
		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
				     xmin=8000., xmax=10000., zmax=1700.,zmin=500.)
		update!(model, bfrac) # adjust bounds just inside the bounds
		=#

	else
		error("invalid attrib")
	end
	if(δ==0.0)
		verbose && Models.print(model,string(attrib))
		return model
	elseif(δ > 0.0)
		mgrid_out=broadcast(x->range(x[1],stop=x[end],step=δ),model.mgrid)
		model_out=Medium_zeros(mgrid_out)
		update!(model_out,model)
		Models.interp_spray!(model, model_out, :interp, :B1)
		verbose && Models.print(model_out,string(attrib))
		return model_out
	else
		error("invalid δ")
	end

end


text:
const NUM_HANDLERS = 4
const QSIZE = 32

IN_QUEUE = Channel{ParamsAndSettings}(QSIZE)

"""
Wait to pull a job off the job queue and sent it
to the calculator. FIXME should also just generate text output.
"""
function calc_one()
	while true
		@debug "calc_one entered"
		params = take!( IN_QUEUE )
		@debug "params taken from IN_QUEUE; got params uuid=$(params.settings.uuid)"
		results = do_run_a( params.cache_key, params.sys, params.settings )
		@debug "model run OK; putting results into STASHED_RESULTS"
		res_text = results_to_html( results.uuid, BASE_RESULTS, results )
		@debug "calc_one; stashing to uuid=$(results.uuid)"
		STASHED_RESULTS[ results.uuid ] = res_text
		@debug "calc_one; caching to to cache_key=$(results.cache_key)"
		CACHED_RESULTS[ results.cache_key ] = res_text
	end
end

function submit_job( cache_key, sys :: TaxBenefitSystem, settings :: Settings )
    uuid = UUIDs.uuid4()
	@debug "submit_job entered uuid=$uuid"
	settings.uuid = uuid
    put!( IN_QUEUE, ParamsAndSettings(uuid, cache_key, sys, settings ))
	@debug "submit exiting queue is now $IN_QUEUE"
    return uuid
end

#
# Set up job queues
#
for i in 1:NUM_HANDLERS # start n tasks to process requests in parallel
    errormonitor(@async calc_one())
 end


text:
"""
    get_quotes(ticker)
    get_quotes(ticker::Array)

Get quote of a symbol:
https://developer.tdameritrade.com/quotes/apis/get/marketdata/{symbol}/quotes

`ticker` can be a single `String` or an `Array` or multiple `String`s

# Examples
```julia
get_quotes("GE")
get_quotes(["BA","GE"])
```
"""
function get_quotes(ticker::AbstractString)
    uri = construct_api("marketdata/$ticker/quotes")
    head = ["Authorization" => "Bearer "*AUTH_KEY.ACCESS_TOKEN]
    # Authorization is a header
    return @pipe HTTP.get(uri, head).body |>
                 JSON3.read |> _[Symbol(ticker)]
end

function get_quotes(tickers::AbstractArray{T,1}) where T<:AbstractString
    kwargs = (; symbol = join(tickers, ",")) # NamedTuple with 1 element
    uri = construct_api("marketdata/quotes", kwargs)
    @show uri
    head = ["Authorization" => "Bearer "*AUTH_KEY.ACCESS_TOKEN]
    # Authorization is a header
    return @pipe HTTP.get(uri, head).body |>
                 JSON3.read
end


text:
euclidex = """
    while a != b
        if a > b
            a = a - b
            p = p - q
            r = r - s
        else
            q = q - p
            b = b - a
            s = s - r
        end
    end
"""

fermat = """
    while r != 0
        if r > 0
            r = r - v
            v = v + 2
        else
            r = r + u
            u = u + 2
        end
    end
"""

wensley = """
    while d>= E
        if P < a+b
            b = b/2
            d = d/2
        else
            a = a+b
            y = y+d/2
            b = b/2
            d = d/2
        end
    end
"""

lcm = """
    while x != y
        if x > y
            x = x - y
            v = v + u
        else
            y = y - x
            u = u + v
        end
    end
"""

knuth = """
    while (s >= d) && (r != 0)
        if 2*r-rp+q < 0
        t  = r
        r  = 2*r-rp+q+d+2
        rp = t
        q  = q+4
        d  = d+2
        elseif (2*r-rp+q >= 0) && (2*r-rp+q < d+2)
        t  = r
        r  = 2*r-rp+q
        rp = t
        d  = d+2
        elseif (2*r-rp+q >= 0) && (2*r-rp+q >= d+2) && (2*r-rp+q < 2*d+4)
        t  = r
        r  = 2*r-rp+q-d-2
        rp = t
        q  = q-4
        d  = d+2
        else # ((2*r-rp+q >= 0) && (2*r-rp+q >= 2*d+4))
        t  = r
        r  = 2*r-rp+q-2*d-4
        rp = t
        q  = q-8
        d  = d+2
        end
    end
"""

mannadiv = """
    while y3 != 0
        if y2 + 1 == x2
            y1 = y1 + 1
            y2 = 0
            y3 = y3 - 1
        else
            y2 = y2 + 1
            y3 = y3 - 1
        end
    end
"""

divbin = """
    while b != B
        x = 2*x
        b = b/2
        if r >= b
            x = x+1
            r = r-b
        end
    end
"""


text:
using MCAnalyzer
using Test
using InteractiveUtils

f() = mark_start()
g() = mark_end()

buf = IOBuffer()
code_native(buf, f, Tuple{})
asm = String(take!(buf))
@test occursin("movl\t\$111, %ebx", asm)
buf = IOBuffer()
code_native(buf, g, Tuple{})
asm = String(take!(buf))
@test occursin("movl\t\$222, %ebx", asm)

function mysum(A)
    acc = zero(eltype(A))
    for i in eachindex(A)
        mark_start()
        @inbounds acc += A[i]
    end
    mark_end()
    return acc
end
code_native(buf, mysum, Tuple{Vector{Float64}})
asm = String(take!(buf))
@test occursin("movl\t\$111, %ebx", asm)
@test occursin("movl\t\$222, %ebx", asm)

@test_nowarn analyze(mysum, Tuple{Vector{Float64}})

using LinearAlgebra

function mynorm(out, X)
    for i in 1:length(X)
        @inbounds out[i] = X[i] / norm(X)
    end
    out
end

fname, io = mktemp()
@test_nowarn MCAnalyzer.code_llvm(io, mynorm, Tuple{Vector{Float64}, Vector{Float64}})
close(io)
rm(fname)


text:
include("../coin-or-common.jl")

name = "Ipopt"
version = Ipopt_version

sources = [
    GitSource("https://github.com/coin-or/Ipopt.git", Ipopt_gitsha),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Ipopt*

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la

LIBASL="-lasl"
if [[ "${target}" == *-linux-* ]]; then
  LIBASL="${LIBASL} -lrt"
fi

export CPPFLAGS="${CPPFLAGS} -I${includedir} -I${includedir}/coin"
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
fi

./configure \
    --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --enable-shared \
    --enable-static \
    --with-pic \
    --disable-dependency-tracking \
    lt_cv_deplibs_check_method=pass_all \
    --with-lapack-lflags=-lopenblas \
    --with-mumps-cflags="-I${includedir}/mumps_seq" \
    --with-mumps-lflags="-ldmumps -lzmumps -lcmumps -lsmumps -lmumps_common -lmpiseq -lpord -lmetis -lopenblas -lgfortran -lpthread" \
    --with-asl-lflags="${LIBASL}"

# parallel build fails
make
make install
"""

platforms = supported_platforms(;experimental=true)
platforms = expand_cxxstring_abis(platforms)
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libipopt", :libipopt),
    LibraryProduct("libipoptamplinterface", :libipoptamplinterface),
    ExecutableProduct("ipopt", :amplexe),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ASL_jll", ASL_version),
    Dependency("MUMPS_seq_jll", compat="=$(MUMPS_seq_version)"),
    Dependency("OpenBLAS32_jll", OpenBLAS32_version),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version = gcc_version,
    julia_compat = "1.6"
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct SAddHandler end

gethandler(::Bool, ::Type{typeof(+)}, ::Type{<:SBitstreamLike}, ::Type{<:SBitstreamLike}) = SAddHandler()
init_state(::SAddHandler) = (id = 0,)

function (handler::SAddHandler)(buffer, netlist, state, inputs, outputs)
    # update netlist with inputs
    setsigned!(netlist, inputs[1], true)
    setsigned!(netlist, inputs[2], true)

    # compute output naming
    lname, rname = handle_broadcast_name(name.(inputs), netsize.(inputs), netsize(outputs[1]))

    # add broadcast signals
    push!(netlist, Net(name = "add$(state.id)_$(name(inputs[1]))_bcast", size = netsize(outputs[1]), signed = true))
    push!(netlist, Net(name = "add$(state.id)_$(name(inputs[2]))_bcast", size = netsize(outputs[1]), signed = true))

    # update netlist with output
    setsigned!(netlist, outputs[1], true)

    outsize = netsize(outputs[1])
    write(buffer, """
        $stdcomment
        // BEGIN add$(state.id)
        """)
    write(buffer, """
        assign add$(state.id)_$(name(inputs[1]))_bcast_p = $(lname("_p"));
        assign add$(state.id)_$(name(inputs[1]))_bcast_m = $(lname("_m"));
        assign add$(state.id)_$(name(inputs[2]))_bcast_p = $(rname("_p"));
        assign add$(state.id)_$(name(inputs[2]))_bcast_m = $(rname("_m"));
        """)
    write_bcast_instantiation(buffer, "add$(state.id)", outsize, """
        stoch_add add$(state.id)_pp (
                .CLK(CLK),
                .nRST(nRST),
                .a(add$(state.id)_$(name(inputs[1]))_bcast_p[add$(state.id)_i]),
                .b(add$(state.id)_$(name(inputs[2]))_bcast_p[add$(state.id)_i]),
                .y($(name(outputs[1]))_p[add$(state.id)_i])
            );
        stoch_add add$(state.id)_mm (
                .CLK(CLK),
                .nRST(nRST),
                .a(add$(state.id)_$(name(inputs[1]))_bcast_m[add$(state.id)_i]),
                .b(add$(state.id)_$(name(inputs[2]))_bcast_m[add$(state.id)_i]),
                .y($(name(outputs[1]))_m[add$(state.id)_i])
            );""")
    write(buffer, """
            // END add$(state.id)
            \n""")

    return buffer, (id = state.id + 1,)
end


text:
## Functions to plot SimpleSDMLayer elements more easily

# Plot layer as a heatmap with worldmap background
"""
    plot_layer(layer::SimpleSDMLayer; shape=true, kw...)

Plots the layer over a grey background based on a WorldClim layer. Note that the
warning about multiple series sharing a color bar is normal. If `shape=true`, it
will also draw a shape countour around the land based on the shapefile return by
`worldmap` (this is slower). Normal plot options can be passed afterwards.
"""
function plot_layer(layer::SimpleSDMLayer; shape=true, kw...)
    ## Arguments
    # layer: SimpleSDMLayer to plot
    # scatter: add observations as points in scatter plot, requires to define occ
    # occ: observations to represent if scatter=true
    # kw: optional plotting arguments

    # Create background layer
    coords = boundingbox(layer)
    baselayer = similar(SimpleSDMPredictor(WorldClim, BioClim, 1; coords...))

    # Create empty plot
    sdm_plot = heatmap(baselayer; c=:lightgrey, lab="", msw=0.0, ms=0.0, frame=:box)

    # Add SDM output as heatmap
    heatmap!(
        sdm_plot,
        layer;
        legend=true,
        clim=extrema(layer), # colorbar limits
        kw..., # additional keyword arguments
    )

    # Add shape countour around land
    if shape
        # Load & clip worldmap background to SimpleSDMLayer (from shp in /assets folder)
        worldmap = clip(worldshape(50), layer)
        # Redraw polygons' outer lines over heatmap values
        for p in worldmap # loop for each polygon
            # Get outer lines coordinates
            xy = map(x -> (x.x, x.y), p.points)
            # Add outer lines to plot
            plot!(sdm_plot, xy; c=:grey, lab="")
        end
    end

    return sdm_plot
end

# Custom plot recipe for layers
@recipe function plot(layer::T) where {T<:SimpleSDMLayer}
    seriestype --> :heatmap
    if get(plotattributes, :seriestype, :heatmap) in [:heatmap, :contour]
        aspect_ratio --> 92.60 / 60.75
        xlims --> (layer.left, layer.right)
        ylims --> (layer.bottom, layer.top)
        xguide --> "Longitude"
        yguide --> "Latitude"
        lg = copy(layer.grid)
        lg[lg .== nothing] .= NaN
        longitudes(layer), latitudes(layer), lg
    elseif get(plotattributes, :seriestype, :histogram) in [:histogram, :density]
        filter(!isnothing, layer.grid)
    end
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
module Theta0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    [                                                                           -12828156665.2727*xi.^10.*(2*l1 - l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 2004399.47894886*xi.^8.*(2*l1 - l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 7990976.98521566*xi.^6.*(2*l1 - l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 84279.8353909465*xi.^4.*(2*l1 - l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 355.555555555555*xi.^2.*(2*l1 - l2 - l3).^2 + 1 12828156665.2727*sqrt(3)*xi.^10.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 405890894.487145*sqrt(3)*xi.^8.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 7990976.98521566*sqrt(3)*xi.^6.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 84279.8353909465*sqrt(3)*xi.^4.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 355.555555555555*sqrt(3)*xi.^2.*(l2 - l3).*(2*l1 - l2 - l3) -2405279374.73864*xi.^9.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 30441817.0865359*xi.^7.*(8*l1 - 4*l2 - 4*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 898984.910836762*xi.^5.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 6320.98765432099*xi.^3.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 26.6666666666667*xi.*(2*l1 - l2 - l3);
12828156665.2727*sqrt(3)*xi.^10.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 405890894.487145*sqrt(3)*xi.^8.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 7990976.98521566*sqrt(3)*xi.^6.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 84279.8353909465*sqrt(3)*xi.^4.*(l2 - l3).*(2*l1 - l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 355.555555555555*sqrt(3)*xi.^2.*(l2 - l3).*(2*l1 - l2 - l3)                                                                                -38484469995.8182*xi.^10.*(l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 1217672683.46143*xi.^8.*(l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 23972930.955647*xi.^6.*(l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 252839.506172839*xi.^4.*(l2 - l3).^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 1066.66666666667*xi.^2.*(l2 - l3).^2 + 1             -4810558749.47727*sqrt(3)*xi.^9.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 + 121767268.346143*sqrt(3)*xi.^7.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 - 1797969.82167352*sqrt(3)*xi.^5.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 + 12641.975308642*sqrt(3)*xi.^3.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) - 26.6666666666667*sqrt(3)*xi.*(l2 - l3);
2405279374.73864*xi.^9.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 - 30441817.0865359*xi.^7.*(8*l1 - 4*l2 - 4*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 + 898984.910836762*xi.^5.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 - 6320.98765432099*xi.^3.*(4*l1 - 2*l2 - 2*l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) + 26.6666666666667*xi.*(2*l1 - l2 - l3)                                                                                4810558749.47727*sqrt(3)*xi.^9.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 - 121767268.346143*sqrt(3)*xi.^7.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 + 1797969.82167352*sqrt(3)*xi.^5.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 - 12641.975308642*sqrt(3)*xi.^3.*(l2 - l3).*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) + 26.6666666666667*sqrt(3)*xi.*(l2 - l3)                                            -51312626661.0909*xi.^10.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^5 + 1623563577.94858*xi.^8.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^4 - 31963907.9408627*xi.^6.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^3 + 337119.341563786*xi.^4.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2).^2 - 1422.22222222222*xi.^2.*(l1.^2 - l1.*l2 - l1.*l3 + l2.^2 - l2.*l3 + l3.^2) + 1]
end
end


text:
using StochasticBlockModel
using Test

@testset "Dataset struct" begin

    dataset = StochasticBlockModel.Dataset(Matrix([0 1; 1 0]), 2)
    @test dataset.n == 2
    @test dataset.m == 1
    @test dataset.k == [1; 1]

    expected_A = [ 2  1  0  1  0  0  0  0
                    1  2  1  0  1  0  0  1
                    0  1  0  1  1  1  0  0
                    1  0  1  2  0  0  1  0
                    0  1  1  0  2  0  1  0
                    0  0  1  0  0  2  1  1
                    0  0  0  1  1  1  0  1
                    0  1  0  0  0  1  1  2]
    dataset = StochasticBlockModel.Dataset("../instances/S1/g01.in")
    @test dataset.A == expected_A
    @test dataset.n == 8
    @test dataset.m == 19
    @test dataset.n_communities == 2
    @test dataset.k == [4; 6; 4; 5; 5; 5; 4; 5]

    # Test if self-edges correctly count as 2 in the adjacency matrix
    @test_throws ArgumentError StochasticBlockModel.Dataset(Matrix([1 0; 0 0]), 2)
    # Test if negative values in the adjacency matrix cause an Exception
    @test_throws DomainError StochasticBlockModel.Dataset(Matrix([-1 0; 0 0]), 2)
    # Test if zero value in the number of clusters cause an Exception
    @test_throws DomainError StochasticBlockModel.Dataset(Matrix([2 0; 0 0]), 0)
    # Test if negative value in the number of clusters cause an Exception
    @test_throws DomainError StochasticBlockModel.Dataset(Matrix([2 1; 1 0]), -2)

    @test_throws ArgumentError StochasticBlockModel.Dataset(Matrix([2 1; 1 0]), 4, 2, 2, [3; 1])
    @test_throws ArgumentError StochasticBlockModel.Dataset(Matrix([2 1; 1 0]), 2, 4, 2, [3; 1])
    @test_throws DomainError StochasticBlockModel.Dataset(Matrix([-2 1; 1 0]), 2, 2, 2, [3; 1])
    @test_throws DomainError StochasticBlockModel.Dataset(Matrix([2 1; 1 0]), 2, 2, -2, [3; 1])

    A = [0 1 1;
         0 0 0;
         0 0 0]
    @test_throws ArgumentError StochasticBlockModel.Dataset(A, 3, 1, 1, [0; 1; 1])

    # Get a list of all files
    path = "../instances/S1"
    for file in readdir(path)
        filepath = string(path, "/", file)
        dataset = StochasticBlockModel.Dataset(filepath)
        StochasticBlockModel.saveDataset(dataset, string(filepath,".test"))
        dataset_test = StochasticBlockModel.Dataset(string(filepath,".test"))
        rm(string(filepath,".test"))
        @test dataset.n == dataset_test.n
        @test dataset.m == dataset_test.m
        @test dataset.n_communities == dataset_test.n_communities
        @test dataset.k == dataset_test.k
        @test dataset.A == dataset_test.A
    end
end


text:
SaveAdaptive(sol::AbstractODESolution, N::Int=500; kwargs...) = SaveAdaptive(sol, (sol.t[1],sol.t[end]); N=N, kwargs...)
function SaveAdaptive(sol::Union{Function,AbstractODESolution}, Tspan::Tuple{Real,Real}; N::Int=500, curvature::Real=0.003, Ntol::Real=0.08, maxiter::Int=30)
    @assert Tspan[1] < Tspan[2]
    test = sol(Tspan[1])
    for _ in 1:maxiter
        T = reduce(vcat, [PlotUtils.adapted_grid(x->sol(x)[i], Tspan; max_curvature=curvature)[1] for i in 1:length(test)]) |> unique! |> sort!
        if length(T) > N
            curvature *= 1.2
        elseif length(T) < (1-Ntol) * N
            curvature *= 0.5
        else
            return Homogenize(T,N)
        end
    end
    throw("SaveAdaptive: DNF in $maxiter iterations.")
end

Homogenize(sol::AbstractODESolution, N::Int=500) = Homogenize(sol.t, N)
function Homogenize(V::AbstractVector, N::Int=500)
    Ts = unique(V)
    for i in 1:(N-length(Ts))
        s = findmax(diff(Ts))[2]
        insert!(Ts,s+1,Ts[s] + (Ts[s+1]-Ts[s])/2)
    end;    Ts
end

Dehomogenize(sol::AbstractODESolution, N::Int=500) = Dehomogenize(sol.t, N)
function Dehomogenize(V::AbstractVector, N::Int=500)
    Ts = unique(V)
    for i in 1:(length(Ts)-N)
        s = findmin(diff(Ts))[2]
        deleteat!(Ts,s+1)
    end;    Ts
end

"""
    SaveConfidence(sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true) -> Matrix
    SaveConfidence(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true) -> Matrix
Returns a `Matrix` of with `N` rows corresponding to the number of evaluations of each `ODESolution` in `sols`.
The colums correspond to the various components of the evaluated solutions.
E.g. for an `ODESolution` with 3 components, the 4. column in the `Matrix` corresponds to the evaluated first components of `sols[2]`.
"""
function SaveConfidence(sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    mapreduce(sol->SaveConfidence(sol,N; sigdigits=sigdigits,adaptive=adaptive), hcat, sols)
end
function SaveConfidence(sol::AbstractODESolution, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    Ts = adaptive ? SaveAdaptive(sol, N) : range(sol.t[1], sol.t[end]; length=N)
    round.(Unpack(map(sol,Ts)); sigdigits=sigdigits)
end

function SaveConfidence(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    length(Planes) != length(sols) && throw("Dimensional Mismatch: length(Planes)=$(length(Planes)), length(sols)=$(length(sols)).")
    mapreduce((PL,sol)->SaveConfidence(PL, sol, N; sigdigits=sigdigits,adaptive=adaptive), hcat, Planes, sols)
end
function SaveConfidence(PL::Plane, sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    mapreduce(sol->SaveConfidence(PL,sol,N; sigdigits=sigdigits,adaptive=adaptive), hcat, sols)
end
function SaveConfidence(PL::Plane, sol::AbstractODESolution, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    Ts = adaptive ? SaveAdaptive(sol, N) : range(sol.t[1], sol.t[end]; length=N)
    round.(Deplanarize(PL,sol,Ts); sigdigits=sigdigits)
end
# function SaveConfidence(sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
#     d = length(sols[1].u[1]);    Res = Array{Float64}(undef,N,d*length(sols))
#     for (i,sol) in enumerate(sols)
#         Res[:,((i-1)*d+1):(d*i)] = SaveConfidence(sol, N; sigdigits=sigdigits, adaptive=adaptive)
#     end;    Res
# end
# function SaveConfidence(sol::AbstractODESolution, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
#     T = adaptive ? SaveAdaptive(sol, N) : range(sol.t[1], sol.t[end]; length=N)
#     round.(Unpack(sol.(T)); sigdigits=sigdigits)
# end

"""
    SaveGeodesics(sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true) -> Matrix
Returns a `Matrix` of with `N` rows corresponding to the number of evaluations of each `ODESolution` in `sols`.
The colums correspond to the various components of the evaluated solutions.
Since the solution objects for geodesics contain the velocity as the second half of the components, only the first half of the components is saved.
"""
function SaveGeodesics(sols::AbstractVector{<:AbstractODESolution}, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    d = length(sols[1].u[1])/2 |> Int;      Res = Array{Float64}(undef,N,d*length(sols))
    for (i,sol) in enumerate(sols)
        Res[:,((i-1)*d+1):(d*i)] = SaveGeodesics(sol, N; sigdigits=sigdigits, adaptive=adaptive)
    end;    Res
end
function SaveGeodesics(sol::AbstractODESolution, N::Int=500; sigdigits::Int=7, adaptive::Bool=true)
    SaveConfidence(sol, N; sigdigits=sigdigits, adaptive=adaptive)[:, 1:Int(length(sol.u[1])/2)]
end

"""
    SaveDataSet(DS::DataSet; sigdigits::Int=0)
Returns a `DataFrame` whose columns respectively constitute the x-values, y-values and standard distributions associated with the data points.
For `sigdigits > 0` the values are rounded to the specified number of significant digits.
"""
function SaveDataSet(DS::AbstractDataSet; sigdigits::Int=0)
    !(xdim(DS) == ydim(DS) == Int(size(ysigma(DS),1)/Npoints(DS))) && throw("Not programmed yet.")
    sig = ysigma(DS)
    !(typeof(sig) <: AbstractVector) && throw("Sigma not a vector, but instead $(typeof(sig)).")
    if sigdigits < 1
        return DataFrame([xdata(DS) ydata(DS) ysigma(DS)], :auto)
    else
        return DataFrame(round.([xdata(DS) ydata(DS) ysigma(DS)]; sigdigits=sigdigits), :auto)
    end
end
SaveDataSet(DM::AbstractDataModel; sigdigits::Int=0) = SaveDataSet(Data(DM); sigdigits=sigdigits)


text:
using bad

# construct the usual prior
pnull  = .2 # response under TAU
pmcr   = pnull + .1 # minimal clinically relevant
prior1 = Beta(mean = .35, sd = .1)  # start with an subjective prior
prior2 = update(prior1, 4, 10)      # update with phase I data
prior3 = .8*prior2 + .2*Beta(1, 1)  # robustify
prior  = prior3 <= min(2*pmcr, 1.0) # restrict to plausible range

mtoer = Power(prior  | pnull)
power = Power(prior >= pmcr)
ess   = SampleSize(prior)

α, β = .05, .2
problem = Problem(
        minimise(ess),
        mtoer <= α,
        power >= 1 - β
    )
design = optimise(problem; verbosity = 0)

prior_ = (0.8*update(prior2, 2, 10) + 0.2*Beta(1, 1)) <= 0.7
x1 = 6
partial_stage_two = 2, 6
design_ = deepcopy(design)
n2_original, c2_original = n2(design, x1), c2(design, x1)
toer = design.problem.toer.score
pow  = design.problem.power.score
α_   = toer(design, x1, partial_stage_two = partial_stage_two)
β_   = 1 - pow(design, x1, partial_stage_two = partial_stage_two)
function ctoer(n2_new, c2_new)
    design_.n2[x1 + 1] = n2_new
    design_.c2[x1 + 1] = c2_new
    toer(design_, x1, partial_stage_two = partial_stage_two)
end
pow = update(pow, prior_)
function cep(n2_new, c2_new)
    design_.n2[x1 + 1] = n2_new
    design_.c2[x1 + 1] = c2_new
    pow(design_, x1, partial_stage_two = partial_stage_two)
end

function f()
    n2_new, c2_new = partial_stage_two[2], -Inf
    a, b = ctoer(n2_new, c2_new), cep(n2_new, c2_new)
    while (a > α_) | (b < 1 - β_)
        for c2_new in vcat([-Inf], 0:(n2_new - 1), [Inf])
            a, b = ctoer(n2_new, c2_new), cep(n2_new, c2_new)
            if (a <= α_) & (b >= 1 - β_)
                return n2_new, c2_new
            end
        end
        n2_new += 1
    end
end

α_, 1 - β_
n2_new, c2_new = 25, 7.0
ctoer(n2_new, c2_new), cep(n2_new, c2_new)

xx1, nn2, cc2 = adapt(design, prior_, (2, 7), verbosity = 3)

xx1, nn2, cc2 = adapt(design, prior_, (5, 7), verbosity = 3)


text:
import LinearAlgebra
using StatsFuns
using SpecialFunctions
using OffsetArrays

### MvUniform
struct MvUniform <: Gen.Distribution{Vector{Float64}} end

const mvuniform = MvUniform()

function Gen.random(::MvUniform, lows::AbstractVector{T},
            highs::AbstractVector{U}) where {T<:Real, U<:Real}
    return [Gen.random(uniform, lows[i], highs[i]) for i in 1:length(lows)]
end

function Gen.logpdf(::MvUniform, x::Array{Float64, 1}, lows::Array{Float64, 1},
            highs::Array{Float64,1})
    if any(x < lows || x > highs)
        return -Inf
    else
        return sum([Gen.logpdf(uniform, x[i], lows[i], highs[i]) for i in 1:length(lows)])
    end
end

function Gen.logpdf_grad(::MvUniform, x::AbstractVector{T},
            lows::AbstractVector{U},
            highs::AbstractVector{V}) where {T<:Real, U<:Real, V<:Real}
    (nothing, nothing, nothing)
end

(::MvUniform)(lows, highs) = Gen.random(MvUniform(), lows, highs)

has_output_grad(::MvUniform) = false
has_argument_grads(::MvUniform) = (false, false)

export mvuniform

### Direction
struct ChooseDirection <: Gen.Distribution{Vector{Float64}} end

const choose_direction = ChooseDirection()

function Gen.random(::ChooseDirection, camera_position::AbstractVector{T},
            location::AbstractVector{U}, detection_sd::V) where {T<:Real,
            U<:Real, V<:Real}
    detection_location = mvnormal(location, LinearAlgebra.Diagonal([detection_sd for _ in 1:length(location)]))

    direction = (detection_location-camera_position)/LinearAlgebra.norm(detection_location-camera_position)

    return direction
end

function Gen.logpdf(::ChooseDirection, x::AbstractVector{T}, camera_position::AbstractVector{U},
            location::AbstractVector{V}, detection_sd::W) where {T<:Real,
            U<:Real, V<:Real, W<:Real}
    k = length(x)
    constant = -0.5 * (k-1) * log(2 * pi * detection_sd)

    term_1 = 0.5 / detection_sd * (LinearAlgebra.dot(location-camera_position, x)/2 - LinearAlgebra.norm(location-camera_position)^2)

    term_2 = log(SpecialFunctions.erfc(-LinearAlgebra.dot(location - camera_position, x) / (2 * sqrt(2*detection_sd))))

    return constant + term_1 + term_2
end

function Gen.logpdf_grad(::ChooseDirection, x::AbstractVector{T}, camera_position::AbstractVector{U},
            location::AbstractVector{V}, detection_sd::W) where {T<:Real,
            U<:Real, V<:Real, W<:Real}
    return (nothing, nothing, nothing, nothing)
end

(::ChooseDirection)(position, location, sd) = Gen.random(ChooseDirection(), position, location, sd)

has_output_grad(::ChooseDirection) = false
has_argument_grads(::ChooseDirection) = false

export choose_direction

###ZeroCategorical
struct ZeroCategorical <: Gen.Distribution{Int} end

"""
    zero_categorical(probs::AbstractArray{U, 1}) where {U <: Real}
Given a vector of probabilities `probs` where `sum(probs) = 1`, sample an `Int` `i` from the set {0, 1, .., `length(probs)-1`} with probability `probs[i]`.
"""
const zero_categorical = ZeroCategorical()

function Gen.logpdf(::ZeroCategorical, x::Int, probs::AbstractArray{U,1}) where {U <: Real}
    (x >= 0 && x < length(probs)) ? log(probs[x]) : -Inf
end

function Gen.logpdf_grad(::ZeroCategorical, x::Int, probs::AbstractArray{U,1})  where {U <: Real}
    (nothing, nothing)
end

function Gen.random(::ZeroCategorical, probs::AbstractArray{U,1}) where {U <: Real}
    one_indexed_probs = OffsetArray(probs, 1:length(probs))
    rand(Distributions.Categorical(one_indexed_probs))-1
end
is_discrete(::ZeroCategorical) = true

(::ZeroCategorical)(probs) = random(ZeroCategorical(), probs)

has_output_grad(::ZeroCategorical) = false
has_argument_grads(::ZeroCategorical) = false

export zero_categorical

max_stars_repo_name: zhangir-azerbayev/MetaGen-New

max_issues_repo_name: zhangir-azerbayev/MetaGen-New

max_forks_repo_name: zhangir-azerbayev/MetaGen-New


text:
baremodule PyPreferences

function use_system end
function use_conda end
# function use_jll end
function use_inprocess end

# API to be used from PyCall
function assert_configured end
function instruction_message end

# function diagnose end
function status end

module Implementations
include("core.jl")
end

let prefs = Implementations.setup_non_failing()
    global const python = prefs.python
    global const inprocess = prefs.inprocess
    global const conda = prefs.conda
    global const python_fullpath = prefs.python_fullpath
    global const libpython = prefs.libpython
    global const python_version = prefs.python_version
    global const PYTHONHOME = prefs.PYTHONHOME
end

const pyprogramname = python_fullpath
const pyversion_build = python_version

end  # baremodule PyPreferences


text:
export CLLData
"""
    CLLData

A dataset of different omic measurements for Chronic lymphocytic leukaemia (CLL)
patient samples. The data can be loaded via:

```
cll_data = CLLData.load()
```

`cll_data` is a named tuple with fields:
* `X`: The features.
* `y`: The response, namely Ibrutinib sensitivity.
* `ngr`: The number of features in each of the three feature groupings, namely
drug sensitivity, methylation and RNAseq data.
* `foldid`: A `Vector{Int}` with values in 1,..,10 that assign each of the rows of `X` to
a fold to be used in cross-validation.

## References

The dataset was originally published in:

Dietrich, Sascha, et al. "Drug-perturbation-based stratification of blood cancer."
The Journal of clinical investigation 128.1 (2018): 427-445.

It was used in the context of side-information by:

Velten, Britta, and Wolfgang Huber.
"Adaptive penalization in high-dimensional regression and classification
 with external covariates using variational Bayes."
Biostatistics (2019).

The `foldid` assignment into folds is the same as the one used by the above publication.

The dataset was copied from the Bioconductor MOFAdata package, available at:
https://bioconductor.org/packages/release/data/experiment/html/MOFAdata.html
"""
module CLLData

using JLD2

const DATA = joinpath(@__DIR__, "cll_data.jld2")

function load()
    JLD2.@load DATA cll_data
    cll_data
end

end


text:
function Base.show(io::IO, fid::File)
    if isvalid(fid)
        intent = h5f_get_intent(fid)
        RW_MASK   = H5F_ACC_RDONLY | H5F_ACC_RDWR
        SWMR_MASK = H5F_ACC_SWMR_READ | H5F_ACC_SWMR_WRITE
        rw = (intent & RW_MASK) == H5F_ACC_RDONLY ? "(read-only" : "(read-write"
        swmr = (intent & SWMR_MASK) != 0 ? ", swmr) " : ") "
        print(io, "HDF5.File: ", rw, swmr, fid.filename)
    else
        print(io, "HDF5.File: (closed) ", fid.filename)
    end
end

function Base.show(io::IO, g::Group)
    if isvalid(g)
        print(io, "HDF5.Group: ", name(g), " (file: ", g.file.filename, ")")
    else
        print(io, "HDF5.Group: (invalid)")
    end
end

function Base.show(io::IO, prop::Properties)
    if prop.class == H5P_DEFAULT
        print(io, "HDF5.Properties: default class")
    elseif isvalid(prop)
        print(io, "HDF5.Properties: ", h5p_get_class_name(prop.class), " class")
    else
        print(io, "HDF5.Properties: (invalid)")
    end
end

function Base.show(io::IO, dset::Dataset)
    if isvalid(dset)
        print(io, "HDF5.Dataset: ", name(dset), " (file: ", dset.file.filename, " xfer_mode: ", dset.xfer.id, ")")
    else
        print(io, "HDF5.Dataset: (invalid)")
    end
end

function Base.show(io::IO, attr::Attribute)
    if isvalid(attr)
        print(io, "HDF5.Attribute: ", name(attr))
    else
        print(io, "HDF5.Attribute: (invalid)")
    end
end
function Base.show(io::IO, attr::Attributes)
    print(io, "Attributes of ", attr.parent)
end

function Base.show(io::IO, dtype::Datatype)
    print(io, "HDF5.Datatype: ")
    if isvalid(dtype)
        h5t_committed(dtype) && print(io, name(dtype), " ")
        print(io, h5lt_dtype_to_text(dtype))
    else
        # Note that h5i_is_valid returns `false` on the built-in datatypes (e.g. H5T_NATIVE_INT),
        # apparently because they have refcounts of 0 yet are always valid. Just temporarily turn
        # off error printing and try the call to probe if dtype is valid since H5LTdtype_to_text
        # special-cases all of the built-in types internally.
        local text
        try
            text = silence_errors(() -> h5lt_dtype_to_text(dtype))
        catch
            text = "(invalid)"
        end
        print(io, text)
    end
end

function Base.show(io::IO, dspace::Dataspace)
    if !isvalid(dspace)
        print(io, "HDF5.Dataspace: (invalid)")
        return
    end
    print(io, "HDF5.Dataspace: ")
    type = h5s_get_simple_extent_type(dspace)
    if type == H5S_NULL
        print(io, "H5S_NULL")
        return
    elseif type == H5S_SCALAR
        print(io, "H5S_SCALAR")
        return
    end
    # otherwise type == H5S_SIMPLE
    sz, maxsz = get_extent_dims(dspace)
    sel = h5s_get_select_type(dspace)
    if sel == H5S_SEL_HYPERSLABS && h5s_is_regular_hyperslab(dspace)
        start, stride, count, _ = get_regular_hyperslab(dspace)
        ndims = length(start)
        print(io, "(")
        for ii in 1:ndims
            s, d, l = start[ii], stride[ii], count[ii]
            print(io, range(s + 1, length = l, step = d == 1 ? nothing : d))
            ii != ndims && print(io, ", ")
        end
        print(io, ") / (")
        for ii in 1:ndims
            print(io, 1:maxsz[ii])
            ii != ndims && print(io, ", ")
        end
        print(io, ")")
    else
        print(io, sz)
        if maxsz != sz
            print(io, " / ", maxsz)
        end
        if sel != H5S_SEL_ALL
            print(io, " [irregular selection]")
        end
    end
end

"""
    SHOW_TREE_ICONS = Ref{Bool}(true)

Configurable option to control whether emoji icons (`true`) or a plain-text annotation
(`false`) is used to indicate the object type by `show_tree`.
"""
const SHOW_TREE_ICONS = Ref{Bool}(true)

"""
    SHOW_TREE_MAX_DEPTH = Ref{Int}(5)

Maximum recursive depth to descend during printing.
"""
const SHOW_TREE_MAX_DEPTH = Ref{Int}(5)

"""
    SHOW_TREE_MAX_CHILDREN = Ref{Int}(50)

Maximum number of children to show at each node.
"""
const SHOW_TREE_MAX_CHILDREN = Ref{Int}(50)

function Base.show(io::IO, ::MIME"text/plain", obj::Union{File,Group,Dataset,Attributes,Attribute})
    if get(io, :compact, false)::Bool
        show(io, obj)
    else
        show_tree(io, obj)
    end
end

_tree_icon(::Type{Attribute}) = SHOW_TREE_ICONS[] ? "🏷" : "[A]"
_tree_icon(::Type{Group})     = SHOW_TREE_ICONS[] ? "📂" : "[G]"
_tree_icon(::Type{Dataset})   = SHOW_TREE_ICONS[] ? "🔢" : "[D]"
_tree_icon(::Type{Datatype})  = SHOW_TREE_ICONS[] ? "📄" : "[T]"
_tree_icon(::Type{File})      = SHOW_TREE_ICONS[] ? "🗂" : "[F]"
_tree_icon(::Type)            = SHOW_TREE_ICONS[] ? "❓" : "[?]"
_tree_icon(obj) = _tree_icon(typeof(obj))
_tree_icon(obj::Attributes) = _tree_icon(obj.parent)

_tree_head(io::IO, obj) = print(io, _tree_icon(obj), " ", obj)
_tree_head(io::IO, obj::Datatype) = print(io, _tree_icon(obj), " HDF5.Datatype: ", name(obj))

_tree_count(parent::Union{File,Group}, attributes::Bool) =
    length(parent) + (attributes ? length(HDF5.attributes(parent)) : 0)
_tree_count(parent::Dataset, attributes::Bool) =
    attributes ? length(HDF5.attributes(parent)) : 0
_tree_count(parent::Attributes, _::Bool) = length(parent)
_tree_count(parent::Union{Attribute,Datatype}, _::Bool) = 0

function _show_tree(io::IO, obj::Union{File,Group,Dataset,Datatype,Attributes,Attribute}, indent::String="";
                    attributes::Bool = true, depth::Int = 1)
    isempty(indent) && _tree_head(io, obj)
    isvalid(obj) || return

    INDENT = "   "
    PIPE   = "│  "
    TEE    = "├─ "
    ELBOW  = "└─ "

    limit = get(io, :limit, false)::Bool
    counter = 0
    nchildren = _tree_count(obj, attributes)

    @inline function childstr(io, n, more=" ")
        print(io, "\n", indent, ELBOW * "(", n, more, n == 1 ? "child" : "children", ")")
    end
    @inline function depth_check()
        counter += 1
        if limit && counter > max(2, SHOW_TREE_MAX_CHILDREN[] ÷ depth)
            childstr(io, nchildren - counter + 1, " more ")
            return true
        end
        return false
    end

    if limit && nchildren > 0 && depth > SHOW_TREE_MAX_DEPTH[]
        childstr(io, nchildren)
        return nothing
    end

    if attributes && !isa(obj, Attribute)
        obj′ = obj isa Attributes ? obj.parent : obj
        h5a_iterate(obj′, H5_INDEX_NAME, H5_ITER_INC) do _, cname, _, _
            depth_check() && return herr_t(1)

            name = unsafe_string(cname)
            icon = _tree_icon(Attribute)
            islast = counter == nchildren
            print(io, "\n", indent, islast ? ELBOW : TEE, icon, " ", name)
            return herr_t(0)
        end
    end

    typeof(obj) <: Union{File, Group} || return nothing

    h5l_iterate(obj, H5_INDEX_NAME, H5_ITER_INC) do loc_id, cname, _, _
        depth_check() && return herr_t(1)

        name = unsafe_string(cname)
        child = obj[name]
        icon = _tree_icon(child)

        islast = counter == nchildren
        print(io, "\n", indent, islast ? ELBOW : TEE, icon, " ", name)
        nextindent = indent * (islast ? INDENT : PIPE)
        _show_tree(io, child, nextindent; attributes = attributes, depth = depth + 1)

        close(child)
        return herr_t(0)
    end
    return nothing
end

show_tree(obj; kws...) = show_tree(stdout, obj; kws...)
function show_tree(io::IO, obj; kws...)
    buf = IOBuffer()
    _show_tree(IOContext(buf, io), obj; kws...)
    print(io, String(take!(buf)))
end


text:
# generate samples for detTest and stabilityTest
using Distributed;
addprocs(30);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

Δt = 0.25;
N = 5;
NN = 1000;
NN2 = 5000;
NNTest = 600;

caseList = [13,33,123];
dataDict = Dict();
dataDict2 = Dict();
dataDictTest = Dict();

for ci in 1:length(caseList)
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathListData = pmap(i -> simuPath(T,pDistr), 1:NN);
        pathDict[T] = Dict();
        for i in 1:length(pathListData)
            pathDict[T][i] = pathListData[i];
        end
    end
    dataDict[ci] = pathDict;
end
# reverse the old generated sample
for ci in 1:2
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathDict[T] = reverseScen(dataDict[ci][T],τ,pDistr);
    end
    dataDict[ci] = pathDict;
end
save("pathHist_1000.jld","pathDict",dataDict);

for ci in 1:length(caseList)
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathListData = pmap(i -> simuPath(T,pDistr), 1:NNTest);
        pathDict[T] = Dict();
        for i in 1:length(pathListData)
            pathDict[T][i] = pathListData[i];
        end
    end
    dataDictTest[ci] = pathDict;
end
# reverse the old generated sample
for ci in 1:2
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathDict[T] = reverseScen(dataDictTest[ci][T],τ,pDistr);
    end
    dataDictTest[ci] = pathDict;
end

save("pathHist_600.jld","pathDict",dataDictTest);

for ci in 1:length(caseList)
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathListData = pmap(i -> simuPath(T,pDistr), 1:NN2);
        pathDict[T] = Dict();
        for i in 1:length(pathListData)
            pathDict[T][i] = pathListData[i];
        end
    end
    dataDict2[ci] = pathDict;
end
for ci in 1:2
    pathDict = Dict();
    TList = [24,36,48,72,96];
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end

        pathDict[T] = reverseScen(dataDict2[ci][T],τ,pDistr);
    end
    dataDict2[ci] = pathDict;
end

save("pathHist_5000.jld","pathDict",dataDict2);


text:
@testset "383.ransom-note.jl" begin
    @test can_construct("a", "b") == false
    @test can_construct("aa", "ab") == false
    @test can_construct("aa", "aba") == true
end


text:
using Documenter, AdaMeanShift

makedocs(;
    modules=[AdaMeanShift],
    format = Documenter.HTML(assets=["assets/invenia.css"]),
    pages=[
        "API" => "index.md",
    ],
    sitename="AdaMeanShift.jl",
    authors="Francesco Alemanno",
)

deploydocs(;
    repo="github.com/francescoalemanno/AdaMeanShift.jl.git",
    target="build"
)


text:
#=
  https://edu.swi-prolog.org/mod/assign/view.php?id=254&forceview=1
  """
  Coin tosses

  http://cplint.eu/p/coin_tosses.swinb

  Coin tosses

  Consider a process where you repeatedly toss coins and record the results.

  After each toss, you continue tossing with probability 0.8.

  Write a predicate tosses(Coin,L) that, given an initial coin id Coin, returns in L the list of
  results of coin tosses obtained using the process above.

  Moreover, write a predicate length_tosses(N) that returns the number of coin tosses.

  Compute the probability of the sequence [h,h,h] using MCINTYRE.

  Compute the probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] was observed
  using rejection sampling and Metropolis Hastings.

  Compute the probability of the sequences of 10 coin tosses using MCINTYRE.

  Compute the expected length of the sequences of coin tosses using MCINTYRE.
  """

  Note: The cplint version (corrected by  Fabrizio Riguzzi) at http://cplint.eu/p/coin_tosses_hakank_rzf.swinb
  give another solutions:
   - first problem: probability of the sequence [h,h,h] -> 0.015
   - second problem: probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] -> about 0.03
   - third problem: probability of the sequences of 10 coin tosses -> about 0.026
   - fourth problem: expected length of the sequences of coin tosses -> about 4.

  Cf ~/blog/coin_tosses.blog
     ~/webppl/coin_tosses.wppl

  Note: compared to BLOG and WebPPL models, this is a for loop solutions.
        (To check: It seems that using recursion as in those models don't
                   work in Turing.
        )

=#

#=
ulia> @time include("coin_tosses.jl")
Problem 1
Probability of problem 1:

Distributions of variable (num:20)
0.00000 =>   35055  (0.876375)
1.00000 =>    4945  (0.123625)

Problem 2
Probability of problem 2:

Distributions of variable (num:20)
0.00000 =>   29875  (0.746875)
1.00000 =>   10125  (0.253125)

Problem 3
Probability of problem 3:

Distributions of variable (num:20)
0.00000 =>   38991  (0.974775)
1.00000 =>    1009  (0.025225)

Problem 4
chns MCMC chain (10000×14×4 Array{Float64,3}):

Iterations        = 1:10000
Thinning interval = 1
chns            = 1, 2, 3, 4
Samples per chain = 10000
parameters        = i, tosses[1], tosses[2], tosses[3], tosses[4], tosses[5], tosses[6], tosses[7], tosses[8], tosses[9], tosses[10]
internals         = le, lp, weight

Summary Statistics
  parameters      mean       std   naive_se      mcse          ess      ⋯
      Symbol   Float64   Float64    Float64   Float64      Float64   Fl ⋯

           i    5.4932    2.8637     0.0143    0.0150   39858.1687    1 ⋯
   tosses[1]    1.5024    0.5000     0.0025    0.0024   38948.5470    0 ⋯
   tosses[2]    1.5983    0.6639     0.0033    0.0032   38547.3980    1 ⋯
   tosses[3]    1.6799    0.7610     0.0038    0.0039   39447.3220    0 ⋯
   tosses[4]    1.7387    0.8228     0.0041    0.0041   39791.3739    1 ⋯
   tosses[5]    1.8036    0.8677     0.0043    0.0040   39909.2109    1 ⋯
   tosses[6]    1.8379    0.8995     0.0045    0.0044   40658.2973    0 ⋯
   tosses[7]    1.8668    0.9222     0.0046    0.0045   40137.0842    1 ⋯
   tosses[8]    1.8955    0.9403     0.0047    0.0047   39153.5236    1 ⋯
   tosses[9]    1.9117    0.9531     0.0048    0.0049   38279.3118    1 ⋯
  tosses[10]    1.9292    0.9628     0.0048    0.0051   39562.1308    0 ⋯
                                                         1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           i    1.0000    3.0000    5.0000    8.0000   10.0000
   tosses[1]    1.0000    1.0000    2.0000    2.0000    2.0000
   tosses[2]    1.0000    1.0000    1.0000    2.0000    3.0000
   tosses[3]    1.0000    1.0000    1.0000    2.0000    3.0000
   tosses[4]    1.0000    1.0000    1.0000    2.0000    3.0000
   tosses[5]    1.0000    1.0000    2.0000    3.0000    3.0000
   tosses[6]    1.0000    1.0000    2.0000    3.0000    3.0000
   tosses[7]    1.0000    1.0000    1.0000    3.0000    3.0000
   tosses[8]    1.0000    1.0000    2.0000    3.0000    3.0000
   tosses[9]    1.0000    1.0000    1.0000    3.0000    3.0000
  tosses[10]    1.0000    1.0000    1.0000    3.0000    3.0000

Probability of problem 4:

Distributions of variable (num:20)
1.00000 =>    8045  (0.201125)
2.00000 =>    6381  (0.159525)
3.00000 =>    5062  (0.126550)
11.00000 =>    4328  (0.108200)
4.00000 =>    4067  (0.101675)
5.00000 =>    3231  (0.080775)
6.00000 =>    2676  (0.066900)
7.00000 =>    2057  (0.051425)
8.00000 =>    1718  (0.042950)
9.00000 =>    1334  (0.033350)
10.00000 =>    1101  (0.027525)
Mean:4.580975

 24.121595 seconds (193.14 M allocations: 12.141 GiB, 35.09% gc time)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function coin_tosses(problem)
    head = 1
    tail = 2

    if problem == 3 || problem == 4
        n = 10
    else
        n = 5
    end
    i ~ DiscreteUniform(1,n)
    # We create a fixed length list which seems to be much easier
    # than dynamically adding elements to tosses, especially
    # since for problem 2 we always access to tosses[1:5]
    tosses = tzeros(n)

    # We always do the first toss
    i = 1
    tosses[i] ~ Categorical([0.5,0.5])
    while rand(flip(0.8))
        i += 1
        if i > n
            break
        end
        tosses[i] ~ Categorical([0.5,0.5])
    end

    # We must add dummy value (here 3) for the unhandled cases in tosses[1:n]
    for t in 1:n
        if tosses[t] == 0
            tosses[t] ~ Categorical([0,0,1]) # dummy 3
        end
    end

    # * first problem: probability of the sequence [h,h,h] -> 0.015
    # * second problem: probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] -> about 0.03
    # * third problem: probability of the sequences of 10 coin tosses -> about 0.026
    # * fourth problem: expected length of the sequences of coin tosses -> about 4.
    if problem == 1
        t ~ Dirac(tosses[1] == head && tosses[2] == head && tosses[3] == head)
    elseif problem == 2
        true ~ Dirac(tosses[1] == head)
        true ~ Dirac(tosses[2] == head)
        true ~ Dirac(tosses[3] == head)
        t ~ Dirac(tosses[1] == head && tosses[2] == head && tosses[3] == head &&
                  tosses[4] == head && tosses[5] == head)
    elseif problem == 3
        t ~ Dirac(i == 10)
    else
        return i
    end
end

function run_problem(problem)
    println("Problem $problem")
    model = coin_tosses(problem)
    num_chns = 4

    # HH has problem with this!
    # chns = sample(model, MH(), MCMCThreads(), 100_000, num_chns)
    # chns = sample(model, MH(), 100_000)

    # chns = sample(model, PG(15), MCMCThreads(), 1_000, num_chns)

    # chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, SMC(1000), 10_000)

    chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)

    display(chns)
    println("Probability of problem $problem:")

    show_var_dist_pct(chns,:i)
    if problem != 4
        show_var_dist_pct(chns,:t)
    end
end

for problem in 1:4
    run_problem(problem)
    println("\n")
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text: include("FEMShapeFunction.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
module DemoCards

import Base: basename
using Dates
using Mustache
using Literate
using ImageCore
using FileIO, JSON, YAML
using Suppressor # suppress log generated by 3rd party tools, e.g., Literate
import HTTP
using Documenter

const JULIA_COMPAT = let regex=r"julia\s*=\s*\"([\d\.]*)\""
    lines = filter(readlines(normpath(@__DIR__, "..", "Project.toml"))) do line
        occursin(regex, line)
    end
    VersionNumber(match(regex, lines[1]).captures[1])
end
const config_filename = "config.json"
const template_filename = "index.md"
# directly copy these folders without processing
const ignored_dirnames = ["assets"]

function verbose_mode()
    if haskey(ENV, "DOCUMENTER_DEBUG")
        rst = ENV["DOCUMENTER_DEBUG"]
        return lowercase(strip(rst)) == "true"
    else
        return false
    end
end

include("compat.jl")

include("types/card.jl")
include("types/section.jl")
include("types/page.jl")

include("utils.jl")
include("show.jl")

include("Themes/Themes.jl")
using .CardThemes

include("generate.jl")
include("preview.jl")

export makedemos, cardtheme, preview_demos

"""

This package breaks the rendering of the whole demo page into three types:

* `DemoPage` contains serveral `DemoSection`s;
* `DemoSection` contains either serveral `DemoSection`s or serveral `DemoCard`s;
* `DemoCard` consists of cover image, title and other necessary information.
"""
DemoCards

end # module


text:
export fill_bit!

function Base.convert(::Type{ST}, x::T2) where {T1 <: Number, T2 <: Number, ST <: BinarySite{T1}}
    round(ST, T1(x))
end

for ST1 in [:Bit, :Spin, :Half]
    @eval Base.convert(::Type{$ST1{T}}, x::$ST1) where T = $ST1{T}(T(x.value))
    @eval Base.convert(::Type{T1}, x::$ST1{T2}) where {T1, T2} = T1(x.value)
    @eval Base.getindex(::Type{$ST1}, args::T...) where T = $ST1{T}[args...]
end

Base.convert(::Type{ST1}, x::ST2) where {ST1 <: BinarySite, ST2 <: BinarySite} =
    x == up(ST2) ? up(ST1) : down(ST1)

Base.convert(::Type{SA}, x::Int) where {S, T, N, SA <: StaticArray{S, Bit{T}, N}} =
    _convert(Size(SA), eltype(SA), SA, x)

get_bit(::Type{T}, x::Int, i::Int) where T = T((x >> (i - 1)) & 1)

function fill_bit!(A::Array{T}, n::Int) where T
    @inbounds for i in 1:length(A)
        A[i] = get_bit(T, n, i)
    end
    A
end

@generated function _convert(::Size{s}, ::Type{Bit{T}}, ::Type{SA}, x::Int) where {s, T, SA}
    v = [:(get_bit($T, x, $i)) for i = 1:prod(s)]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

# Arrays
(::Type{T})(x::Bit) where {T <: Number} = T(x.value)

for IntType in (:Int8, :Int16, :Int32, :Int64, :Int128, :BigInt)

@eval Base.$IntType(x::Bit) = convert($IntType, x)
@eval Base.convert(::Type{$IntType}, x::Bit{T}) where T = convert($IntType, x.value)

@eval begin
        Base.$IntType(x::AbstractArray{Bit{T}, N}) where {T, N} = convert($IntType, x)

        function Base.convert(::Type{$IntType}, x::AbstractArray{Bit{T}, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(InexactError(:convert, $IntType, x))
            end

            sum(convert($IntType, each) << (i-1) for (i, each) in enumerate(x))
        end

        function Base.convert(::Type{$IntType}, x::AbstractArray{Spin{T}, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(InexactError(:convert, $IntType, x))
            end

            $IntType(sum($IntType(div(each.value+1, 2)) << (i-1) for (i, each) in enumerate(x)))
        end

        function Base.convert(::Type{$IntType}, x::AbstractArray{Half{T}, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(InexactError(:convert, $IntType, x))
            end

            $IntType(sum($IntType(each.value+0.5) << (i-1) for (i, each) in enumerate(x)))
        end
    end
end


text:
const _pickle = PyNULL()

pickle() = _pickle.o == C_NULL ? copy!(_pickle, pyimport(PyReCall.pyversion.major ≥ 3 ? "pickle" : "cPickle")) : _pickle

function Base.serialize(s::AbstractSerializer, pyo::PyObject)
    Base.serialize_type(s, PyObject)
    if pyo.o == C_NULL
        serialize(s, pyo.o)
    else
        b = PyBuffer(PyReCall(pickle()["dumps"], PyObject, pyo))
        serialize(s, unsafe_wrap(Array, Ptr{UInt8}(pointer(b)), sizeof(b)))
    end
end

pybytes(b) = PyObject(@pycheckn ccall(@pysym(PyString_FromStringAndSize),
                                      PyPtr, (Ptr{UInt8}, Int),
                                      b, sizeof(b)))

function Base.deserialize(s::AbstractSerializer, t::Type{PyObject})
    b = deserialize(s)
    if isa(b, PyPtr)
        @assert b == C_NULL
        return PyNULL()
    else
        return PyReCall(pickle()["loads"], PyObject, pybytes(b))
    end
end


text:
@testset "MPPI (PointMass)" begin
    tseed!(1)
    etype = LyceumMuJoCo.PointMass
    env = etype()
    T = 300
    K = 8
    H = 10

    mppi = MPPI(
        env_tconstructor = n -> tconstruct(etype, n),
        covar = Diagonal(0.1^2*I, size(actionspace(env), 1)),
        lambda = 0.01,
        K =  K,
        H = H,
        gamma = 0.99
    )
    f = (a, s, o) -> getaction!(a, s, mppi)
    experiment = ControllerIterator(f, env; T=T, plotiter=T+1)
    for x in experiment
    end
    @test isapprox(getreward(env), 1, atol=0.001)
end


text:
using ChemEquations
using Documenter

DocMeta.setdocmeta!(ChemEquations, :DocTestSetup, :(using ChemEquations))

makedocs(;
    modules=[ChemEquations],
    authors="zlatanvasovic <zlatanvasovic@gmail.com> and contributors",
    repo="https://github.com/zlatanvasovic/ChemEquations.jl/blob/{commit}{path}#L{line}",
    sitename="ChemEquations.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://zlatanvasovic.github.io/ChemEquations.jl",
        assets=String[],
        analytics="UA-33643623-2",
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Library" => [
            "lib/public.md",
            "lib/internals.md",
        ]
    ]
)

deploydocs(;
    repo="github.com/zlatanvasovic/ChemEquations.jl",
)


text:
module ApproximateGPs

using Reexport
@reexport using AbstractGPs
@reexport using GPLikelihoods
using Distributions
using LinearAlgebra
using Statistics
using StatsBase
using FastGaussQuadrature
using SpecialFunctions
using ChainRulesCore
using FillArrays
using PDMats: chol_lower
using IrrationalConstants: sqrt2, invsqrtπ

using AbstractGPs: AbstractGP, FiniteGP, LatentFiniteGP, ApproxPosteriorGP, At_A, diag_At_A

include("utils.jl")

export DefaultQuadrature, Analytic, GaussHermite, MonteCarlo
include("expected_loglik.jl")

export SparseVariationalApproximation, Centered, NonCentered
include("sparse_variational.jl")

using ForwardDiff

export LaplaceApproximation
export build_laplace_objective, build_laplace_objective!
export approx_lml  # TODO move to AbstractGPs, see https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/issues/221
include("laplace.jl")

include("deprecations.jl")

end


text:
include("helper_functions.jl")
using GeometryTypes

abstract type obstacle end

struct line_obstacle <: obstacle
    id::Int64
    p1::Vector{Float64}
    p2::Vector{Float64}
    points::Vector{Vector{Float64}}
    occupied_cells::Vector{Vector{Int64}}
    lines::Vector{LineSegment{Point{2, Float32}}}

    function line_obstacle(id::Int64, p1::Vector{Float64}, p2::Vector{Float64}, params::parameters)
        occupied_cells = intersected_cells([p1, p2], params)
        lines = LineSegment{Point{2, Float32}}[]
        points = Vector{Vector{Float64}}()
        push!(points, p1)
        push!(points, p2)
        push!(lines, LineSegment(Point2f0(p1[1], p1[2]), Point2f0(p2[1], p2[2])))
        new(id, p1, p2, points, occupied_cells, lines)
    end
end

struct convex_polygon <: obstacle
    id::Int64
    # Points are defined such that clockwise connection defines the polygon
    points::Vector{Vector{Float64}}
    occupied_cells::Vector{Vector{Int64}}
    lines::Vector{LineSegment{Point{2, Float32}}}

    function convex_polygon(id::Int64, points::Vector{Vector{Float64}}, params::parameters)
        length(points) < 3 ? throw(DomainError("Polygon must have at least 3 points")) : nothing
        occupied_cells = intersected_cells(points, params)
        lines = LineSegment{Point{2, Float32}}[]
        for i = 1:(length(points)-1)
            push!(lines, LineSegment(Point2f0(points[i][1], points[i][2]), Point2f0(points[i+1][1], points[i+1][2])))
        end
        push!(lines, LineSegment(Point2f0(points[end][1], points[end][2]), Point2f0(points[1][1], points[1][2])))
        new(id, points, occupied_cells, lines)
    end
end

function collision_check(p1::Vector{Float64}, p2::Vector{Float64}, this::T) where {T <: obstacle}
    edge = LineSegment(Point2f0(p1[1], p1[2]), Point2f0(p2[1], p2[2]))
    for line in this.lines
        if intersects(edge, line)[1]
            return true
        end
    end
    return false
end

function plot_obstacle(this::T) where {T <: obstacle}
    x_data = []
    y_data = []
    for pt in this.points
        push!(x_data, pt[1])
        push!(y_data, pt[2])
    end
    push!(x_data, this.points[1][1])
    push!(y_data, this.points[1][2])

    plot!(x_data, y_data, legend=false, color="blue")
end


text:
# Licensed under MIT License, see LICENSE.md
# Copyright 2016 Gandalf Software, Inc., Scott P. Jones

"""
String macro with more Swift-like syntax, plus support for emojis and LaTex names
"""
macro u_str(str) ; s_interp_parse(str) ; end

"""
Interpolates one or more strings using more Swift-like syntax
julia> x = "World"; @sinterpolate "Hello \\(x)"
"Hello World"
"""
macro sinterpolate(args...) ; s_interp_parse(args...) ; end

"""
Handle Unicode character constant, of form \\u{<hexdigits>}
"""
function s_parse_unicode(io, s,  i)
    done(s,i) && throw(ArgumentError("Incomplete \\u{...} in $(repr(s))"))
    c, i = next(s, i)
    c != '{' && throw(ArgumentError("\\u missing opening { in $(repr(s))"))
    done(s,i) && throw(ArgumentError("Incomplete \\u{...} in $(repr(s))"))
    c, i = next(s, i)
    n::UInt32 = 0
    k = 0
    while c != '}'
        done(s, i) && throw(ArgumentError("\\u{ missing closing } in $(repr(s))"))
        (k += 1) > 6 && throw(ArgumentError("Unicode constant too long in $(repr(s))"))
        n = n<<4 + c - ('0' <= c <= '9' ? '0' :
                        'a' <= c <= 'f' ? 'a' - 10 :
                        'A' <= c <= 'F' ? 'A' - 10 :
                        throw(ArgumentError("\\u missing closing } in $(repr(s))")))
        c, i = next(s,i)
    end
    k == 0 && throw(ArgumentError("\\u{} has no hex digits in $(repr(s))"))
    ((0x0d800 <= n <= 0x0dfff) || n > 0x10ffff) &&
    throw(ArgumentError("Invalid Unicode character constant $(repr(s))"))
    print(io, Char(n))
    i
end

"""
Handle Emoji character, of form \\:<name>:
"""
function s_parse_emoji(io, s,  i)
    beg = i-2 # start location
    c, i = next(s, i)
    while c != ':'
        done(s, i) && throw(ArgumentError("\\: missing closing : in $(repr(s))"))
        c, i = next(s, i)
    end
    emojistr = get(Base.REPLCompletions.emoji_symbols, s[beg:i-1], "")
    emojistr == "" && throw(ArgumentError("Invalid Emoji name in $(repr(s))"))
    print(io, emojistr)
    i
end

"""
Handle LaTex character/string, of form \\<name>
"""
function s_parse_latex(io, s,  i)
    beg = i # start location
    c, i = next(s, i)
    while c != '>'
        done(s, i) && throw(ArgumentError("\\< missing closing > in $(repr(s))"))
        c, i = next(s, i)
    end
    latexstr = get(Base.REPLCompletions.latex_symbols, string("\\", s[beg:i-2]), "")
    latexstr == "" && throw(ArgumentError("Invalid LaTex name in $(repr(s))"))
    print(io, latexstr)
    i
end

"""
Handle HTML character/string, of form \\&<name>;
"""
function s_parse_html(io, s,  i)
    beg = i # start location
    c, i = next(s, i)
    while c != ';'
        done(s, i) && throw(ArgumentError("\\& missing ending ; in $(repr(s))"))
        c, i = next(s, i)
    end
    htmlstr = HTMLNames.lookupname(s[beg:i-2])
    htmlstr == "" && throw(ArgumentError("Invalid HTML name in $(repr(s))"))
    print(io, htmlstr)
    i
end

"""
Handle Unicode name, of form \\N{<name>}, from Python
"""
function s_parse_uniname(io, s,  i)
    done(s, i) && throw(ArgumentError("\\N incomplete in $(repr(s))"))
    c, i = next(s, i)
    c != '{' && throw(ArgumentError("\\N missing initial { in $(repr(s))"))
    done(s, i) && throw(ArgumentError("\\N{ incomplete in $(repr(s))"))
    beg = i # start location
    c, i = next(s, i)
    while c != '}'
        done(s, i) && throw(ArgumentError("\\N{ missing closing } in $(repr(s))"))
        c, i = next(s, i)
    end
    unichar = UnicodeNames.lookupchar(uppercase(s[beg:i-2]))
    unichar == typemax(UInt32) && throw(ArgumentError("Invalid Unicode name in $(repr(s))"))
    print(io, Char(unichar))
    i
end

"""
String interpolation parsing
Based on code resurrected from Julia base:
https://github.com/JuliaLang/julia/blob/deab8eabd7089e2699a8f3a9598177b62cbb1733/base/string.jl
"""
function s_print_unescaped(io, s::AbstractString)
    i = start(s)
    while !done(s,i)
        c, i = next(s,i)
        if !done(s,i) && c == '\\'
            c, i = next(s,i)
            if c == 'u'
                i = s_parse_unicode(io, s, i)
            elseif c == ':'	# Emoji
                i = s_parse_emoji(io, s, i)
            elseif c == '&'	# HTML
                i = s_parse_html(io, s, i)
            elseif c == '<'	# LaTex
                i = s_parse_latex(io, s, i)
            elseif c == 'N'	# Unicode name
                i = s_parse_uniname(io, s, i)
            else
                c = (c == '0' ? '\0' :
                     c == '"' ? '"'  :
                     c == '\'' ? '\'' :
                     c == 'a' ? '\a' :
                     c == 'b' ? '\b' :
                     c == 't' ? '\t' :
                     c == 'n' ? '\n' :
                     c == 'v' ? '\v' :
                     c == 'f' ? '\f' :
                     c == 'r' ? '\r' :
                     c == 'e' ? '\e' :
                     throw(ArgumentError("Invalid \\$c sequence in $(repr(s))")))
                write(io, UInt8(c))
            end
        else
            print(io, c)
        end
    end
end

s_unescape_string(s::AbstractString) = sprint(endof(s), s_print_unescaped, s)

function s_print_escaped(io, s::AbstractString, esc::AbstractString)
    i = start(s)
    while !done(s,i)
        c, i = next(s, i)
        c == '\0'       ? print(io, "\\0") :
        c == '\e'       ? print(io, "\\e") :
        c == '\\'       ? print(io, "\\\\") :
        c in esc        ? print(io, '\\', c) :
        '\a' <= c <= '\r' ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
        isprint(c)      ? print(io, c) :
        		  print(io, "\\u{", hex(c), "}")
    end
end

s_escape_string(s::AbstractString) = sprint(endof(s), s_print_escaped, s, "\"")

function s_interp_parse(s::AbstractString, unescape::Function, p::Function)
    sx = []
    i = j = start(s)
    while !done(s, j)
        c, k = next(s, j)
        if c == '\\' && !done(s, k)
            if s[k] == '('
                # Handle interpolation
                if !isempty(s[i:j-1])
                    push!(sx, unescape(s[i:j-1]))
                end
                ex, j = parse(s, k, greedy=false)
                if isa(ex, Expr) && is(ex.head, :continue)
                    throw(ParseError("Incomplete expression"))
                end
                push!(sx, esc(ex))
                i = j
            elseif s[k] == '%'
                # Move past \\, c should point to '%'
                c, k = next(s, k)
                done(s, k) && throw(ParseError("Incomplete % expression"))
                # Handle interpolation
                if !isempty(s[i:j-1])
                    push!(sx, unescape(s[i:j-1]))
                end
                if s[k] == '('
                    # Need to find end to parse to
                    _, j = parse(s, k, greedy=false)
                    # This is a bit hacky, and probably doesn't perform as well as it could,
                    # but it works! Same below.
                    str = "(fmt" * s[k:j-1] * ")"
                else
                    # Move past %, c should point to letter
                    beg = k
                    while true
                        c, k = next(s, k)
                        done(s, k) && throw(ParseError("Incomplete % expression"))
                        s[k] == '(' && break
                    end
                    _, j = parse(s, k, greedy=false)
                    str = string("(cfmt(\"", s[beg-1:k-1], "\",", s[k+1:j-1], ")")
                end
                ex, _ = parse(str, 1, greedy=false)
                if isa(ex, Expr) && is(ex.head, :continue)
                    throw(ParseError("Incomplete expression"))
                end
                push!(sx, esc(ex))
                i = j
            elseif s[k] == '{'
                # Move past \\, c should point to '{'
                c, k = next(s, k)
                done(s, k) && throw(ParseError("Incomplete {...} Python format expression"))
                # Handle interpolation
                if !isempty(s[i:j-1])
                    push!(sx, unescape(s[i:j-1]))
                end
                beg = k # start location
                c, k = next(s, k)
                while c != '}'
                    done(s, k) && throw(ArgumentError("\\{ missing closing } in $(repr(s))"))
                    c, k = next(s, k)
                end
                c != '(' && throw(ParseError("Missing (expr) in Python format expression"))
                # Need to find end to parse to
                _, j = parse(s, k, greedy=false)
                # This is a bit hacky, and probably doesn't perform as well as it could,
                # but it works! Same below.
                str = "(pyfmt" * s[k:j-1] * ")"
                ex, _ = parse(str, 1, greedy=false)
                if isa(ex, Expr) && is(ex.head, :continue)
                    throw(ParseError("Incomplete expression"))
                end
                push!(sx, esc(ex))
                i = j
            else
                j = k
            end
        else
            j = k
        end
    end
    if !isempty(s[i:end])
        push!(sx, unescape(s[i:j-1]))
    end
    length(sx) == 1 && isa(sx[1], ByteString) ? sx[1] : Expr(:call, :sprint, p, sx...)
end

s_interp_parse(s::AbstractString, u::Function) = s_interp_parse(s, u, print)
s_interp_parse(s::AbstractString) =
    s_interp_parse(s, x -> isvalid(UTF8String, s_unescape_string(x))
                   ? s_unescape_string(x) : throw(ArgumentError("Invalid UTF-8 sequence")))


text:
"""
    prob_dde_qs

Delay differential equation model of Quorum Sensing (QS) of Pseudomonas putida IsoF in
continuous cultures.

# References

Buddrus-Schiemann et al. (2014). Analysis of N-Acylhomoserine Lactone Dynamics in Continuous
Cultures of Pseudomonas Putida IsoF By Use of ELISA and UHPLC/qTOF-MS-derived Measurements
and Mathematical Models, Analytical and Bioanalytical Chemistry.
"""
const prob_dde_qs =
  let S₀ = 1, τ = 2, D = 0.1, γₛ = 1.3e-12, Kₘ = 0.38, nₛ = 1.3, a = 0.66, αₐ = 2.3e-19,
    βₐ = 2.3e-18, C₁ = 70e-9, n₁ = 2.3, γₐ = 0.05, αC = 4e4, R = 5e-7, γ₃ = 0.08,
    Kₑ = 1.5e-4, αₗ = 1.1e-8, C₂ = 70e-9, n₂ = 2.5, γₗ = 0.005
    function hill_growth(x, K, n)
      if x > K
        return 1/(1 + (K / x)^n)
      end

      if K == 0
        return one(x)
      end

      y = (x / K)^n
      return y / (1 + y)
    end

    global function f_dde_qs!(du, u, h, p, t)
      if u[1] < 0
        # eq. 1 and 2 not defined for u[1] < 0
        du[1] = 0
        du[2] = 0
      else
        tmp = u[2] * hill_growth(u[1], Kₘ, nₛ)
        du[1] = D * (S₀ - u[1]) - γₛ * tmp
        du[2] = a * tmp - D * u[2]
      end

      if u[4] < 0
        # eq. 3 not defined for u[4] < 0
        du[3] = 0
        du[4] = max(0, αC * (R - u[4]) * u[3] - γ₃ * u[4])
      else
        du[3] = (αₐ + βₐ * hill_growth(u[4], C₁, n₁)) * u[2] -
          (γₐ + D) * u[3] - αC * (R - u[4]) * u[3] + γ₃ * u[4] -
          Kₑ * u[3] * u[5]
        du[4] = αC * (R - u[4]) * u[3] - γ₃ * u[4]
      end

      Cₜ = h(p, t - τ; idxs = 4)
      if Cₜ < 0
        # eq. 5 not defined for Cₜ < 0
        du[5] = 0
      else
        du[5] = αₗ * hill_growth(Cₜ, C₂, n₂) * u[2] - (γₗ + D) * u[5]
      end

      nothing
    end

    # history function for the fourth component
    function h_dde_qs(p, t; idxs = 4)
      idxs == 4 || error("history function is only implemented for the fourth component")
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      7.6e-8
    end

    DDEProblem(f_dde_qs!, [1, 8.4e8, 2.5e-9, 7.6e-8, 5e-15], h_dde_qs, (0.0, 45.0);
               constant_lags = [τ])
  end


text:
# TODO: add a type for kinetic energy

struct Hamiltonian{M<:AbstractMetric, Tlogπ, T∂logπ∂θ}
    metric::M
    ℓπ::Tlogπ
    ∂ℓπ∂θ::T∂logπ∂θ
end
Base.show(io::IO, h::Hamiltonian) = print(io, "Hamiltonian(metric=$(h.metric))")

struct DualValue{V<:AbstractScalarOrVec{<:AbstractFloat}, G<:AbstractVecOrMat{<:AbstractFloat}}
    value::V    # cached value, e.g. logπ(θ)
    gradient::G # cached gradient, e.g. ∇logπ(θ)
    function DualValue(value::V, gradient::G) where {V, G}
        # Check consistence
        if value isa AbstractFloat
            # If `value` is a scalar, `gradient` is a vector
            @assert gradient isa AbstractVector "`typeof(gradient)`: $(typeof(gradient))"
        else
            # If `value` is a vector, `gradient` is a matrix
            @assert gradient isa AbstractMatrix "`typeof(gradient)`: $(typeof(gradient))"
        end
        return new{V,G}(value, gradient)
    end
end

Base.similar(dv::DualValue{<:AbstractVecOrMat{T}}) where {T<:AbstractFloat} =
    DualValue(zeros(T, size(dv.value)...), zeros(T, size(dv.gradient)...))

# `∂H∂θ` now returns `(logprob, -∂ℓπ∂θ)`
function ∂H∂θ(h::Hamiltonian, θ::AbstractVecOrMat)
    res = h.∂ℓπ∂θ(θ)
    return DualValue(res[1], -res[2])
end

∂H∂r(h::Hamiltonian{<:UnitEuclideanMetric}, r::AbstractVecOrMat) = copy(r)
∂H∂r(h::Hamiltonian{<:DiagEuclideanMetric}, r::AbstractVecOrMat) = h.metric.M⁻¹ .* r
∂H∂r(h::Hamiltonian{<:DenseEuclideanMetric}, r::AbstractVecOrMat) = h.metric.M⁻¹ * r

struct PhasePoint{T<:AbstractVecOrMat{<:AbstractFloat}, V<:DualValue}
    θ::T  # Position variables / model parameters.
    r::T  # Momentum variables
    ℓπ::V # Cached neg potential energy for the current θ.
    ℓκ::V # Cached neg kinect energy for the current r.
    function PhasePoint(θ::T, r::T, ℓπ::V, ℓκ::V) where {T, V}
        @argcheck length(θ) == length(r) == length(ℓπ.gradient) == length(ℓπ.gradient)
        if any(isfinite.((θ, r, ℓπ, ℓκ)) .== false)
            @warn "The current proposal will be rejected due to numerical error(s)." isfinite.((θ, r, ℓπ, ℓκ))
            # NOTE eltype has to be inlined to avoid type stability issue; see #267
            ℓπ = DualValue(
                map(v -> isfinite(v) ? v : -eltype(T)(Inf), ℓπ.value),
                ℓπ.gradient
            )
            ℓκ = DualValue(
                map(v -> isfinite(v) ? v : -eltype(T)(Inf), ℓκ.value),
                ℓκ.gradient
            )
        end
        new{T,V}(θ, r, ℓπ, ℓκ)
    end
end

Base.similar(z::PhasePoint{<:AbstractVecOrMat{T}}) where {T<:AbstractFloat} =
    PhasePoint(
        zeros(T, size(z.θ)...),
        zeros(T, size(z.r)...),
        similar(z.ℓπ),
        similar(z.ℓκ),
    )

phasepoint(
    h::Hamiltonian,
    θ::T,
    r::T;
    ℓπ=∂H∂θ(h, θ),
    ℓκ=DualValue(neg_energy(h, r, θ), ∂H∂r(h, r))
) where {T<:AbstractVecOrMat} = PhasePoint(θ, r, ℓπ, ℓκ)

# If position variable and momentum variable are in different containers,
# move the momentum variable to that of the position variable.
# This is needed for AHMC to work with CuArrays and other Arrays (without depending on it).
phasepoint(
    h::Hamiltonian,
    θ::T1,
    _r::T2;
    r=safe_rsimilar(θ, _r),
    ℓπ=∂H∂θ(h, θ),
    ℓκ=DualValue(neg_energy(h, r, θ), ∂H∂r(h, r))
) where {T1<:AbstractVecOrMat,T2<:AbstractVecOrMat} = PhasePoint(θ, r, ℓπ, ℓκ)
# ensures compatibility with ComponentArrays
function safe_rsimilar(θ, _r)
    r = similar(θ)
    copyto!(r, _r)
    r
end

Base.isfinite(v::DualValue) = all(isfinite, v.value) && all(isfinite, v.gradient)
Base.isfinite(v::AbstractVecOrMat) = all(isfinite, v)
Base.isfinite(z::PhasePoint) = isfinite(z.ℓπ) && isfinite(z.ℓκ)

###
### Negative energy (or log probability) functions.
### NOTE: the general form (i.e. non-Euclidean) of K depends on both θ and r.
###

neg_energy(z::PhasePoint) = z.ℓπ.value + z.ℓκ.value

neg_energy(h::Hamiltonian, θ::AbstractVecOrMat) = h.ℓπ(θ)

neg_energy(
    h::Hamiltonian{<:UnitEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVector} = -sum(abs2, r) / 2

neg_energy(
    h::Hamiltonian{<:UnitEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractMatrix} = -vec(sum(abs2, r; dims=1)) / 2

neg_energy(
    h::Hamiltonian{<:DiagEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVector} = -sum(abs2.(r) .* h.metric.M⁻¹) / 2

neg_energy(
    h::Hamiltonian{<:DiagEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractMatrix} = -vec(sum(abs2.(r) .* h.metric.M⁻¹; dims=1) ) / 2

function neg_energy(
    h::Hamiltonian{<:DenseEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVecOrMat}
    mul!(h.metric._temp, h.metric.M⁻¹, r)
    return -dot(r, h.metric._temp) / 2
end

energy(args...) = -neg_energy(args...)

####
#### Momentum refreshment
####

phasepoint(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    θ::AbstractVecOrMat{T},
    h::Hamiltonian
) where {T<:Real} = phasepoint(h, θ, rand(rng, h.metric))

abstract type AbstractMomentumRefreshment end

"Completly resample new momentum."
struct FullMomentumRefreshment <: AbstractMomentumRefreshment end

refresh(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    ::FullMomentumRefreshment,
    h::Hamiltonian,
    z::PhasePoint,
) = phasepoint(h, z.θ, rand(rng, h.metric))

"""
Partial momentum refreshment with refresh rate `α`.

See equation (5.19) [1]

    r' = α⋅r + sqrt(1-α²)⋅G

where r is the momentum and G is a Gaussian random variable.

## References

1. Neal, Radford M. "MCMC using Hamiltonian dynamics." Handbook of markov chain monte carlo 2.11 (2011): 2.
"""
struct PartialMomentumRefreshment{F<:AbstractFloat} <: AbstractMomentumRefreshment
    α::F
end

refresh(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    ref::PartialMomentumRefreshment,
    h::Hamiltonian,
    z::PhasePoint,
) = phasepoint(h, z.θ, ref.α * z.r + sqrt(1 - ref.α^2) * rand(rng, h.metric))


text:
using Pkg
Pkg.add("BenchmarkTools")
Pkg.activate(joinpath(@__DIR__, ".."))
#Pkg.add("https://github.com/blenessy/AMCL.jl")

using BenchmarkTools
using MilagroCrypto

curves = isempty(Base.ARGS) ? ["BLS381"] : copy(Base.ARGS)

rng = csprng(octet("insecure"))
msg = octet(ones(UInt8, 200)) # constant 200 byte message

for c in curves
    println("== BLS_$(c)_KEY_PAIR_GENERATE ==")
    sk1 = octet(eval(Meta.parse("AMCL.BGS_$(c)")))
    sk2 = octet(eval(Meta.parse("AMCL.BGS_$(c)")))
    pk1 = octet(16 * eval(Meta.parse("AMCL.BFS_$(c)"))) # should fit BLS48
    pk2 = octet(16 * eval(Meta.parse("AMCL.BFS_$(c)"))) # should fit BLS48
    api = eval(Meta.parse("AMCL.BLS_$(c)_KEY_PAIR_GENERATE"))
    @assert api(rng, sk1, pk1) == AMCL.BLS_OK
    @btime $api(rng, $sk2, $pk2);

    println("== BLS_$(c)_SIGN ==")
    sig1 = octet(eval(Meta.parse("AMCL.BFS_$(c)")) + 1)
    sig2 = octet(eval(Meta.parse("AMCL.BFS_$(c)")) + 1)
    api = eval(Meta.parse("AMCL.BLS_$(c)_SIGN"))
    @assert api(sig1, msg, sk1) == AMCL.BLS_OK
    @btime $api($sig2, msg, $sk2);

    println("== BLS_$(c)_ADD_G1 ==")
    sig = octet(eval(Meta.parse("AMCL.BFS_$(c)")) + 1)
    api = eval(Meta.parse("AMCL.BLS_$(c)_ADD_G1"))
    @assert api(sig1, sig2, sig) == AMCL.BLS_OK
    @btime $api($sig1, $sig2, $sig)

    println("== BLS_$(c)_ADD_G2 ==")
    pk = octet(16 * eval(Meta.parse("AMCL.BFS_$(c)"))) # should fit BLS48
    api = eval(Meta.parse("AMCL.BLS_$(c)_ADD_G2"))
    @assert api(pk1, pk2, pk) == AMCL.BLS_OK
    @btime $api($pk1, $pk2, $pk)

    println("== BLS_$(c)_VERIFY ==")
    api = eval(Meta.parse("AMCL.BLS_$(c)_VERIFY"))
    @assert api(sig, msg, pk) == AMCL.BLS_OK
    @btime $api($sig, msg, $pk);
end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние прямоугольные перегородки (все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)
#РЕЗУЛЬТАТ: Робот - в исходном положении, и в 4-х приграничных клетках, две из которых имеют ту же широту, а две - ту же долготу, что и Робот, стоят маркеры.
function mark_angles(r)
    for i = 0:3
        moves(r,i)
    end
end
function moves(r,HS)
    y = 0
    flag = true
    while isborder(r,HorizonSide(HS)) || flag==0
        if isborder(r,HorizonSide(HS))==0
            move!(r,HorizonSide(HS))
            y += 1
        end
        if isborder(r,HorizonSide(HS))
            flag, y = round(r,HS,y)
        end
        if flag == false
            break
        end
    end
    putmarker!(r)
    return(r,HS,y)
end
function round(r,HS,y)
    x = 0
    flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 1) % 4))==1
        move!(r,HorizonSide((HS + 1) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))
        flag = false
    end
    if  isborder(r,HorizonSide(HS))==0
        move!(r,HorizonSide(HS))
        y += 1
    end
    while isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))==0
        move!(r,HorizonSide(HS))
        y += 1
    end
    for i = 1:x
        move!(r,HorizonSide((HS+3) % 4))
    end
    return flag,y
end
function return(r,HS,y)
    while y > 0
        if isborder(r,HorizonSide((HS + 2) % 4))==0
            move!(r,HorizonSide(((HS + 2) % 4)))
        else
            x = 0
            while isborder(r,HorizonSide((HS + 2) % 4)) && !isborder(r,HorizonSide((HS + 1) % 4))==1
                move!(r,HorizonSide((HS + 1) % 4))
                x += 1
            end
            move!(r,HorizonSide((HS + 2) % 4))
            while isborder(r,HorizonSide((HS + 3) % 4))==1
                move!(r,HorizonSide((HS + 2) % 4))
                y -= 1
            end
            for i = 1:x
                move!(r,HorizonSide((HS + 3) % 4))
            end
        end
        y -= 1
    end
end
function all(r)
    for i = 0:3
        moves(r,i)
    end
end
mark_angles(r)