text:
using JuMP, MosekTools, COSMO
using SketchyCGAL
using LinearAlgebra, SparseArrays
using BSON

include("utils.jl")

G = graph_from_file(joinpath(@__DIR__, "data/gset/G1"))
n = size(G, 1)

function solve_with_jump(C; optimizer=Mosek.Optimizer)
    sdp = Model(optimizer)
    @variable(sdp, X[1:n, 1:n] in PSDCone())
    @constraint(sdp, diag(X) .== 1)
    @objective(sdp, Min, sum(C.*X))
    optimize!(sdp)
    Xopt = value.(X)

    # Optionally, we can apply a correction to enforce feasibility
    # Xopt[diagind(X)] .= ones(n)
    return Xopt
end

# Builds the LMI of the dual problem
function build_affine_exp(C, y)
    n = size(C, 1)
    aff_exp = spzeros(GenericAffExpr{Float64, VariableRef}, n, n)

    # This is a bit of a hack
    for (j, k, v) ∈ zip(findnz(C)...)
        aff_exp[j,k] = 1
    end

    for i in 1:n
        add_to_expression!(aff_exp[i,i], -y[i])
    end
    for (j, k, v) ∈ zip(findnz(C)...)
        add_to_expression!(aff_exp[j,k], C[j,k])
    end

    # hack part 2
    for (j, k, v) ∈ zip(findnz(C)...)
        aff_exp[j,k] -= 1
    end

    return aff_exp
end

function solve_dual_with_jump(C; optimizer=Mosek.Optimizer)
    sdp = Model(optimizer)
    @variable(sdp, y[1:n])
    @objective(sdp, Max, sum(y))
    A = Symmetric(build_affine_exp(C, y))
    psd_constraint = @constraint(sdp, A in PSDCone())
    optimize!(sdp)
    Xopt = value.(dual(psd_constraint))
    return Xopt
end

## ---------- Problem setup ----------
# Parameters
n = size(G, 1)
d = n

# Data
# objective = min -1/4( ∑_ij w_ij*(1-yi*yj) ) = -1/4⟨diag(∑ᵢ w_ij) - W), Y⟩
# note that this reformulation uses the fact that Tr(Y) = 1
C = -0.25*(Diagonal(G*ones(n)) - G)
b = ones(n)

# Scaling variables -- so trace is bounded by 1
scale_C = 1 / norm(C)
scale_X = 1 / n

# Linear map
# AX = diag(X)
function A!(y, X)
    @views y .= X[diagind(X)]
    return nothing
end

# Adjoint
# A*z = Diagonal(z)
function A_adj!(S::SparseMatrixCSC, z)
    @views S[diagind(S)] .+= z
    return nothing
 end

# primative 1: u -> C*u
function C_u!(v, u)
    mul!(v, C_const, u)
end

#primative 2: (u, z) -> (A*z)u
function Aadj_zu!(v, u, z)
    v .= u .* z
end

#primative 3: u -> A(uu^T)
function A_uu!(z, u)
    z .= u.*u
end

## Solve with JuMP + COSMO
Xopt = solve_with_jump(C)

# We can also try the dual version of the problem + chordal decomposition,
#   but in this case, this approach seems much slower
# Xopt = solve_dual_with_jump(C;
#     optimizer=with_optimizer(
#         COSMO.Optimizer,
#         decompose = true,
#         merge_strategy = COSMO.CliqueGraphMerge
# ))

# Uncomment the line below to save the optimal soln
# BSON.bson("G67_Xopt", Dict("Xopt" => Xopt))

pstar_jump = sum(C .* Xopt)
@show pstar_jump

## Solve with CGAL
@time XT, yT = SketchyCGAL.cgal_full(
    C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
    max_iters=250,
    print_iter=25
)

pstar_cgal = sum(C .* XT * 1/scale_X)
@show pstar_cgal

## Solve with SketchyCGAL
R = 20
ηt(t) = 2.0/(t + 1.0)
δt(t) = 1.0
tt = @timed scgal_full(
    C, b, A!, A_adj!, A_uu!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
    max_iters=1_000,
    print_iter=100,
    R=R,
    # logging=true,
    # logging_primal=true,
    ηt=ηt,
    # δt=δt
)

soln = tt.value
trial_time = tt.time
trial_alloc = tt.bytes/1e6 #(MB)

Xhat = SketchyCGAL.construct_Xhat(soln)
pstar_scgal = sum(C .* Xhat) * 1/scale_X
@show pstar_scgal

## Compute true objective (MAXCUT after rounding scheme)
cache = zeros(n, R)
@time true_obj_scgal = SketchyCGAL.compute_objective(C, soln.UT, soln.ΛT, cache=cache)# * 1/scale_X
@show true_obj_scgal

## Compute primal infeasibility
@time p_infeas_scgal = SketchyCGAL.compute_primal_infeas_mc(soln.UT, soln.ΛT, cache=cache)
@show p_infeas_scgal


text:
using Test
using MPI
using Random

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
sz = MPI.Comm_size(comm)
filename = MPI.bcast(tempname(), 0, comm)

MPI.Barrier(comm)

# Collective write
fh = MPI.File.open(comm, filename, read=true, write=true, create=true)
MPI.File.set_view!(fh, 0, MPI.Datatype(Int64), MPI.Datatype(Int64))
MPI.File.write_at_all(fh, rank*2, ArrayType([Int64(rank+1) for i = 1:2]))

MPI.File.sync(fh)

# Noncollective read
if rank == 0
    data = zeros(Int64, (2,sz))
    MPI.File.read_at!(fh, 0, data)
    @test data == [j for i = 1:2, j=1:sz]
end

MPI.File.sync(fh)
MPI.Barrier(comm)

if rank == sz-1
    MPI.File.write_at(fh, 0, ArrayType([Int64(-1) for i = 1:2]))
end

MPI.File.sync(fh)

# Collective read
data = zeros(Int64, 1)
MPI.File.read_at_all!(fh, rank*2, data)
@test data == [rank == 0 ? -1 : rank+1]
close(fh)

ext: jl


text:
# src/qccdevcontrol.jl
# Created by Anabel Ovide and Alejandro Villoria, 23 April, 2021
# MIT license
# Sub-module QCCDevCtrl

module QCCDev_Feasible
export load_checks, OperationNotAllowedException, isallowed_load, isallowed_swap
export  isallowed_linear_transport

using ..QCCDevDes_Types
using ..QCCDevControl_Types
using ..QCCDev_Utils

"""
Default error message for QCCD operations.
"""
struct OperationNotAllowedException <: Exception
  msg ::String
end

"Nicer way to throw error"
opError(x) = throw(OperationNotAllowedException(x))

"""
Function `time_check()` — checks if given time is correct

# Arguments
* `qdc:: Time_t` — Actual qdc device's time.
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.

The function throws an error if time is not correct.
"""
function _time_check(t_qdc:: Time_t, t::Time_t, id::Symbol)
  0 ≤ t_qdc ≤ t  || opError("Time must be higher than $t_qdc")
  haskey(OperationTimes, id) || opError("Time model for $id not defined.")
end

"""
Function `isallowed_load()` — checks if load operation is posisble

# Arguments
* `qdc::QCCDevControl` — Actual device's status.
* `loading_zone::Symbol` — Desired place to load an ion
* `t::Time_t` — Time at which the operation commences.  Must be no earlier than the latest time
                given to previous function calls.
# checks
* Check time — Call _time_check function.
* Check maximum capacity — Check device's maximum capacity is not exeeded.
* Check trap exist — Check current loading_zone exist.
* Check hole is avialable — Checl loading hole is not busy.
"""
function isallowed_load(qdc::QCCDevControl, loading_zone::Symbol, t::Time_t)
    _time_check(qdc.t_now, t, :load)
    haskey(qdc.loadingZones, loading_zone) || opError("Loading zone with id $loading_zone doesn't exist.")
    qdc.loadingZones[loading_zone].hole != nothing && opError("Loading hole is busy.")
end

"""
Function `isallowed_swap()` — checks if load operation is posisble

# Arguments
* `qdc::QCCDevControl` — Actual device's status.
* `ion𝑖_idx`, 𝑖=1,2, the (1-based) indices of the two ions.  Must be in the same gate zone.
* `t::Time_t` — Time at which the operation commences.  Must be no earlier than the latest time
                given to previous function calls.
# checks
* Check time — Call _time_check function.
* Check if two ions exists
* Check if two ions are in the same zone
* check if ions are in same chain and are adjacents
* Check if chain is in a gate zone
"""
function isallowed_swap(qdc::QCCDevControl, ion1_idx:: Int, ion2_idx:: Int , t::Time_t)
    _time_check(qdc.t_now, t, :load)
    haskey(qdc.qubits, ion1_idx) || opError("Qubit with id $ion1_idx doesn't exist.")
    haskey(qdc.qubits, ion2_idx) || opError("Qubit with id $ion2_idx doesn't exist.")
    qdc.qubits[ion1_idx].position == qdc.qubits[ion2_idx].position ||
                    opError("Qubits with ids $ion1_idx  and $ion2_idx are not in the same zone.")

    zone = giveZone(qdc, qdc.qubits[ion1_idx].position)
    zone.zoneType == :gateZone || opError("Swap can only be done in Gate Zones.")

    pos1 = map( y -> findall(x->x==ion1_idx, y), zone.chain)
    pos2 = map( y -> findall(x->x==ion2_idx, y), zone.chain)
    check = x ->  isempty(pos1[x]) && isempty(pos2[x]) || !isempty(pos1[x]) && !isempty(pos2[x]) ||
                  opError("Qubits with ids $ion1_idx and $ion2_idx are not in the same chain.")
    map(x -> check(x) , 1:length(pos1))

    pos1 = collect(Iterators.flatten(pos1))[1]
    pos2 = collect(Iterators.flatten(pos2))[1]
    pos1 == pos2 + 1 || pos1 == pos2 - 1 ||
                    opError("Qubits with ids $ion1_idx and $ion2_idx are not adjacents.")
end

"""
Function `isallowed_linear_transport()` — checks feasibility of `linear_transport`.
# Arguments: See `linear_transport`
# checks
* Check time — Call _time_check function.
* Check OperationTimes - Checks if time model is defined for `linear_transport`
* Check ion_idx - Checks if ion is in the device.
* Check destination_idx - Cehcks if destination zone is in the device.
* Check ion position - Check if ion's position exists in the device.
* Check ion chain - Checks if the ion is in the correct chain position to leave the device.
* Check ends - Checks if ion position and destination are adjacent.
* Check capacity - Checks if destination zone is not currenly full.
"""
function isallowed_linear_transport(qdc           :: QCCDevControl,
                                    t             :: Time_t,
                                    ion_idx       :: Int,
                                    destination_idx      :: Symbol)

  _time_check(qdc.t_now, t, :linear_transport)

  if ion_idx ∉ keys(qdc.qubits)
    opError("Ion with ID $ion_idx is not in device")
  end

  destination = giveZone(qdc, destination_idx)
  if isnothing(destination)
    opError("Zone with ID $destination_idx is not in device")
  end

  currentPosition = giveZone(qdc, qdc.qubits[ion_idx].position)
  if isnothing(currentPosition)
    opError("Ion with ID $ion_idx is nowhere (?)")
  end

  chain = nothing
  if currentPosition.end0 == destination.id

    chain = currentPosition.zoneType == :loadingZone ?
        currentPosition.hole : first(currentPosition.chain)
  elseif currentPosition.end1 == destination.id

    chain = currentPosition.zoneType == :loadingZone ?
      currentPosition.hole : last(currentPosition.chain)
  else
    opError("Can't do linear transport to a non-adjacent zone.")
  end

  if chain != ion_idx && chain != [ion_idx]
    opError("Ion $ion_idx can't leave the trap since it's not in the correct end position.")
  end

  destinationIsLoadingZone = destination.zoneType == :loadingZone
  if (destinationIsLoadingZone && !isnothing(destination.hole)) ||
    (!destinationIsLoadingZone && !isempty(destination.chain) &&
    sum(length, destination.chain) == destination.capacity)
    opError("Destination zone with ID $destination_idx cannot hold more ions.")
  end
end

end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
Base.@kwdef struct CTSimulation <: AbstractSimulation
    nphotons::Int = 10000
end

"""
    simulate(sim::CTSimulation, sinog::CTSinogram; <keyword arguments>)

Simulate a CT scan by generating `sim.nphotons` Poisson-distributed photons for
each angle, which are then detected with a probability given by
``p_{ϕ_j}(x_i) = exp(-s_{i,j})``, where ``s_{i,j}`` is `sinog[i,j]`.

The keyword arguments are the same as those for [`simulate_ct`](@ref).

See also: [`simulate_ct`](@ref)
"""
@inline function simulate(sim::CTSimulation, sinog::CTSinogram; kwargs...)
    simulate_ct(sinog; sim.nphotons, kwargs...)
end

"""
    StatsBase.sample(data::CTSinogram, xs::AbstractVector[; nsamples=1000, nblks=1, nbins=nothing])

Compute `nsamples` from `xs` according to the distribution given by each column
of `data`. Optionally, if `nblks > 1`, `nblks × nsamples` samples are computed.
The sampled data are collected into a histogram of length `nbins` for each
column in `data`. The resulting data have size `(nbins, size(data, 2), nblks)`.
"""
function StatsBase.sample(
    data::CTSinogram,
    xs::AbstractVector;
    nsamples::Integer = 1000,
    nblks::Integer = 1,
    nbins::Optional{Integer} = nothing,
    progress::Bool = false,
)
    nd, nϕ = size(data)
    @assert length(xs) == nd "Dimension mismatch: `xs` vector should match first dimension of `data`"
    nbins = maybe(nd, nbins)
    ζ = half(xs)
    bins = linspace(-ζ..ζ, nbins+1)
	sampled = similar(data, (nbins, nϕ, nblks))
    pxs = fill(similar(sampled, nsamples), 1)
    Threads.resize_nthreads!(pxs)
    p = Progress(
        size(sampled, 3); dt=0.2, desc="Computing samples...", enabled=progress)
	Threads.@threads for k ∈ axes(sampled, 3)
        id = Threads.threadid()
        px = pxs[id]
        @inbounds for j ∈ axes(sampled, 2)
            ws = StatsBase.weights(view(data, :, j))
            StatsBase.sample!(xs, ws, px)
            h = normalize(
                StatsBase.fit(StatsBase.Histogram, px, bins);
                mode = :probability,
            )
            sampled[:,j,k] .= h.weights
        end
        next!(p)
    end
	sampled
end

@inline function StatsBase.sample(
    data::CTSinogram,
    xsi::Interval;
    kwargs...
)
    xs = linspace(xsi, size(data, 1))
    sample(data, xs; kwargs...)
end

"""generate_photons(n::Integer, nx::Integer, nϕ::Integer)

Generate ``⟨n⟩`` Poisson distributed photons per projection
angle over an array of `nx` detectors.

Return a `nd × nϕ` matrix with the generated photons.
```
"""
@inline function generate_photons(n::Integer, nx::Integer, nϕ::Integer)
    d = Distributions.Poisson(n)
    emit = Distributions.DiscreteUniform(1,nx)
    rngs_gen = [Random.MersenneTwister() for _ ∈ 1:Threads.nthreads()]
    rngs_measure = [Random.MersenneTwister() for _ ∈ 1:Threads.nthreads()]
    photons = zeros(Int, nx, nϕ)
    p = Progress(nϕ, 0.2, "Generating photons: ")
    Threads.@threads for i ∈ axes(photons, 2)
        id = Threads.threadid()
        @inbounds g = rngs_gen[id]
        @inbounds m = rngs_measure[id]
        num = rand(g, d)
        foreach(1:num) do _
            x = rand(m, emit)
            @inbounds photons[x, i] += 1
        end
        next!(p)
    end
    photons
end

"""simulate_ct(sinog::AbstractMatrix{T}; <keyword arguments>) where {T <: Real}

Simulate a low dose CT scan. This samples `sinog` with ``⟨n⟩``
random photons per projection angle.

# Arguments
- `sinog`: sinogram data.
- `nphotons::Integer=10000`: mean number of photons.
- `ϵ::Real=1`: detectors quantum efficiency.
- `take_log::Bool=false`: whether to take the logarithm of the
  resampled intensities to obtain the corresponding
  sinogram.
"""
function simulate_ct(
    sinog::AbstractMatrix{T};
    nphotons::Integer = 10000,
    ϵ::Real = 1,
    take_log::Bool = false,
) where {T <: Real}
    nd, nϕ = size(sinog)
    photons = generate_photons(nphotons, nd, nϕ)
    rngs_absorp = [Random.MersenneTwister() for _ ∈ 1:Threads.nthreads()]
    rngs_clicks = [Random.MersenneTwister() for _ ∈ 1:Threads.nthreads()]
    measure = Distributions.Uniform()
    detector = Distributions.Uniform()
    low_sample = similar(sinog)
    p = Progress(length(low_sample), 0.2, "Tracking photons: ")
    Threads.@threads for k ∈ LinearIndices(low_sample)
        id = Threads.threadid()
        @inbounds a = rngs_absorp[id]
        @inbounds c = rngs_clicks[id]
        @inbounds u = exp(-sinog[k])
        @inbounds n = photons[k]
        @inbounds low_sample[k] = count(1:n) do _
            rand(a, measure) <= u && rand(c, detector) <= ϵ
        end
        next!(p)
    end
    max_photons = maximum(low_sample)
    if take_log
        smax = maximum(sinog)
        @. low_sample = -log(low_sample / max_photons)
    else
        @. low_sample = low_sample / max_photons
    end
    low_sample
end

simulate_ct(nphotons::Integer; kwargs...) =
    x -> simulate_ct(x; nphotons, kwargs...)
simulate_ct(; kwargs...) = x -> simulate_ct(x; kwargs...)

"""
    sample_sinogram_external(sinog::AbstractMatrix{T}; <keyword arguments>) where {T <: Real}

Simulate a low dose CT scan. This Resample `sinog` with `n` photons.

# Arguments
- `sinog`: sinogram data.
- `sinog_path='_tmp_sinog.dat'`: path to a file where to write the sinogram to.
- `resampled_path='_tmp_resampled.dat'`: path to a file where to read the resampled sinogram from.
- `n=10000`: mean number of photons.
- `ϵ=1`: detectors quantum efficiency.
- `take_log=true`: whether to take the logarithm of the resampled intensities to obtain the corresponding sinogram.
- `verbosity=0`: set verbosity level.
- `progress=false`: show progress bar.
- `options=[]`: Additional options to be passed to the program.
"""
function sample_sinogram_external(
    sinog::AbstractMatrix{T};
    sinog_path = "_tmp_sinog.dat",
    resampled_path = "_tmp_resampled.dat",
    log_path = "_tmp_log.txt",
    n = 10000,
    ϵ = 1,
    eps = nothing,
    take_log = true,
    verbosity = nothing,
    progress = false,
    options = [],
) where {T<:Real}
    ϵ = isnothing(eps) ? ϵ : eps
    tmp_path = mk_temp_dir()
    sinog_path = joinpath(tmp_path, sinog_path)
    resampled_path = joinpath(tmp_path, resampled_path)
    log_path = joinpath(tmp_path, log_path)
    vlevel = verbosity_level(verbosity)
    write_ct_image(sinog_path, sinog)
    nd, nϕ = size(sinog)
    cmd_args = []
    if !isnothing(vlevel)
        append!(cmd_args, ["$vlevel"])
    end
    if !progress
        append!(cmd_args, ["--xpbar"])
    end
    append!(cmd_args, options)
    cmd = @cmd "ctrun resample -ILT1 -D \$nd -P \$nϕ -n \$n -e \$ϵ \$cmd_args"
    run(pipeline(
        cmd,
        stdin = sinog_path,
        stdout = resampled_path,
        stderr = log_path,
    ))
    read_ct_image(resampled_path; rows = nd, cols = nϕ)
end

@deprecate resample_sinogram sample_sinogram_external false


text:
using ImageMagick
using CairoMakie, AbstractPlotting, MakieGallery

tmp = joinpath(@__DIR__, "test.png")
AbstractPlotting.current_backend[] = CairoMakie.CairoBackend(tmp)
database = MakieGallery.load_database()
filter!(database) do entry
    "2d" in entry.tags &&
    "Text rotation" != entry.title &&
    "fem polygon 2d" != lowercase(entry.title) &&
    "Hbox" != entry.title &&
    !("heatmap" in entry.tags) && # why though, they worked -.-
    !("image" in entry.tags)
end

ref_path = MakieGallery.download_reference(v"0.0.9")

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)
rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)

MakieGallery.record_examples(test_record_path)
MakieGallery.run_comparison(test_record_path, ref_path, tested_diff_path)
rm(tmp, force = true)


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6710a438-60b3-11eb-3f4d-1f20bf10632b
begin
	import Pkg
	Pkg.activate("..")
	using Plots   # pour les graphiques
	using PlutoUI # pour la partie interactive
end

# ╔═╡ b8c02e5c-60ae-11eb-3406-71a8af749eb1
md"""
# Lois Continues
"""

# ╔═╡ e64dfd22-60ae-11eb-0a06-05c3bd1198ff
md"""
# Lois continues

Afin de faciliter la recherche et compréhension des concepts, je vais utiliser parfois des abbreviations en anglais. Je veux simplement présenter des différentes loix, enlister leurs paramètres et essayer de les expliquer tout d'abord pour confirmer ma propre compréhension de celles-ci et puis pour faciliter l'apprentissage de ceux qui me suivront.
"""

# ╔═╡ 8e20343c-657f-11eb-07c6-c90cd368c3c2
md"""
* poisson
* bernoulli / geometric ?
* exponential
* student
* F
* X^2
"""

# ╔═╡ 090dfb62-658e-11eb-1635-1d021fed9d2e
plotly()

# ╔═╡ e4886b96-60af-11eb-3446-3781f9e0a830
md"""
## Loi Gaussienne (Normale)

La loi des lois.

\begin{equation}
\mathcal{N}(\mu, \sigma) = f(x)
= \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x - \mu}{\sigma})^2}
\end{equation}
"""

# ╔═╡ 87df4316-624f-11eb-3c51-4719053a20a8
function Normale(x, μ, σ)
	exp(-0.5 * ((x - μ)/σ)^2) / (σ*(2*π)^2)
end

# ╔═╡ dd4d901a-6582-11eb-2060-61d87809df7f
md"""
$(@bind go_normal Button("Reset"))
"""

# ╔═╡ eabbaf28-6250-11eb-086b-5bd8739b310d
let
	go_normal
	md"""
	μ $(@bind μ Slider(-10:0.1:10, default=0, show_value=true))

	σ $(@bind σ Slider(0.1:0.1:10, default=1, show_value=true))

	lower $(@bind n_lower NumberField(-3:0, default=(μ - 5*σ))) upper $(@bind n_upper
	NumberField(1:4, default=(μ + 5*σ)))
	"""
end

# ╔═╡ 948f84e0-6583-11eb-0a94-b7cff41c8a8d
x = n_lower:0.1:n_upper;

# ╔═╡ 6c8090ac-6583-11eb-1e5d-5fc9d3bcfe76
begin
	n_pdf = plot(x, Normale.(x, μ, σ), legend=false)
	title!("Densité de Probabilité")
end

# ╔═╡ 19afa73c-6588-11eb-2487-9da0cdeb4612
begin
	plot(x, cumsum(Normale.(x, μ, σ)), legend=false)
	title!("Fonction de répartition")
end

# ╔═╡ d09ee836-658d-11eb-082c-e79a6e76d453

# ╔═╡ Cell order:
# ╟─b8c02e5c-60ae-11eb-3406-71a8af749eb1
# ╟─e64dfd22-60ae-11eb-0a06-05c3bd1198ff
# ╠═8e20343c-657f-11eb-07c6-c90cd368c3c2
# ╠═6710a438-60b3-11eb-3f4d-1f20bf10632b
# ╠═090dfb62-658e-11eb-1635-1d021fed9d2e
# ╟─e4886b96-60af-11eb-3446-3781f9e0a830
# ╠═87df4316-624f-11eb-3c51-4719053a20a8
# ╟─6c8090ac-6583-11eb-1e5d-5fc9d3bcfe76
# ╟─948f84e0-6583-11eb-0a94-b7cff41c8a8d
# ╟─dd4d901a-6582-11eb-2060-61d87809df7f
# ╟─eabbaf28-6250-11eb-086b-5bd8739b310d
# ╟─19afa73c-6588-11eb-2487-9da0cdeb4612
# ╠═d09ee836-658d-11eb-082c-e79a6e76d453


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

JLLWrappers.@generate_main_file_header("CHiPPS_BLIS")
JLLWrappers.@generate_main_file("CHiPPS_BLIS", UUID("33c101cd-e3cd-5787-bb96-d137b8b1ecc4"))
end  # module CHiPPS_BLIS_jll


text:
module VI
export value_iteration, bellman_backup

function state_backup(sas, T, C, V_old)
    return get(T, sas, 0.)*(get(C, sas, 0.) + V_old[sas[end]])
end

function bellman_backup(s, S, A, T, C, V_old)
    Q = ones(length(A)).*typemax(Float64)
    for a in  eachindex(A)
        for s_next in eachindex(S)
            # Intutively, (sparse) matrix operations will be a lot faster than dict access
            if (s, a, s_next) in keys(T)
                if Q[a]==typemax(Float64)
                    Q[a]=0.
                end
                Q[a] += state_backup((s, a, s_next), T, C, V_old)
            end
        end
    end
    return minimum(Q)
end

function value_iteration(S::AbstractArray{Int64,1},
    A::AbstractArray{Int64,1},
    T,C; ϵ::Float64=0.001, V_old = zeros(length(S)))
    V_new = copy(V_old)
    residuals = ones(length(S)).*typemax(Float64)
    n = 0
    while maximum(residuals)>ϵ
        n+=1
        for s in eachindex(S)
            V_new[s] = bellman_backup(s, S, A, T, C, V_old)
        end
        residuals = abs.(V_new.-V_old)
        V_old, V_new = V_new, V_old
    end
    # println("Finished after $n iterations")
    return V_new
end
end


text:
function resnet50(;pretrained = true, keephead = false)
    model = Metalhead.resnet50()
    if !pretrained
        return model
    end

    weights = BSON.load(datadep"weights-resnet50/resnet50.bson")[:weights]
    Flux.loadparams!(model, weights)

    if !keephead
        model = model[1:end-3]
    end

    return model
end


text:
#math_graphics.jl:
#-------------------------------------------------------------------------------

#==Aliases
===============================================================================#
#Real data type used for display functionnality (low-level):
const DReal = Float64

#==Constants
===============================================================================#
const DInf = convert(DReal, Inf)
const DNaN = convert(DReal, NaN)

const COLOR_TRANSPARENT = ARGB32(0, 0, 0, 0)

const COLOR_BLACK = RGB24(0, 0, 0)
const COLOR_WHITE = RGB24(1, 1, 1)

const COLOR_RED = RGB24(1, 0, 0)
const COLOR_GREEN = RGB24(0, 1, 0)
const COLOR_BLUE = RGB24(0, 0, 1)

#==Data structures.
===============================================================================#
#TODO: Use library version?
abstract type Point end #So we can use as a function
struct Point2D <: Point
	x::DReal
	y::DReal
end

abstract type DirectionalVector end
struct Vector2D <: DirectionalVector
	x::DReal
	y::DReal
end
Vector2D(p::Point2D) = Vector2D(p.x, p.y)
Point2D(p::Vector2D) = Point2D(p.x, p.y)

#Plot extents along one dimension:
struct PExtents1D
	min::DReal
	max::DReal
end
PExtents1D(;min::Real=DNaN, max::Real=DNaN) = PExtents1D(min, max)

#Could use BoundingBox, but cannot control type (and is technically not a BB).
#TODO: PExtents2D from PExtents1D?
struct PExtents2D
	xmin::DReal
	xmax::DReal
	ymin::DReal
	ymax::DReal
end
PExtents2D(;xmin::Real=DNaN, xmax::Real=DNaN, ymin::Real=DNaN, ymax::Real=DNaN) =
	PExtents2D(xmin, xmax, ymin, ymax)
PExtents2D(xext::PExtents1D, yext::PExtents1D) =
	PExtents2D(xext.min, xext.max, yext.min, yext.max)

#Transform used for Canvas2D/CanvasF1 (only scale/offset - no rotation).
#xd = (xs + x0) * xu
#yd = (ys + y0) * yu
#where: xu/yu: user coordinates & xd/yd: device coordinates
#NOTE: Format used instead of xs*xu + x0 for improved numeric stability?
#TODO: Does struct cause copy on function calls instead of using pointers??
struct Transform2D
	xs::DReal
	x0::DReal
	ys::DReal
	y0::DReal
end

struct LineStyle
	style::Symbol
	width::Float64 #Device units (~pixels)
	color::Colorant
end

#==Constructor-like functions
===============================================================================#
#=Inputs
 -ext: Visible plot extents
 -inputb: Bounding box of display area.
NOTE:
 -Assumes device y-coordinate increases as we descend on the screen.
=#
function Transform2D(ext::PExtents2D, inputb::BoundingBox)
	xs = (inputb.xmax-inputb.xmin) / (ext.xmax-ext.xmin)
	ys = -(inputb.ymax-inputb.ymin) / (ext.ymax-ext.ymin)
	x0 = + inputb.xmin/xs - ext.xmin
	y0 = + inputb.ymax/ys - ext.ymin

	return Transform2D(xs, x0, ys, y0)
end

#==Accessors
===============================================================================#
xvalues(ext::PExtents2D) = PExtents1D(ext.xmin, ext.xmax)
yvalues(ext::PExtents2D) = PExtents1D(ext.ymin, ext.ymax)

#==Basic operations
===============================================================================#
#Compute vector norm:
function vecnorm(v::Vector2D)
	return sqrt(v.x*v.x+v.y*v.y)
end

#Simple scaling - no commutativity:
Base.:*(s::DReal, pt::Point2D) = Point2D(s*pt.x, s*pt.y)
Base.:*(s::Real, pt::Point2D) = DReal(s)*pt
Base.:+(p1::Point2D, p2::Point2D) = Point2D(p1.x+p2.x, p1.y+p2.y)
Base.:-(p1::Point2D, p2::Point2D) = Point2D(p1.x-p2.x, p1.y-p2.y)

Base.:*(s::DReal, v::Vector2D) = Vector2D(s*v.x, s*v.y)
Base.:*(s::Real, v::Vector2D) = DReal(s)*v
Base.:+(v1::Vector2D, v2::Vector2D) = Vector2D(v1.x+v2.x, v1.y+v2.y)

Base.:+(p::Point2D, v::Vector2D) = Point2D(p.x+v.x, p.y+v.y)
Base.:-(p::Point2D, v::Vector2D) = Point2D(p.x-v.x, p.y-v.y)
Base.:+(v::Vector2D, p::Point2D) = p+v

function union(e1::PExtents1D, e2::PExtents1D)
	umin(v1, v2) = isnan(v1) ? v2 : min(v1, v2)
	umax(v1, v2) = isnan(v1) ? v2 : max(v1, v2)
	return PExtents1D(
		umin(e1.min, e2.min),
		umax(e1.max, e2.max),
	)
end
function union(elist::Vector{PExtents1D})
	result = PExtents1D(DNaN, DNaN)
	for ext in elist
		result = union(result, ext)
	end
	return result
end

function union(e1::PExtents2D, e2::PExtents2D)
	umin(v1, v2) = isnan(v1) ? v2 : min(v1, v2)
	umax(v1, v2) = isnan(v1) ? v2 : max(v1, v2)
	return PExtents2D(
		umin(e1.xmin, e2.xmin),
		umax(e1.xmax, e2.xmax),
		umin(e1.ymin, e2.ymin),
		umax(e1.ymax, e2.ymax),
	)
end
function union(elist::Vector{PExtents2D})
	result = PExtents2D(DNaN, DNaN, DNaN, DNaN)
	for ext in elist
		result = union(result, ext)
	end
	return result
end

#Overwrite with new extents, if defined
#-------------------------------------------------------------------------------
function Base.merge(base::PExtents1D, new::PExtents1D)
	baseifnan(bv, nv) = isnan(nv) ? bv : nv
	emin = baseifnan(base.min, new.min)
	emax = baseifnan(base.max, new.max)
	return PExtents1D(emin, emax)
end

function Base.merge(base::PExtents2D, new::PExtents2D)
	baseifnan(bv, nv) = isnan(nv) ? bv : nv
	xmin = baseifnan(base.xmin, new.xmin)
	xmax = baseifnan(base.xmax, new.xmax)
	ymin = baseifnan(base.ymin, new.ymin)
	ymax = baseifnan(base.ymax, new.ymax)
	return PExtents2D(xmin, xmax, ymin, ymax)
end

Base.isfinite(ext::PExtents1D) = isfinite(ext.min) && isfinite(ext.max)

Base.isfinite(ext::PExtents2D) =
	isfinite(ext.xmin) && isfinite(ext.xmax) &&
	isfinite(ext.ymin) && isfinite(ext.ymax)

#==Mapping/interpolation functions
===============================================================================#

#Map axis -> device coordinates:
function map2dev(xf::Transform2D, pt::Point2D)
	x = (pt.x + xf.x0)*xf.xs
	y = (pt.y + xf.y0)*xf.ys
	return Point2D(x, y)
end
#Map device -> axis coordinates (reverse mapping):
function map2axis(xf::Transform2D, pt::Point2D)
	x = pt.x/xf.xs - xf.x0
	y = pt.y/xf.ys - xf.y0
	return Point2D(x, y)
end

#Map vectors between axis -> device coordinate systems:
function map2dev(xf::Transform2D, pt::Vector2D)
	x = pt.x*xf.xs
	y = pt.y*xf.ys
	return Vector2D(x, y)
end
#Map vectors between device -> axis coordinate systems (reverse mapping):
function map2axis(xf::Transform2D, pt::Vector2D)
	x = pt.x/xf.xs
	y = pt.y/xf.ys
	return Vector2D(x, y)
end

#Interpolate between two points.
#-------------------------------------------------------------------------------
function interpolate(p1::Point2D, p2::Point2D, x::DReal)
	m = (p2.y-p1.y) / (p2.x-p1.x)
	return m*(x-p1.x)+p1.y
end

#Last line


text:
module NIntegration

### Dependencies
using DataStructures
using StaticArrays

### Implementation
path = dirname(realpath(@__FILE__)) # works even for symlinks
include(joinpath(path, "types.jl"))
include(joinpath(path, "rules.jl"))
include(joinpath(path, "utils.jl"))
include(joinpath(path, "integration.jl"))
include(joinpath(path, "3d.jl"))

### Exports
export nintegrate, weightedpoints, weights, points

end # module


text:
# indicator of a PSD

export IndPSD

"""
**Indicator of the set of positive semi-definite cone**

    IndPSD(;scaling=false)

Returns the indicator of the set
```math
C = \\{ X : X \\succeq 0 \\}.
```
The argument to the function can be either a `Symmetric` or `Hermitian` object,
or an object of type `AbstractVector{Float64}` holding a symmetric matrix in (lower triangular) packed storage.

If `scaling = true` then the vectors `y` and `x` in
`prox!(y::AbstractVector{Float64}, f::IndPSD, x::AbstractVector{Float64}, args...)`
have the off-diagonal elements multiplied with `√2` to preserve inner products,
see Vandenberghe 2010: http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf .

I.e. when when `scaling=true`, let `X,Y` be matrices and

`x = (X_{1,1}, √2⋅X_{2,1}, ... ,√2⋅X_{n,1}, X_{2,2}, √2⋅X_{3,2}, ..., X_{n,n})`,

`y = (Y_{1,1}, √2⋅Y_{2,1}, ... ,√2⋅Y_{n,1}, Y_{2,2}, √2⋅Y_{3,2}, ..., Y_{n,n})`

then `prox!(Y, f, X)` is equivalent to `prox!(y, f, x)`.
"""
struct IndPSD <: ProximableFunction
    scaling::Bool
end

IndPSD(;scaling=false) = IndPSD(scaling)

function (f::IndPSD)(X::HermOrSym{T}) where T <: RealOrComplex
  F = eigen(X)
  for i in eachindex(F.values)
    #Do we allow for some tolerance here?
    if F.values[i] <= -1e-14
      return +Inf
    end
  end
  return 0.0
end

is_convex(f::IndPSD) = true
is_cone(f::IndPSD) = true

function prox!(Y::HermOrSym{T}, f::IndPSD, X::HermOrSym{T}, gamma::Real=1.0) where T <: RealOrComplex
  n = size(X, 1)
  F = eigen(X)
  for i in eachindex(F.values)
    F.values[i] = max.(0.0, F.values[i])
  end
  for i = 1:n
    for j = 1:n
      Y.data[i,j] = 0.0
      for k = 1:n
        Y.data[i,j] += F.vectors[i,k]*F.values[k]*F.vectors[j,k]
      end
    end
  end
  return 0.0
end

fun_name(f::IndPSD) = "indicator of positive semidefinite cone"
fun_dom(f::IndPSD) = "Symmetric, Hermitian, AbstractArray{Float64}"
fun_expr(f::IndPSD) = "x ↦ 0 if A ⪰ 0, +∞ otherwise"
fun_params(f::IndPSD) = "none"

function prox_naive(f::IndPSD, X::HermOrSym{T}, gamma::Real=1.0) where T <: RealOrComplex
  F = eigen(X)
  return F.vectors * Diagonal(max.(0.0, F.values)) * F.vectors', 0.0
end

"""
Scales the diagonal of `x` with `val`, where `x` is the lower triangualar part
of a matrix, stored column by column.
"""
function scale_diagonal!(x, val)
    n = Int(sqrt(1/4+2*length(x))-1/2)
    k = -n
    for i = 1:n
         k += n - i + 2    #Calculate indices of diagonal elements recursively (paralell faster?)
         x[k] *= val      #Scale diagonal
    end
end

### Below: with AbstractVector argument

function (f::IndPSD)(x::AbstractVector{T}) where T <: Float64
  y = copy(x)
  f.scaling && scale_diagonal!(y, sqrt(2)) #If scaling, scale diagonal (eigenvalues scaled by sqrt(2))

  Z = dspev!(:N, :L, y)
  for i in 1:length(Z)
    #Do we allow for some tolerance here?
    if Z[i] <= -1e-14
      return +Inf
    end
  end
  return 0.0
end

function prox!(y::AbstractVector{Float64}, f::IndPSD, x::AbstractVector{Float64}, gamma::Real=1.0)
  y .= x              # Copy x since dspev! corrupts input

  f.scaling && scale_diagonal!(y, sqrt(2)) #If scaling, scale diagonal

  (W, Z) = dspevV!(:L, y)
  W = max.(W, 0.0)         # NonNeg eigenvalues
  M = Z.*W'             # Equivalent to Z*diagm(W) without constructing W matrix
  M = M*Z'              # Now let M = Z*diagm(W)*Z'
  n = length(W)
  k = 1
  for j in 1:n, i in j:n  # Store lower diagonal of M in y
    y[k] = M[i,j]
    k = k+1
  end

  f.scaling && scale_diagonal!(y, 1/sqrt(2))  #If scaling, un-scale diagonal

  return 0.0
end

function prox_naive(f::IndPSD, x::AbstractVector{T}, gamma::Real=1.0) where T<:Float64
  n = Int(sqrt(1/4+2*length(x))-1/2)  # Formula for size of matrix
  X = Array{T,2}(undef, n, n)
  k = 1
  for j = 1:n, i = j:n                # Store y in M
    X[i,j] = x[k]                   # Lower half
    if i != j
      X[j,i] = x[k]               # Strictly upper half
    end
    k = k+1
  end
  # Scale diagonal elements by sqrt(2)
  # See Vandenberghe 2010 http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf
  # It's equivalent to scaling off-diagonal by 1/sqrt(2) and working with sqrt(2)*X
  if f.scaling
    for i = 1:n
      X[i,i] *= sqrt(2)
    end
  end
  X, v = prox_naive(f, Symmetric(X), gamma)

  if f.scaling  #Scale diagonal elements back
    for i = 1:n
      X[i,i] /= sqrt(2)
    end
  end

  y = similar(x)
  k = 1
  for j = 1:n, i = j:n                # Store Lower half of X in y
    y[k] = X[i,j]
    k = k+1
  end
  return y, 0.0
end


text:
###################
# plate generator #
###################

"""
Generator that makes many independent application of a kernel generator,
similar to 'map'.  The arguments are a tuple of vectors, each of of length N,
where N is the nubmer of applications of the kernel.
"""
struct PlateOfDist{T} <: Generator{PersistentVector{T},VectorDistTrace{T}}
    kernel::Distribution{T}
end

accepts_output_grad(plate::PlateOfDist) = false # TODO
has_argument_grads(plate::PlateOfDist) = has_argument_grads(plate.kernel)

function plate(kernel::Distribution{T}) where {T}
    PlateOfDist{T}(kernel)
end

function get_static_argument_types(plate::PlateOfDist)
    [Vector{typ} for typ in get_static_argument_types(plate.kernel)]
end

function simulate(gen::PlateOfDist{T}, args) where {T}
    len = length(args[1])
    values = Vector{T}(undef, len)
    score = 0.
    for key=1:len
        kernel_args = get_args_for_key(args, key)
        values[key] = random(gen.kernel, kernel_args...)
        score += logpdf(gen.kernel, values[key], kernel_args...)
    end
    persist_vec = Gen.PersistentVector{T}(values)
    call = CallRecord(score, persist_vec, args)
    Gen.VectorDistTrace{T}(persist_vec, call)
end

function generate(gen::PlateOfDist{T}, args, constraints) where {T}
    len = length(args[1])
    values = Vector{T}(undef, len)
    score = 0.
    weight = 0.
    for key=1:len
        kernel_args = get_args_for_key(args, key)
        is_constrained = has_leaf_node(constraints, key)
        if is_constrained
            value = get_leaf_node(constraints, key)
        else
            value = random(gen.kernel, kernel_args...)
        end
        values[key] = value
        lpdf = logpdf(gen.kernel, value, kernel_args...)
        if is_constrained
            weight += lpdf
        end
        score += lpdf
    end
    # TODO also check that there are no extra constraints
    persist_vec = Gen.PersistentVector{T}(values)
    call = CallRecord(score, persist_vec, args)
    trace = Gen.VectorDistTrace{T}(persist_vec, call)
    (trace, weight)
end

function update(gen::PlateOfDist{T}, new_args, args_change::Nothing, trace::VectorDistTrace{T}, constraints) where {T}
    (new_length, prev_length) = get_prev_and_new_lengths(new_args, trace)
    @assert new_length == prev_length
    # TODO handle increases or decreaeses int e length

    prev_args = get_call_record(trace).args
    score = get_call_record(trace).score
    weight = 0.
    values = trace.values
    for key=1:new_length
        prev_value = values[key]
        prev_kernel_args = get_args_for_key(prev_args, key)
        kernel_args = get_args_for_key(new_args, key)
        is_constrained = has_leaf_node(constraints, key)
        if is_constrained
            @assert false
            value = get_leaf_node(constraints, key)
            values = assoc(values, key, value)
        else
            value = prev_value
        end
        lpdf = logpdf(gen.kernel, value, kernel_args...)
        prev_lpdf = logpdf(gen.kernel, prev_value, prev_kernel_args...)
        weight += lpdf - prev_lpdf
        score += lpdf - prev_lpdf
    end
    # TODO also check that there are no extra constraints
    # TODO handle retchange
    call = CallRecord(score, values, new_args)
    new_trace = Gen.VectorDistTrace{T}(values, call)
    (new_trace, weight, EmptyAssignment(), nothing)
end

# TODO handle other selection types
function backprop_trace(gen::PlateOfDist{T}, trace::VectorDistTrace{T}, selection::EmptyAddressSet, retval_grad::Nothing) where {T}
    call = get_call_record(trace)
    args = call.args
    n_args = length(args)
    len = length(args[1])

    has_grads = has_argument_grads(gen.kernel)
    arg_grad = Vector(undef, n_args)
    for (i, has_grad) in enumerate(has_grads)
        if has_grad
            arg_grad[i] = Vector(undef, len)
        else
            arg_grad[i] = nothing
        end
    end

    #value_trie = DynamicAssignment()
    #gradient_trie = DynamicAssignment()

    for key=1:len
        value = trace.values[key]
        kernel_args = get_args_for_key(args, key)
        #subtrace = trace.subtraces[key]
        #if has_internal_node(selection, key)
            #sub_selection = get_internal_node(selection, key)
        #else
            #sub_selection = EmptyAddressSet()
        #end
        #kernel_retval_grad = (retval_grad == nothing) ? nothing : retval_grad[key]
        kernel_grads = logpdf_grad(gen.kernel, value, kernel_args...)
        value_grad = kernel_grads[1]
        kernel_arg_grad = (kernel_grads[2:end]...,)
        #(kernel_arg_grad::Tuple, kernel_value_trie, kernel_gradient_trie) = backprop_trace(
            #gen.kernel, subtrace, sub_selection, kernel_retval_grad)
        #set_internal_node!(value_trie, key, kernel_value_trie)
        #set_internal_node!(gradient_trie, key, kernel_gradient_trie)
        for (i, grad, has_grad) in zip(1:n_args, kernel_arg_grad, has_grads)
            if has_grad
                arg_grad[i][key] = grad
            end
        end
    end
    value_trie = EmptyAssignment()
    gradient_trie = EmptyAssignment()
    ((arg_grad...,), value_trie, gradient_trie)
end

export plate

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
#=
Case folding for Unicode Str types

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

const _wide_upper       = ChrBase._wide_upper
const _wide_lower_latin = ChrBase._wide_lower_latin

function uppercase_first(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _is_lower_a(ch) || return str
        buf, out = _allocate(UInt8, len)
        unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, ch - 0x20)
        Str(C, buf)
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _is_upper_a(ch) || return str
        buf, out = _allocate(UInt8, len)
        unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, ch + 0x20)
        Str(C, buf)
    end
end

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:ASCIICSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _is_lower_a(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function _lower(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:ASCIICSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _is_upper_a(ch) && set_codeunit!(out, ch + 0x20)
        out += 1
    end
    Str(C, buf)
end

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:LatinCSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _can_upper_al(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _is_lower_a(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function lowercase(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _is_upper_a(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function uppercase_first(str::MaybeSub{S}) where {C<:LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _can_upper_al(ch) || return str
        buf, out = _allocate(UInt8, len)
        set_codeunit!(out, ch - 0x20)
        len > 1 && unsafe_copyto!(out + 1, pnt+1, len-1)
        Str(C, buf)
    end
end

# Special handling for characters that can't map into Latin1
function uppercase_first(str::MaybeSub{S}) where {C<:_LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if _wide_lower_latin(ch)
            buf, out = _allocate(UInt16, len)
            _widen!(out, pnt, pnt + len)
            set_codeunit!(out, _wide_upper(ch))
            Str(_UCS2CSE, buf)
        elseif _can_upper_al(ch)
            buf8, out8 = _allocate(UInt8, len)
            len > 1 && unsafe_copyto!(out8, pnt, len)
            set_codeunit!(out8, ch - 0x20)
            Str(_LatinCSE, buf8)
        else
            str
        end
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:Latin_CSEs,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _is_upper_al(ch) || return str
        buf, out = _allocate(UInt8, len)
        set_codeunit!(out, ch + 0x20)
        len > 1 && unsafe_copyto!(out+1, pnt+1, len-1)
        Str(C, buf)
    end
end

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:_LatinCSE}
    fin = beg + len
    cur = beg + off
    # Need to scan the rest of the string to see if _widenupper needs to be called
    while cur < fin
        _wide_lower_latin(get_codeunit(cur)) && return _widenupper(beg, off, len)
        cur += 1
    end
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _can_upper_al(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function _widen!(dst::Ptr{T}, src::Ptr{S}, fin::Ptr{S}) where {T<:CodeUnitTypes, S<:CodeUnitTypes}
    while src < fin
        set_codeunit!(dst, get_codeunit(src)%T)
        dst += sizeof(T)
        src += sizeof(S)
    end
    nothing
end
const _narrow! = _widen!  # When this is optimized for SSE/AVX/etc. instructions, will be different

function _widenupper(beg::Ptr{UInt8}, off, len)
    buf, out = _allocate(UInt16, len)
    fin = bytoff(out, len)
    cur = beg + off
    _widen!(out, beg, cur)
    out = bytoff(out, off)
    while out < fin
        ch = get_codeunit(cur)
        set_codeunit!(out, _can_upper_al(ch) ? (ch - 0x20) : _wide_upper(ch))
        cur += 1
        out += 2
    end
    Str(_UCS2CSE, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _can_upper_al(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function uppercase(str::MaybeSub{S}) where {C<:_LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            ch = get_codeunit(pnt)
            _wide_lower_latin(ch) && return _widenupper(beg, pnt-beg, len)
            _can_upper_al(ch) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function _lower(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:Latin_CSEs}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _is_upper_al(ch) && set_codeunit!(out, ch + 0x20)
        out += 1
    end
    Str(C, buf)
end

function lowercase(str::MaybeSub{S}) where {C<:Latin_CSEs,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _is_upper_al(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

_is_latin_ucs2(len, pnt) = _check_mask_ul(pnt, len, _latin_mask(UInt16))

# result must have at least one character > 0xff, so if the only character(s)
# > 0xff became <= 0xff, then the result may need to be narrowed and returned as _LatinStr

function _lower(::Type{C}, beg, off, len) where {C<:_UCS2CSE}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    lenw = len*sizeof(CU)
    fin = out + lenw
    out += off
    flg = false
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _is_upper_a(ch) && set_codeunit!(out, ch += 0x20)
        elseif ch <= 0xff
            _is_upper_l(ch) && set_codeunit!(out, ch += 0x20)
        elseif ch <= 0xffff
            if _can_lower_bmp(ch)
                ch = _lower_bmp(ch)
                flg = ch <= 0xff
                set_codeunit!(out, ch)
            end
        end
        out += sizeof(CU)
    end
    if flg && (src = reinterpret(Ptr{UInt16}, pointer(buf)); _is_latin_ucs2(lenw, src))
        buf8 = _allocate(len)
        _narrow!(pointer(buf8), src, src + lenw)
        Str(_LatinCSE, buf8)
    else
        Str(C, buf)
    end
end

function _lower(::Type{C}, beg, off, len) where {C<:Union{UCS2CSE,UTF32_CSEs}}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len*sizeof(CU))
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0xff
            _is_upper_al(ch) && set_codeunit!(out, ch += 0x20)
        elseif ch <= 0xffff
            _can_lower_bmp(ch) && set_codeunit!(out, _lower_bmp(ch))
        elseif ch <= 0x1ffff
            _can_lower_slp(ch) && set_codeunit!(out, _lower_slp(ch))
        end
        out += sizeof(CU)
    end
    Str(C, buf)
end

function lowercase_first(str::MaybeSub{S}) where {C<:_UCS2CSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        (ch <= 0xff ? _is_upper_al(ch) : ch <= 0xffff ? _can_lower_bmp(ch) :
         (ch <= 0x1ffff && _can_lower_slp(ch))) ||
         return str
        cl = _lowercase(ch)
        if ch > 0xff && cl <= 0xff && _check_mask_ul(pnt+1, len-1, _latin_mask(UInt16))
            buf8, out8 = _allocate(UInt8, len)
            len > 1 && _narrow!(out8, pnt, pnt + len)
            set_codeunit!(out8, cl)
            Str(_LatinCSE, buf8)
        else
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
            set_codeunit!(out, cl)
            Str(C, buf)
        end
    end
end

function uppercase_first(str::MaybeSub{S}) where {C<:Union{UCS2_CSEs,UTF32_CSEs},S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        cp = _titlecase(ch)
        ch == cp && return str
        buf, out = _allocate(codeunit(C), len)
        len > 1 && unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, cp)
        Str(C, buf)
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:Union{UCS2CSE,UTF32_CSEs},S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _can_lower_ch(ch) || return str
        buf, out = _allocate(codeunit(C), len)
        len > 1 && unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, _lowercase(ch))
        Str(C, buf)
    end
end

function lowercase(str::MaybeSub{S}) where {C<:Union{UCS2_CSEs,UTF32_CSEs},S<:Str{C}}
    @preserve str begin
        CU = codeunit(C)
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            _can_lower_ch(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, ncodeunits(str))
            pnt += sizeof(CU)
        end
    end
    str
end

function _upper(::Type{C}, beg, off, len) where {C<:Union{UCS2_CSEs,UTF32_CSEs}}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len*sizeof(CU))
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _is_lower_a(ch) && set_codeunit!(out, ch -= 0x20)
        elseif ch <= 0xff
            set_codeunit!(out, _uppercase_l(ch))
        elseif ch <= 0xffff
            _can_upper_bmp(ch) && set_codeunit!(out, _upper_bmp(ch))
        elseif ch <= 0x1ffff
            _can_upper_slp(ch) && set_codeunit!(out, _upper_slp(ch))
        end
        out += sizeof(CU)
    end
    Str(C, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:Union{UCS2_CSEs,UTF32_CSEs},S<:Str{C}}
    @preserve str begin
        CU = codeunit(C)
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            _can_upper_ch(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, ncodeunits(str))
            pnt += sizeof(CU)
        end
        str
    end
end


text:
module MaxHelpingHandSetFlagOnSpawnController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/SetFlagOnSpawnController" SetFlagOnSpawnController(x::Integer, y::Integer, flag::String="flag_name", enable::Bool=false, onlyOnRespawn::Bool=false)

const placements = Ahorn.PlacementDict(
    "Set Flag On Spawn Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SetFlagOnSpawnController
    )
)

function Ahorn.selection(entity::SetFlagOnSpawnController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SetFlagOnSpawnController, room::Maple.Room) = Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/set_flag_on_spawn", 0, 0)

end

max_stars_repo_name: JaThePlayer/MaxHelpingHand

max_issues_repo_name: JaThePlayer/MaxHelpingHand

max_forks_repo_name: JaThePlayer/MaxHelpingHand


text:
# Get All Cookies
"""
    cookies(session::Session)::Vector{Cookie}

The Get Page Source command returns a string serialization of the DOM of the current browsing context active document.
"""
function cookies(session::Session)::Vector{Cookie}
    @unpack addr, id = session
    response =
        HTTP.get("$addr/session/$id/cookie", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    convert(Vector{Cookie}, Cookie.(JSON3.read(response.body).value))
end


text:
"""
    alexnet(; nclasses = 1000)

Create an AlexNet model
([reference](https://papers.nips.cc/paper/2012/file/c399862d3b9d6b76c8436e924a68c45b-Paper.pdf)).

# Arguments
- `nclasses`: the number of output classes
"""
function alexnet(; nclasses = 1000)
  layers = Chain(Chain(Conv((11, 11), 3 => 64, stride = (4, 4), relu, pad = (2, 2)),
                       MaxPool((3, 3), stride = (2, 2)),
                       Conv((5, 5), 64 => 192, relu, pad = (2, 2)),
                       MaxPool((3, 3), stride = (2, 2)),
                       Conv((3, 3), 192 => 384, relu, pad = (1, 1)),
                       Conv((3, 3), 384 => 256, relu, pad = (1, 1)),
                       Conv((3, 3), 256 => 256, relu, pad = (1, 1)),
                       MaxPool((3, 3), stride = (2, 2)),
                       AdaptiveMeanPool((6,6))),
                 Chain(flatten,
                       Dropout(0.5),
                       Dense(256 * 6 * 6, 4096, relu),
                       Dropout(0.5),
                       Dense(4096, 4096, relu),
                       Dense(4096, nclasses)))

  return layers
end

"""
    AlexNet(; pretrain = false, nclasses = 1000)

Create a `AlexNet`.
See also [`alexnet`](#).

!!! warning
    `AlexNet` does not currently support pretrained weights.

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes
"""
struct AlexNet{T}
  layers::T
end

function AlexNet(; pretrain = false, nclasses = 1000)
  layers = alexnet(nclasses = nclasses)
  pretrain && pretrain_error("AlexNet")

  AlexNet(layers)
end

@functor AlexNet

(m::AlexNet)(x) = m.layers(x)


text:
@testset "func: symmetrise" begin
    g = σ3 ⊗ one(σ3) - one(σ3) ⊗ σ3
    vg = (g ⊗ one(g) - one(g) ⊗ g)
    basis = U1Space(-1=>1, 1=>1)
    SG = PEPOKit.SymmetryGroup(vg, basis)
    H = σ1 ⊗ σ1 - σ2 ⊗ σ2
    L = PEPOKit.StandardLiouvillian(PEPOKit.coherent(H) + PEPOKit.dissipator(2), nothing)
    symL = PEPOKit.symmetrise(SG,L)

    @test PEPOKit.unsymmetrise(symL) ≈ L
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


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

mutable struct VirtualNetworkGatewayListConnectionsResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{VirtualNetworkGatewayConnectionListEntity} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function VirtualNetworkGatewayListConnectionsResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(VirtualNetworkGatewayListConnectionsResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(VirtualNetworkGatewayListConnectionsResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type VirtualNetworkGatewayListConnectionsResult

const _property_map_VirtualNetworkGatewayListConnectionsResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_VirtualNetworkGatewayListConnectionsResult = Dict{Symbol,String}(Symbol("value")=>"Vector{VirtualNetworkGatewayConnectionListEntity}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ VirtualNetworkGatewayListConnectionsResult }) = collect(keys(_property_map_VirtualNetworkGatewayListConnectionsResult))
Swagger.property_type(::Type{ VirtualNetworkGatewayListConnectionsResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkGatewayListConnectionsResult[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayListConnectionsResult }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayListConnectionsResult[property_name]

function check_required(o::VirtualNetworkGatewayListConnectionsResult)
    true
end

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


text:
using BinaryBuilder, Pkg

name = "RAxML"
version = v"8.2.12"
hash = "a33ff40640b4a76abd5ea3a9e2f57b7dd8d854f6"
sources = [
    GitSource("https://github.com/stamatak/standard-RAxML.git", hash),
]

script = raw"""
cd ${WORKSPACE}/srcdir/standard-RAxML

# see section IV in manual for license spec
# https://github.com/stamatak/standard-RAxML/blob/master/manual/NewManual.pdf
install_license /usr/share/licenses/GPL-3.0+

# Sequential
make -f Makefile.gcc -j${nproc}
rm *.o
make -f Makefile.SSE3.gcc -j${nproc}
rm *.o
make -f Makefile.AVX.gcc -j${nproc}
rm *.o

# PTHREADS
make -f Makefile.PTHREADS.gcc -j${nproc}
rm *.o
make -f Makefile.SSE3.PTHREADS.gcc -j${nproc}
rm *.o
make -f Makefile.AVX.PTHREADS.gcc -j${nproc}
rm *.o

# MPI version
if [[ "$target" == *-mingw* ]]; then
    GLOBAL_DEPS="axml.h globalVariables.h rmq.h rmqs.h ${includedir}/mpi.h"
    LIBRARIES="-lm -lmsmpi -I${includedir} -L${libdir}"

    CFLAGS="-D_WAYNE_MPI -D_GNU_SOURCE -fomit-frame-pointer -funroll-loops -O2 -msse -I${includedir} -L${libdir} -lmsmpi"
    make -f Makefile.MPI.gcc -j${nproc} \
        CC=${target}-gcc \
        CFLAGS="$CFLAGS" \
        GLOBAL_DEPS="${GLOBAL_DEPS}" \
        LIBRARIES="${LIBRARIES}"
    rm *.o

    CFLAGS="-D_WAYNE_MPI -D__SIM_SSE3 -O2 -D_GNU_SOURCE -msse3 -fomit-frame-pointer -funroll-loops -I${includedir} -L${libdir} -lmsmpi"
    make -f Makefile.SSE3.MPI.gcc -j${nproc} \
        CC=${target}-gcc \
        CFLAGS="$CFLAGS" \
        GLOBAL_DEPS="${GLOBAL_DEPS}" \
        LIBRARIES="${LIBRARIES}"
    rm *.o

    CFLAGS="-D_WAYNE_MPI -D__SIM_SSE3 -D__AVX -O2 -D_GNU_SOURCE -msse3 -fomit-frame-pointer -funroll-loops -I${includedir} -L${libdir} -lmsmpi"
    make -f Makefile.AVX.MPI.gcc -j${nproc} \
        CC=${target}-gcc \
        CFLAGS="$CFLAGS" \
        GLOBAL_DEPS="${GLOBAL_DEPS}" \
        LIBRARIES="${LIBRARIES}"
    rm *.o
else
    make -f Makefile.MPI.gcc -j${nproc}
    rm *.o
    make -f Makefile.SSE3.MPI.gcc -j${nproc}
    rm *.o
    make -f Makefile.AVX.MPI.gcc -j${nproc}
    rm *.o
fi

for f in raxmlHPC*; do
    install -Dvm 755 $f ${bindir}/${f}${exeext};
done
"""

# currently RAxML assumes intel cpus with at least SSE instructions
platforms = filter(p -> BinaryBuilder.proc_family(p) == "intel", supported_platforms())

products = [
    ExecutableProduct("raxmlHPC", :raxmlHPC),
    ExecutableProduct("raxmlHPC-SSE3", :raxmlHPC_SSE3),
    ExecutableProduct("raxmlHPC-AVX", :raxmlHPC_AVX),
    ExecutableProduct("raxmlHPC-PTHREADS", :raxmlHPC_PTHREADS),
    ExecutableProduct("raxmlHPC-PTHREADS-SSE3", :raxmlHPC_PTHREADS_SSE3),
    ExecutableProduct("raxmlHPC-PTHREADS-AVX", :raxmlHPC_PTHREADS_AVX),
    ExecutableProduct("raxmlHPC-MPI", :raxmlHPC_MPI),
    ExecutableProduct("raxmlHPC-MPI-SSE3", :raxmlHPC_MPI_SSE3),
    ExecutableProduct("raxmlHPC-MPI-AVX", :raxmlHPC_MPI_AVX),
]

dependencies = [
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4"); platforms=filter(!Sys.iswindows, platforms)),
    Dependency(PackageSpec(name="MicrosoftMPI_jll", uuid="9237b28f-5490-5468-be7b-bb81f5f5e6cf"); platforms=filter(Sys.iswindows, platforms)),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")


text:
ntime = 10
Δt = 0.1

ts = Timeseries{eltype(Δt)}(ntime, Δt, 1)
@test typeof(ts) <: AbstractArray
@test length(ts) == ntime+1
@test size(ts) == (ntime+1,)
@test indices(ts, 1) == 0:ntime

compute_timeseries!(ts)
t = collect(0:Δt:ntime*Δt)
@test_approx_eq_eps(ts, t, eps())

ts1 = Timeseries(ntime, Δt)
compute_timeseries!(ts1)
@test ts1 == ts


text:
"""
    set(x::Coloring)
    set(x::Blending)
These methods set the color or blend defined in `x`.
"""
set(color::Coloring) = Luxor.setcolor(Luxor.sethue(color.hue)..., color.alpha)
set(color::Blending) = Luxor.setblend(Luxor.blend(color.direction..., color.hue...))

"""
    draw(plot::Plot{T}) where T<:Geometry
These methods draw a waterfall plot using all of the styling and values stored in `plot`.
"""
function draw(shape::T) where T <: Shape
    set(shape.color)
    _draw(shape)
    return nothing
end

function draw(lab::T) where T <: Label
    Luxor.fontsize(lab.scale * FONTSIZE)
    Luxor.setcolor(Luxor.sethue("black")...)
    _draw(lab)

    Luxor.fontsize(FONTSIZE)
    return nothing
end

# Draw some type values.
draw(x::Cascade{T}) where T<:Geometry = draw(collect_data(x))
draw(x::T) where T<:Geometry = draw(x.shape)

# Draw all Type values
draw(x::T; kwargs...) where T<:Plot = _draw(x)
draw(x::T) where T<:Axis = _draw(x)
draw(x::Annotation) = _draw(x)
draw(x::Handle) = _draw(x)
draw(x::String) = nothing
draw(x::Missing) = nothing

draw(lst::AbstractArray) = [draw(x) for x in lst]
draw(pair::Pair) = draw(collect(pair))

"""
    _draw(lab::T) where T <: Label
    _draw(shape::T) where T <: Shape
These are helper methods for drawing different types of labels and shapes. They do not
include steps like setting colors and font size, since these processes are consistent across
`Label` and `Shape`.
"""
function _draw(lab::Label{String})
    Luxor.text(lab.text, lab.position; halign=lab.halign, valign=lab.valign, angle=lab.angle)
    return nothing
end

function _draw(lab::Label{Vector{String}})
    Luxor.textbox(lab.text, lab.position;
        alignment = lab.halign,
        leading = lab.leading*lab.scale*FONTSIZE,
    )
    return nothing
end

# _draw(lab::Label{Missing}) = nothing

# Draw shapes.
_draw(shape::Arrow) = Luxor.arrow(shape.position..., linewidth=2.0)
_draw(shape::Box) = Luxor.box(shape.position..., shape.style)
_draw(shape::Line) = Luxor.line(shape.position..., shape.style)
_draw(shape::Poly) = Luxor.poly(shape.position, close=true, shape.style)

# Draw all Type values.
_draw(x::Any) = draw(collect(values(x)))

"""
    _draw_legend()
"""
function _draw_legend()
    Luxor.text("mean", Luxor.Point(WIDTH-5*SEP,5*SEP); halign=:left, valign=:middle)
    draw(Box(
        (Luxor.Point(WIDTH-7.5*SEP,4.5*SEP),Luxor.Point(WIDTH-5.5*SEP,5.5*SEP)),
        _define_from(Coloring, "black"),
        :stroke,
    ))
    return nothing
end

"""
"""
function padding(ax::XAxis)
    lab = ax.ticksublabels
    lines = maximum(length.(getfield.(lab,:text)))
    return lines * FONTSIZE
end


text:
# 4-box biosphere impulse response, see Joos (1996) appendix A.3

let
	# to get weights for each box, take bioBoxCoeff.*bioboxtime
	biodata = [
		#bioBoxCoeff	bioboxtime
		-0.71846		2.181818				# detritus
		0.70211			2.85714					# ground vegetation
		0.013414		20						# wood
		0.0029323		100						# soil organic carbon
	]

	bioboxweight	= biodata[:,1] .* biodata[:,2]
	bioboxtime		= biodata[:,2]

	baseNPP			= 60.0		# GtC/year

	global const bioboxlength = length(bioboxweight)

	global function biocarbon!(s::ClimateState, p::ClimateParams)
		@unpack Concentration, BioReservoir, Temp, Temp_land = s
		@unpack timestep, fertilization, bioQ10factor = p
		# BioNPP [GtC/year]: Net Primary Production of biosphere, increases with CO2 concentration (fertilization)
		BioNPP = baseNPP * (1 + fertilization*log(Concentration[:CO2]/conc_CO2_preind))

		# BioReservoir [GtC]: carbon content of 4 biosphere boxes, temperature feedback controlled by bioQ10factor
		dt_BioReservoir = BioNPP * bioboxweight - BioReservoir * bioQ10factor^(Temp_land/10) ./ bioboxtime
		BioReservoir[:] += timestep * dt_BioReservoir

		# Net increase of terrestrial biomass
		NetFluxBiosphere = sum(dt_BioReservoir)

		@pack! s = NetFluxBiosphere, BioReservoir
	end

	global function init_biocarbon!(s::ClimateState, p::ClimateParams)
		@unpack fertilization, equilibriumCO2 = p
		startNPP = baseNPP * (1 + fertilization*log(equilibriumCO2/conc_CO2_preind))	# GtC/year
		NetFluxBiosphere = 0.0
		BioReservoir = startNPP * bioboxweight .* bioboxtime

		@pack! s = NetFluxBiosphere, BioReservoir
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
push!(LOAD_PATH, "./jl")

using Base.Test
using BQL

include("../citydb_json.jl")

C = Const(42)
@test C(citydb) == 42

I = This()
@test I(42) == 42

Departments = Field(:departments)
departments = Departments(citydb)
@test isa(departments, Array)

O = Select(:x => Const(42))
@test O(citydb) == Dict("x" => 42)

Name = Field(:name)
Department_Names = Departments >> Name
department_names = Department_Names(citydb)
@test "POLICE" in department_names

Num_Dept = Count(Departments)
@test Num_Dept(citydb) == 35

Employees = Field(:employees)
Departments_With_Num_Empl = Departments >> Select(:name => Name, :num_empl => Count(Employees))
departments_with_num_empl = Departments_With_Num_Empl(citydb)
@test Dict("name" => "POLICE", "num_empl" => 13570) in departments_with_num_empl

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module Null_models_I_and_II
using LinearAlgebra
using DelimitedFiles

export sample_from_Null_model_I, sample_from_Null_model_II, PairwiseHammingDist,translate_fasta_to_num_matrix, transform_MSA_fasta, export_fasta_file
include("utils.jl")

##Dirac delta
function deltak(i::Int64,j::Int64)
  i==j ? 1.0 : 0.0
    end

###compute hamming distance matrix D
 function PairwiseHammingDist(sample::Array{Int64,2})
    Y=sample[1:end,1:end]
    Y = Y'
    (N,M) = size(Y)
    hav = zeros(M,M)

  @inbounds  for m in 1:M
               for n in (m+1):M
                hd = 0.
                for i = 1:N
                 hd += Y[i,m]!=Y[i,n]
                 end
               hd = hd / N

              hav[m,n] = hav[n,m]=hd
               end
              end
    return hav
end

function PairwiseHammingDist(sample::Any)
   if typeof(sample)==Array{Char,2}
     (M,N)=size(sample)
	 Y=Array{Int}(undef,M,N);
	 for i in 1:M
		for j in 1:N
         Y[i,j]=Di21[sample[i,j]]
        end
     end
    elseif typeof(sample)==Array{Int64,2}
          Y=sample[1:end,1:end]
    elseif typeof(sample)==String
        sample=translate_fasta_to_num_matrix(sample)
        Y=sample[1:end,1:end]
    else
     error("Null_models_I_and_II.jl - PairwiseHammingDist: input parameter format is wrong.\n")
   end
    Y = Y'
    (N,M) = size(Y)
    hav = zeros(M,M)

  @inbounds  for m in 1:M
               for n in (m+1):M
                hd = 0.
                for i = 1:N
                 hd += Y[i,m]!=Y[i,n]
                 end
               hd = hd / N

              hav[m,n] = hav[n,m]=hd
               end
              end
    return hav
end

#given column k and two rows m,n an exchange of the align entries a_km and a_kn is attempted.
function swap_align!(align::Array{Int64,2},m1::Int64,m2::Int64,k::Int64)
    alpha=align[m1,k]
    bet=align[m2,k]
    align[m1,k]=bet
    align[m2,k]=alpha
end

### cost function
function Energy(dtarget::Array{Float64,2},dnew::Array{Float64,2})
 return sum((dtarget .- dnew) .^2)
end

## frobenious distance between matrices
function compute_l2_error(d1::Array{Float64,2},d2::Array{Float64,2})

     return sqrt(sum((d1 .- d2).^2))/size(d1)[1]
end

# return delta matrix which store the increment in the hamming distance entries
function delta_d(align::Array{Int64,2},delta::Array{Float64,2},m1::Int64,m2::Int64,k::Int64)
    (M,L)=size(align)
    #delta .=0
 @inbounds @simd for j in 1:M
        if j!=m1&&j!=m2
        delta[j,m1]=-(deltak(align[m2,k],align[j,k])-deltak(align[m1,k],align[j,k]))/L
        delta[j,m2]=-delta[j,m1]
        delta[m1,j]=delta[j,m1]
        delta[m2,j]=delta[j,m2]
            end

        end

   return delta
end

function delta_d!(align::Array{Int64,2},delta::Array{Float64,2},m1::Int64,m2::Int64,k::Int64)
    (M,L)=size(align)
 @fastmath @inbounds @simd for j in 1:M
        if j!=m1&&j!=m2
        delta[j,m1]=-1.0*(deltak(align[m2,k],align[j,k])-deltak(align[m1,k],align[j,k]))/L
        delta[j,m2]=-1.0*delta[j,m1]
        delta[m1,j]=delta[j,m1]
        delta[m2,j]=delta[j,m2]
            end

        end

end

## change produced at d_current as consequence of swap entries m1 and m2
 function change_d!(d_current::Array{Float64,2},delta::Array{Float64,2},m1::Int64,m2::Int64)
        M=size(d_current)[1]

    @fastmath @inbounds @simd for j in 1:M
                     if j!=m1&&j!=m2
                     d_current[j,m1]=d_current[j,m1]+delta[j,m1]
                     d_current[j,m2]=d_current[j,m2]+delta[j,m2]
                     d_current[m1,j]=d_current[j,m1]
                     d_current[m2,j]=d_current[j,m2]
                     end
                    end

        end

##Frobenious norm betwenn hamming distances matrices at times t+1 (dinic) respect to dtarget
## given an alignment permutation
 function dE(dtarget::Array{Float64,2},dinic::Array{Float64,2},delta::Array{Float64,2},m1::Int64,m2::Int64)
  output=0.0
  s=size(dtarget)[1]
  @fastmath @inbounds for j in 1:s
  if j!=m1&&j!=m2
  output+=2*(dinic[j,m1]-dtarget[j,m1])*delta[j,m1]+delta[j,m1]^2+2*(dinic[j,m2]-dtarget[j,m2])*delta[j,m2]+delta[j,m2]^2
              end
       end
    return 2.0*output
 end

   function do_swap(align::Array{Int64,2},delta::Array{Float64,2},d_target::Array{Float64,2},d_current::Array{Float64,2},temp::Float64,energ::Float64)

    (M,L)=size(align)
    m1=rand(1:M)
    m2=rand(1:M)
    k=rand(1:L)
    de=0.0
    while align[m1,k]==align[m2,k]
           m1=rand(1:M)
           m2=rand(1:M)
           k=rand(1:L)
    end

     delta_d!(align,delta,m1,m2,k)###
     de=dE(d_target,d_current,delta,m1,m2)####
     if de<0.0 || rand()< exp(-de/temp)
            swap_align!(align,m1,m2,k)
            energ=energ+de
            change_d!(d_current,delta,m1,m2)###
               return 1,energ,d_current
           else
            return 0,energ,d_current
           end

    end

function update_param_SA(alignment::Array{Int64,2},alignment_best::Array{Int64,2},n_it_tot::Int64,n_it_current::Int64,flag::Int64,acceptance_rate::Float64, energy_min ::Float64,energy_min_prev ::Float64,energy::Float64,T_init::Float64,T_min::Float64,n_it_max::Int64,T_factor_slow::Float64,T_factor_fast::Float64)
 (M,L)=size(alignment)
 stop_condition = 0
 n_it_max =n_it_max#N_IT_MAX;
 flag_max = 100;
 n_swap_min = convert(Int64,round(M*L/5))
 acceptance_threshold_low = 0.0001
 acceptance_threshold_high = 0.95
 T=T_init
 T_min = T_min#2/(10*M)
 #T_factor_slow =T_factor_slow#T_FACTOR_SLOW,
 #T_factor_fast =T_factor_slow#T_FACTOR_FAST;

 #Update number of iterations for next MCMC
    if acceptance_rate!=0
      n_it_current = convert(Int64,round(n_swap_min / acceptance_rate))
    else
        n_it_current = 2*n_it_current
    end

    if n_it_current > M*L
        n_it_current = M*L
    end
##If a better solution has been found and acceptance rate is high enough, decrease temperature
    if energy_min_prev > energy_min && acceptance_rate > acceptance_threshold_low

#linear cooling scheme with two slopes
        if acceptance_rate > acceptance_threshold_high
        T = T*T_factor_fast;
        else
        T = T*T_factor_slow
        end
 #
        flag -=1;
        if flag<0
        flag=0
        end
    ## Else, give it another chance at this temperature, starting from the best solution
    else
        alignment .= alignment_best
        energy = energy_min;
        flag+=1
    end
    ## Test various stop conditions
    if T<T_min
    stop_condition=1
     println("SA terminated -- Temperature reached its minimum")
    elseif n_it_tot>n_it_max
        stop_condition = 1
     println("SA terminated -- Maximum number of iterations")

    elseif flag>flag_max
        stop_condition=1
    println("SA terminated - Could not find better solution OR acceptance rate too low $(flag).")
        end
    return (stop_condition,alignment,n_it_current,flag,energy,T);
end

function shuffle_alignment(alignment::Array{Int64,2},T::Float64,d_curr::Array{Float64,2}, d_target::Array{Float64,2})

 (M,L)=size(alignment)
 n_it_tot=0
 n_it_curr = 10000
 n_it_max=2*M*L
 stop_condition = 0
 stop_condition_max = 100
 err_l2 = compute_l2_error(d_curr, d_target)
 d_pos=zeros(M,M)

    while stop_condition<stop_condition_max
        n_swap_done = 0
        for i in 1:n_it_curr
         (add,energy,d_curr) = do_swap(alignment,d_pos,d_target,d_curr,T,err_l2)
             n_swap_done = n_swap_done + add;
        end
        acceptance_rate = n_swap_done/(n_it_curr-1)
        n_it_tot = n_it_tot + n_it_curr;
        err_l2 = compute_l2_error(d_curr,d_target);
        if acceptance_rate < 0.99
            T = 2*T
        else
            T = T / 1.1
        end
        println("n_it_tot:$(n_it_tot)")
        if n_it_tot>n_it_max
            stop_condition=stop_condition_max
        end
    end
    println("$(n_it_tot) iterations for shuffling\n")
    return T,d_curr
end

### run simmulating annealing with a linear cooling scheme with two slopes
function run_SA(align::Array{Int64,2},T_init::Float64,T_min::Float64,T_factor_slow::Float64,T_factor_fast::Float64,shuffle::Int64,num_iter_max::Int64,d_curr::Array{Float64,2},d_target::Array{Float64,2})

    (M,L)=size(align)
    T = T_init
    energy=0.0
    energy_min=0.0
    energy_min_prev=100.0
    acceptance_rate = 0.0;
    n_it_max = num_iter_max
    n_it_current = 1
    n_it_tot = 0
    stop_condition = 0
    flag=0
    ## Initial state
    println("runMC\n")
    alignment_best=align[1:end,1:end]
    err_l2 = compute_l2_error(d_curr,d_target)
    err_l2_min = err_l2
    energy = Energy(d_curr, d_target)

    energy_min = energy
    #println("Initial energy : $(energy) -- Initial l2 error: $(err_l2)")

    d_i=d_curr[1:end,1:end]

    if shuffle!=0

        println("Shuffling system - T = $(T) ...")
        (T,d_curr)=shuffle_alignment(align,T,d_i,d_target)
        err_l2 = compute_l2_error(d_curr, d_target)
        err_l2_min = err_l2
        energy = Energy(d_curr, d_target)
        energy_min = energy;
        println("Done.");
        println("After shuffling : energy = $(energy) -- l2error=$(err_l2) -- T =$(T)");
    end

    ## Sampling
    println("\nSampling ... \n")
        d_pos=zeros(M,M)

    while stop_condition==0
     #Calculation
    println("n_it_current: $(n_it_current)")
    println(" T : $(T)");
    println("energy : $(energy)");
    #println("energy 2 : $(Energy(d_target,d_curr))")
    println("acceptance_rate : $(acceptance_rate)");
      n_swap_done = 0;
           for i in 1:n_it_current

         (add,energy,d_curr) = do_swap(align,d_pos,d_target,d_curr,T,energy)

        n_swap_done = n_swap_done + add;
             if energy < energy_min
             energy_min_prev = energy_min
             energy_min = energy
             acceptance_rate = n_swap_done/i;
                if acceptance_rate < 0.001
             alignment_best .= align
                end
             end
          end

        acceptance_rate =  n_swap_done/n_it_current
        n_it_tot = n_it_tot + n_it_current
        err_l2 = compute_l2_error(d_curr,d_target)

        if err_l2<err_l2_min
          err_l2_min = err_l2
        end
        ## Just to have alignment_best initialized
        if acceptance_rate > 0.001
            alignment_best .= align
         end
  println("n_it_tot:$(n_it_tot)")

 ## Updating parameters
 (stop_condition,align,n_it_current,flag,energy,T)=update_param_SA(align,alignment_best,n_it_tot,n_it_current,flag,acceptance_rate, energy_min,energy_min_prev,energy,T,T_min,n_it_max,T_factor_slow,T_factor_fast)

    end

    println(" \ndone\n");
    println("n_it_tot:$(n_it_tot)")
    println("Final l2 error = $(err_l2_min)")

  return align

end

####### given a MSA in numeric format this function provide a randomized MSA according to Null-model I
function sample_from_Null_model_I(infile::String;outfile="",shuffle_temp=10.0)
  msa=translate_fasta_to_num_matrix(infile)
  println("Shuffling alignment")

  dtarget=PairwiseHammingDist(msa);
  shuffle_alignment(msa,shuffle_temp,dtarget,dtarget)
  msa_fasta=transform_MSA_fasta(msa)
  if !isempty(outfile)
  export_fasta_file(outfile,msa_fasta)
  println("randomized MSA has been exported with fasta format!!")
  end

 return msa_fasta
end

####### given a MSA in numeric format this function provide a randomized MSA according to Null-model II
function sample_from_Null_model_II(infile::String;outfile="",target_dist="",shuffle=1,shuffle_temp=10.0,T_factor_slow=0.8,T_factor_fast=0.1,min_temp="default",num_iter_max=20000000)
  msa=translate_fasta_to_num_matrix(infile)
   if !isempty(target_dist)
     dtarget=readdlm(target_dist)
    else
    dtarget=PairwiseHammingDist(msa);
   end

  (M,L)=size(msa)
  if min_temp=="default"
  	t_min=2/(10*M)
  else
  	t_min=min_temp
  end
  println("Parameters:")
  if shuffle==1
     println("Simulating annealing start from shuffled alignment")
  end
  println("minimal temperature  T_min=$(t_min)")
  println("slow cooling rate = $(T_factor_slow), fast cooling rate = $(T_factor_fast)")
  println("maximum number of iterations $(num_iter_max)")

  result=run_SA(msa,shuffle_temp,t_min,T_factor_slow,T_factor_fast,shuffle,num_iter_max,dtarget,dtarget)
  msa_fasta=transform_MSA_fasta(result)
  if !isempty(outfile)
  export_fasta_file(outfile,msa_fasta)
  println("randomized MSA has been exported with fasta format!!")
  end
 return msa_fasta
end

end


text:
#############################################################################
# lambdamin_max.jl
# Handles maximum and minimum eigenvalue of a symmetric positive definite matrix
# (and imposes the constraint that its argument be PSD)
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
export lambdamax, lambdamin

### Lambda max

struct LambdaMaxAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function LambdaMaxAtom(x::AbstractExpr)
        children = (x,)
        m, n = size(x)
        if m == n
            return new(:lambdamax, hash(children), children, (1, 1))
        else
            error("lambdamax can only be applied to a square matrix.")
        end
    end
end

function sign(x::LambdaMaxAtom)
    return Positive()
end

function monotonicity(x::LambdaMaxAtom)
    return (Nondecreasing(),)
end

function curvature(x::LambdaMaxAtom)
    return ConvexVexity()
end

function evaluate(x::LambdaMaxAtom)
    eigvals(evaluate(x.children[1]))[end]
end

lambdamax(x::AbstractExpr) = LambdaMaxAtom(x)

# Create the equivalent conic problem:
#   minimize t
#   subject to
#            tI - A is positive semidefinite
#            A      is positive semidefinite
function conic_form!(x::LambdaMaxAtom, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        m, n = size(A)
        t = Variable()
        p = minimize(t, t*Matrix(1.0I, n, n) - A ⪰ 0)
        cache_conic_form!(unique_conic_forms, x, p)
    end
    return get_conic_form(unique_conic_forms, x)
end

### Lambda min

struct LambdaMinAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function LambdaMinAtom(x::AbstractExpr)
        children = (x,)
        m, n = size(x)
        if m == n
            return new(:lambdamin, hash(children), children, (1,1))
        else
            error("lambdamin can only be applied to a square matrix.")
        end
    end
end

function sign(x::LambdaMinAtom)
    return Positive()
end

function monotonicity(x::LambdaMinAtom)
    return (Nondecreasing(),)
end

function curvature(x::LambdaMinAtom)
    return ConcaveVexity()
end

function evaluate(x::LambdaMinAtom)
    eigvals(evaluate(x.children[1]))[1]
end

lambdamin(x::AbstractExpr) = LambdaMinAtom(x)

# Create the equivalent conic problem:
#   maximize t
#   subject to
#            A - tI is positive semidefinite
#            A      is positive semidefinite
function conic_form!(x::LambdaMinAtom, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        m, n = size(A)
        t = Variable()
        p = maximize(t, A - t*Matrix(1.0I, n, n) ⪰ 0)
        cache_conic_form!(unique_conic_forms, x, p)
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
using LoopVectorization
using MuladdMacro

export jacobi, jacobi_series, jacobi_sum, jacobi_weight
"""
    jacobi_weight(α, β, x)

Compute the weighting function (1-x)^α * (1-x)^β .
"""
function jacobi_weight(α, β, x)
	return (1 - x)^α * (1 - x)^β
end

# calcac_startb computes the a and c terms of the recurrence relation for the
# jacobi polynomials.  It also computes three pieces of b which do not depend on
# x.  This is not a formal docstring because documeter complains if I document
# the function and do not export it.

function calcac_startb(n, α, β)
    a = (2n) * (n + α + β) * (2n + α + β - 2)
    c = 2 * (n + α - 1) * (n + β - 1) * (2n + α + β)
    b1 = (2n + α + β - 1)
    b2 = (2n + α + β)
    b2 = b2 * (b2 - 2)
    b3 = α^2 - β^2
    return a, c, b1, b2, b3
end

"""
    jacobi(n, α, β, x)

Compute the Jacobi polynomial of order n and weights α,β at point x.

This function uses a recurrence relation and is numerical stable to very high
order.  The computation time is linear w.r.t. the order n.  jacobi(n, a, b, x)
is also linear w.r.t. the size of argument x.  x should be passed as an array,
and jacobi should not be called with dot syntax, for best performance.

See also: [`jacobi_series`](@ref), [`jacobi_sum`](@ref)
"""
function jacobi(n, α, β, x)
    # the body of this function has largely been copy/pasted to jacobi_sequence
    # and jacobi_sum.  Changes should be propagated.
	if n == 0
        return ones(size(x))
    elseif n == 1
        return (α + 1) .+ (α + β + 2) .* ((x .- 1) ./ 2)
    end
	Pnm1 = (α + 1) .+ (α + β + 2) .* ((x .- 1) ./ 2)
	a, c, b1, b2, b3 = calcac_startb(2, α, β)
	inva = 1 / a
	# .- c .* Pnm2, but Pnm2 == 1, so drop it
	Pn = similar(Pnm1)
    Pn .= @muladd (b1 .* (b2 .* x .+ b3) .* Pnm1 .- c) .* inva
    if n == 2
        return Pn
    end
	Pnm2 = similar(Pn)
	for i = 3:n
		# assign Pnm2 to Pn to reuse a buffer
		Pnm2, Pnm1, Pn = Pnm1, Pn, Pnm2
		a, c, b1, b2, b3 = calcac_startb(i, α, β)
        inva = 1 / a
		@avx unroll = 4 for i in eachindex(Pn)
			Pn[i] = (b1 * (b2 * x[i] + b3) * Pnm1[i] - c * Pnm2[i]) * inva
		end
    end
	return Pn
end

"""
    jacobi_series(ns, α, β, x)

Compute a series of Jacobi polynomials of orders n.  Returns an array with shape
(size(x)..., length(ns)).  That is, the _final_ dimension contains the modes
and the first dimension(s) are spatial.

See also: [`jacobi`](@ref), [`jacobi_sum`](@ref)
"""
function jacobi_series(ns, α, β, x)
    # the body of this function is largely copied from jacobi.  Changes to
    # either should be propagated.
    # allocate the output buffer
    out = Array{eltype(x), ndims(x)+1}(undef, size(x)..., length(ns))
    plane_idx = ndims(x)+1
    lowest_plane = 1
    if ns[1] == 0
        # set the first element of out to ones, if the first n the user wants
        # is n=0
        fill!(selectdim(out, plane_idx, lowest_plane), 1.)
        lowest_plane += 1
    end

    if lowest_plane > length(ns)
        return out
    end
    # seed the recurrence relation
	Pnm1 = (α + 1) .+ (α + β + 2) .* ((x .- 1) ./ 2)
	a, c, b1, b2, b3 = calcac_startb(2, α, β)
	inva = 1 / a

    # unlike the regular jacobi function, we do not need to allocate
    # a scratch buffer, since we can use the lowest plane of the output
    # bufffer with impunity.  This only works because there is no weighting
    # so the lowest plane is guaranteed to contain a "pure" Jacobi polynomial
    # the exception is one buffer for Pnm2 to start up since we don't know if
    # the user opted to keep it

    # .- c .* Pnm2, but Pnm2 == 1, so drop it
    Pn = @muladd (b1 .* (b2 .* x .+ b3) .* Pnm1 .- c) .* inva
    Pnm2 = similar(Pn)

    if ns[lowest_plane] == 1
        view = selectdim(out, plane_idx, lowest_plane)
        view .= Pnm1
        lowest_plane += 1
    end
    if lowest_plane > length(ns)
        return out
    end
	if ns[lowest_plane] == 2
        view = selectdim(out, plane_idx, lowest_plane)
        view .= Pn
        lowest_plane += 1
    end
    if lowest_plane > length(ns)
        return out
    end

    # now do the loop that avoids recurrence
	for i = 3:ns[end]
        # Pnm2, Pnm1, Pn are purposfully not part of out
		Pnm2, Pnm1, Pn = Pnm1, Pn, Pnm2
		a, c, b1, b2, b3 = calcac_startb(i, α, β)
        inva = 1 / a
		@avx unroll = 4 for i in eachindex(Pn)
			Pn[i] = (b1 * (b2 * x[i] + b3) * Pnm1[i] - c * Pnm2[i]) * inva
        end
        for j in eachindex(ns)
			if ns[j] == i
                # this order is relevant to the sum
                view = selectdim(out, plane_idx, lowest_plane)
                view .= Pn
                lowest_plane += 1
                break
			end
		end
    end
	return out
end

"""
    jacobi_sum(ns, weights, α, β, x)

Compute a sum of Jacobi polynomial of order n weighted by weights.

See also: [`jacobi`](@ref), [`jacobi_series`](@ref)
"""
function jacobi_sum(ns, weights, α, β, x)
    # the body of this function is largely copied from jacobi.  Changes to
    # either should be propagated.

	# start by allocating the output array, be a little clever
	# around the initialization value
	out = similar(x)
	if ns[1] == 0
		fill!(out, weights[1])
	else
		fill!(out, 0)
	end

	# Pnm1 = P_1, Pnm2= P_0, == 1
	# this block of setup is the same as before
	Pnm1 = (α + 1) .+ (α + β + 2) .* ((x.-1)./2)
	a, c, b1, b2, b3 = calcac_startb(2, α, β)
	inva = 1 / a
	# .- c .* Pnm2, but Pnm2 == 1, so drop it
    Pn = @muladd (b1 .* (b2 .* x .+ b3) .* Pnm1 .- c) .* inva

	# add P_2 if we can.  Two iffs is faster than
	# a searchsorted
	if ns[1] == 1
		out .+= weights[1] .* Pnm1
	elseif ns[2] == 1
		out .+= weights[2] .* Pnm1
	end
	if ns[1] == 2
		out .+= weights[1] .* Pn
	elseif ns[2] == 2
		out .+= weights[2] .* Pn
	end

	Pnm2 = similar(Pn);
	for i = 3:ns[end]
		# assign Pnm2 to Pn to reuse a buffer
		Pnm2, Pnm1, Pn = Pnm1, Pn, Pnm2
		a, c, b1, b2, b3 = calcac_startb(i, α, β)
        inva = 1 / a
		@avx unroll=4 for i in eachindex(Pn)
			Pn[i] = (b1 * (b2 * x[i] + b3) * Pnm1[i] - c * Pnm2[i]) * inva
		end
		# i is the order that was just computed
		for j in eachindex(ns)
			if ns[j] == i
				# this order is relevant to the sum
				out .+= weights[j] .* Pn
				break
			end
		end
    end
	return out
end


text:
using Tokens
using Test

function testpacked31()
    @testset "Packed31" begin
        p0_0 = Packed31(0%UInt32)
        p15_255 = Packed31(15%UInt8, 255%UInt32)
        p15_256 = Packed31(15%UInt8, 256%UInt32)
        p15_257 = Packed31(15%UInt8, 257%UInt32)
        @test p0_0.nibble == Nibble(0)
        @test p15_255.nibble == Nibble(15)
        @test p15_256.nibble == Nibble(15)
        @test p15_257.nibble == Nibble(15)
        @test p15_255.len == 255
        @test p15_256.len == 256
        @test p15_257.len == 257
        p15_max = Packed31(15%UInt8,UInt32(2^27 - 1))
        @test UInt32(p15_max) == UInt32(2^31 - 1)

        @test p0_0.nibble == Nibble(0)
        @test p0_0.nibble == Nibble(0)

    end
end


text:
function Pair(arg0::Object, arg1::Object)
    return Pair((Object, Object), arg0, arg1)
end

function Pair(arg0::Pair)
    return Pair((Pair,), arg0)
end

function create(::Type{Pair}, arg0::Object, arg1::Object)
    return jcall(Pair, "create", Pair, (Object, Object), arg0, arg1)
end

function equals(obj::Pair, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_first(obj::Pair)
    return jcall(obj, "getFirst", Object, ())
end

function get_key(obj::Pair)
    return jcall(obj, "getKey", Object, ())
end

function get_second(obj::Pair)
    return jcall(obj, "getSecond", Object, ())
end

function get_value(obj::Pair)
    return jcall(obj, "getValue", Object, ())
end

function hash_code(obj::Pair)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::Pair)
    return jcall(obj, "toString", JString, ())
end


text:
using Luna
Luna.set_fftw_mode(:estimate)
import Logging: @warn, disable_logging, Info
import PyPlot: plt

disable_logging(Info) # disable progress logging to avoid filling the shell window

function run(item)
    if isfile(item) && endswith(item, ".jl")
        println("============ Running example $item")
        try
            include(item)
        catch e
            bt = catch_backtrace()
            @warn sprint(showerror, e, bt)
        finally
            plt.close("all")
        end
    elseif isdir(item)
        for iitem in readdir(item)
            run(joinpath(item, iitem))
        end
    end
end

run(joinpath(Utils.lunadir(), "examples"))


text:
include("get_div_grad.jl")

function test_cgs()
  cgs_tol = 1.0e-6

  # Symmetric and positive definite systems (cubic spline matrix).
  n = 10
  A = spdiagm(-1 => ones(n-1), 0 => 4*ones(n), 1 => ones(n-1))
  b = A * [1:n;]
  (x, stats) = cgs(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)

  # Symmetric indefinite variant.
  A = A - 3 * I
  b = A * [1:n;]
  (x, stats) = cgs(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)

  # Nonsymmetric and positive definite systems.
  A = [i == j ? 10.0 : i < j ? 1.0 : -1.0 for i=1:n, j=1:n]
  b = A * [1:n;]
  (x, stats) = cgs(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)

  # Nonsymmetric indefinite variant.
  A = [i == j ? 10*(-1.0)^(i*j) : i < j ? 1.0 : -1.0 for i=1:n, j=1:n]
  b = A * [1:n;]
  (x, stats) = cgs(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)

  # Code coverage.
  (x, stats) = cgs(sparse(A), b)
  show(stats)

  # Sparse Laplacian.
  A = get_div_grad(16, 16, 16)
  b = ones(size(A, 1))
  (x, stats) = cgs(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)

  # Test b == 0
  (x, stats) = cgs(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = spdiagm(-1 => ones(Int, n-1), 0 => 4*ones(Int, n), 1 => ones(Int, n-1))
  b = A * [1:n;]
  (x, stats) = cgs(A, b)
  @test stats.solved

  # Test with Jacobi (or diagonal) preconditioner
  A = ones(10,10) + 9 * I
  b = 10 * [1:10;]
  M = 1/10 * opEye(10)
  (x, stats) = cgs(A, b, M=M)
  show(stats)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CGS: Relative residual: %8.1e\n", resid)
  @test(resid ≤ cgs_tol)
  @test(stats.solved)
end

test_cgs()


text:
using Mayday
using JuMP
using Base.Test

include("../examples/simple_sos_problem.jl")

include("../examples/spotopt_van_der_pol.jl")
test_spotopt_van_der_pol(Mayday.chebyshev_polynomial_first_kind)

include("degrees.jl")
test_variable_degrees()

include("chebyshev.jl")
test_univariate_chebyshev_polynomials_range()


text:
# Test comparison for Kepler-11 example
const DAYS_IN_YEAR = 365.2425
const SR_TO_AU = 0.0046491;   #  /* http://en.wikipedia.org/wiki/Solar_radius */
if VERSION >= v"0.7"
  using Pkg, Libdl
  #global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".",joinpath(Pkg.devdir(),"CORBITS")])
  global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[joinpath(dirname(pathof(CORBITS)),"..")])
else
  global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".","..","../..","/usr/local/lib"])
end

# Radius and Mass of Kepler-11 from kepler.nasa.gov
R11 = 1.10;
M11 = 0.95;
# Periods hardwired rather than from file
P = [ 10.30375, 13.02502, 22.68719, 31.99590, 46.68876, 74.34319, 118.37774 ]
use = ones(Cint,length(P))
use[6] = 0
a = cbrt.((P./DAYS_IN_YEAR).^2 .* M11)
r_star = R11 * SR_TO_AU
r = zeros(length(P))
e = zeros(length(P))
omega = zeros(length(P))
inc = zeros(length(P))
Omega = zeros(length(P))

function prob_of_transits_approx(a::Vector{Cdouble},r_star::Cdouble,r::Vector{Cdouble}, e::Vector{Cdouble},
                                 Omega::Vector{Cdouble}, omega::Vector{Cdouble}, inc::Vector{Cdouble}, use::Vector{Cint})
  @assert(length(a) == length(r) == length(e) == length(Omega) == length(omega) == length(inc) == length(use) )
  @assert(length(a) >=1 )
  local num_pl = length(a)
  return ccall( (:prob_of_transits_approx_arrays, LIB_CORBITS), Cdouble,
      (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint),
       a, r_star, r, e, Omega, omega, inc, use, num_pl)
end

function rand_Rayleigh(sigma::Real)
    # // inverse CDF based on http://en.wikipedia.org/wiki/Rayleigh_distribution
    return sqrt(-2.0*log(rand()) * sigma^2);
end

function prob( sigma_i::Cdouble )
  local sum = 0.0
  local NTRIALS = 10000
  for i in 1:NTRIALS
    #sum_k = 0.0
    for k in 1:length(P)
      for j in 1:length(P)
        inc[j] = rand_Rayleigh(sigma_i) * pi/180.0
        Omega[j] = 2pi*rand()
        use[j] = (k!=j) ? 1 : 0
      end # j
      val = prob_of_transits_approx(a,r_star,r,e,Omega,omega,inc,use)
      sum += val
    end # k
  end # i
  sum /= NTRIALS
  return sum
end # prob

# Print outputs for comparison to C++ version on Kepler-11 example
#for i in linspace(0.0,6.0,61)
for i in range(0.0,stop=6.0,length=61)
   curp = prob(i *sqrt(2.0/pi) )
   print(i, ' ', curp, '\n')
end


text:
# Load DataFrame package
require("DataFrames")
using DataFrames

# Load a very simple lm() function
require(Pkg.dir("DataFrames", "demo", "lm.jl"))

# Load a CSV file into a DataFrame
df = read_table(Pkg.dir("DataFrames", "demo", "toy_example.csv"))

# Run a simple linear model that predicts A using B and C
lm_fit = lm(:(A ~ B + C), df)

# Print out a summary of the results
print(lm_fit)

#
# Behind The Scenes
#

# Generate a Formula object
f = Formula(:(A ~ B + C))

# Generate a ModelFrame object
mf = model_frame(f, df)

# Generate a ModelMatrix object
mm = model_matrix(mf)