text:
using Documenter,DiffEqBase,DiffEqPDEBase,DiffEqProblemLibrary,DiffEqBiological

makedocs(modules=[DiffEqBase,DiffEqPDEBase,DiffEqProblemLibrary,DiffEqBiological],
         doctest=false, clean=true,
         format = :html,
         sitename="DifferentialEquations.jl",
         authors="Chris Rackauckas",
         pages = Any[
         "Home" => "index.md",
         "Tutorials" => Any[
           "tutorials/ode_example.md",
           "tutorials/sde_example.md",
           "tutorials/rode_example.md",
           "tutorials/dde_example.md",
           "tutorials/dae_example.md",
           "tutorials/discrete_stochastic_example.md",
           "tutorials/jump_diffusion.md",
           "tutorials/bvp_example.md",
           "tutorials/additional.md"
         ],
         "Basics" => Any[
           "basics/overview.md",
           "basics/common_solver_opts.md",
           "basics/solution.md",
           "basics/plot.md",
           "basics/integrator.md",
           "basics/problem.md",
           "basics/faq.md",
           "basics/compatibility_chart.md"
         ],
         "Problem Types" => Any[
           "types/discrete_types.md",
           "types/ode_types.md",
           "types/dynamical_types.md",
           "types/split_ode_types.md",
           "types/steady_state_types.md",
           "types/bvp_types.md",
           "types/sde_types.md",
           "types/rode_types.md",
           "types/dde_types.md",
           "types/dae_types.md",
           "types/jump_types.md",
         ],
         "Solver Algorithms" => Any[
           "solvers/discrete_solve.md",
           "solvers/ode_solve.md",
           "solvers/dynamical_solve.md",
           "solvers/split_ode_solve.md",
           "solvers/steady_state_solve.md",
           "solvers/bvp_solve.md",
           "solvers/jump_solve.md",
           "solvers/sde_solve.md",
           "solvers/rode_solve.md",
           "solvers/dde_solve.md",
           "solvers/dae_solve.md",
           "solvers/benchmarks.md"
         ],
         "Additional Features" => Any[
           "features/performance_overloads.md",
           "features/diffeq_arrays.md",
           "features/diffeq_operator.md",
           "features/noise_process.md",
           "features/linear_nonlinear.md",
           "features/callback_functions.md",
           "features/callback_library.md",
           "features/monte_carlo.md",
           "features/io.md",
           "features/low_dep.md",
           "features/progress_bar.md"
         ],
         "Analysis Tools" => Any[
           "analysis/parameterized_functions.md",
           "analysis/parameter_estimation.md",
           "analysis/bifurcation.md",
           "analysis/sensitivity.md",
           "analysis/global_sensitivity.md",
           "analysis/uncertainty_quantification.md",
           "analysis/neural_networks.md",
           "analysis/dev_and_test.md"
         ],
         "Domain Modeling Tools" => Any[
             "models/multiscale.md",
             "models/physical.md",
             "models/financial.md",
             "models/biological.md",
             "models/external_modeling.md"
         ],
         "APIs" => Any[
             "apis/diffeqbio.md"
         ],
         "Extra Details" => Any[
             "extras/timestepping.md",
         ]
         ])

deploydocs(
   repo = "github.com/JuliaDiffEq/DiffEqDocs.jl.git",
   target = "build",
   osname = "linux",
   julia = "1.1",
   deps = nothing,
   make = nothing)


text:
@testset "Testing Transformation Submodule" begin

N = 256
I = reshape(shepp_logan(N),N,N,1)

IInterp = interpolateToGrid(I, [1.0,1.0,1.0]./256, [0.0,0.0,0.0], [0.0,0.0,1.0],
                            [2.0,2.0,2.0], [0.4,0.4,0.0], [400,400,1])
exportImage("img/trafo1.png", IInterp[:,:,1])
@testImg("trafo1.png")

IA = makeAxisArray(reshape(I,N,N,1), [1.0,1.0,1.0]./256, -[0.5,0.5,0.0].+[1.0,1.0,0.0]./256 )
IInterp2 = interpolateToGrid(IA, [0.0,0.0,1.0], [2.0,2.0,2.0], [0.4,0.4,0.0], [400,400,1])
exportImage("img/trafo2.png", IInterp2[:,:,1])
@testImg("trafo2.png")

IM = ImageMeta(IA, Dict{Symbol,Any}(:rotation=>[0.0,0.0,1.0]))
IInterp3 = interpolateToGrid(IM, [2.0,2.0,2.0], [0.4,0.4,0.0], [400,400,1])
exportImage("img/trafo3.png", IInterp3[:,:,1])
@testImg("trafo3.png")

end

ext: jl


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [3969*x1^2*x7^2-630*x1^2*x7*x9+25*x1^2*x9^2-4158*x1*x3^2*x7+330*x1*x3^2*x9+7308*x1*x3*x5*x7-580*x1*x3*x5*x9+4662*x1*x7*x9^2-370*x1*x9^3+64*x2^2*x8^2-848*x2*x3*x5*x8-784*x2*x8*x9^2+1089*x3^4-3828*x3^3*x5+6173*x3^2*x5^2-2442*x3^2*x9^2+9486*x3*x5*x9^2+3770*x9^4-128*x2^2*x8+848*x2*x3*x5-128*x2*x5*x8+784*x2*x9^2+848*x3*x5^2+784*x5*x9^2-10584*x1*x7+840*x1*x9+64*x2^2+128*x2*x5+1424*x2*x8+5544*x3^2-19178*x3*x5+64*x5^2-14938*x9^2-1424*x2-1424*x5+14977, 128*x2*x8^2-848*x3*x5*x8-784*x8*x9^2-256*x2*x8+848*x3*x5-128*x5*x8+784*x9^2+128*x2+128*x5+1424*x8-1424, -8316*x1*x3*x7+660*x1*x3*x9+7308*x1*x5*x7-580*x1*x5*x9-848*x2*x5*x8+4356*x3^3-11484*x3^2*x5+12346*x3*x5^2-4884*x3*x9^2+9486*x5*x9^2+848*x2*x5+848*x5^2+11088*x3-19178*x5, 0, 7308*x1*x3*x7-580*x1*x3*x9-848*x2*x3*x8-3828*x3^3+12346*x3^2*x5+9486*x3*x9^2+848*x2*x3-128*x2*x8+1696*x3*x5+784*x9^2+128*x2-19178*x3+128*x5-1424, 0, 7938*x1^2*x7-630*x1^2*x9-4158*x1*x3^2+7308*x1*x3*x5+4662*x1*x9^2-10584*x1, 128*x2^2*x8-848*x2*x3*x5-784*x2*x9^2-128*x2^2-128*x2*x5+1424*x2, -630*x1^2*x7+50*x1^2*x9+330*x1*x3^2-580*x1*x3*x5+9324*x1*x7*x9-1110*x1*x9^2-1568*x2*x8*x9-4884*x3^2*x9+18972*x3*x5*x9+15080*x9^3+1568*x2*x9+1568*x5*x9+840*x1-29876*x9];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module GRAPE

include("workspace.jl")
include("result.jl")
include("optimize.jl")
include("backend_optim.jl")
include("backend_lbfgsb.jl")

import QuantumControlBase

"""
```julia
opt_result = optimize(problem; method=:GRAPE, kwargs...)
```

optimizes [`problem`](@ref QuantumControlBase.ControlProblem)
using GRadident Ascent Pulse Engineering (GRAPE), see
[`GRAPE.optimize_grape`](@ref).
"""
QuantumControlBase.optimize(problem, method::Val{:GRAPE}) = optimize_grape(problem)
QuantumControlBase.optimize(problem, method::Val{:grape}) = optimize_grape(problem)

end


text:
__precompile__(true)

module IntervalContractors

export plus_rev, minus_rev, inv_rev,
       mul_rev, div_rev, power_rev,
       sign_rev, max_rev, min_rev,
        sqr_rev, sqrt_rev, abs_rev,
        exp_rev, exp2_rev, exp10_rev, expm1_rev,
        log_rev, log2_rev, log10_rev, log1p_rev,
        sin_rev, cos_rev, tan_rev,
        asin_rev, acos_rev, atan_rev,
        sinh_rev, cosh_rev, tanh_rev,
        asinh_rev, acosh_rev, atanh_rev,
        mul_rev_IEEE1788

using IntervalArithmetic

const half_pi = IntervalArithmetic.half_pi(Float64)  # interval
const two_pi = IntervalArithmetic.two_pi(Float64)  # interval

#
# Base.:∪(f::Function, g::Function) = X -> ( f(X) ∪ g(X) )
# Base.:∩(f::Function, g::Function) = X -> ( f(X) ∩ g(X) )  # or f(g(X)) for contractors

include("arithmetic.jl")
include("transformations.jl")
include("powers.jl")
include("exponential.jl")
include("trig.jl")
include("inverse_trig.jl")
include("hyperbolic.jl")
include("inverse_hyperbolic.jl")
include("extrema.jl")

"""
Dictionary mapping functions to their reverse functions.
"""
const reverse_operations = Dict(
                    :+     => :plus_rev,
                    :-     => :minus_rev,
                    :*     => :mul_rev,
                    :/     => :div_rev,
                    :^     => :power_rev,
                    :()    => :()   # empty operation
                    )

for f in (  :sqrt, :sqr, :abs,
            :exp, :exp2, :exp10, :expm1,
            :log, :log2, :log10, :log1p,
            :sin, :cos, :tan,
            :asin, :acos, :atan,
            :sinh, :cosh, :tanh,
            :asinh, :acosh, :atanh,
            :inv, :sign, :max, :min)
    reverse_operations[f] = Symbol(f, "_rev")
end

end


text:
# TEST:
function priorsample(y=range(-8, 4, length=200))
  tau = rand(Gamma(.5, 1))
  omega = rand(InverseGamma(2.5, tau))
  psi = rand(Normal(-1, .5))
  mu = rand(Normal(1, 1))
  sigma = CDE.Util.scalefromaltskewt(sqrt(omega), psi)
  phi = CDE.Util.skewfromaltskewt(sqrt(omega), psi)
  nu = rand(LogNormal(3, .5))
  println("(tau,omega,psi): $((tau,omega,psi))")
  println("(sigma,phi): $((sigma,phi))")
  return y, pdf.(CDE.Util.SkewT(mu, sigma, nu, phi), y)
end
dopriorsample() = let
  _y, ypdf = priorsample()
  plot(_y, ypdf, label=nothing); savefig("results/tmp.pdf"); closeall()
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module Replay

using Crayons

include("FakePTYs.jl")
using .FakePTYs: open_fake_pty

const CTRL_C = '\x03'
const UP_ARROW = "\e[A"
const DOWN_ARROW = "\e[B"
const RIGHT_ARROW = "\e[C"
const LEFT_ARROW = "\e[D"

export CTRL_C, UP_ARROW, DOWN_ARROW, RIGHT_ARROW, LEFT_ARROW
export replay

function clearline(; move_up::Bool=false)
    buf = IOBuffer()
    print(buf, "\x1b[2K") # clear line
    print(buf, "\x1b[?25l") # hide cursor
    print(buf, "\x1b[999D") # rollback the cursor
    move_up && print(buf, "\x1b[1A") # move up
    print(buf |> take! |> String)
end

function clearlines(H::Integer)
    for _ in 1:H
        clearline(move_up=true)
    end
end

function type_with_ghost_core(line::AbstractString, mode; display_prompt=false)
    if !display_prompt
        # we assume we're in julian mode
        spacestring = " "
        dummy = repeat(spacestring, length("julia> "))
        line = dummy * line
    end
    clearline()
    for index in collect(eachindex(line))
        if display_prompt
            if Base.get_have_color()
                print(mode.prompt_prefix, mode.prompt, mode.prompt_suffix)
            else
                print(mode.prompt)
            end
        end
        println(join(line[begin:index]))
        clearline(move_up=true)
        duration = if 30 < length(line)
            0.0125
        elseif 15 < length(line) < 30
            0.05
        else
            0.1
        end
        sleep(duration)
    end
end

function type_with_ghost(repl_script::AbstractString, mode)
    lines = split(String(repl_script), '\n'; keepempty=false)
    H = length(lines)
    for (i, line) in enumerate(lines)
        display_prompt = (i == 1)
        type_with_ghost_core(line, mode; display_prompt)
        println()
    end
    clearlines(H)
end

function setup_pty(julia_project="@."::AbstractString, cmd=String = "--color=yes")
    pts, ptm = open_fake_pty()
    blackhole = Sys.isunix() ? "/dev/null" : "nul"
    julia_exepath = joinpath(Sys.BINDIR, Base.julia_exename())
    replproc = withenv(
        "JULIA_HISTORY" => blackhole,
        "JULIA_PROJECT" => "$julia_project",
        "TERM" => "",
    ) do
        # Install packages
        run(`$(julia_exepath) -e 'using Pkg; Pkg.instantiate()'`)
        # Initialize REPL
        run(
            ```$(julia_exepath) -i -e 'print("\x1b[?25l")' $(split(cmd))```,
            pts,
            pts,
            pts;
            wait=false,
        )
    end
    Base.close_stdio(pts)
    return replproc, ptm
end

function replay(
    instructions::Vector{<:AbstractString},
    buf::IO=stdout;
    use_ghostwriter=false,
    julia_project="@.",
    cmd=String = "--color=yes",
)
    print("\x1b[?25l") # hide cursor
    replproc, ptm = setup_pty(julia_project, cmd)
    # Prepare a background process to copy output from process until `pts` is closed
    output_copy = Base.BufferStream()
    tee = @async try
        while !eof(ptm)
            l = readavailable(ptm)
            write(buf, l)
            Sys.iswindows() && (sleep(0.1); yield(); yield()) # workaround hang - probably a libuv issue?
            write(output_copy, l)
        end
        close(output_copy)
        close(ptm)
    catch ex
        close(output_copy)
        close(ptm)
        if !(ex isa Base.IOError && ex.code == Base.UV_EIO)
            rethrow() # ignore EIO on ptm after pts dies
        end
    end
    # wait for the definitive prompt before start writing to the TTY
    readuntil(output_copy, "julia>")
    sleep(0.1)
    readavailable(output_copy)

    name = :julian
    prompt = "julia> "
    prompt_prefix = Base.text_colors[:bold] * Base.text_colors[:green]
    prompt_suffix = Base.color_normal
    mode = (; name, prompt, prompt_prefix, prompt_suffix)
    current_mode_name = :julian

    # let's replay!
    for cell in instructions
        sleep(1)

        ghost_script = cell
        if current_mode_name == :julian
            mode = if startswith(cell, ';')
                ghost_script = cell[begin+1:end] # remove ';'
                # shell mode
                name = :shell
                prompt = "shell> "
                prompt_prefix = Base.text_colors[:bold] * Base.text_color[:red]
                prompt_suffix = Base.color_normal
                current_mode_name = name
                (; name, prompt, prompt_prefix, prompt_suffix)
            elseif startswith(cell, ']')
                ghost_script = cell[begin+1:end] # remove ']'
                # pkg repl
                name = :repl
                prompt_prefix = Base.text_colors[:bold] * Base.text_colors[:blue]
                active_project_dir, _ =
                    splitext(Base.active_project() |> dirname |> basename)
                if occursin(r"v[0-9].[0-9]", active_project_dir)
                    prompt = "(@v$(VERSION.major).$(VERSION.minor)) pkg> "
                else
                    prompt = "($(active_project_dir)) pkg> "
                end
                prompt_suffix = Base.color_normal
                current_mode_name = name
                (; name, prompt, prompt_prefix, prompt_suffix)
            elseif startswith(cell, '?')
                ghost_script = cell[begin+1:end] # remove '?'
                # help mode
                name = :help
                prompt = "help?> "
                prompt_prefix = Base.text_colors[:bold] * Base.text_colors[:yellow]
                prompt_suffix = Base.color_normal
                # help mode should back to julian mode
                current_mode_name = :julian
                (; name, prompt, prompt_prefix, prompt_suffix)
            else
                # julian mode
                name = :julian
                prompt = "julia> "
                prompt_prefix = Base.text_colors[:bold] * Base.text_colors[:green]
                prompt_suffix = Base.color_normal
                current_mode_name = :julian
                (; name, prompt, prompt_prefix, prompt_suffix)
            end
        else
            if endswith(cell, CTRL_C)
                # julian mode
                name = :julian
                prompt = "julia> "
                prompt_prefix = Base.text_colors[:bold] * Base.text_colors[:green]
                prompt_suffix = Base.color_normal
                current_mode_name = :julian
                mode = (; name, prompt, prompt_prefix, prompt_suffix)
            end
        end

        bytesavailable(output_copy) > 0 && readavailable(output_copy)

        use_ghostwriter && type_with_ghost(ghost_script, mode)

        if endswith(cell, CTRL_C)
            write(ptm, cell)
        else
            if length(split(string(cell), '\n'; keepempty=false)) == 1
                write(ptm, cell, "\n")
            else
                write(ptm, cell)
            end
        end
        readuntil(output_copy, "\n")
        # wait for the next prompt-like to appear
        # NOTE: this is rather inaccurate because the Pkg REPL mode is a special flower
        readuntil(output_copy, "\n")
        readuntil(output_copy, "> ")
    end

    sleep(1)
    use_ghostwriter && type_with_ghost("exit()", mode)
    write(ptm, "exit()\n")
    sleep(1)
    wait(tee)
    success(replproc) || Base.pipeline_error(replproc)
    close(ptm)
    print("\x1b[?25h") # unhide
    return buf
end

function replay(repl_script::AbstractString, args...; kwargs...)
    replay(split(String(repl_script), '\n'; keepempty=false), args...; kwargs...)
end

end # module


text:
struct EchelleSpectrum{W <: Number,F <: Number} <: AbstractSpectrum{W,F}
    wave::Matrix{<:Real}
    flux::Matrix{<:Real}
    meta::Dict{Symbol,Any}
end

EchelleSpectrum(wave, flux, meta::Dict{Symbol,Any}) = EchelleSpectrum(collect(wave), collect(flux), meta)

function Base.show(io::IO, spec::EchelleSpectrum)
    println(io, "EchelleSpectrum($(eltype(spec.wave)), $(eltype(spec.flux)))")
    print(io, "  # orders: $(size(spec, 1))")
    for (key, val) in spec.meta
        print(io, "\n  $key: $val")
    end
end

function Base.getindex(spec::EchelleSpectrum, i::Integer)
    wave = spec.wave[i, :]
    flux = spec.flux[i, :]
    meta = merge(Dict(:Order => i), spec.meta)
    return Spectrum(wave, flux, meta)
end

function Base.getindex(spec::EchelleSpectrum, I::AbstractVector)
    waves = spec.wave[I, :]
    flux = spec.flux[I, :]
    meta = merge(Dict(:Orders => (first(I), last(I))), spec.meta)
    return EchelleSpectrum(waves, flux, meta)
end

Base.lastindex(spec::EchelleSpectrum) = lastindex(spec.flux, 1)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "WORLD"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mmorise/World.git", "1f8815f7a85e344ea9f98ca6b1383578a161d4e2"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd World/
mkdir -p ${libdir}
$CXX -O3 -Wall -shared -fPIC -o ${libdir}/libworld.${dlext} -Isrc src/*.cpp
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libworld", :libworld)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using LightGraphs, MetaGraphs, LinearAlgebra, SparseArrays, SuiteSparse

function load_matpower(filename)
  instance_name = split(filename, '.')[1]

  touch(instance_name*".temp")
  f = open(filename)
  out = open(instance_name*".temp", "w")

  # removing all ';' at end of lines
  while !eof(f)
    line = readline(f)
    if length(line) > 0 && line[1] != '%' && line[1] != 'f'
      s = split(line, ";")
      println(out, s[1])
    end
  end
  close(f)
  close(out)

  data = DelimitedFiles.readdlm(instance_name*".temp")
  rm(instance_name*".temp")
  return data
end

function find_numarray(i_start, data)
  i_debut = i_start
  while !isa(data[i_debut, 1], Int)
    i_debut+=1
  end
  i_fin=i_debut
  while !isa(data[i_fin,1], SubString)
    i_fin += 1
  end
  i_debut, i_fin-1
end

checkfor(data, line_ind, name) = (data[line_ind, 1] == name) || error("Expected ", name, " at line ", line_ind, ", got ", data[line_ind,1], " instead.")

function read_sparsity_pattern(instance_path::String)

    data = load_matpower(instance_path)

    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut-1, "mpc.bus")
    nb_bus = i_fin-i_debut+1
    index_bus = Dict( data[i+i_debut-1,1] => i for i=1:nb_bus)

    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.gen")

    #initialize network graph G
    sp = spzeros(nb_bus,nb_bus)

    ## Link information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.branch")
    for i=i_debut:i_fin
      if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
      else
          orig = index_bus[data[i,1]]
          dest = index_bus[data[i,2]]
          sp[orig,dest] = 1
      end
    end

    sp_sym = sp + sp'

    diag = zeros(nb_bus)

    for i=1:nb_bus
        sum_col = sum(sp_sym[i,:])
        diag[i] = Int(sum_col + 1)
    end
    return sp_sym + Diagonal(diag)
    # return sp+sp'+nb_bus*sparse(I, nb_bus, nb_bus)
end

using MAT

function read_sparsity_pattern_MAT(instance_path::String)

    file = matopen(instance_path)
    data = read(file, "mpc")
    close(file)

    nb_bus = size(data["bus"],1)
    index_bus = Dict(data["bus"][i,1] => i for i=1:nb_bus)

    #initialize network graph G
    sp = spzeros(nb_bus,nb_bus)

    data_branch = data["branch"]
    for i=1:size(data_branch,1)
      if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
      else
          orig = index_bus[data_branch[i,1]]
          dest = index_bus[data_branch[i,2]]
          sp[orig,dest] = 1
      end
    end

    sp_sym = sp + sp'

    diag = zeros(nb_bus)

    for i=1:nb_bus
        sum_col = sum(sp_sym[i,:])
        diag[i] = Int(sum_col + 1)
    end
    return sp_sym + Diagonal(diag)
    # return sp+sp'+nb_bus*sparse(I, nb_bus, nb_bus)
end

function chordal_ext_cholesky(sparsity_pattern)
    A = sparsity_pattern
    nb_edges_A = (nnz(A) - size(A,1))/2
    #computing cholesky factorisation of A
    F = cholesky(A) #NOTE: order = F.p
    # println(F.p)
    # computing L + LT
    L = sparse(F.L)
    nb_edges_L = nnz(L) - size(A,1)
    nb_added_edges = nb_edges_L - nb_edges_A
    SP = L + L'
    #inverting permutation to get chordal extension of sparsity_pattern
    H = SP[invperm(F.p), invperm(F.p)]
    return H, F.p, nb_added_edges
end

function construct_graph_from_matrix(L)
    n = size(L,1)
    H = MetaGraph(n)
    for i in 1:n
        set_props!(H, i, Dict(:name => "node$i"))
    end
    for i in 1:n
        for j in 1:i
            if L[i,j] != 0
                MetaGraphs.add_edge!(H,i,j)
            end
        end
    end
    return H
end

##########################################################################################################
function cliques_max(G, order)
    n = nv(G)
    C0 = []
    C = C0
    cliques_dict = Dict{Int64,Array{Any,1}}()
    k = 0
    for i=1:n
        v = order[i]
        Nv = neighbors(G,v)
        # println(v)
        C = [Nv[j] for j=1:length(Nv) if Nv[j] ∉ order[1:i]]
        C = [C ; v]
        # println(C)
        if C ⊈ C0
            k +=1
            cliques_dict[k] = C
            C0 = C
        end
    end
    return cliques_dict
end

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

function weighted_graph(cliques_dict)
    nb_cl = length(cliques_dict)
    G = MetaGraph(nb_cl)
    for i in 1:nb_cl
        for j in 1:i
            inter = length(intersect(cliques_dict[i], cliques_dict[j]))
            if inter > 0
                LightGraphs.add_edge!(G,i,j)
                set_prop!(G, i, j, :weight, inter)
            end
        end
    end
    return G
end

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

function Prim_algo(G)
    A = MetaGraph(nv(G))
    for i = 1:nv(G)
        set_prop!(A,i, :name, "C$i")
    end
    V = Set([ v for v in LightGraphs.vertices(G)])
    Vprim = Set([1])
    while Vprim != V
        max_weight = 0
        edge_to_add = LightGraphs.Edge(0,0)
        node_to_add = 0
        for edge in LightGraphs.edges(G)
            u = src(edge)
            v = dst(edge)
            if (u ∈ Vprim && v ∉ Vprim)
                weight_uv = get_prop(G,LightGraphs.Edge(u,v), :weight)
                if weight_uv > max_weight
                    max_weight = weight_uv
                    edge_to_add = LightGraphs.Edge(u,v)
                    node_to_add = v
                end
            elseif (u ∉ Vprim && v ∈ Vprim)
                weight_uv = get_prop(G,LightGraphs.Edge(u,v), :weight)
                if weight_uv > max_weight
                    max_weight = weight_uv
                    edge_to_add = LightGraphs.Edge(u,v)
                    node_to_add = u
                end
            end
        end
        LightGraphs.add_edge!(A, edge_to_add)
        set_prop!(A, edge_to_add, :weight, max_weight)
        push!(Vprim, node_to_add)
    end
    return A
end


text:
__precompile__()

"""
A package for the quick creation of customizable exceptions.

Links:
- Repo: https://github.com/paveloom-j/Exceptions.jl
- Docs: https://paveloom-j.github.io/Exceptions.jl
"""
baremodule Exceptions

"This module contains all inner parts of this package."
baremodule Internal

export @exception

using Base

# Load exceptions
Base.include(Internal, "include/exceptions.jl")

# Load macros
Base.include(Internal, "include/exception.jl")

end

# Export the second-level components
using .Internal

# Export the first-level components
export @exception

end


text:
import Radio
import DSP

M          = 4
modulation = Radio.PSK(M)
modem      = Radio.Modem( modulation, samplesPerSymbol=2 )
txData     = rand( 0:M-1, 1000 )
symbols    = Radio.modulate( modem, txData )
rxData     = Radio.demodulate( modem, symbols )
taps       = modem.txFilter.h

p = DSP.welch_pgram( symbols, window=DSP.hanning )
freqs = fftshift(DSP.freq( p ))
power = log10( fftshift(DSP.power( p )) )

using PyPlot

clf()
subplot( 2,1,1 )
scatter( real(symbols), imag(symbols) )

subplot( 2,1,2 )
plot( freqs, power )


text:
# Take the number 192 and multiply it by each of 1, 2, and 3:
#
# 192 × 1 = 192
# 192 × 2 = 384
# 192 × 3 = 576
#
# By concatenating each product we get the 1 to 9 pandigital, 192384576. We
# will call 192384576 the concatenated product of 192 and (1,2,3)
#
# The same can be achieved by starting with 9 and multiplying by 1, 2, 3, 4,
# and 5, giving the pandigital, 918273645, which is the concatenated product of
# 9 and (1,2,3,4,5).
#
# What is the largest 1 to 9 pandigital 9-digit number that can be formed as
# the concatenated product of an integer with (1,2, ... , n) where n > 1?

using ProjectEulerSolutions

# Simple using Julia's builtin digits function.
function p038solution(n::Integer=10)::Integer
    pandigitals = Integer[]
    # Loop through initial multipliers and grab digits
    for i in 1:n
        digs = digits(i)
        prepend!(digs, digits(i * 2))
        mult = i * 3
        j = 3
        while ndigits(mult) <= 9 - length(digs)
            prepend!(digs, digits(mult))
            j += 1
            mult = i * j
        end
        # Turn into set to make sure there are 9 digits and none are 0
        digitset = Set(digs)
        if length(digitset) == 9 && ! in(0, digitset)
            # Convert back into integer
            append!(pandigitals, sum([digs[k+1] * 10^k for k in 0:8]))
        end
    end
    sort!(pandigitals, rev=true)
    return pandigitals[1]
end

p038 = Problems.Problem(p038solution)

Problems.benchmark(p038, 9999)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
###############
# recreate manual-traitdispatch.jl: f2

@traitfn ttf2{X,Y<:Integer; D1{Y}, D4{X,Y}}(x::X,y::Y) = x + sin(y)
@traitfn ttf2{S,T<:Integer; D1{S}, D1{T}  }(s::S,t::T) = sin(s) - sin(t)
@traitfn ttf2{X,Y<:AbstractFloat; D1{X}, D1{Y}  }(x::X,y::Y) = cos(x) - cos(y)

@test ttf2(4,5)==(4 + sin(5))
@test ttf2(4,5.) == cos(4)-cos(5.)

@test ttf2(MTT1(5),4) == sin(5)-sin(4)
@test ttf2(MTT2(5),3) == sin(3)+5

@test length(methods(f2))==length(methods(ttf2))


text:
function Base.display(cell::Cell)
    print("origin = $(round.(cell.boundary.origin, digits=3)).")
    print("\twidth =  $(round.(cell.boundary.widths[1], digits=3)).")
    print("\tisleaf = $(isnothing(cell.children))\n")
    println("   ids = $(cell.data.ids)")
end

function Base.display(octree::Octree)
    print("origin = $(round.(octree.root.boundary.origin, digits=3)).")
    print("\twidth =  $(round.(octree.root.boundary.widths[1], digits=3)).")
    println("\tisleaf = $(isnothing(octree.root.children)). depth=$(octree.maxdepth)")
end

struct MyRefinery <: AbstractRefinery
    tolerance::Float64
    K::Int
    X::AbstractArray
end

function in(x::Array, r::HyperRectangle{N,T}) where {N,T<:AbstractFloat}
    for i in eachindex(x)
        x[i] < r.origin[i] && return false
        x[i] > r.origin[i] + r.widths[i] && return false
    end
    return true
end

function needs_refinement(r::MyRefinery, cell)
    maximum(cell.boundary.widths) <= r.tolerance && return false
    length(cell.data.ids) > r.K
end

function refine_data(r::MyRefinery, cell::Cell, indices)
    newpoints = Int32[]
    for i in cell.data.ids
        if r.X[:, i] in child_boundary(cell, indices)
            push!(newpoints, i)
        end
    end
    return LeafData(ids = newpoints,
                    code = cell.data.code * (indices .- 1),
                    depth = cell.data.depth + 1)
end

"""
    Octree(X::Array{S}, δ::T; K::Int=4) where {S<:M, T<:M} where M <: AbstractFloat

Create a quadtree, octree, or higher D-dimensional analog for the given dataset X with
dimensions D x N. The minimum number of points per cell is given by K and the minimum
side length of each cell is given by δ.
"""
function Octree(X::Array{S}, δ::T, K::Int) where {S<:M,T<:M} where {M<:AbstractFloat}
    r = MyRefinery(δ, K, X)
    # Make bounding box
    origin = SVector(minimum(X, dims = 2) .- δ/2...)
    widths = SVector(fill(
        maximum(maximum(X, dims = 2) - minimum(X, dims = 2)) .+ δ/2,
        size(X, 1),
    )...)
    root = Cell(origin, widths, LeafData(collect(1:size(X, 2)), LocationalCode(), Int16(0)))
    print("Making tree..")
    adaptivesampling!(root, r)
    println(" finished.")

    # Find max depth
    print("Finding max depth..")
    # TODO the following hangs when there is a Makie window open.
    maxdepth = maximum(map(x -> x.data.depth, allleaves(root)))
    println(" finished.")

    return Octree(X, root, maxdepth)
end

function Octree(X::Array{S}) where S<:AbstractFloat
    """
    Use a heuristic to find appropriate values of δ and K.
    """
    δ = Inf
    for i in rand(1:size(X, 2), floor(Int, size(X, 2) / 2))
        for j in rand(1:size(X, 2), floor(Int, size(X, 2) / 2))
            i == j && continue
            δ = min(δ, norm(X[:, i] - X[:, j]))
        end
    end
    K = ceil(Int, log10(size(X, 2)))
    println("Creating Octree with minimum side length δ = $(round(δ, digits=5)) and K=$K.")
    return Octree(X, δ, K)
end

const NEIGHBOR_INDICES = (
    (-1, -1, -1),
    (-1, -1, 0),
    (-1, -1, 1),
    (-1, 0, -1),
    (-1, 0, 0),
    (-1, 0, 1),
    (-1, 1, -1),
    (-1, 1, 0),
    (-1, 1, 1),
    (0, -1, -1),
    (0, -1, 0),
    (0, -1, 1),
    (0, 0, -1),
    (0, 0, 0),
    (0, 0, 1),
    (0, 1, -1),
    (0, 1, 0),
    (0, 1, 1),
    (1, -1, -1),
    (1, -1, 0),
    (1, -1, 1),
    (1, 0, -1),
    (1, 0, 0),
    (1, 0, 1),
    (1, 1, -1),
    (1, 1, 0),
    (1, 1, 1),
)

"""
    knn(octree::Octree, p, r)

    Find neighbors with fixed radius r of a given point p using an octree.
    """
function knn(octree::Octree, p::Vector, r::T) where T<:AbstractFloat
    ℓ = clamp(floor(Int, log(2, octree.root.boundary.widths[1] / r)), 1, octree.maxdepth)

    code = LocationalCode(octree::Octree, p)
    neighbors = Set{Int}()
    for (δ1, δ2, δ3) in NEIGHBOR_INDICES
        xloc = clamp(parse(Int, code.x[1:ℓ], base = 2) + δ1, 0, 2^(ℓ) - 1)
        yloc = clamp(parse(Int, code.y[1:ℓ], base = 2) + δ2, 0, 2^(ℓ) - 1)
        zloc = clamp(parse(Int, code.z[1:ℓ], base = 2) + δ3, 0, 2^(ℓ) - 1)
        xloc = bitstring(xloc)[end-ℓ+1:end]
        yloc = bitstring(yloc)[end-ℓ+1:end]
        zloc = bitstring(zloc)[end-ℓ+1:end]
        for i in code_to_cell(octree, LocationalCode(xloc, yloc, zloc)).data.ids
            if norm(octree.X[:, i] - p) < r
                push!(neighbors, i)
            end
        end
    end
    return neighbors
end


text:
"""
	Bipartition

Bit vector representing how a branch separates two exclusive sets of species in a tree. Unrooted tree bipartitions uniquely characterise branches invariantly to the position of the root or tree origin.
"""
struct Bipartition
	v::BitVector

	function Bipartition(v::BitVector)
		v = v[1] ? .!v : v

		return new(v)
	end
end


text:
export crammer_svm

using Flux
using LinearAlgebra
import SCS

"""
Implement multi-class kernel-based SVM

# Inputs

- `embed_support`: (`feat_dim`, `n_ways`*`k_shots`, `num_tasks`)
- `batch :: MetaDataSample`

# References

- http://jmlr.csail.mit.edu/papers/volume2/crammer01a/crammer01a.pdf
"""
function crammer_svm(embed_support, support_onehot, batch::MetaDataSample; C_reg = Float32(0.1))
    # reshape (m, n, c, num_samples, num_tasks) to (m, n, c, num_samples * num_tasks)
    n_support = batch.support_n_ways * batch.support_k_shots
    # construct dual QP problem: finding Q, p, G, h, A, b
    # (`n_ways`*`k_shots`, `n_ways`*`k_shots`, `num_tasks`)
    kernel_mat = Utils.gram_matrix(embed_support, embed_support)
    Q = Utils.batched_kronecker(
        kernel_mat,
        repeat(
            Array{Float32}(I, batch.support_n_ways, batch.support_n_ways),
            outer = (1, 1, size(batch)),
        ),
    )
    p = Float32.(-1 * support_onehot)
    G = repeat(
        Array{Float32}(
            I,
            n_support * batch.support_n_ways,
            n_support * batch.support_n_ways,
        ),
        outer = (1, 1, size(batch)),
    )
    h = C_reg * support_onehot

    A = Utils.batched_kronecker(
        repeat(Array{Float32}(I, n_support, n_support), outer = (1, 1, size(batch))),
        ones(Float32, 1, batch.support_n_ways, size(batch)),
    )
    b = zeros(Float32, n_support, size(batch))
    return Q, p, G, h, A, b
end


text:
include("../../../datastructs/src/datastructs_module.jl")

# module parseutil

# export is_left_bracket,
#     is_right_bracket,
#     is_left_curly_brace,
#     is_right_curly_brace,
#     is_comma,
#     is_colon,
#     json_syntax_chars,
#     json_quote,
#     is_json_whitespace,
#     is_json_syntax

# function is_json_syntax(char_str)
#     bool = char_str .== json_syntax_chars()
#     return any(bool)
# end

# end # parseutil

module lexer
# using ..parseutil
using ..datastructs

export  lex_json,
        TokenType,
        TokenElement,
        TokenList,
        is_left_bracket,
        is_right_bracket,
        is_left_curly_brace,
        is_right_curly_brace,
        is_comma,
        is_colon

function is_left_bracket(t::T) where {T <: Unsigned}
    return t == UInt8('[')
end
function is_right_bracket(t::T) where {T <: Unsigned}
    return t == UInt8(']')
end

function is_left_curly_brace(t::T) where {T <: Unsigned}
    return t == UInt8('{')
end
function is_right_curly_brace(t::T) where {T <: Unsigned}
    return t == UInt8('}')
end

function is_comma(t::T) where {T <: Unsigned}
    return t == UInt8(',')
end
function is_colon(t::T) where {T <: Unsigned}
    return t == UInt8(':')
end

# function json_syntax_chars()
#     return ['{', '}', '[', ']', ':', ',']
# end

function json_quote()
    return UInt8('"')
end

function is_json_whitespace(char_str::T) where {T <: Unsigned}
    # space " " 32
    # newline "\n" 10
    # carriage return "\r" 13
    # horztab "\t" 9
    return char_str == 32 || char_str == 10 || char_str == 9 || char_str == 13
end

@enum TokenType begin
    left_curly_t
    right_curly_t
    left_brace_t
    right_brace_t
    colon_t
    comma_t

    null_t
    bool_false_t
    bool_true_t

    number_t
    string_t

    unknown_t
end

struct TokenElement
    type::TokenType
    array_index::Int
end
function TokenElement()
    return TokenElement(unknown_t::TokenType, -1)
end

struct TokenList
    token_stack::datastructs.Stack{TokenElement}
    number_stack::datastructs.Stack{Float64}
    string_stack::datastructs.Stack{String}
end
function TokenList(init_size::Int = 0)
    token_stack = datastructs.Stack{TokenElement}(init_size)
    number_stack = datastructs.Stack{Float64}()
    string_stack = datastructs.Stack{String}()

    return TokenList(token_stack, number_stack, string_stack)
end

function Base.getindex(token_list::TokenList, idx::Int)
    return token_list.token_stack[idx]
end

function Base.getindex(token_list::TokenList, I::Vararg{Int, N}) where {N}
    return token_list.token_stack[I]
end

function Base.show(io::IO, token_list::TokenList)
    println("TokenList object with:")
    println("\ttoken_stack with $(token_list.token_stack.ptr) elements")
    println("\tnumber_stack with $(token_list.number_stack.ptr) elements")
    println("\tstring_stack with $(token_list.string_stack.ptr) elements")

    return io
end

function is_left_bracket(t::TokenElement)::Bool
    return t.type == left_brace_t::TokenType
end

function is_right_bracket(t::TokenElement)::Bool
    return t.type == right_brace_t::TokenType
end

function is_left_curly_brace(t::TokenElement)::Bool
    return t.type == left_curly_t::TokenType
end

function is_right_curly_brace(t::TokenElement)::Bool
    return t.type == right_curly_t::TokenType
end

function is_comma(t::TokenElement)::Bool
    return t.type == comma_t::TokenType
end

function is_colon(t::TokenElement)::Bool
    return t.type == colon_t::TokenType
end

function lex_string(json_text, start_idx::Int)

    str_ptr_1::Int = 0
    str_ptr_2::Int = 0

    json_string = nothing
    if json_text[start_idx] != json_quote()
        return (start_idx, str_ptr_1, str_ptr_2)
    end
    L = length(json_text)
    firstidx = start_idx + 1
    @inbounds for ii = firstidx:L
        if json_text[ii] == json_quote()
            json_string = json_text.s[firstidx:ii-1]
            # json_text = json_text[(ii+1):end]
            start_idx = ii+1

            str_ptr_1 = firstidx
            str_ptr_2 = ii-1
            break
        end
    end
    return (start_idx, str_ptr_1, str_ptr_2)
end

function lex_number(json_text, start_idx::Int, valid_chars; number_type = Float64)
    json_num = 0.0
    # valid_chars = ('-', 'e', 'E', '.', '0', '1', '2', '3', '4', '5', '6', '7', '8', '9' )

    success = false
    L = length(json_text)
    idx = 0
    @inbounds for cc in start_idx:L
        if any(json_text[cc] .== valid_chars) # cc in valid_chars
            idx += 1
        else
            break
        end
    end
    if idx > 0
        success = true
        # json_string = json_text[start_idx:idx+start_idx-1]
        json_string = SubString(json_text.s, start_idx:idx+start_idx-1)
        json_num = parse(number_type, json_string)
        start_idx += idx
        # json_text = json_text[(idx+1):end]
    end
    return (start_idx, json_num, success)

end

function lex_bool(json_text, start_idx::Int, str_true, str_false)
    success = false

    Lt = 4 # length("true")
    Lf = 5 # length("false")

    L = length(json_text)
    Ltidx = min( Lt+start_idx-1, L )
    Lfidx = min( Lf+start_idx-1, L )

    # substr_t = SubString(json_text.s, start_idx:Ltidx)
    # substr_f = SubString(json_text.s, start_idx:Lfidx)

    substr_t = json_text[ start_idx:Ltidx ]
    substr_f = json_text[ start_idx:Lfidx ]

    # if lowercase(json_text[start_idx:Ltidx]) == "true"
    json_bool = false
    if substr_t == str_true # "true"
        json_bool = true
        #json_text = json_text[(Lt+1):end]
        start_idx += Lt
        success = true
    elseif substr_f == str_false# "false"
    # elseif lowercase(json_text[start_idx:Lfidx]) == "false"
        json_bool = false
        #json_text = json_text[(Lf+1):end]
        start_idx += Lf
        success = true
    end

    return (start_idx, json_bool, success)
end

function lex_null(json_text, start_idx::Int, null_str)
    success = false
    Ln = 4 # length("null")
    L = length(json_text)
    idx = start_idx:min(start_idx+Ln,L)

    substr = SubString(json_text.s, idx)
    # if lowercase(json_text[idx]) == "null"
    if substr == null_str # "null"
        success = true
        start_idx += Ln
        # json_text = json_text[(Ln+1):end]
    end
    return (start_idx, success)
end

function lex_json(json_text::T) where {T <: AbstractString}

    # tokens = Vector{Any}(undef, 0)
    null_str = codeunits("null")
    str_true = codeunits("true")
    str_false = codeunits("false")
    valid_chars = ('-', 'e', 'E', '.', '0', '1', '2', '3', '4', '5', '6', '7', '8', '9' ) .|> UInt8

    Ls = length(json_text)
    if Ls < 1
        return tokens
    end

    start_idx::Int = 1
    # tokens = datastructs.Stack{Any}(div(Ls,10)) # heuristics

    token_list = TokenList()
    utf8 = codeunits(json_text)

    @inbounds while start_idx <= Ls

        # println(start_idx)
        # println(tokens)

        (start_idx, str_ptr_1, str_ptr_2) = lex_string(utf8, start_idx)
        if str_ptr_1 != 0
            json_string = utf8.s[str_ptr_1:str_ptr_2]
            push!(token_list.string_stack, json_string)
            idx = token_list.string_stack.ptr
            push!(token_list.token_stack, TokenElement(string_t::TokenType, idx))
            continue
        end

        (start_idx, json_num, success) = lex_number(utf8, start_idx, valid_chars)
        if success
            push!(token_list.number_stack, json_num)
            idx = token_list.number_stack.ptr
            push!(token_list.token_stack, TokenElement(number_t::TokenType, idx))
            continue
        end

        (start_idx, json_bool, success) = lex_bool(utf8, start_idx, str_true, str_false)
        if success
            if json_bool
                push!(token_list.token_stack, TokenElement(bool_true_t::TokenType, 0))
            else
                push!(token_list.token_stack, TokenElement(bool_false_t::TokenType, 0))
            end
            continue
        end

        (start_idx, success) = lex_null(utf8, start_idx, null_str)
        if success
            # push!(tokens, json_string)
            push!(token_list.token_stack, TokenElement(null_t::TokenType, 0))
            continue
        end

        if is_json_whitespace(utf8[start_idx])
            # json_text = json_text[2:end]
            start_idx += 1
        # elseif  is_json_syntax(json_text[start_idx])

        elseif is_left_curly_brace(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(left_curly_t::TokenType, 0))
            start_idx += 1
        elseif is_right_curly_brace(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(right_curly_t::TokenType, 0))
            start_idx += 1
        elseif is_left_bracket(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(left_brace_t::TokenType, 0))
            start_idx += 1
        elseif is_right_bracket(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(right_brace_t::TokenType, 0))
            start_idx += 1
        elseif is_colon(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(colon_t::TokenType, 0))
            start_idx += 1
        elseif is_comma(utf8[start_idx])
            push!(token_list.token_stack, TokenElement(comma_t::TokenType, 0))
            # json_text = json_text[2:end]
            start_idx += 1
        else
            fi = max(0, start_idx-10)
            ei = min(Ls, start_idx + 10)
            println("$(utf8.s[fi:ei])")
            error("Unexpected character: $(utf8.s[start_idx])")
        end

    end

    return token_list
end

end # lexer

module parser
using ..datastructs
# using ..parseutil
using ..lexer

function parse_json(text)
    Lt = length(text)
    if Lt < 0
        return nothing
    end

    # transform everything into ascii
    # text = ascii(text)
    token_list = lex_json(text) # @time
    L = length(token_list.token_stack)
    if L < 1
        return nothing
    end

    start_idx = 1
    (start_idx, json_obj, tokens) = composite_object_parse(token_list, start_idx) # @time
    return json_obj
end

# key = strip_key(tokens, key)
function strip_key(token_list::TokenList, token_element::TokenElement)

    if token_element.type != lexer.string_t
        error("The key in an object must be a string")
    end
    idx = token_element.array_index
    str = token_list.string_stack[idx]
    if str[1] == '"'
        str = chop(str; head = 1, tail = 0)
    end
    if str[end] == '"'
        str = chop(str; head = 0, tail = 1)
    end
    return String(str)
end

function is_type_of(token_element::TokenElement, typedecl::TokenType)
    return token_element.type == typedecl
end

function is_string_type(token_element)
    return is_type_of(token_element, lexer.string_t)
end

function is_number_type(token_element)
    return is_type_of(token_element, lexer.number_t)
end

function composite_object_parse(tokens::TokenList, start_idx::Int)

    token_element = tokens[start_idx]

    if is_left_bracket(token_element) # array begin
        start_idx += 1
        return parse_array(tokens, start_idx)

    elseif is_left_curly_brace(token_element) # object begin
        start_idx += 1
        return parse_object(tokens, start_idx)

    else

        if is_string_type(token_element)

            start_idx += 1
            val = tokens.string_stack[ token_element.array_index ]
            return (start_idx, val, tokens)

        elseif is_number_type(token_element)

            start_idx += 1
            val = tokens.number_stack[ token_element.array_index ]
            return (start_idx, val, tokens)

        elseif is_type_of(token_element, lexer.null_t)
            start_idx += 1
            return (start_idx, NaN, tokens)

        elseif is_type_of(token_element, lexer.bool_true_t)
            start_idx += 1
            return (start_idx, true, tokens)

        elseif is_type_of(token_element, lexer.bool_false_t)
            start_idx += 1
            return (start_idx, false, tokens)

        end
        # default
        start_idx += 1
        return (start_idx, token_element, tokens)
    end
end

function parse_array(tokens, start_idx)

    array = Vector{Any}(undef, 0)

    t = tokens[1]
    if is_right_bracket(t)
        # return (array, tokens[2:end])
        start_idx += 1
        return (start_idx, array, tokens)
    end

    # Lt = length(tokens)
    while true # avoid infinite loop

        (start_idx, element_value, tokens) = composite_object_parse(tokens, start_idx)
        push!(array, element_value)

        t = tokens[start_idx]
        if is_right_bracket(t)
            start_idx += 1
            return (start_idx, array, tokens)
        elseif !is_comma(t)
            error("Expected a comma after object in the array")
        else
            start_idx += 1
            # tokens = tokens[2:end]
        end

    end

    return (start_idx, array, tokens)
end

function parse_object(tokens, start_idx)

    json_obj = Dict{String, Any}()

    t = tokens[start_idx]
    if is_right_curly_brace(t)
        return (start_idx, json_obj, tokens)
    end

    while true

        key = tokens[start_idx]
        key = strip_key(tokens, key)
        start_idx += 1 # start index is increased if it is a valid key == doesnt throw an error
        # if isa(key, String)
        #     start_idx += 1
        #     # tokens = tokens[2:end]
        # else
        #     error("Invalid key: $(key)")
        # end

        if !is_colon(tokens[start_idx])
            error("Expected a colon, but got: $(token[start_idx])")
        end

        start_idx += 1
        (start_idx, value, tokens) = composite_object_parse(tokens, start_idx)
        json_obj[key] = value

        t = tokens[start_idx]
        if is_right_curly_brace(t)
            start_idx += 1
            return (start_idx, json_obj, tokens)
        elseif !is_comma(t)
            error("Expected comma after object but got: $(t)")
        end

        start_idx += 1
        # tokens = tokens[2:end]
    end
    return (start_idx, json_obj, tokens)
end

end # parser

module json
# using ..parseutil
using ..lexer
using ..parser

function parse(json_text)
    return parser.parse_json(json_text)
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using MPCCmod

MPCCmod_success=true
println("On teste le module MPCCmod")
f(x)=x[1]-x[2]

if false
#On transforme la sortie NL en type MPCC :
#Constructeur 1:
mpcc1 = MPCCmod.MPCC(f,ones(2),G,H,1,-Inf*ones(2),Inf*ones(2),c,ones(1),lcon,Inf*ones(1))

#Constructeur 2:
nlp=ADNLPModel(f,ones(2), lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1,1e-6)

println("")
println("KDB test")
println("")
xkdb,fkdb,statkdb = MPCCsolve.solve(mpcc,0.1,0.01,0.1,0.01,0.0,0.01,"KDB")

#Test le solve
xss,fss,statss = MPCCsolve.solve(mpcc)
xb,fb,statb = MPCCsolve.solve(mpcc,0.1,0.01,0.1,0.01,0.1,0.01,"Butterfly")
println("")
println("KDB test")
println("")
xkdb,fkdb,statkdb = MPCCsolve.solve(mpcc,0.1,0.01,0.1,0.01,0.0,0.01,"KDB")
xks,fks,statks = MPCCsolve.solve(mpcc,0.1,0.01,0.1,0.01,0.1,0.01,"KS")

if !(norm(-1-fss)<1e-6 || norm(-1-fb)<1e-6 || norm(-1-fkdb)<1e-6 || norm(-1-fks)<1e-6)
 println("Distance à l'objectif")
 println(norm(-1-fss)<1e-6, norm(-1-fb)<1e-6 , norm(-1-fkdb)<1e-6, norm(-1-fks)<1e-6)
 MPCCmod_success=false
end
if !(MPCCmod.viol_comp(mpcc,xss)<1e-6 || MPCCmod.viol_comp(mpcc,xb)<1e-6 || MPCCmod.viol_comp(mpcc,xkdb)<1e-6 || MPCCmod.viol_comp(mpcc,xks)<1e-6)
 println("Satisfaction de la complémentarité")
 println(MPCCmod.viol_comp(mpcc,xss)<1e-6,MPCCmod.viol_comp(mpcc,xb)<1e-6,MPCCmod.viol_comp(mpcc,xkdb)<1e-6,MPCCmod.viol_comp(mpcc,xks)<1e-6)
 MPCCmod_success=false
end
end

#Bilan : donner une sortie bilan
if MPCCmod_success==true
 println("MPCCmod.jl passes the test !")
else
 println("MPCCmod.jl contains some error")
end


text:
# This file contains a number of additions to BlockArrays.jl. These are completely
# independent of Stheno.jl, and will (hopefully) move over to BlockArrays.jl at some point.

using FillArrays, LinearAlgebra, BlockArrays
using FillArrays: Fill

import Base: +, *, size, getindex, eltype, copy, \, vec, getproperty
import LinearAlgebra: UpperTriangular, LowerTriangular, logdet, Symmetric, transpose,
    adjoint, AdjOrTrans, AdjOrTransAbsMat, cholesky
import BlockArrays: BlockArray, BlockVector, BlockMatrix, BlockVecOrMat, getblock,
    blocksize, setblock!, nblocks

# Do some character saving.
const BV{T} = BlockVector{T}
const BM{T} = BlockMatrix{T}
const ABV{T} = AbstractBlockVector{T}
const ABM{T} = AbstractBlockMatrix{T}
const ABVM{T} = AbstractBlockVecOrMat{T}
const LUABM{T} = Union{ABM, LowerTriangular{T, <:ABM}, UpperTriangular{T}, <:ABM}

const AdjOrTransABM{T} = Union{ABM{T}, AdjOrTrans{T, <:ABM{T}}}
const AdjTransTriABM{T} = Union{AbstractTriangular{T, <:ABM{T}}, AdjOrTransABM{T}}
const TriABM{T} = Union{AbstractTriangular{T, <:ABM{T}}, ABM{T}}
const BlockTri{T} = AbstractTriangular{T, <:ABM{T}}
const BlockAdjOrTrans{T} = AdjOrTrans{T, <:ABM{T}}

unbox(X::AbstractBlockMatrix) = X
unbox(X::Symmetric) = unbox(X.data)
unbox(X::AbstractMatrix) = X

####################################### Various util #######################################

import BlockArrays: BlockVector, BlockMatrix
export BlockVector, BlockMatrix, blocksizes, blocklengths

"""
    BlockVector(xs::Vector{<:AbstractVector{T}}) where T

Construct a `BlockVector` from a collection of `AbstractVector`s.
"""
function BlockArrays.BlockVector(xs::Vector{<:AbstractVector{T}}) where T
    x = BlockVector{T}(undef_blocks, length.(xs))
    for (n, x_) in enumerate(xs)
        setblock!(x, x_, n)
    end
    return x
end

"""
    BlockMatrix(Xs::Matrix{<:AbstractVecOrMat{T}}) where T

Construct a `BlockMatrix` from a matrix of `AbstractVecOrMat`s.
"""
function BlockMatrix(Xs::Matrix{<:AbstractVecOrMat{T}}) where T
    X = BlockMatrix{T}(undef_blocks, size.(Xs[:, 1], Ref(1)), size.(Xs[1, :], Ref(2)))
    for q in 1:nblocks(X, 2), p in 1:nblocks(X, 1)
        setblock!(X, Xs[p, q], p, q)
    end
    return X
end
BlockMatrix(Xs::Vector{<:AbstractVecOrMat}) = BlockMatrix(reshape(Xs, length(Xs), 1))
BlockMatrix(x::AbstractVecOrMat) = BlockMatrix([x])

"""
    BlockMatrix(xs::Vector{<:AM}, P::Int, Q::Int)

Construct a block matrix with `P` rows and `Q` columns of blocks.
"""
BlockMatrix(xs::Vector{<:AM}, P::Int, Q::Int) = BlockMatrix(reshape(xs, P, Q))

"""
    blocksizes(X::AbstractBlockArray, d::Int)

Get a vector containing the block sizes over the `d`th dimension of `X`.
"""
function blocksizes(X::AbstractBlockArray, d::Int)
    @assert d > 0 && d <= ndims(X)
    idxs = [1 for _ in 1:length(size(X))]
    block_sizes = Vector{Int}(undef, nblocks(X, d))
    for p in eachindex(block_sizes)
        idxs[d] = p
        block_sizes[p] = blocksize(X, idxs...)[d]
    end
    return block_sizes
end
blocksizes(X::AbstractBlockArray) = ([blocksizes(X, d) for d in 1:length(size(X))]...,)
blocklengths(x::BlockVector) = blocksizes(x, 1)

################################# Symmetric BlockMatrices ##############################

const BS{T} = Symmetric{T, <:AbstractBlockMatrix{T}}
unbox(X::BS) = X.data
nblocks(X::BS) = nblocks(unbox(X))
nblocks(X::BS, i::Int) = nblocks(unbox(X), i)
blocksize(X::BS, N::Int...) = blocksize(unbox(X), N...)
blocksizes(X::BS, d::Int...) = blocksizes(unbox(X), d...)

function getblock(X::BS, p::Int, q::Int)
    @assert blocksizes(X, 1) == blocksizes(X, 2)
    X_, uplo = unbox(X), X.uplo
    if p < q
        return uplo == 'U' ? getblock(X_, p, q) : transpose(getblock(X_, q, p))
    elseif p == q
        return Symmetric(getblock(X_, p, q))
    else
        return uplo == 'U' ? transpose(getblock(X_, q, p)) : getblock(X_, p, q)
    end
end

######################## Util for triangular block matrices ######################

@inline unbox(U::UpperTriangular{T, <:ABM{T}} where T) = U.data
function blocksize(U::UpperTriangular{T, <:ABM{T}} where T, p::Int, q::Int)
    return blocksize(unbox(U), p, q)
end
blocksizes(U::UpperTriangular{T, <:ABM{T}} where T, d::Int) = blocksizes(unbox(U), d)
blocksizes(U::UpperTriangular{T, <:ABM{T}} where T) = blocksizes(unbox(U))
function getblock(U::UpperTriangular{T, <:ABM{T}}, p::Int, q::Int) where T
    @assert blocksizes(U, 1) == blocksizes(U, 2)
    if p > q
        return Zeros{T}(blocksize(U, p, q)...)
    elseif p == q
        return UpperTriangular(getblock(unbox(U), p, q))
    else
        return getblock(unbox(U), p, q)
    end
end
nblocks(U::UpperTriangular{T, <:ABM{T}} where T, d::Int...) = nblocks(unbox(U), d...)
function BlockMatrix(U::UpperTriangular{T, <:ABM{T}}) where T
    B = BlockMatrix{T}(undef_blocks, blocksizes(U)...)
    for q in 1:nblocks(U, 2)
        for p in 1:q-1
            setblock!(B, getblock(U, p, q), p, q)
        end
        setblock!(B, UpperTriangular(getblock(U, q, q)), q, q)
        for p in q+1:nblocks(U, 1)
            setblock!(B, Zeros{T}(blocksize(U, p, q)), p, q)
        end
    end
    return B
end

@inline unbox(L::LowerTriangular{T, <:ABM{T}} where T) = L.data
function blocksize(L::LowerTriangular{T, <:ABM{T}} where T, p::Int, q::Int)
    return blocksize(unbox(L), p, q)
end
blocksizes(L::LowerTriangular{T, <:ABM{T}} where T, d::Int) = blocksizes(unbox(L), d)
blocksizes(L::LowerTriangular{T, <:ABM{T}} where T) = blocksizes(unbox(L))
function getblock(L::LowerTriangular{T, <:ABM{T}}, p::Int, q::Int) where T
    @assert blocksizes(L, 1) == blocksizes(L, 2)
    if p > q
        return getblock(unbox(L), p, q)
    elseif p == q
        return LowerTriangular(getblock(unbox(L), p, q))
    else
        return Zeros{T}(blocksize(L, p, q)...)
    end
end
nblocks(L::LowerTriangular{T, <:ABM{T}} where T, d::Int...) = nblocks(unbox(L), d...)
function BlockMatrix(L::LowerTriangular{T, <:ABM{T}}) where T
    B = BlockMatrix{T}(undef_blocks, blocksizes(L)...)
    for q in 1:nblocks(L, 2)
        for p in 1:q-1
            setblock!(B, Zeros{T}(blocksize(L, p, q)), p, q)
        end
        setblock!(B, LowerTriangular(getblock(L, q, q)), q, q)
        for p in q+1:nblocks(L, 1)
            setblock!(B, getblock(L, p, q), p, q)
        end
    end
    return B
end

####################################### Copying ######################################

function copy(a::BlockVector{T}) where T
    b = BlockVector{T}(undef_blocks, blocksizes(a, 1))
    for p in 1:nblocks(b, 1)
        setblock!(b, copy(getblock(a, p)), p)
    end
    return b
end

function copy(A::BlockMatrix{T, V}) where {T, V<:AbstractMatrix{T}}
    B = BlockMatrix{T, V}(undef_blocks, blocksizes(A, 1), blocksizes(A, 2))
    for q in 1:nblocks(B, 2), p in 1:nblocks(B, 1)
        setblock!(B, copy(getblock(A, p, q)), p, q)
    end
    return B
end

copy(B::BS{T, <:ABM{T}} where T) = Symmetric(copy(unbox(B)))
copy(L::LowerTriangular{T, <:BS{T}} where T) = LowerTriangular(copy(unbox(L)))
copy(U::UpperTriangular{T, <:BS{T}} where T) = UpperTriangular(copy(unbox(U)))

####################################### Transposition ######################################

nblocks(A::AdjOrTrans{T, <:TriABM{T}} where T) = reverse(nblocks(A'))
function nblocks(A::AdjOrTrans{T, <:TriABM{T}} where T, d::Int) where {N}
    @assert d ∈ (1, 2)
    return nblocks(A)[d]
end
function blocksizes(A::AdjOrTrans{T, <:TriABM{T}} where T, d::Int)
    @assert d ∈ (1, 2)
    return blocksizes(A', d == 1 ? 2 : 1)
end
function getblock(A::Adjoint{T, <:TriABM{T}} where T<:Real, p::Int, q::Int)
    return getblock(A.parent, q, p)'
end
function getblock(A::Transpose{T, <:TriABM{T}} where T<:Real, p::Int, q::Int)
    return transpose(getblock(A.parent, q, p))
end

####################################### Multiplication #####################################

"""
    are_conformal(A::BlockVecOrMat, B::BlockVecOrMat)

Test whether two block matrices (or vectors) are conformal. This criterion is stricter than
that for general matrices / vectors as we additionally require that each block be conformal
with block of the other matrix with which it will be multiplied. This ensures that the
result is itself straightforwardly representable as `BlockVecOrMat`.
"""
are_conformal(A::AVM, B::AVM) = blocksizes(A, 2) == blocksizes(B, 1)

"""
    *(A::AdjTransTriABM, x::BlockVector)

Matrix-vector multiplication between `BlockArray`s. Fails if block are not conformal.
"""
function *(
    A::Union{ABM{T}, BlockTri{T}, BlockAdjOrTrans{T}, AdjOrTrans{T, <:BlockTri{T}}},
    x::ABV{T},
) where T
    @assert are_conformal(A, x)
    y = BlockVector{T}(undef_blocks, blocksizes(A, 1))
    P, Q = nblocks(A)
    for p in 1:P
        setblock!(y, getblock(A, p, 1) * getblock(x, 1), p)
        for q in 2:Q
            setblock!(y, getblock(y, p) + getblock(A, p, q) * getblock(x, q), p)
        end
    end
    return y
end
function *(A::AdjTransTriABM{T}, b::AV{T}) where {T}
    @assert nblocks(A, 2) == 1
    return A * BlockVector([b])
end
function *(A::AdjTransTriABM{T}, b::FillArrays.Zeros{T,1}) where T
    @show size(A), size(b)
    return invoke(*, Tuple{AbstractMatrix, typeof(b)}, A, b)
end

"""
    *(A::AdjTransTriABM{T}, B::AdjTransTriABM{T}) where T

Matrix-matrix multiplication between `BlockArray`s. Fails if blocks are not conformal.
"""
function *(
    A::Union{ABM{T}, BlockTri{T}, BlockAdjOrTrans{T}, AdjOrTrans{T, <:BlockTri{T}}},
    B::Union{ABM{T}, BlockTri{T}, BlockAdjOrTrans{T}, AdjOrTrans{T, <:BlockTri{T}}},
) where T
    @assert are_conformal(A, B)
    C = BlockMatrix{T}(undef_blocks, blocksizes(A, 1), blocksizes(B, 2))
    P, Q, R = nblocks(A, 1), nblocks(A, 2), nblocks(B, 2)
    for p in 1:P, r in 1:R
        setblock!(C, getblock(A, p, 1) * getblock(B, 1, r), p, r)
        for q in 2:Q
            setblock!(C, getblock(C, p, r) + getblock(A, p, q) * getblock(B, q, r), p, r)
        end
    end
    return C
end
*(A::AdjTransTriABM{T}, B::AM{T}) where T = A * BlockMatrix([B])
*(A::AM{T}, B::AdjTransTriABM{T}) where T = BlockMatrix([A]) * B

function *(A::LazyPDMat{T, <:Symmetric{T, <:ABM{T}}}, B::AdjTransTriABM{T}) where T
    return unbox(unbox(A)) * B
end

const UpperOrAdjLower{T} = Union{
    UpperTriangular{T, <:ABM{T}},
    AdjOrTrans{T, <:LowerTriangular{T, <:ABM{T}}},
}
const LowerOrAdjUpper{T} = Union{
    LowerTriangular{T, <:ABM{T}},
    AdjOrTrans{T, <:UpperTriangular{T, <:ABM{T}}},
}

function \(U::UpperOrAdjLower{T}, x::ABV{T}) where T<:Real
    @assert are_conformal(unbox(U), x)
    y = BlockVector{T}(undef_blocks, blocksizes(U, 1))
    for p in reverse(1:nblocks(y, 1))
        setblock!(y, getblock(x, p), p)
        for p′ in p+1:nblocks(y, 1)
            setblock!(y, getblock(y, p) - getblock(U, p, p′) * getblock(y, p′), p)
        end
        setblock!(y, getblock(U, p, p) \ getblock(y, p), p)
    end
    return y
end

function _block_ldiv_mat_upper(U, X)
    @assert are_conformal(unbox(U), X)
    Y = BlockMatrix{eltype(U)}(undef_blocks, blocksizes(U, 1), blocksizes(X, 2))
    for q in 1:nblocks(Y, 2), p in reverse(1:nblocks(Y, 1))
        setblock!(Y, getblock(X, p, q), p, q)
        for p′ in p+1:nblocks(Y, 1)
            setblock!(Y, getblock(Y, p, q) - getblock(U, p, p′) * getblock(Y, p′, q), p, q)
        end
        setblock!(Y, getblock(U, p, p) \ getblock(Y, p, q), p, q)
    end
    return Y
end
function \(U::UpperOrAdjLower{T}, X::Union{ABM{T}, BlockTri{T}}) where T<:Real
    return _block_ldiv_mat_upper(U, X)
end
function \(U::UpperOrAdjLower{T}, X::Adjoint{T, <:UpperTriangular{T, <:ABM{T}}}) where T<:Real
    return _block_ldiv_mat_upper(U, X)
end

function \(L::LowerOrAdjUpper{T}, x::ABV{T}) where T<:Real
    @assert are_conformal(unbox(L), x)
    y = BlockVector{T}(undef_blocks, blocksizes(L, 1))
    for p in 1:nblocks(y, 1)
        setblock!(y, getblock(x, p), p)
        for p′ in 1:p-1
            setblock!(y, getblock(y, p) - getblock(L, p, p′) * getblock(y, p′), p)
        end
        setblock!(y, getblock(L, p, p) \ getblock(y, p), p)
    end
    return y
end

function _block_ldiv_mat_lower(L, X)
    @assert are_conformal(unbox(L), X)
    Y = BlockMatrix{eltype(L)}(undef_blocks, blocksizes(L, 1), blocksizes(X, 2))
    for q in 1:nblocks(Y, 2), p in 1:nblocks(Y, 1)
        setblock!(Y, getblock(X, p, q), p, q)
        for p′ in 1:p-1
            setblock!(Y, getblock(Y, p, q) - getblock(L, p, p′) * getblock(Y, p′, q), p, q)
        end
        setblock!(Y, getblock(L, p, p) \ getblock(Y, p, q), p, q)
    end
    return Y
end
function \(L::LowerOrAdjUpper{T}, X::Union{ABM{T}, BlockTri{T}}) where T<:Real
    return _block_ldiv_mat_lower(L, X)
end
function \(L::LowerOrAdjUpper{T}, X::Adjoint{T, <:UpperTriangular{T, <:ABM{T}}}) where T<:Real
    return _block_ldiv_mat_lower(L, X)
end

"""
    _chol(A::Symmetric{T, <:BlockMatrix{T}}) where T<:Real

Get the Cholesky decomposition of `A` in the form of a `BlockMatrix`.

Only works for `A` where `is_block_symmetric(A) == true`. Assumes that we want the
upper triangular version.
"""
function _chol(A::Symmetric{T, <:BlockMatrix{T, V}}) where {T<:Real, V<:AbstractMatrix{T}}

    U = BlockMatrix{T, V}(undef_blocks, blocksizes(A, 1), blocksizes(A, 1))

    # Do an initial pass to fill each of the blocks with Zeros. This is cheap
    for q in 1:nblocks(U, 2), p in 1:nblocks(U, 1)
        setblock!(U, Zeros{T}(blocksize(A, p, q)...), p, q)
    end

    # Fill out the upper triangle with the Cholesky
    for j in 1:nblocks(A, 2)

        # Update off-diagonals.
        for i in 1:j-1
            setblock!(U, getblock(A, i, j), i, j)
            for k in 1:i-1
                Uki, Ukj = getblock(U, k, i), getblock(U, k, j)
                setblock!(U, getblock(U, i, j) - Uki' * Ukj, i, j)
            end
            setblock!(U, getblock(U, i, i)' \ getblock(U, i, j), i, j)
        end

        # Update diagonal.
        setblock!(U, getblock(A, j, j), j, j)
        for k in 1:j-1
            Ukk, Ukj = getblock(U, k, k), getblock(U, k, j)
            setblock!(U, getblock(U, j, j) - Ukj' * Ukj, j, j)
        end
        setblock!(U, cholesky(getblock(U, j, j)).U, j, j)
    end
    return U
end

# Get only implement one of these methods for now.
function cholesky(A::Symmetric{T, <:BlockMatrix{T, <:AbstractMatrix{T}}} where T<:Real)
    return Cholesky(_chol(A), :U, 0)
end

logdet(A::Cholesky{T, <:AbstractBlockMatrix{T}}) where T = 2 * logdet(A.U)
function logdet(U::UpperTriangular{T, <:AbstractBlockMatrix{T}}) where T
    val = zero(T)
    for n in 1:nblocks(U, 1)
        val += logdet(UpperTriangular(getblock(unbox(U), n, n)))
    end
    return val
end

function +(u::UniformScaling, X::AbstractBlockMatrix)
    @assert blocksizes(X, 1) == blocksizes(X, 2)
    Y = copy(X)
    for p in 1:nblocks(Y, 1)
        setblock!(Y, getblock(Y, p, p) + u, p, p)
    end
    return Y
end
+(u::UniformScaling, X::Symmetric{T, <:ABM{T}} where T) = Symmetric(u + unbox(X))
function +(X::AbstractBlockMatrix, u::UniformScaling)
    @assert blocksizes(X, 1) == blocksizes(X, 2)
    Y = copy(X)
    for p in 1:nblocks(Y, 1)
        setblock!(Y, u + getblock(Y, p, p), p, p)
    end
    return Y
end
+(X::Symmetric{T, <:ABM{T}} where T, u::UniformScaling) = Symmetric(unbox(X) + u)

# Define addition and subtraction for compatible block matrices and vectors.
import Base: +, -
for foo in [:+, :-]
    @eval function $foo(A::BV{T}, B::BV{T}) where T
        @assert blocksizes(A, 1) == blocksizes(B, 1)
        C = BlockVector{T}(undef_blocks, blocksizes(A, 1))
        for p in 1:nblocks(C, 1)
            setblock!(C, $foo(getblock(A, p), getblock(B, p)), p)
        end
        return C
    end
    @eval function $foo(A::BM{T}, B::BM{T}) where T
        @assert blocksizes(A, 1) == blocksizes(B, 1)
        @assert blocksizes(A, 2) == blocksizes(B, 2)
        C = BlockMatrix{T}(undef_blocks, blocksizes(A, 1), blocksizes(A, 2))
        for q in 1:nblocks(C, 2), p in 1:nblocks(C, 1)
            setblock!(C, $foo(getblock(A, p, q), getblock(B, p, q)), p, q)
        end
        return C
    end
end

#################################### Broadcasting ##########################################
# Override the usual broadcasting machinery for BlockArrays. This is a pretty disgusting
# hack. At the time of writing it, I was also writing my first year report, so was short of
# time. This is generally an open problem for AbstractBlockArrays which really does need to
# be resolved at some point.

# Very specific `broadcast` method for particular use case. Needs to be generalised.
function broadcast(f, A::AbstractBlockVector, b::Real)
    return BlockVector([broadcast(f, getblock(A, p), b) for p in 1:nblocks(A, 1)])
end
function broadcast(f, A::AbstractBlockMatrix, b::Real)
    return BlockMatrix([broadcast(f, getblock(A, p, q), b)
        for p in 1:nblocks(A, 1), q in 1:nblocks(A, 2)])
end


text:
@testset "94.binary-tree-inorder-traversal.jl" begin
    @test postorder_traversal(TreeNode{Int}([1])) == [1]
    @test postorder_traversal(nothing) == Int[]
    @test postorder_traversal(TreeNode{Int}([1, nothing, 2, 3])) == [3, 2, 1]
    @test postorder_traversal(TreeNode{Int}([1, 2, 3, 4, 5, 6])) == [4, 5, 2, 6, 3, 1]
    @test postorder_traversal(TreeNode{Int}([3, 9, 4, nothing, nothing, 5, 7])) ==
        [9, 5, 7, 4, 3]
end


text:
using StatsBase
using PyPlot

function linreg2(x, y, w::WeightVec)
    size(x) == size(y) || throw(DimensionMismatch("x and y must be the same size"))
    my = mean(y, w)
    mx, vx = mean(x, w)
    # don't need to worry about the scaling (n vs n - 1) since they cancel in the ratio
    b = cov(x, y)/vx
    a = my - b*mx
    return (a, b)
end

x = collect(linspace(0, 10, 20))
y = randn(20)
w = WeightVec(ones(20))
cov(x, y, w)
linreg(x, y)
linreg(x, y, ones(20))

linreg2(x, y, w)

max_stars_repo_name: gabrielgellner/LeastSquares

max_issues_repo_name: gabrielgellner/LeastSquares

max_forks_repo_name: gabrielgellner/LeastSquares


text:
export Payment, CashFlow, npv, irr

struct Payment{T, S<:Real}
    date::T
    amount::S
end

struct CashFlow{T, S<:Real, R<:DayCountConvention}
    flow::Array{Payment{T,S},1}
    daycount::R
    basedate::Date
end

Base.length(cf::CashFlow) = length(cf.flows)

Base.convert(::Type{Payment{T}}, p::Payment{Date}, daycount, base = today()) where {T<:Real} =
    Payment(yearfraction(daycount, base, p.date, T), p.amount)

Base.convert(::Type{CashFlow{T}}, cf::CashFlow{Date}) where {T<:Real} =
    CashFlow([convert(Payment{T}, p, cf.daycount, cf.basedate) for p in cf.flow], cf.daycount, cf.basedate)

Base.show(io::IO, cf::CashFlow) = begin
    println(io, rpad("Day count convention:", 23, ' ') * "$(cf.daycount)")
    println(io, rpad("Base date:", 23, ' ') * "$(cf.basedate)")
    println(io, "Cash flow:")
    for p in cf.flow
        println(io, rpad("  $(p.date)", 23, ' ') * "$(p.amount)")
    end
end

npv(flow, rate) = sum(p.amount / (1 + rate)^p.date for p in flow)
npv(cf::CashFlow{T}, rate::T) where {T<:Real} = npv(cf.flow, rate)
npv(cf::CashFlow, rate::T) where {T<:Real} = npv(convert(CashFlow{T}, cf), rate)

irr(flow, r0) = fzero(r -> npv(flow, r), r0)
irr(cf::CashFlow{T}, r0::T = 0.1) where {T<:Real} = irr(cf.flow, r0)
irr(cf::CashFlow, r0::T = 0.1) where {T<:Real} = irr(convert(CashFlow{T}, cf), r0)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EigenAnalysis(proj)

The eigenanalysis of the covariance with a given projection `proj`.

## Projections

* `:V` - Uncorrelated variables (PCA transform)
* `:VD` - Uncorrelated variables and variance one (DRS transform)
* `:VDV` - Uncorrelated variables and variance one (SDS transformation)

The `:V` projection used in the PCA transform projects the data on the eigenvectors
V of the covariance matrix.

The `:VD` projection used in the DRS transform. Similar to the `:V` projection,
but the eigenvectors are multiplied by the squared inverse of the eigenvalues D.

The `:VDV` projection used in the SDS transform. Similar to the `:VD` transform,
but the data is projected back to the basis of the original variables using the Vᵀ matrix.

See [https://geostatisticslessons.com/lessons/sphereingmaf](https://geostatisticslessons.com/lessons/sphereingmaf)
for more details about these three variants of eigenanalysis.

# Examples

```julia
EigenAnalysis(:V)
EigenAnalysis(:VD)
EigenAnalysis(:VDV)
```
"""
struct EigenAnalysis <: Transform
  proj::Symbol

  function EigenAnalysis(proj)
    @assert proj ∈ (:V, :VD, :VDV) "invalid projection"
    new(proj)
  end
end

assertions(::Type{EigenAnalysis}) = [assert_continuous]

isrevertible(::Type{EigenAnalysis}) = true

function apply(transform::EigenAnalysis, table)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # original columns names
  names = Tables.columnnames(table)

  # table as matrix
  X = Tables.matrix(table)

  # center the data
  μ = mean(X, dims=1)
  Y = X .- μ

  # eigenanalysis of covariance
  S, S⁻¹ = eigenmatrices(transform, Y)

  # project the data
  Z = Y * S

  # table with transformed columns
  𝒯 = (; zip(names, eachcol(Z))...)
  newtable = 𝒯 |> Tables.materializer(table)

  newtable, (μ, S, S⁻¹)
end

function revert(::EigenAnalysis, newtable, cache)
  # transformed column names
  names = Tables.columnnames(newtable)

  # table as matrix
  Z = Tables.matrix(newtable)

  # retrieve cache
  μ, S, S⁻¹ = cache

  # undo projection
  Y = Z * S⁻¹

  # undo centering
  X = Y .+ μ

  # table with original columns
  𝒯 = (; zip(names, eachcol(X))...)
  𝒯 |> Tables.materializer(newtable)
end

function reapply(transform::EigenAnalysis, table, cache)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # original columns names
  names = Tables.columnnames(table)

  # table as matrix
  X = Tables.matrix(table)

  # retrieve cache
  μ, S, S⁻¹ = cache

  # center the data
  Y = X .- μ

  # project the data
  Z = Y * S

  # table with transformed columns
  𝒯 = (; zip(names, eachcol(Z))...)
  𝒯 |> Tables.materializer(table)
end

function eigenmatrices(transform, Y)
  proj = transform.proj

  Σ = cov(Y)
  λ, V = eigen(Σ)

  if proj == :V
    S   = V
    S⁻¹ = transpose(V)
  elseif proj == :VD
    Λ   = Diagonal(sqrt.(λ))
    S   = V * inv(Λ)
    S⁻¹ = Λ * transpose(V)
  elseif proj == :VDV
    Λ   = Diagonal(sqrt.(λ))
    S   = V * inv(Λ) * transpose(V)
    S⁻¹ = V * Λ * transpose(V)
  end

  S, S⁻¹
end

"""
    PCA()

The PCA transform is a shortcut for
`ZScore() → EigenAnalysis(:V)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).
"""
PCA() = ZScore() → EigenAnalysis(:V)

"""
    DRS()

The DRS transform is a shortcut for
`ZScore() → EigenAnalysis(:VD)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).
"""
DRS() = ZScore() → EigenAnalysis(:VD)

"""
    SDS()

The SDS transform is a shortcut for
`ZScore() → EigenAnalysis(:VDV)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).
"""
SDS() = ZScore() → EigenAnalysis(:VDV)


text:
struct MapsTo{X,Y}
    x::X
    y::Y
end

export ↦, mapsto

mapsto(x, y) = x ↦ y

↦(x::X, y::Y) where {X,Y} = MapsTo{X,Y}(x, y)

Base.first(t::MapsTo) = t.x
Base.last(t::MapsTo) = t.y

Base.Pair(t::MapsTo) = t.x => t.y

Base.show(io::IO, t::MapsTo) = print(t.x, " ↦ ", t.y)

logdensity(d, t::MapsTo) = logdensity(d, t.y)


text:
#######################################################################
########################### SNAP ######################################
#######################################################################

function virial(c:: Configuration, p::SNAP)
    A = get_snap(c, p)
    virial = A[end-5:end, 1:(end-1)] * p.β[1:(end-1)]
    return sum(virial[1:3])
end

function virial(r:: Vector{Configuration}, p::SNAP)
    n = length(r)
    v = Vector{Real}(undef, n)
    for j = 1:n
        v[j] = virial(r[j], p)
    end
    return v
end

function virial_stress(c:: Configuration, p::SNAP)
    A = get_snap(c, p)
    virial = A[end-5:end, 1:(end-1)] * p.β[1:(end-1)]
    return virial
end

function virial_stress(r:: Vector{Configuration}, p::SNAP)
    n = length(r)
    v = [Vector{Real}(undef, 6) for i = 1:n]
    for i = 1:n
        v[i] = virial_stress(r[i], p)
    end
    return v
end


text:
# simple iid bootstrap
function bootstrap(data)
    n = size(data,1)
    data[rand(1:n,n),:]
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module EuclideanWrapper

using JavaCall

include("OnedWrapper/OnedWrapper.jl")
include("ThreedWrapper/ThreedWrapper.jl")
include("TwodWrapper/TwodWrapper.jl")

end


text:
function warp_and_intersect(moving, fixed, tfm::IdentityTransformation)
    if axes(moving) == axes(fixed)
        return moving, fixed
    end
    inds = intersect.(axes(moving), axes(fixed))
    return view(moving, inds...), view(fixed, inds...)
end

function warp_and_intersect(moving, fixed, tfm)
    moving = warp(moving, tfm)
    inds = intersect.(axes(moving), axes(fixed))
    return view(moving, inds...), view(fixed, inds...)
end

"""
    I, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation())

Find the best shift `I` (expressed as a tuple) aligning `moving` to `fixed`,
possibly after an initial transformation `initial_tfm` applied to `moving`.
`mm` is the sum-of-square errors at shift `I`.
`mxshift` represents the maximum allowed shift, in pixels.
`thresh` is a threshold on the minimum allowed overlap between `fixed` and the shifted `moving`,
expressed in pixels if `normalization=:pixels` and in sum-of-squared-intensity if `normalization=:intensity`.

One can compute an overall transformation by composing `initial_tfm` with the returned shift `I`.
"""
function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation())
    moving, fixed = warp_and_intersect(moving, fixed, initial_tfm)
    mms = mismatch(fixed, moving, mxshift; normalization=normalization)
    best_i = indmin_mismatch(mms, thresh)
    mm = mms[best_i]
    return best_i.I, ratio(mm, thresh, typemax(eltype(mm)))
end

"""
    itfm = arrayscale(ptfm, SD)

Convert the physical-space transformation `ptfm` into one, `itfm`, that operates on index-space
for arrays.
For example, suppose `ptfm` is a pure 3d rotation, but you want to apply it to an array
in which the sampling along the three axes is (0.5mm, 0.5mm, 2mm). Then by setting
`SD = Diagonal(SVector(1, 1, 4))` (the ratio of scales along each axis),
one obtains an `itfm` suitable for warping the array.

Any translational component of `ptfm` is interpreted in physical-space units, not index-units.

`SD` does not even have to be diagonal, if the array sampling is skewed.
The columns of `SD` should correspond to the physical-coordinate displacement
achieved by shifting by one array element along each axis of the array.
Specifically, `SD[:,j]` should be the physical displacement of one array voxel along dimension `j`.
"""
arrayscale(ptfm::AbstractAffineMap, SD::AbstractMatrix) =
    arrayscale(ptfm, LinearMap(SD))

arrayscale(ptfm::AbstractAffineMap, scale::LinearMap) =
    inv(scale) ∘ ptfm ∘ scale

arrayscale(ptfm::AbstractAffineMap, scale::UniformScaling) = ptfm

"""
    pt = pscale(it, SD)

Convert an index-scaled translation `it` into a physical-space translation `pt`.
See [`arrayscale`](@ref) for more information.
"""
pscale(t::Translation, SD::AbstractMatrix) = pscale(t, LinearMap(SD))
pscale(t::Translation, scale::LinearMap) = Translation(scale(t.translation))
pscale(t::Translation, scale::UniformScaling) = Translation(scale.λ*t.translation)

#returns new minbounds and maxbounds with range sizes change by fac
function scalebounds(minb, maxb, fac::Real)
    orng = maxb.-minb
    newradius = fac.*orng./2
    ctrs = minb.+(orng./2)
    return ctrs.-newradius, ctrs.+newradius
end

"""
minb, maxb = default_linmap_bounds(img::AbstractArray{T,N}; dmax=0.05, ndmax=0.05) where {T,N}

Returns two matrices describing a search space of linear transformation matrices.
(Linear transformation matrices can encode rotations, scaling, shear, etc)
`minb` and  `maxb` contain the minimum and maximum acceptable values of an
NxN transformation matrix.  The center of the space is the identity matrix.
The size of the space can be specified with `dmax` and `ndmax` kwargs.
These represent the maximum (absolute) difference from the identity matrix for elements
along the diagonal and off the diagnonal, respectively.
e.g. `dmax=0.05` implies that diagonal elements can range from 0.95 to 1.05.
The space is centered on the identity matrix
"""
function default_linmap_bounds(img::AbstractArray{T,N}; dmax=0.05, ndmax=0.05) where {T, N}
    deltas = fill(abs(ndmax), N,N)
    for i=1:N
        deltas[i,i] = abs(dmax)
    end
    return Matrix(1.0*I,N,N).-deltas, Matrix(1.0*I,N,N).+deltas
end

"""
m = default_lin_minwidths(img::AbstractArray{T,N}; dmin=1e-3, ndmin=1e-3) where {T,N}

Returns a NxN matrix describing granularity of a search space of linear transformation matrices.
This can be useful for setting the `minwidth` parameter of QuadDIRECT when performing a
full affine registration. `dmin` and `ndmin` set the tolerances for diagonal and
off-diagonal elements of the linear transformation matrix, respectively.
"""
function default_lin_minwidths(img::AbstractArray{T,N}; dmin=1e-5, ndmin=1e-5) where {T, N}
    mat = fill(abs(ndmin), N,N)
    for i=1:N
        mat[i,i] = abs(dmin)
    end
    return mat[:]
end

"""
    θ = default_minrot(ci::CartesianIndices, SD=I; Δc=0.1)

Compute the rotation `θ` that results in largest change in coordinates
of size `Δc` (in pixels) for any index in `ci`.
"""
function default_minrot(ci::CartesianIndices, SD=I; Δc=0.1)
    L = -Inf
    for x in CornerIterator(ci)
        x′ = SD*SVector(Tuple(x))  # position of corner point in physical space
        L = max(L, norm(x′))
    end
    S2 = SD'*SD
    if SD == I
        λ = 1
    else
        F = eigen(S2)
        λ = minimum(F.values)
    end
    ℓ = sqrt(λ)*Δc
    return 2*asin(ℓ/(2*L))
end

default_minwidth_rot(ci::CartesianIndices{2}, SD=I; kwargs...) =
    [default_minrot(ci, SD; kwargs...)]
function default_minwidth_rot(ci::CartesianIndices{3}, SD=I; kwargs...)
    θ = default_minrot(ci, SD; kwargs...)
    return [θ, θ, θ]
end

"""
    getSD(A::AbstractArray)

If your image is not uniformily sampled, use this to get the `SD` matrix, which represents spacing along all axes of an image.

# Examples
```julia-repl
julia> myimage
Normed ImageMeta with:
  data: 3-dimensional AxisArray{N2f14,3,...} with axes:
    :x, 0.0 μm:0.71 μm:6.39 μm
    :l, 0.0 μm:0.71 μm:6.39 μm
    :z, 0.0 μm:6.2 μm:55.8 μm
And data, a 10×10×10 Array{N2f14,3} with eltype Normed{UInt16,14}
  properties:
    imagineheader: <suppressed>

julia> getSD(myimage)
3×3 SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
 0.71  0.0   0.0
 0.0   0.71  0.0
 0.0   0.0   6.2
```
"""
getSD(A::AbstractArray) = getSD(spacedirections(A)) #takes advantage of how spacedirections automatically strips out the time component
function getSD(sd::NTuple{N,Tuple}) where N
    SD = zeros(N,N)
    denom = oneunit(sd[1][1]) #dividing by a unit should remove unit inconsistancies
    for i = 1:N
        for j = 1:N
            SD[i,j] = sd[i][j]/denom
        end
    end
    return SMatrix{N,N}(SD) #returns static matrix for faster type inferrence
end

"""
    imgq = qsmooth(img)
    imgq = qsmooth(T, img)

Create a smoothed version of `img`, smoothing with the kernel of a quadratic B-spline.
Use this on your `fixed` image in preparation for registration, and pass `presmoothed`
as an option. (Do not smooth `moving`.)

`T` allows you to specify the output eltype (default `Float32`).
"""
function qsmooth(::Type{T}, img::AbstractArray{T2,N}) where {T,N,T2}
    kern1 = centered(T[1/8, 3/4, 1/8])   # quadratic B-spline kernel
    return imfilter(img, kernelfactors(ntuple(i->kern1, Val(N))))
end
qsmooth(img::AbstractArray) = qsmooth(Float32, img)

function qinterp(::Type{T}, moving) where T
    widen1(ax) = first(ax)-1:last(ax)+1

    # This sets up an interpolant object *without* using Interpolations.prefilter
    # Prefiltering is a form of "anti-smoothing," meaning that for quadratic interpolation
    # it defines a coefficient array so that, when you interpolate (aka, smooth),
    # it reconstructs the original values. In this case, we *want* interpolation to
    # perform smoothing, because the user has applied the same smoothing to `fixed`
    # by calling `qsmooth`.
    # Build the coefficient array. Quadratic interpolation requires one beyond-the-edge
    # element, which we set to NaN
    nanT = convert(T, NaN)
    mmpad = PaddedView(nanT, of_eltype(T, moving), widen1.(axes(moving)))
    # Create the interpolant, bypassing prefiltering
    return extrapolate(Interpolations.BSplineInterpolation(T, T.(mmpad), BSpline(Quadratic(Free(OnCell()))), axes(moving)), nanT)
end

#     θ = Inf
#     # R = I + [0 -Δθ; Δθ 0] is a rotmtrx for infinitesimal Δθ
#     # `dRdθ` is the "slope" of the rotation, which allows us to compute
#     # the
#     dRdθ = SD \ (@SMatrix([0 -1; 1 0]) * SD)
#     for c in CornerIterator(ci)
#         dcdθs = dRdθ*SVector(Tuple(c))
#         for dcdθ in dcdθs
#             θ = min(θ, abs(Δc/dcdθ))
#         end
#     end
#     return (θ,)
# end

# function default_minwidth_rot(I::CartesianIndices{3}, SD; d=0.1)
#     slice2d(I, i1, i2) = CartesianIndices((I.indices[i1], I.indices[i2]))
#     return (default_minwidth_rot(slice2d(I, 2, 3), SD[2:3, 2:3]; d=d)...,
#             default_minwidth_rot(slice2d(I, 1, 3), SD[[1,3], [1,3]]; d=d)...,
#             default_minwidth_rot(slice2d(I, 1, 2), SD[1:2, 1:2]; d=d)...)
# end

#sets splits based on lower and upper bounds
function _analyze(f, lower, upper; kwargs...)
    splits = ([[lower[i]; lower[i]+(upper[i]-lower[i])/2; upper[i]] for i=1:length(lower)]...,)
    QuadDIRECT.analyze(f, splits, lower, upper; kwargs...)
end


text:
using ViewerGL
using Distributed
if nprocs() == 1
    addprocs(3)
	@info "adding 3 more CPUs"
end
using Test
using Logging
using SparseArrays
@everywhere using LarSurf

block_size = [2, 2, 2]
segmentation = LarSurf.tetris_brick()
LarSurf.lsp_setup(block_size)
larmodel = LarSurf.lsp_get_surface(segmentation)
V, FV = larmodel
FVtri = LarSurf.triangulate_quads(FV)
objlines = LarSurf.Lar.lar2obj(V, FVtri, "tetris_tri.obj")

ViewerGL.VIEW([
    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

Vs = LarSurf.Smoothing.smoothing_FV_taubin(V, FV, 0.4, -0.2, 2)

ViewerGL.VIEW([
    ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

objlines = LarSurf.Lar.lar2obj(Vs, FVtri, "tetris_tri_taubin.obj")


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# based on Oleg Mazurov's Java Implementation and Jeremy Zerfas' C implementation
# transliterated by Hamza Yusuf Çakır

global const preferred_num_blocks = 12

struct Fannkuch
    n::Int64
    blocksz::Int64
    maxflips::Vector{Int32}
    chksums::Vector{Int32}

    function Fannkuch(n, nthreads)
        nfact = factorial(n)

        blocksz = nfact ÷ (nfact < preferred_num_blocks ? 1 : preferred_num_blocks)
        maxflips = zeros(Int32,nthreads)
        chksums = zeros(Int32, nthreads)

        new(n, blocksz, maxflips, chksums)
    end
end

struct Perm
    p::Vector{Int8}
    pp::Vector{Int8}
    count::Vector{Int32}

    function Perm(n)
        p = zeros(Int8, n)
        pp = zeros(Int8, n)
        count = zeros(Int32, n)

        new(p, pp, count)
    end
end

Base.@propagate_inbounds function first_permutation(perm::Perm, idx)
    p = perm.p
    pp = perm.pp

    for i = 2:length(p)
        p[i] = i - 1
    end

    for i = length(p):-1:2
        ifact = factorial(i-1)
        d = idx ÷ ifact
        perm.count[i] = d
        idx = idx % ifact

        for j = 1:i
            pp[j] = p[j]
        end

        for j = 1:i
            p[j] = j+d <= i ? pp[j+d] : pp[j+d-i]
        end
    end
end

Base.@propagate_inbounds function next_permutation(perm::Perm)
    p = perm.p
    count = perm.count

    first = p[2]
    p[2]  = p[1]
    p[1]  = first

    i = 2
    while (count[i] + 1) >= i
        count[i] = 0
        i += 1

        next = p[1] = p[2]

        for j = 1:i-1
            p[j] = p[j+1]
        end

        p[i] = first
        first = next
    end
    count[i] += 1
end

Base.@propagate_inbounds @inline function count_flips(perm::Perm)
    p = perm.p
    pp = perm.pp

    flips = 1
    first = p[1] + 1
    if p[first] != 0

        copyto!(pp, 2, p, 2, length(p) - 1)

        while true
            flips += 1

            new_first = pp[first]
            pp[first] = first - 1

            if first > 3
                lo = 2; hi = first-1

                # see the note in Jeremy Zerfas' C implementation
                for k = 1:15
                    t = pp[lo]
                    pp[lo] = pp[hi]
                    pp[hi] = t

                    !((lo + 3) <= hi) && break
                    lo += 1
                    hi -= 1
                end
            end

            first = new_first + 1
            pp[first] == 0 && break
        end
    end

    return flips
end

Base.@propagate_inbounds function run_task(f::Fannkuch, perm::Perm, idxmin, idxmax)

    maxflips = 1
    chksum = 0

    let
        i = idxmin
        while true
            if perm.p[1] != 0
                flips = count_flips(perm)
                maxflips = max(maxflips, flips)
                chksum += iseven(i) ? flips : -flips
            end
            if i == idxmax
                break
            end
            i += 1
            next_permutation(perm)
        end
    end
    id = Threads.threadid()
    f.maxflips[id] = max(f.maxflips[id], maxflips)
    f.chksums[id] += chksum
end

function runf(f::Fannkuch)
    factn = factorial(f.n)

    Threads.@threads for idxmin=0:f.blocksz:factn-1
        perm = Perm(f.n)
        first_permutation(perm, idxmin)
        idxmax = idxmin + f.blocksz - 1
        @inbounds run_task(f, perm, idxmin, idxmax)
    end
end

function fannkuchredux(n)
    f = Fannkuch(n,Threads.nthreads())

    runf(f)

    # reduce results
    chk = sum(f.chksums)
    res = maximum(f.maxflips)

    println(chk, "\nPfannkuchen(", n, ") = ", res)
end

n = parse(Int, ARGS[1])
fannkuchredux(n)


text:
# Tests for the LazyDict type.
# Author: Noah C. Benson <n@nben.net>

@testset "LazyDict" begin
    syms = Symbol[:a, :b, :c, :d, :e, :f, :g, :h, :i, :j]
    nums = Real[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
    n = 100
    # First just test them like normal dicts (they must work the same);
    # then we check delays. We use these functions, where the order they
    # are called will produce different results.
    # (Note that it's not a good idea to modify values like this inside
    # anonymous functions.
    fnval = 0
    delayfn1 = () -> begin
        fnval += 1
        return fnval
    end
    delayfn2 = () -> begin
        fnval *= 2
        return fnval
    end
    @testset "LazyIdDict" begin
        compare_test(Air.LazyIdDict{Symbol,Real}(),
                     Base.IdDict{Symbol,Real}(),
                     syms, nums, n)
        compare_test(Air.LazyIdDict{Symbol,Real}(:b=>20, :d=>40, :e=>50),
                     Base.IdDict{Symbol,Real}(:b=>20, :d=>40, :e=>50),
                     syms, nums, n)
        d = Air.LazyIdDict{Symbol,Int}(:a => Delay{Int}(delayfn1),
                                       :b => Delay{Int}(delayfn2))
        fnval = 0
        @test (d[:a], d[:b]) == (1, 2)
        @test (d[:b], d[:a]) == (2, 1)
        d = Air.LazyIdDict{Symbol,Int}(:a => Delay{Int}(delayfn1),
                                       :b => Delay{Int}(delayfn2))
        fnval = 0
        @test (d[:b], d[:a]) == (0, 1)
        @test (d[:a], d[:b]) == (1, 0)
    end
    @testset "LazyDict" begin
        compare_test(Air.LazyDict{Symbol,Real}(),
                     Base.Dict{Symbol,Real}(),
                     syms, nums, n)
        compare_test(Air.LazyDict{Symbol,Real}(:b=>20, :d=>40, :e=>50),
                     Base.Dict{Symbol,Real}(:b=>20, :d=>40, :e=>50),
                     syms, nums, n)
        d = Air.LazyDict{Symbol,Int}(:a => Delay{Int}(delayfn1),
                                     :b => Delay{Int}(delayfn2))
        fnval = 0
        @test (d[:a], d[:b]) == (1, 2)
        @test (d[:b], d[:a]) == (2, 1)
        d = Air.LazyDict{Symbol,Int}(:a => Delay{Int}(delayfn1),
                                     :b => Delay{Int}(delayfn2))
        fnval = 0
        @test (d[:b], d[:a]) == (0, 1)
        @test (d[:a], d[:b]) == (1, 0)
    end
end


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

mutable struct DataLakeStoreAccountUpdateParameters <: SwaggerModel
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags
    properties::Any # spec type: Union{ Nothing, UpdateDataLakeStoreAccountProperties } # spec name: properties

    function DataLakeStoreAccountUpdateParameters(;tags=nothing, properties=nothing)
        o = new()
        validate_property(DataLakeStoreAccountUpdateParameters, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        validate_property(DataLakeStoreAccountUpdateParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type DataLakeStoreAccountUpdateParameters

const _property_map_DataLakeStoreAccountUpdateParameters = Dict{Symbol,Symbol}(Symbol("tags")=>Symbol("tags"), Symbol("properties")=>Symbol("properties"))
const _property_types_DataLakeStoreAccountUpdateParameters = Dict{Symbol,String}(Symbol("tags")=>"Dict{String, String}", Symbol("properties")=>"UpdateDataLakeStoreAccountProperties")
Base.propertynames(::Type{ DataLakeStoreAccountUpdateParameters }) = collect(keys(_property_map_DataLakeStoreAccountUpdateParameters))
Swagger.property_type(::Type{ DataLakeStoreAccountUpdateParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DataLakeStoreAccountUpdateParameters[name]))}
Swagger.field_name(::Type{ DataLakeStoreAccountUpdateParameters }, property_name::Symbol) =  _property_map_DataLakeStoreAccountUpdateParameters[property_name]

function check_required(o::DataLakeStoreAccountUpdateParameters)
    true
end

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


text:
using Distributions
using DifferentialEquations
using DataFrames
using CSV
using PyPlot

"""
	function mcmc(data,cols,p,Population,total::Int)

Metropolis-Hastings MCMC of the SIR ODE model assuming a Poisson likelihood
for the cumulative number of cases and deaths, and a noninformative prior
data = DataFrame of
cols = [active cases, deaths, recovered cases]
p[1:end-1] = array of parameters used in the model
p[end] = day of first infection relative to day 0 of data (Jan 22, 2020)
Population = population of region in data
total = number of MCMC steps
output = (loglikelihoods, parameters, maximum loglikelihood, max loglikelihood parameters)

"""
function mcmc(data::DataFrame,cols::Array,p::Array,Population::Float64,total::Int)

	chi = ll(data,cols,p,Population)
	accept = 0
	pt = copy(p)
	chiml = chi
	pml = copy(p)

	chiout = Array{Float64,1}(undef,total)
	pout = Array{Float64,2}(undef,total,length(p))

	for step in 1:total
		sampleparams!(pt,p)
		chit = ll(data,cols,pt,Population)
		if rand() < exp(chi - chit)
			chi = chit
			p = copy(pt)
			accept += 1
			if chi < chiml
				chiml = chi
				pml = copy(p)
			end
		end
		chiout[step] = chi
		pout[step,:] = p
	end
	return chiout, pout, accept/total, chiml, pml
end

# generate MCMC step guess from current parameter value r using logNormal distirbution for all
# parameters except start day with scale parameter manually adjusted so approximately a third of guesses are accepted
function sampleparams!(pt,p)
	for i in 1:length(p) - 1
		d = Distributions.LogNormal(log(p[i]),.002)
		pt[i] = rand(d)
	end
	pt[end] = p[end] + .05*randn()
end

# loglikelihood of data given parameters at r using a Poisson likelihood function
function ll(data,cols,p,Population)
	total = length(data[!,cols[1]])-1
	prediction = model(p,total,Population)
	chi = 0
	for set in 1:3
		for i in 1:total
			d = Distributions.Poisson(max(prediction[i,set],1.))
			chi -= loglikelihood(d,[data[i,cols[set]]])
		end
	end
	return chi
end

# produce arrays of cases and deaths
function model(pin,total,Population)
	u0 = [0.,0.,0.,1/Population,0.,0.]
	p = pin[1:end-1]
	if pin[end] < total - 1
		start = ceil(Int,pin[end])
		tspan = (pin[end],total-1)
		prob = ODEProblem(sirc!,u0,tspan,p)
		sol = solve(prob,saveat=collect(start:total-1))
		# align solutions to time 0 of data and reorder rows and cols
		prediction = zeros(total,3)
		for set in 1:3
			prediction[max(0,start)+1:end,set] .= sol[set,max(0,-start)+1:end]
		end
		return prediction*Population
	else
		return zeros(total,3)
	end
end

# SIRC ODE model using Differential Equations
function sirc!(du,u,p,t)
	C,D,R,Y,Z,U = u
	beta,sigmac,sigmar,gammad,gammar = p
	du[1] = dC = sigmac*Y - gammar*C - gammad*C
	du[2] = dD = gammad*C
	du[3] = dR = gammar*C
	du[4] = dY = beta*Y*(1-Z) - sigmac*Y - sigmar*Y
	du[5] = dZ = beta*Y*(1-Z)
	du[6] = dU = sigmar*Y
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
# time-set.jl
# timing tests for setindex!

# times reported below are all with julia 1.5
# on 2017 iMac (4.2GHz Intel Core i7) with Mojave 10.14.6

using ConstantArrays
#using FillArrays: Ones
using FillArrays: Ones, Trues # for 0.9
using BenchmarkTools

dim = (2^8,2^8,2^4)
x = randn(dim)
y = similar(x)

mask0 = trues(dim)
maskc = ConstantArray(dim)
maskf = Ones{Bool}(dim)
#const Trues = Ones{Bool, N, Axes} where {N, Axes} # lazy trues()
maskt = Trues(dim)
@assert mask0 == maskt

@assert mask0 == maskc == maskf == maskt

# function wrappers for timing tests

function setter!(y, mask, x) # setindex!
	y[mask] = x
end

function set_checker()
	@assert copyto!(y, x) ==
		setter!(y, mask0, x) ==
		setter!(y, maskc, x) ==
		setter!(y, maskf, x) ==
		setter!(y, maskt, x) == x
end

# time y[mask] = x when mask is trues()
# default setindex! operations are much slower than copyto!
if true
	@info "setter"
	set_checker()
	@btime copyto!($y, $x) # 0.65 ms (baseline)
	@btime setter!($y, $mask0, $x) # 1.6 ms
	@btime setter!($y, $maskc, $x) # 2.5 ms :(
	@btime setter!($y, $maskf, $x) # 2.3 ms :( before code update, # 0.6 after
	@btime setter!($y, $maskt, $x) # 2.3 ms :( before code update, # 0.6 after
end

# here is the solution for a ConstantArray:
# y[mask] = x
function Base.setindex!(y::AbstractArray{T,D}, x, mask::ConstantArray{Bool,D}) where {T,D}
	if mask.value == true
		copyto!(y, x)
	else
		y[mask] = x
	end
end

#=
# here is the essence of the solution for Ones in FillArrays:
# y[mask] = x when mask isa Trues
function Base.setindex!(y::AbstractArray{T,D}, x, mask::Trues{D}) where {T,D}
	copyto!(y, x)
end
=#

if true
	@info "updated setter!"
	set_checker()
	@btime setter!($y, $maskc, $x) # 0.65 ms :)
	@btime setter!($y, $maskf, $x) # 0.65 ms :)
	@btime setter!($y, $maskt, $x) # 0.65 ms :)
end

nothing