text:
module Pixels

using Images
using ArgParse

include("ui.jl")
include("bag.jl")
include("treap.jl")
include("sss.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--file", "-f"
        help = "filename for image to use as input"
        arg_type = String
        default = "images/Bachalpseeflowers.jpg"
    "--type", "-t"
        help = "function type to use: min, avg, or evo"
        arg_type = String
        default = "min"
end
parsed_args = parse_args(ARGS, s)
const target_file = parsed_args["file"]

const fntype = symbol(parsed_args["type"])
(fntype in [:min,:avg,:evo]) || error("Invalid function type provide. Got $fntype; expected min, avg, or evo.")

immutable NeighborStats
    n::Int
    avg::Rgb{Uint8}
    sumvar::Uint
    sqrt_sumvar::Uint8
    NeighborStats(n, avg, sumvar) = new(n, avg, sumvar, uint8(sqrt(sumvar)))
    NeighborStats() = new(0, Rgb{Uint8}(), 0, 0)
end

type Vertex
    x::Int
    y::Int
    color::GLPixel
    filled::Bool
    nominated::Bool
    bagindex::Int
    stats::NeighborStats
    neighbors::Vector{Vertex}
    t::Float64 # Used to preserve key uniqueness and break ties in <(Vertex, Vertex)
    function Vertex(x, y)
        neighbors = Vertex[]
        sizehint(neighbors, 8)
        new(x, y, GLPixel(), false, false, 0, NeighborStats(), neighbors, rand(Float64))
    end
end

Base.start(p::Rgb) = 1
Base.next{T}(p::Rgb{T}, state::Int) =
    state == 1 ? (p.r, 2) :
    state == 2 ? (p.g, 3) :
    state == 3 ? (p.b, 4) : throw("Iteration state error.")
Base.done(p::Rgb, state::Int) = state == 4

type Canvas
    grid::Array{Vertex}
    frontier
    function Canvas(width, height, frontier)
        grid = [Vertex(x, y) for y in 1:height, x in 1:width]

        for x in 1:width, y in 1:height
            for nx in max(x - 1, 1):min(x + 1, width), ny in max(y - 1, 1):min(y + 1, height)
                !(nx == x && ny == y) && push!(grid[y, x].neighbors, grid[ny, nx])
            end
        end
        new(grid, frontier)
    end
end
frontier(canvas) = canvas.frontier
Base.getindex(c::Canvas, a, b) = c.grid[a, b]
neighbors(vertex::Vertex) = vertex.neighbors
isnominated(vertex::Vertex) = vertex.bagindex != 0 || vertex.nominated # TODO: Fix.

function nominate!(vertex, frontier::Bag{Vertex})
    @assert !isnominated(vertex)

    # Bag add
    push!(frontier, vertex)
    vertex.bagindex = length(frontier)
end

function withdraw!(vertex, frontier::Bag{Vertex})
    @assert isnominated(vertex)

    # Bag remove
    index = vertex.bagindex
    deleteat!(frontier, index)
    length(frontier) >= index && (frontier[index].bagindex = index)
    vertex.bagindex = 0
end

function nominate!(vertex, frontier::Treap)
    @assert !vertex.nominated
    add!(frontier, vertex)
    vertex.nominated = true
end

function withdraw!(vertex, frontier::Treap)
    @assert vertex.nominated
    remove!(frontier, vertex)
    vertex.nominated = false
end

hasemptyneighbors(vertex) = vertex.stats.n < length(neighbors(vertex))
emptyneighbors(vertex) = filter((v) -> !v.filled, neighbors(vertex))
function randemptyneighbor(vertex)
    prob = 1.0
    local neighbor
    for v in neighbors(vertex)
        if !v.filled
            if rand() < prob
                neighbor = v
            end
            prob /= 2
        end
    end
    neighbor
end

function fill!(pixels, frontier, vertex, color, outer::Bool)
    # TODO: Don't pass pixels in.
    pixels[vertex.y, vertex.x] = color

    vertex.color = color
    vertex.filled = true
    if outer
        # Outer frontier
        for v in neighbors(vertex)
            # Withdraw in order to re-insert neighbors already on the frontier,
            # as their stats have changed in a way that may influence their
            # storage in the frontier
            isnominated(v) && withdraw!(v, frontier)
            stats!(v)
            !isnominated(v) && !v.filled && nominate!(v, frontier)
        end
        isnominated(vertex) && withdraw!(vertex, frontier)
    else
        # Inner frontier
        for v in neighbors(vertex)
            # NOTE: Does something like the re-insertion above need to happen here as well?
            stats!(v)
            isnominated(v) && !hasemptyneighbors(v) && withdraw!(v, frontier)
        end
        hasemptyneighbors(vertex) && nominate!(vertex, frontier)
    end
end

# TODO: To calculate average without overflow: (R - L)/2 + L
function stats!(vertex)
    n = zero(vertex.stats.n)
    avg = var = Rgb{Float64}() #zero(vertex.stats.avg)

    for v in neighbors(vertex)
        if v.filled
            n += 1
            avg += v.color
            # Promote types before squaring, since Rgb{Uint8} * Rgb{Uint8} => Rgb{Uint8}
            var += oftype(var, v.color) * v.color
        end
    end

    if n == 0
        vertex.stats = NeighborStats()
    else
        avg = avg / n
        var = var / n - avg * avg
        vertex.stats = NeighborStats(n, Rgb{Uint8}(uint8(avg.r), uint8(avg.g), uint8(avg.b)), uint(sum(var)))
    end
end

r() = rand() - 0.5
mutate(color::GLPixel, scale) = GLPixel(
    clamp(color.r + int(scale * r()), typemin(color.r), typemax(color.r)),
    clamp(color.g + int(scale * r()), typemin(color.g), typemax(color.g)),
    clamp(color.b + int(scale * r()), typemin(color.b), typemax(color.b))
)

function nearest(frontier::Bag{Vertex}, color)
    best = frontier[1]
    best_dist::Float64 = Inf
    for v in frontier
        d = sqdiff(v.stats.avg, color) + v.stats.sumvar
        if d < best_dist
            best = v
            best_dist = d
        end
    end
    best
end

function loadimage(imgpath::String)
    im = convert(Array, imread(imgpath))
    w, h = size_spatial(im)
    colors = Pixel[]
    for i in 1:w, j in 1:h
        push!(colors, Pixel(uint8(im[i, j].r * 255), uint8(im[i, j].g * 255), uint8(im[i, j].b * 255)))
    end
    colors
end

SPEED = 1000

function evolve(pixels)
    (height, width) = size(pixels)
    canvas = Canvas(width, height, Bag{Vertex}(height * width))

    fill!(pixels, frontier(canvas), canvas[div(height, 2), div(width, 2)], GLPixel(255, 255, 255), true)
    for i in 1:length(canvas.grid) - 1
        vertex = frontier(canvas)[rand(1:length(frontier(canvas)))]
        color = filter(n -> n.filled, neighbors(vertex))[1].color
        fill!(pixels, frontier(canvas), vertex, mutate(color, 15), true)
        i % SPEED == 0 && produce(pixels)
    end
    pixels
end

function gencolors(imgpath::String)
    colors = shuffle!(loadimage(imgpath))
    sort!(colors, alg=QuickSort, lt=(p::Pixel, q::Pixel) -> hue(p) < hue(q))
    colors
end

function gencolors(n::Int)
    colors = [rand(GLPixel) for i in 1:n]
    # sort!(colors, alg=QuickSort, lt=(p::Pixel, q::Pixel) -> hue(p) < hue(q))
    sort!(colors, rev=true, alg=QuickSort, lt=(p::Pixel, q::Pixel) -> p.r < q.r && p.g < q.g && p.b < q.b)
    # sort!(colors, rev=true, alg=QuickSort, by=(p::Pixel) -> p.b + p.g^2)
    colors
end

function placeavg(pixels)
    (height, width) = size(pixels)
    canvas = Canvas(width, height, Treap{Vertex}())
    # colors = gencolors(length(pixels))
    colors = gencolors(target_file)
    # colors = gencolors(1000)
    i = 1
    fill!(pixels, frontier(canvas), canvas[div(height, 2), div(width, 2)], colors[i += 1], true)
    while i < length(pixels)
        color = colors[mod1(i, length(colors))]
        vertex = nearest(frontier(canvas), color)
        fill!(pixels, frontier(canvas), vertex, color, true)
        i % SPEED == 0 && produce(pixels)
        i += 1
    end
    pixels
end

function placemin(pixels)
    (height, width) = size(pixels)
    canvas = Canvas(width, height, Treap{Vertex}())
    # colors = gencolors(length(pixels))
    # colors = gencolors("Bachalpseeflowers.jpg")
    colors = gencolors(target_file)
    i = 1

    fill!(pixels, frontier(canvas), canvas[div(height, 2), div(width, 2)], colors[i += 1], false)
    # fill!(pixels, frontier(canvas), canvas[200, 200], colors[i += 1], false)
    # fill!(pixels, frontier(canvas), canvas[1, 1], colors[i += 1], false)

    while i < length(pixels)
        color = colors[mod1(i, length(colors))]
        candidate = nearest(frontier(canvas), color)
        vertex = randemptyneighbor(candidate)
        fill!(pixels, frontier(canvas), vertex, color, false)
        i % SPEED == 0 && produce(pixels)
        i += 1
    end
    pixels
end

fn = fntype == :min ? placemin :
     fntype == :avg ? placeavg :
     fntype == :evo ? evolve : throw("fntype must be :min, :avg, or :evo.")

if fntype == :min || fntype == :evo
    # Minimum selection
    Base.size(p::Vertex, n) = n == 1 ? 3 : error("Vertices are one-dimensional.")
    Base.getindex(p::Vertex, n) = p.color[n]
    <(p::Vertex, q::Vertex) = p.color == q.color ? p.t < q.t : shuffless(p.color, q.color)
    <(p::Vertex, q::Rgb) = shuffless(p.color, q)
    <(p::Rgb, q::Vertex) = shuffless(p, q.color)
    dist_sq{T}(p::Rgb{T}, q::Rgb{T}) = sqdiff(p, q)
    dist_sq{T}(p::Vertex, q::Rgb{T}) = sqdiff(p.color, q)
    dist_sq{T}(p::Rgb{T}, q::Vertex) = sqdiff(p, q.color)
else
    # Average selection
    Base.size(p::Vertex, n) = n == 1 ? 4 : error("n must be 1.")
    Base.getindex(p::Vertex, n) = n < 4 ? uint8(p.stats.avg[n]) : n == 4 ? uint8(p.stats.sqrt_sumvar) : error("Index out of bounds: $n")
    # If the colors are equal, compare t values; we need this because the treap enforces uniqueness.
    <(p::Vertex, q::Vertex) = shuffless(p, q) ? true : shuffmore(p, q) ? false : p.t < q.t
    <(p::Vertex, q::Rgb) = shuffless(p, q)
    <(p::Rgb, q::Vertex) = shuffless(p, q)
    dist_sq{T}(p::Rgb{T}, q::Rgb{T}) = sqdiff(p, q)
    dist_sq{T}(p::Vertex, q::Rgb{T}) = sqdiff(p.stats.avg, q) + p.stats.sumvar
    dist_sq{T}(p::Rgb{T}, q::Vertex) = sqdiff(p, q.stats.avg) + q.stats.sumvar

    # Avg color: Hacked version with p[4] == 0
    Base.start(p::Rgb) = 1
    Base.next{T}(p::Rgb{T}, state::Int) =
        state == 1 ? (p.r, 2) :
        state == 2 ? (p.g, 3) :
        state == 3 ? (p.b, 4) :
        state == 4 ? (zero(T), 5) : throw("Iteration state error.")
    Base.done(p::Rgb, state::Int) = state == 5
    Base.getindex{T}(p::Rgb{T}, n::Int) = n == 1 ? p.r : n == 2 ? p.g : n == 3 ? p.b : n == 4 ? zero(T) : throw("Color indexing error.")
    Base.size(p::Rgb) = (4,)
    Base.size(p::Rgb, n) = n == 1 ? 4 : throw("Invalid dimension.")
end

# placemin(zeros(GLPixel, 100, 100))
display(1000, 750, fn; title="Pixels")

end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
using BulletsMaggots, ProgressMeter, Statistics

build_table()
function jonah_experiment(reps::Int = 1000)
    println("Playing $reps games")
    t = zeros(Int, reps)
    PM = Progress(reps)
    for i = 1:reps
        t[i] = jonah_solver(random_code(), false, 10)
        next!(PM)
    end
    println("Average steps to solve:   ", mean(t))
    println("Median steps to solve:    ", median(t))
    println("Standard deviation:       ", std(t))
    println("Range in number of steps: ", minimum(t), "-", maximum(t))
end

ext: jl


text:
export
    WignerFunction,
    WignerSurface,
    wigner

#=
    Wigner function by Laguerre Polynominal
=#

function calc_wigner(
    x_range::AbstractRange,
    p_range::AbstractRange,
    dim::Integer
)
    𝐰 = Array{ComplexF64,4}(undef, dim, dim, length(x_range), length(p_range))
    xs = collect(x_range)
    ps = collect(p_range)
    @sync for m in 1:dim
        Threads.@spawn for n in 1:dim
            𝐰[m, n, :, :] .=
                gaussian_function(xs, ps) .*
                coefficient_of_wave_function(m, n) .*
                z_to_power(m, n, xs, ps) .*
                laguerre(m, n, xs, ps)
        end
    end

    return 𝐰
end

function load_wigner(
    bin_path::String,
    x_range::AbstractRange,
    p_range::AbstractRange,
    dim::Integer
)
    return load_𝐰(bin_path, x_range, p_range, dim)
end

function create_wigner(
    x_range::AbstractRange,
    p_range::AbstractRange,
    dim::Integer
)
    bin_name = gen_wigner_bin_name(x_range, p_range, dim)
    hash = my_artifact_hash(bin_name, my_artifacts[])

    if !isnothing(hash) && my_artifact_exists(hash)
        return load_wigner(my_artifact_path(hash), x_range, p_range, dim)
    end

    𝐰 = calc_wigner(x_range, p_range, dim)
    save_𝐰(bin_name, 𝐰)

    return 𝐰
end

"""
    WignerFunction

`WignerFunction` to calculate Wigner function for a quantum state in gevin `x` and `p` range.

## Arguments

* `x_range`: First physical domain range in phase space.
* `p_range`: Second physical domain range in phase space.
* `dim`: Photon number truncation of the quantum state.

## Example

```julia-repl
julia> xs = -1:0.1:1; ps = -1:0.1:1; dim=35;

julia> wf = WignerFunction(xs, ps, dim=dim);

julia> wf(VacuumState(Matrix, dim=dim));
```
"""
mutable struct WignerFunction{T<:Integer, U<:AbstractRange}
    𝐰::Array{ComplexF64,4}
    x_range::U
    p_range::U
    dim::T

    function WignerFunction(
        x_range::U,
        p_range::U,
        dim::T
    ) where {T<:Integer, U<:AbstractRange}
        𝐰 = create_wigner(x_range, p_range, dim)

        return new{T, U}(𝐰, x_range, p_range, dim)
    end
end

function WignerFunction(x_range::AbstractRange, p_range::AbstractRange; dim)
    return WignerFunction(x_range, p_range, dim)
end

struct WignerSurface{T<:AbstractRange}
    x_range::T
    p_range::T
    𝐰_surface::Matrix{Float64}

    function WignerSurface(
        x_range::T,
        p_range::T,
        𝐰_surface::Matrix{Float64}
    ) where {T<:AbstractRange}
        return new{T}(x_range, p_range, 𝐰_surface)
    end
end

function (wf::WignerFunction)(ρ::AbstractMatrix{T}) where {T}
    𝐰_surface = Matrix{Float64}(undef, length(wf.x_range), length(wf.p_range))
    @sync for i in 1:length(wf.x_range)
        Threads.@spawn for j in 1:length(wf.p_range)
            𝐰_surface[i, j] = real(sum(ρ .* wf.𝐰[:, :, i, j]))
        end
    end

    return WignerSurface(wf.x_range, wf.p_range, 𝐰_surface)
end

# TODO: (wf::WignerFunction)(state::StateVector) = wf(ρ(state))

"""
    wigner(ρ::AbstractMatrix{T}, x_range::AbstractRange, p_range::AbstractRange)

Calculate the Wigner fuction of the quantum state in the gevin `x` and `p` range.

## Arguments

* `ρ`: Density matrix of a quantum state.
* `x_range`: First physical domain range in phase space.
* `p_range`: Second physical domain range in phase space.

## Example

```julia-repl
julia> xs = -1:0.1:1; ps = -1:0.1:1;

julia> wigner(VacuumState(Matrix, dim=35), xs, ps);
```
"""
function wigner(ρ::AbstractMatrix{T}, x_range::AbstractRange, p_range::AbstractRange) where {T}
    dim = size(ρ, 1)
    wf = WignerFunction(x_range, p_range, dim=dim)

    return wf(ρ)
end

# TODO: wigner(v::AbstractVector{T}, x_range::AbstractRange, p_range::AbstractRange) where {T}


text:
abstract type MomentSettings end

"Moment calculation settings for histogram methods"
struct HistogramSettings{S<:Integer,T<:AbstractVector{S},S1<:Real,T1<:AbstractVector{S1}} <: MomentSettings
    time_shift_sample_points::T
    bin_edges::T1
    ti_grid
    n_time_points
    n_evaluation_points
end

function HistogramSettings(time_shift_sample_points,bin_edges)
    n_time_points = length(time_shift_sample_points)
    n_evaluation_points = length(bin_edges)
    ti_grid = make_grid(time_shift_sample_points,
                        bin_edges)
    return HistogramSettings(time_shift_sample_points,
                             bin_edges,
                             ti_grid,
                             n_time_points,
                             n_evaluation_points)
end

struct Moments{T<:AbstractFloat}
    M1::Matrix{T}
    M2::Matrix{T}
end

"Moments structure"
struct MomentSolution{T<:AbstractFloat}
    moments::Moments{T}
    errors::Union{Nothing,Matrix{T}}
    observation::Union{Observation,EnsembleObservation}
    evaluation_points
    settings::MomentSettings
end
#TODO: Add types
#NOTE: Make τ vector a field?

"Ensemble moments structure"
struct EnsembleMoments{T<:AbstractFloat}
    M1::Vector{Matrix{T}}
    M2::Vector{Matrix{T}}
end

function MomentSolution(observation,settings::MomentSettings)

    moments,errors = build_moments(observation, settings)
    return MomentSolution(moments,
                           errors,
                           observation,
                           center(settings.bin_edges),
                           settings)
end

function make_grid(tau_vector::AbstractVector,edge_vector::LinRange)
    return broadcast((x,y)->(x,y),tau_vector,(1:edge_vector.lendiv)')
end

center(x::LinRange) = (x[1:end-1]+x[2:end])/2

"Calculating moments with histograms"
function build_moments(observation::Observation, settings::MomentSettings)
    moments, errors = moments_map(observation.X,
                                   settings.time_shift_sample_points,
                                   settings.bin_edges,
                                   settings.ti_grid)
    errors = nothing
    return moments, errors
end

"Calculating ensemble moments with histograms"
function build_moments(observation::EnsembleObservation, settings::MomentSettings)
    ensembleMoments = moments_map(observation.X,
                                 settings.time_shift_sample_points,
                                 settings.bin_edges,
                                 settings.ti_grid)

    errors = nothing

    M1mean = mean(ensembleMoments.M1)
    M2mean = mean(ensembleMoments.M2)

    return Moments(M1mean, M2mean), errors
end

# BUG: When bins have no entries they give NaN
function get_moments(X, iX, tau)
    inc = X[iX .+ tau] .- X[iX]
    M1 = mean(inc)
    residuals = inc .- M1
    M2 = mean(residuals.^2)
    #KSerror = KS(residuals/(M2*sqrt(2)))
    KSerror = nothing
    return M1, M2, KSerror
end
function KS(Y::AbstractArray{T}) where T
    Ysort = sort(Y)
    D = zero(T)
    N = length(Y)
    for (i,v) in enumerate(Ysort)
        Dtmp = max(NCDF(v) - convert(T,(i-one(T))/N),
                   convert(T,i/N) - NCDF(v))
        D = max(D,Dtmp)
    end
    return D
end
function NCDF(x::T) where T<:Real
    return oftype(x,0.5)*(one(T)+erf(x))
end
function erf(x::T) where T<:Real
    if x < zero(T)
        return -erfpos(-x)
    else
        return erfpos(x)
    end
end
function erfpos(x::T) where T<:Real
    a1 = oftype(x,0.278393)
    a2 = oftype(x,0.230389)
    a3 = oftype(x,0.000972)
    a4 = oftype(x,0.078108)
    x1 = muladd(x,a4,a3)
    x2 = muladd(x,x1,a2)
    x3 = muladd(x,x2,a1)
    x4 = muladd(x,x3,one(T))
    x5 = x4*x4
    return one(T) - inv(x5*x5)
end
find_bin(r::LinRange,x::Real) = floor(Int, ((x-r.start)/(r.stop-r.start)*r.lendiv + 1))
in_range(r::LinRange,x::Real) = (x >= r.start) && (x <= r.stop)
function moments_map(X::AbstractVector,tau_vector,edge_vector::LinRange,ti_grid)

    bin_index = map(X) do y
        in_range(edge_vector,y) ? find_bin(edge_vector,y) : 0
    end

    bin_groups = map(1:edge_vector.lendiv) do i
        findall(bin_index[1:end-maximum(tau_vector)].==i)
    end

    M1M2 = map(ti_grid) do ti
        get_moments(X,bin_groups[ti[2]],ti[1])
    end

    M1 = broadcast(first,M1M2)
    M2 = broadcast(x->x[2],M1M2)
    errors = broadcast(last,M1M2)

    return Moments(M1, M2), errors
end
function moments_map(X::AbstractMatrix{T},tau_vector,edge_vector::LinRange,ti_grid) where T
    #HACK: Do I need this explicit type declaration?
    vecMoments::Vector{Moments{T}} = map(eachcol(X)) do y
        moments_map(y,tau_vector,edge_vector,ti_grid)[1]
    end
    l = length(vecMoments)

    return EnsembleMoments(map(x->vecMoments[x].M1,1:l),
                           map(x->vecMoments[x].M2,1:l))
end


text:
close("all")
transD_GP.getchi2forall(opt)
ax = gcf().axes;
ax[2].set_ylim(15, 40)
χ² = aem.ndatalow + aem.ndatahigh
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
opt.xall[:] .= zall
transD_GP.plot_posterior(aem, opt, burninfrac=0.5, figsize=(10,6),
    qp1=0.01, qp2=0.99, nbins=50, vmaxpc=1.0)
ax = gcf().axes
ax[1].invert_xaxis()
ax[1].step(log10.(ρ[2:end]), z[2:end], color="k", linewidth=3)
ax[1].step(log10.(ρ[2:end]), z[2:end], color="y", linewidth=1.5)
ax[2].set_ylim(280,0)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using Test

include("isogram.jl")

@testset "empty string" begin
    @test isisogram("")
end

@testset "isogram with only lower case characters" begin
    @test isisogram("isogram")
end

@testset "word with one duplicated character" begin
    @test !isisogram("eleven")
end

@testset "longest reported english isogram" begin
    @test isisogram("subdermatoglyphic")
end

@testset "word with duplicated character in mixed case" begin
    @test !isisogram("Alphabet")
end

@testset "hypothetical isogrammic word with hyphen" begin
    @test isisogram("thumbscrew-japingly")
end

@testset "isogram with duplicated non letter character" begin
    @test isisogram("Hjelmqvist-Gryb-Zock-Pfund-Wax")
end

@testset "made-up name that is an isogram" begin
    @test isisogram("Emily Jung Schwartzkopf")
end

@testset "word with one duplicated character from the end of the alphabet" begin
    @test !isisogram("zzyzx")
end

@testset "word with duplicated character in mixed case, lowercase first" begin
    @test !isisogram("alphAbet")
end

@testset "hypothetical word with duplicated character following hyphen" begin
    @test !isisogram("thumbscrew-jappingly")
end

@testset "isogram with duplicated hyphen" begin
    @test isisogram("six-year-old")
end

@testset "duplicated character in the middle" begin
    @test !isisogram("accentor")
end

@testset "same first and last characters" begin
    @test !isisogram("angola")
end


text:
# ---
# title: 1638. Count Substrings That Differ by One Character
# id: problem1638
# author: Tian Jun
# date: 2021-01-26
# difficulty: Medium
# categories: Hash Table, String, Trie, Rolling Hash
# link: <https://leetcode.com/problems/count-substrings-that-differ-by-one-character/description/>
# hidden: true
# ---
#
# Given two strings `s` and `t`, find the number of ways you can choose a non-
# empty substring of `s` and replace a **single character** by a different
# character such that the resulting substring is a substring of `t`. In other
# words, find the number of substrings in `s` that differ from some substring in
# `t` by **exactly** one character.
#
# For example, the underlined substrings in `" _compute_ r"` and `" _computa_
# tion"` only differ by the `'e'`/`'a'`, so this is a valid way.
#
# Return _the number of substrings that satisfy the condition above._
#
# A **substring** is a contiguous sequence of characters within a string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aba", t = "baba"
#     Output: 6
#     Explanation: The following are the pairs of substrings from s and t that differ by exactly 1 character:
#     (" _a_ ba", " _b_ aba")
#     (" _a_ ba", "ba _b_ a")
#     ("ab _a_ ", " _b_ aba")
#     ("ab _a_ ", "ba _b_ a")
#     ("a _b_ a", "b _a_ ba")
#     ("a _b_ a", "bab _a_ ")
#     The underlined portions are the substrings that are chosen from s and t.
#
#
#  **Example 2:**
#
#
#
#     Input: s = "ab", t = "bb"
#     Output: 3
#     Explanation: The following are the pairs of substrings from s and t that differ by 1 character:
#     (" _a_ b", " _b_ b")
#     (" _a_ b", "b _b_ ")
#     (" _ab_ ", " _bb_ ")
#     The underlined portions are the substrings that are chosen from s and t.
#
#
# **Example 3:**
#
#
#
#     Input: s = "a", t = "a"
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: s = "abe", t = "bbc"
#     Output: 10
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length, t.length <= 100`
#   * `s` and `t` consist of lowercase English letters only.
#
#
## @lc code=start
using LeetCode

function count_substrings(s::String, t::String)
    function count_substr(s::SubString{String}, t::String)
        res, len_s, len_t = 0, length(s), length(t)
        for i in 1:len_t
            j, k, diff = 1, i, 0
            while diff < 2 && j <= len_s && k <= len_t
                res += diff
                diff += (s[j] != t[k])
                j += 1; k += 1;
            end
            res += (diff == 1)
        end
        res
    end
    res = 0
    for i in 1:length(s)
        res += count_substr(SubString(s, i), t)
    end
    return res
end
## @lc code=end


text:
include("../base.jl")

f(x, y) = (y / 5.5) + x^2

println("hola")
println("hola")

x = collect(-5:0.5:5)
y = collect(0:0.5:10)

nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i = 1:nvalues
    objs[i] = f(x[i], y[i])
end

"""
N=["expr","op", "pre-op", "var"]
T=["x", "sin", "cos", "tan", "log", "+", "-", "/", "*", "(", ")"]

# <expr> <op> <expr>

R=["<expr> -> (<expr> <op> <expr>)
         | <pre-op> ( <expr> )
         | <var>",

   "<op> -> +
        | -
        | /
        | *",

    "<pre-op> -> sin
                | cos
                | tan
                | log",

    "<var> -> x"
]

S="expr"
"""

N = ["expr", "op", "var", "cons", "int", "digit"]
T = [
    "x",
    "y",
    "+",
    "-",
    "/",
    "*",
    "(",
    ")",
    "0",
    "1",
    "2",
    "3",
    "4",
    "5",
    "6",
    "7",
    "8",
    "9",
    ".",

]

#"\"\"\"\nInterpolate using `\\\$`\n\"\"\"\n"
# <expr> <op> <expr>

R = [
    "<expr> -> (<expr> <op> <expr>)
             | <var>
             | <cons>",
    "<op> -> +
           | -
           | /
           | *",
    "<var> -> x
            | y",
    "<cons> -> <int> | <int>.<int>",
    "<int> -> <int><digit> | <digit>",
    "<digit> -> 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9",
]

R2 = [
    "<expr> -> (<expr> <op> <expr>)
             | <var>",
    "<op> -> +
           | -
           | /
           | *",
    "<var> -> x
            | y"
]

S = "expr"

println("Hola")
println("Hola")

setRandomSeed(10)

myrandom() = rand(GenJ._experimentInfo._rng, collect(0:9))
setGEInfo(N, T, R, S, ["x", "y"], maxProductions = 50, maxDepth = 8)

setIndividualType(GEGenotype)

setStopCondition(maxIterations=50)
setEvaluator([FitnessFunction(compareFunctions, objs, x, y, weight=-1)])
setGenerator(rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
setSelector(tournamentSelector, 4)
setCrossoverOperator(proteinCross, probability = 0.9)
setNullMutation()
setMutationOperator(
    proteinMutation,
    uniformMutation,
    0,
    255,
    10,
    probability = 0.5,
)
setReplacementOperator(replaceAllPopulation)

@time runGenJ(verbose = false)

#genPopulation!(GenJ)
#evaluate!(GenJ, GenJ._population)

"""

Ind1 = fullGenerator(GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng)
Ind2 = fullGenerator(GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng)

#children = onePointCrossDepth(Ind1, Ind2, GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng)
#Ind2 = generateTree(GenJ._experimentInfo._GPExperimentInfo._grammar, Ind1._ind._representation)
#Ind2 = growGenerator(GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng, Ind1._ind._representation)
println("Arbol1->(join(getPhenotype(Ind1._tree)))")
#Ind2 = fullGeneratorProductions(GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng)
println("Arbol2->(getPhenotype(Ind2))")

#println(compareFunctions(Ind1, GenJ._experimentInfo._GPExperimentInfo, objs, x, y))
#println(compareFunctions(Ind2, GenJ._experimentInfo._GPExperimentInfo, objs, x, y))

Indcopy = copyGenotype(Ind1)

println(getPhenotype(GenJ._experimentInfo._GPExperimentInfo._grammar, Ind1))
println(getPhenotype(GenJ._experimentInfo._GPExperimentInfo._grammar, Ind2))
#println(getPhenotype(GenJ._experimentInfo._GPExperimentInfo._grammar, children[1]))
#println(getPhenotype(GenJ._experimentInfo._GPExperimentInfo._grammar, children[2]))

"""
#println("Parent1: lenInd = $(length(Ind1._ind._representation))  -> selections = $(length(Ind1._selections)) -> productions = $(Ind1._productions)")
#println("Parent2: lenInd = $(length(Ind2._ind._representation))  -> selections = $(length(Ind2._selections)) -> productions = $(Ind2._productions)")

#children = proteinCross(Ind1, Ind2, GenJ._experimentInfo._GPExperimentInfo, GenJ._experimentInfo._rng, kPointCross, 2)

#println("Child1 lenInd = $(length(children[1]._ind._representation))  -> selections = $(length(children[1]._selections)) -> productions = $(children[1]._productions)")
#println("Child2 lenInd = $(length(children[2]._ind._representation))  -> selections = $(length(children[2]._selections)) -> productions = $(children[2]._productions)")

"""
result = getPhenotype(children[1]._tree)
result2 = getPhenotype(children[2]._tree)
println(join(result))

printTree(Ind1._tree)

println()
println()
printTree(Ind2._tree)

result = getPhenotype(Ind1._tree)
result2 = getPhenotype(Ind2._tree, GenJ._experimentInfo._GPExperimentInfo)

child1 = getPhenotype(children[1]._tree, GenJ._experimentInfo._GPExperimentInfo)
child2 = getPhenotype(children[2]._tree, GenJ._experimentInfo._GPExperimentInfo)

println("Padre1: result")
println("Padre2: result2")
println("Hijo1: child1")
println("Hijo2: child2")
"""

"""utils
println(compareFunctions(Ind1, GenJ._experimentInfo._GPExperimentInfo, objs, x, y))

println(compareFunctions(Ind2, GenJ._experimentInfo._GPExperimentInfo, objs, x, y))
println(compareFunctions(children[1], GenJ._experimentInfo._GPExperimentInfo, objs, x, y))
println(compareFunctions(children[2], GenJ._experimentInfo._GPExperimentInfo, objs, x, y))
"""

#println(evaluate(Ind1, GenJ._experimentInfo._GPExperimentInfo, Dict("x" => 5)))
"""

x=5

println(auxstack)
println(selection)
println(result)

println()

println(auxstack2)
println(selection2)
println(result2)

"""

#println(Meta.parse(result))
#println(Meta.parse(result2))

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using HydrogenDimerMRCC

r0 = 0.776777 # ground state vibrational bond length
R  = 6.0      # centre-of-mass intermolecular distance
dimer1 = PolarCoords(r0, pi/2, pi/4)
dimer2 = PolarCoords(r0, pi/4, 4*pi/3)
dparam = DimerParameters(R, dimer1, dimer2)

filename = "MRCCinputfiles/example1/MINP"
make_mrcc_input_file(dparam, filename)


text:
export timesum
################################################################
"""
    timesum(snoop::Vector{Tuple{Float64, Core.MethodInstance}}, unit = :s)

Calculates the total time measured by a snoop macro. `unit` can be :s or :ms.

# Examples
```julia
using SnoopCompile
data = @snoopi begin
    using MatLang
    MatLang_rootpath = dirname(dirname(pathof("MatLang")))

    include("\$MatLang_rootpath/test/runtests.jl")
end
println(timesum(data, :ms))
```
"""
function timesum(snoop::Vector{Tuple{Float64, Core.MethodInstance}}, unit::Symbol = :s)
    if isempty(snoop)
        t = 0.0
    else
        t = sum(first, snoop)
    end
    if unit == :s
        t = round(t, digits = 6)
    elseif unit == :ms
        t = round(t * 1000, digits = 3)
    else
        throw(ArgumentError("unit should be `:s` or `:ms`"))
    end
    return t
end
################################################################
"""
    os_string, os_func = detectOS()

Returns Operating System of a machine as a string and a function `os_func()` that will return
`true` on the current OS.

# Examples
```julia
julia> detectOS()
("windows", Base.Sys.iswindows)
```
"""
function detectOS()
allos_funs = [
         Sys.iswindows,
         Sys.isapple,
         Sys.islinux,
         Sys.isbsd]
    if VERSION >= v"1.1"
        allos_funs = [allos_funs...,
                  Sys.isdragonfly,
                  Sys.isfreebsd,
                  Sys.isnetbsd,
                  Sys.isopenbsd]
    end
    if VERSION >= v"1.2"
        push!(allos_funs, Sys.isjsvm)
    end

    os = ""
    osfun = allos_funs[1] # temp
    for osfun in allos_funs
        if osfun()
            os = split(string(osfun), '.')[end][3:end]
            return os, osfun
        end
    end
    @error "os is not detected"
end
################################################################
"""
    standardize_osname(input::String)
    standardize_osname(inputs::Vector{String})

Standardize different names from Github actions, Travis, etc

https://help.github.com/en/actions/reference/virtual-environments-for-github-hosted-runners#supported-runners-and-hardware-resources

# Examples
```jldoctest; setup=:(using SnoopCompile)
julia> SnoopCompile.standardize_osname("ubuntu-latest")
"linux"

julia> SnoopCompile.standardize_osname(["ubuntu-latest", "macos-latest"])
2-element Array{String,1}:
 "linux"
 "apple"
```
"""
function standardize_osname(input::String)
    for (key, standard) in OS_MAP
        if occursin(key, input)
            input = standard
        end
    end
    return input
end
standardize_osname(inputs::Vector{String}) = standardize_osname.(inputs)
standardize_osname(input::Nothing) = input

const OS_MAP = Dict(
    "macos" => "apple",
    "mac" => "apple",
    "osx" => "apple",
    "apple" => "apple",
    "ubuntu" => "linux",
    "linux" => "linux",
    "win" => "windows",
    "windows" => "windows",
)
################################################################

using FilePathsBase

export GoodPath
# to convert WindowsPath before they get inside the SnoopCompile code
# GoodPath(inp::String) = inp |> Path |> _GoodPath |> string
# _GoodPath(path::WindowsPath) = PosixPath((path.drive, path.segments...))
# _GoodPath(path) = path
GoodPath(x::String) = replace(x, "\\" => "/") # doesn't remove / from the end of the strings
goodjoinpath(args...) = GoodPath(joinpath(args...))

################################################################
"""
     searchdir(rootpath::String, pattern::AbstractString)
a function to search a directory.
```julia
julia> searchdir("src", "/bot.jl")
1-element Array{String,1}:
 "src/bot.jl"
```
"""
function searchdir(rootpath::String, pattern::AbstractString)
    found_files = String[]
    for file in walkpath(Path(rootpath))
        if occursin(pattern, GoodPath(string(file)))
            push!(found_files, GoodPath(string(file)))
        end
    end
    return found_files
end

"""
Searches for the file rather a path.
"""
searchdirfile(rootpath::String, file::AbstractString) = searchdir(rootpath, string(Path(file).segments[end]))

"""
    searchdirs(rootpaths::Array{String}, pattern::AbstractString)

a function to search multiple directories.

```julia
julia> searchdirs(["./src", "./test"], "/bot.jl")
2-element Array{String,1}:
 "./src/bot.jl"
 "./test/bot/bot.jl"
```
"""
function searchdirs(rootpaths::Array{String}, pattern::AbstractString)
    unique(
        reduce(vcat,
            searchdir.(rootpaths, pattern)
        )
    )
end

"""
Searches for the file rather a path.
"""
searchdirsfile(rootpaths::Array{String}, file::AbstractString) = searchdirs(rootpaths, basename(file))

"""
search for the script! - needed because of confusing paths when referencing pattern_or_file in CI
"""
function searchdirsboth(rootpaths::Array{String}, pattern_or_file::AbstractString)
    if isfile(pattern_or_file)
        return pattern_or_file
    elseif isfile("../../$pattern_or_file") # if in deps/SnoopCompile
        return "../../$pattern_or_file"
    else
        found_files = searchdirs(rootpaths, pattern_or_file)
        if length(found_files) === 0
            found_files = searchdirsfile(rootpaths, pattern_or_file)
            if length(found_files) === 0
                @error "Couldn't find $(pattern_or_file)"
            elseif length(found_files) > 1
                @error "Multiple $(pattern_or_file) exists at the current directory."
            else
                pattern_or_file = found_files[1]
            end
        elseif length(found_files) > 1
            @error "Multiple $(pattern_or_file) exists at the current directory."
        else
            pattern_or_file = found_files[1]
        end
    end
    return pattern_or_file
end
################################################################

export pathof_noload

"""
Returns a package's path without loading the package in the main Julia process.
May launch a separate Julia process to find the package.

# Examples
```julia
pathof_noload("MatLang")
```
"""
function pathof_noload(package_name::String)
    path = Base.find_package(package_name)
    if isnothing(path)
        cmd = "import $package_name; print(pathof($package_name))"
        try
            path = Base.read(`$(Base.julia_cmd()) -e $cmd`, String)
            return GoodPath(path)
        catch
            try
                path = Base.read(`$(Base.julia_cmd()) --project=@. -e $cmd`, String)
                return GoodPath(path)
            catch
                @error "Couldn't find the path of $package_name"
            end
        end
    else
        return GoodPath(path)
    end
end

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

import Pkg
"""
Should be removed once Pkg allows adding test dependencies to the current environment

Used in Github Action workflow yaml file
"""
function addtestdep()

    if isfile("test/Test.toml")
        toml = Pkg.TOML.parsefile("test/Test.toml")
        test_deps = get(toml, "deps", nothing)
    elseif isfile("test/Project.toml")
        toml = Pkg.TOML.parsefile("test/Project.toml")
        test_deps = get(toml, "deps", nothing)
    else
        toml = Pkg.TOML.parsefile("Project.toml")
        test_deps = get(toml, "extras", nothing)
    end

    if !isnothing(test_deps)
        for (name, uuid) in test_deps
            Pkg.add(Pkg.PackageSpec(name = name, uuid = uuid))
        end
    end
end

################################################################
"""
Get the float version from semver version
```
julia> VersionFloat(v"1.5.0")
"1.4"
```
"""
VersionFloat(v::VersionNumber) = join(split(string(v),'.')[1:2],'.')

# https://github.com/JuliaLang/julia/pull/36223:
"""
    JuliaVersionNumber(v::String)
    JuliaVersionNumber(v::VersionNumber)
returns the Julia version number by following the same specifications as VersionNumber.
`JuliaVersionNumber` will fetch the latest nightly version number if `"nightly"` or `"latest"` is given as the input.
# Examples
```julia
julia> SnoopCompile.JuliaVersionNumber("nightly")
v"1.6.0-DEV"
```
```jldoctest; setup=:(using SnoopCompile)
julia> SnoopCompile.JuliaVersionNumber("1.2.3")
v"1.2.3"
julia> SnoopCompile.JuliaVersionNumber(v"1.2.3")
v"1.2.3"
```
"""
JuliaVersionNumber(v::VersionNumber) = v
function JuliaVersionNumber(v::String)
    if in(v, ["nightly", "latest"])
        version_path = download(
            "https://raw.githubusercontent.com/JuliaLang/julia/master/VERSION",
            joinpath(tempdir(), "VERSION.txt"),
        )
        version_str = replace(Base.read(version_path, String), "\n" => "")
        return VersionNumber(version_str)
    else
        return VersionNumber(v)
    end
end

################################################################
# TODO we need to remove begin and end manually
# https://github.com/timholy/SnoopCompile.jl/pull/98#issuecomment-645969610
toplevel_string(expr::Expr) = replace(string(expr), r"^begin([\s\S]*)end$"=>s"\1")

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

"""
Adds a package if it is not defined in the given module.
"""
function addpkg_ifnotfound(package_name::Symbol, env_modul::Module)
    package_name_str = string(package_name)
    if !isdefined(env_modul, package_name) && Core.eval(env_modul, :( Base.find_package($package_name_str) === nothing ))
        @info "Adding $package_name_str to $env_modul"
        Core.eval(env_modul, quote
            using Pkg
            Pkg.add($package_name_str)
        end)
    end
end

"""
Develop a package if it is not defined in the given module.
"""
function devpkg_ifnotfound(package_name::Symbol, package_path::String, env_modul::Module)
    package_name_str = string(package_name)
    if !isdefined(env_modul, package_name) && Core.eval(env_modul, :( Base.find_package($package_name_str) === nothing ))
        if dirname(Base.current_project()) !== package_path
            @info "Developing $package_name_str to $env_modul"
            Core.eval(env_modul, quote
                using Pkg
                Pkg.develop(PackageSpec(path=$package_path))
            end)
        end
    end
end


text:
function CssiSpaceWeatherDataLoader(arg0::TimeScale)
    return CssiSpaceWeatherDataLoader((TimeScale,), arg0)
end

function get_data_set(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "getDataSet", SortedSet, ())
end

function get_last_daily_predicted_date(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "getLastDailyPredictedDate", AbsoluteDate, ())
end

function get_last_observed_date(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "getLastObservedDate", AbsoluteDate, ())
end

function get_max_date(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function load_data(obj::CssiSpaceWeatherDataLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::CssiSpaceWeatherDataLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
using GLMakie
using Makie.GeometryBasics

@views function PlotMakie(mesh, v, xmin, xmax, ymin, ymax, cmap = :viridis, min_v = minimum(v), max_v = maximum(v))
    ``` Create patch plots for either Quadrangles or Triangles
    ```

    # min_v, max_v = minimum(v), maximum(v)

    f = Figure()

    Axis(f[1, 1])

    p = [Polygon( Point2f0[ (mesh.xv[mesh.e2v[i,j]], mesh.yv[mesh.e2v[i,j]]) for j=1:mesh.nf_el] ) for i in 1:mesh.nel]

    poly!(p, color = v, colormap = cmap, strokewidth = 0, strokecolor = :black, markerstrokewidth = 0, markerstrokecolor = (0, 0, 0, 0), aspect_ratio = :equal, colorrange=(min_v,max_v))

    Colorbar(f[1, 2], colormap = cmap, limits=[min_v max_v], flipaxis = true, size = 25, height = Relative(2/3) )

    xlims!(xmin, xmax)
    ylims!(ymin, ymax)

    display(f)
    return
end


text:
using Documenter, PowerModelsReliability

makedocs(
    modules = [PowerModelsReliability],
    format = :html,
    sitename = "PowerModelsReliability",
    authors = "Frederik Geth, Hakan Ergun",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting Started" => "quickguide.md",
            "Results" => "result-data.md",
        ],
        "Library" => [
            "Network Formulations" => "formulations.md",
            "Problem Specifications" => "specifications.md",
            "Modeling Components" => [
                "Objective" => "objective.md",
                "Variables" => "variables.md",
                "Constraints" => "constraints.md"
            ],
            "File IO" => "parser.md"
        ],
    ]
)
deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/frederikgeth/PowerModelsReliability.jl.git",
    julia = "0.6"
)

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
"""
    dnx(n::T, f::Basic, v::Basic) where {T <: Unsigned}

Computes the n-th derivative of the the function f with respect to the variable v.
This function re implements SymEngine's own `diff` function using an early quitting strategy.

# Input
- `n::T`: The order of the derivative.
- `f::Basic`: The expression to compute the derivative.
- `v::Basic`: The variable with respect to which the derivative will be computed.

# Output
A `SymEngine.Basic` object with the derived expression.
"""
function dnx(n::T, f::Basic, v::Basic) where {T <: Unsigned}
    b0 = Basic(0)
    zt = zero(T)
    ot = one(T)

    f == b0 && return b0
    n == zt && return f
    n == ot && return SymEngine.diff(f, v)
    n > 1 && return dnx(n - ot, SymEngine.diff(f, v), v)
end

"""
    dnx(p::AnalyticAIMProblem{N,T}, n::N, f::Function, v::Function) where {N <: Unsigned, T <: Number}

Computes the n-th derivative of the AIM expressions with respect to ODE's variable.
This function is only a thin wrapper around SymEngine's own `diff` function.
It works as a barrier function that produces a type stable `Basic` result.

# Input
- `p::AnalyticAIMProblem`: The problem data with the expressions to compute the derivative.
- `n::Unsigned`: The order of the derivative.
- `f::Function`: The actual expression to compute the derivative. Either λ0 or S0.
- `v::Function`: The variable with respec to which the derivative will be computed. Either get_ODEvar or get_ODEeigen.

# Output
A `SymEngine.Basic` object with the derived expression.
"""
function dnx(p::AnalyticAIMProblem{N,T}, n::N, f::Function, v::Function) where {N <: Unsigned, T <: Number}
    return convert(SymEngine.Basic, dnx(n, f(p), v(p)))
end

"""
    computePolynomialFactors(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function) where {N <: Unsigned, T <: Number}

Create a second order `Polynomial` object in the variable `ω` by computing derivatives of `λ0` or `S0`.

# Input
- `p::QuadraticEigenvalueProblem`: The problem data with the expressions to compute the derivative.
- `n::Unsigned`: The order of the derivative.
- `f::Function`: the function to extract the polynomial from. Either λ0 or S0.

# Output
An object of type `Polynomial{T}` containing the polynomial resulting from the derivation of the expression.
"""
function computePolynomialFactors(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function) where {N <: Unsigned, T <: Number}
    SYB = SymEngine.Basic

    # Compute the derivative at x0 and expand it to obtain a polynomial in ω
    step1::SYB = convert(SYB, subs(dnx(p, n, f, get_ODEvar), get_ODEvar(p) => get_x0(p)))
    step2::SYB = convert(SYB, expand(step1))

    # Extract the factors of the polynomial in ω
    ω = get_ODEeigen(p)

    p0::T = convert(T, SymEngine.N(coeff(step2, ω, Basic(0))))
    p1::T = convert(T, SymEngine.N(coeff(step2, ω, Basic(1))))
    p2::T = convert(T, SymEngine.N(coeff(step2, ω, Basic(2))))

    poly = Polynomial{T}([p0, p1, p2])

    return poly
end

"""
    createPoly(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function)

Compute the n-th coefficient of the Taylor expansion around x0 for the functions `λ0` or `S0`

# Input
- `p::QuadraticEigenvalueProblem`: The problem data with the expressions to compute the derivative.
- `n::Unsigned`: The order of the derivative.
- `f::Function`: the function to extract the polynomial from. Either λ0 or S0.

# Output
An object of type `Polynomial{T}` containing the polynomial Taylor coefficient.
"""
function createPoly(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function) where {N <: Unsigned, T <: Number}
    computePolynomialFactors(p, n, f)/factorial(n)
end

function createPoly(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function) where {N <: Unsigned, T <: BigFloat}
    computePolynomialFactors(p, n, f)/factorial(big(n))
end

function createPoly(p::QuadraticEigenvalueProblem{N,T}, n::N, f::Function) where {N <: Unsigned, T <: Complex{BigFloat}}
    computePolynomialFactors(p, n, f)/factorial(big(n))
end


text:
using OptiMo, Bazinga
using NLPModelsIpopt, CUTEst
using DataFrames, Query
using CSV, Printf

data = DataFrame();

problem = "HIMMELBI"
nlp = CUTEstModel( problem )
prob = NLPOptiModel( nlp )

TOL_OPTIM = 1e-6
TOL_CVIOL = TOL_OPTIM

# warm-up
out = ipopt( nlp, tol=sqrt(TOL_OPTIM) )
solver = Bazinga.ALPX(  verbose=true,
                        tol_optim=sqrt(TOL_OPTIM),
                        tol_cviol=sqrt(TOL_CVIOL),
                        subsolver=:zerofpr,
                        max_iter=50,
                        max_sub_iter=10 )
out = solver( prob )

# solve!
for maxsubiter in [1000; 2500; 5000; 10000; 12500; 15000; 20000]

    @printf "\n\n max sub iters %d \n" maxsubiter
    solver = Bazinga.ALPX(  verbose=true,
                            tol_optim=TOL_OPTIM,
                            tol_cviol=TOL_CVIOL,
                            subsolver=:zerofpr,
                            max_iter=50,
                            max_sub_iter=maxsubiter )
    out = solver( prob )
    @printf "sub iters %d \n" out.solver[:sub_iterations]
    @printf "time      %f \n" out.time

end

@printf "\n\n IPOPT \n"
out = ipopt( nlp, tol=TOL_OPTIM )
print(out)
@printf "time      %f \n\n" out.elapsed_time

@printf "\n\n"
maxsubiter = 5000
mu_est=1.0
solver = Bazinga.ALPX(  verbose=true,
                        tol_optim=TOL_OPTIM,
                        tol_cviol=TOL_CVIOL,
                        mu_est=mu_est,
                        max_sub_iter=maxsubiter )
out = solver( prob )
@printf "sub iters %d \n" out.solver[:sub_iterations]
@printf "time      %f \n" out.time

mu_est=10.0
solver = Bazinga.ALPX(  verbose=true,
                        tol_optim=TOL_OPTIM,
                        tol_cviol=TOL_CVIOL,
                        mu_est=mu_est,
                        max_sub_iter=maxsubiter )
out = solver( prob )
@printf "sub iters %d \n" out.solver[:sub_iterations]
@printf "time      %f \n" out.time

# finalize
finalize( nlp )


text:
export PropagationModel, source!, receiver!, record

abstract type PropagationModel end

Base.@kwdef struct AcousticSource
  pos = []
  signal = nothing
end

Base.@kwdef struct AcousticReceiver
  pos = []
end

source!(model::PropagationModel, pos, signal=nothing) = append!(model.sources, AcousticSource(pos, signal))
receiver!(model::PropagationModel, pos) = append!(model.receivers, AcousticReceiver(pos))

function record(model::PropagationModel, duration)
  nsamples = round(Int, duration*model.fs)
  x = zeros(nsamples, length(model.receivers))
  for j = 1:length(model.sources)
    sig = generate(model.sources[j].signal, nsamples, model.fs)
    for k = 1:length(model.receivers)
      x[:,k] .+= transmit(model, sig, model.sources[j].pos, model.receivers[k].pos)
    end
  end
  for k = 1:length(model.receivers)
    x[:,k] .+= rand(model.noise, nsamples) * √model.fs
  end
  return x
end


text:
using FileIO

function evalObjGrad( pcof0:: Array{Float64, 1}, params:: Juqbox.objparams, wa:: Juqbox.Working_Arrays, refFileName:: String, writeFile:: Bool=false)
    rtol = 1e-10
    atol = 1e-14
    verbose = false

    # return flag
    success=false

    # evaluate fidelity
    objv, grad = Juqbox.traceobjgrad(pcof0, params, wa, verbose, true);

    if writeFile
        save(refFileName, "obj0", objv, "grad0", grad)
        success=true
    else
        # read reference solution
        # pcofRef = vec(readdlm(refFileName)) # change to JLD2 file to get full accuracy
        # objvRef = pcofRef[1]
        # gradRef = pcofRef[2:end]
        dict = load(refFileName)
        objvRef = dict["obj0"]
        gradRef = dict["grad0"]

        objDiff = abs(objv - objvRef)

        refNorm = norm(gradRef)
        aNorm = norm(grad-gradRef)

        pass1 = false
        if objDiff < atol
            pass1 = true
        elseif abs(objvRef) >= atol && objDiff/abs(objvRef) < rtol
            pass1 = true
        end

        pass2 = false
        if aNorm < atol
            pass2 = true
        elseif refNorm >= atol &&  aNorm/refNorm < rtol
            pass2 = true
        end

        success = pass1 && pass2

        if !success
            println("rel objfunc error =", objDiff, " absolute grad error =", aNorm, " relative grad error =", aNorm/refNorm)
        end

    end

    return success
end


text:
# ------------ Dialect alias and APIs ------------ #

const Dialect = MLIR.API.MlirDialect

get_context(d::Dialect) = MLIR.API.mlirDialectGetContext(d)
is_null(d::Dialect) = MLIR.API.mlirDialectIsNull(d)
Base.:(==)(d1::Dialect, d2::Dialect) = MLIR.API.mlirDialectEqual(d1, d2)
get_namespace(d::Dialect) = MLIR.API.mlirDialectGetNamespace(d)

@doc(
"""
const Dialect = MLIR.API.MlirDialect
""", Dialect)


text:
@require "ms-parse"

export pretty

#
# Convert `ms` into a pretty printed time string.
#
# Examples:
#   343564753 => 52d 20h 42m 24s
#   302 => 302ms
#   12379 => 12.4s
#

function pretty (ms)
  ret = String[]
  time = parse(ms)

  if ms < 1000
    return post(time["milliseconds"], "ms")
  end

  add!(ret, time["days"], "d")
  add!(ret, time["hours"], "h")
  add!(ret, time["minutes"], "m")

  # Display `ms` to `s` when `s` is the highest
  # unit of time

  if time["hours"] == 0
    rms = int(round(time["milliseconds"] / 1000, 1) * 10)
    s = string(time["seconds"]) * "." * string(rms)
    add!(ret, s, "s")
  else
    add!(ret, time["seconds"], "s")
  end

  final = join(ret, " ")
  return final
end

#
# push Any `val` to the given Array `ret` only if
# `val` is not 0
#

function add! (ret::Array, val, postfix::String)
  if val != 0
    push!(ret, post(val, postfix))
  end
end

#
# Add `value` together wth `postfix` as a string.
#

post (val, postfix) = string(val) * string(postfix)

max_stars_repo_name: queckezz/ms-pretty

max_issues_repo_name: queckezz/ms-pretty

max_forks_repo_name: queckezz/ms-pretty


text:
# Imaginary Time Evolution
# I don't think this works!
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Imaginary%20Time%20Evolution][Imaginary Time Evolution:1]]
# Fixme! this does not appear to find ground states!

function _MPO_handed_time_evolver(hs::Vector{Matrix{T}}, τ, L, d) where {T}
    tensors = Array{T, 4}[]
    for h in hs
        O = exp(-τ*h)
        @cast P[(σⁱ, σⁱ′), (σⁱ⁺¹, σⁱ⁺¹′)] |= O[(σⁱ, σⁱ⁺¹), (σⁱ′, σⁱ⁺¹′)] (σⁱ:d, σⁱ′:d)
        U, S, V = svd(P)

        @cast U[1, k, σⁱ, σⁱ′]     := U[(σⁱ, σⁱ′), k] * √(S[k])      (σⁱ:d)
        @cast Ū[k, 1, σⁱ⁺¹, σⁱ⁺¹′] := √(S[k]) * V'[k, (σⁱ⁺¹, σⁱ⁺¹′)] (σⁱ⁺¹:d)
        push!(tensors, U, Ū)
    end
    MPO{L, T}(tensors)
end

function MPO_time_evolvers(h1::Matrix, hi::Matrix, hL::Matrix, τ, L, d)
    if iseven(L)
        odd_hs  = [h1, [hi for _ in 3:2:(L-1)]...]
        even_hs = [[hi for i in 2:2:(L-1)]..., hL]
    else
        odd_hs  = [h1, [hi for _ in 3:2:(L-1)]..., hL]
        even_hs = [hi for i in 2:2:(L-1)]
    end

    Uodd  = _MPO_handed_time_evolver(odd_hs, τ, L, d)
    Ueven = _MPO_handed_time_evolver(even_hs, τ, L, d)
    Uodd, Ueven
end

function imag_time_evolution(ψ::MPS{L, T}, h1::Matrix{T}, hi::Matrix{T}, hL::Matrix{T},
                             β, N, Dcut) where {L, T}
    @warn "This probably still doesn't work!"
    τ = β/N
    d = length(ψ[1][1, 1, :])
    ϕ = ψ  # Ground state guess
    dir = left
    Uodd, Ueven = MPO_time_evolvers(h1, hi, hL, τ, L, d)
    for _ in 1:N
        ϕ1, dir = compress(Uodd  * ϕ,  dir, Dcut=Dcut)
        ϕ,  dir = compress(Ueven * ϕ1, dir, Dcut=Dcut)
        #ϕ,  dir = compress(Uodd  * ϕ2, dir, Dcut=Dcut)
    end
    ϕ
end
# Imaginary Time Evolution:1 ends here


text:
#!/usr/bin/env julia

using Luxor, Colors

using Test

using Random
Random.seed!(42)

function drawlogoclip(x, y)
    foregroundcolors = diverging_palette(rand(0:360), rand(0:360), 200, s = 0.99, b=0.8)
    gsave()
    translate(x-100, y)
    julialogo(action=:clip)
    for i in 1:500
        sethue(foregroundcolors[rand(1:end)])
        circle(rand(-50:350), rand(0:300), 15, :fill)
    end
    #  to test clippreserve(), draw a path, and it is intersected with clip region (julia logo)
    # then filled with solid red
    gsave()
        circle(O, 200, :path)
        clippreserve()
        setcolor("red")
        fillpath()
    grestore()
    clipreset()
    grestore()
end

function clipping_test(fname)
    Drawing(500, 500, fname)
    origin()
    setopacity(.4)
    foregroundcolors = diverging_palette(rand(0:360), rand(0:360), 200, s = 0.99, b=0.8)
    drawlogoclip(0, 0)
    @test finish() == true
    println("...finished test: output in $(fname)")
end

clipping_test("clipping-tests.png")


text:
module Testbyrow

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "AsTable on RHS" begin
    df = DataFrame(a = [1, 2], b = [3, 4], d = [100, 200])

    res = DataFrame(a = [1, 2], b = [3, 4], d = [100, 200], c = [4, 6])

    d = @transform df :c = sum(AsTable([:a, :b]))
    @test d == res

    d = @rtransform df :c = sum(AsTable([:a, :b]))
    @test d == res

    d = @transform df @byrow :c = sum(AsTable([:a, :b]))
    @test d == res

    vars = ["a", "b"]

    d = @rtransform df :c = sum(AsTable(vars))
    @test d == res

    d = @rtransform df :c = sum(AsTable(Between(:a, :b)))
    @test d == res

    d = @transform df :c = sum(AsTable(Not(:d)))
    @test d == res
end

end # module


text:
module GLMakie

using ModernGL, FixedPointNumbers, Colors, GeometryBasics, StaticArrays
using AbstractPlotting, FileIO

using AbstractPlotting: @key_str, Key, broadcast_foreach, to_ndim, NativeFont
using AbstractPlotting: Scene, Lines, Text, Image, Heatmap, Scatter
using AbstractPlotting: convert_attribute, @extractvalue, LineSegments
using AbstractPlotting: @get_attribute, to_value, to_colormap, extrema_nan
using AbstractPlotting: ClosedInterval, (..)
using ShaderAbstractions
using FreeTypeAbstraction

using Base: RefValue
import Base: push!, isopen, show
using Base.Iterators: repeated, drop

using LinearAlgebra

for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
    @eval export $(name)
end

struct GLBackend <: AbstractPlotting.AbstractBackend
end

"""
returns path relative to the assets folder
"""
assetpath(folders...) = joinpath(@__DIR__, "GLVisualize", "assets", folders...)

"""
Loads a file from the asset folder
"""
function loadasset(folders...; kw_args...)
    path = assetpath(folders...)
    isfile(path) || isdir(path) || error("Could not locate file at $path")
    load(path; kw_args...)
end

const deps_path = joinpath(@__DIR__, "..", "deps", "deps.jl")

if isfile(deps_path)
    include(deps_path)
else
    error("""
        The file $(deps_path) does not exist.
        This file is generated during the build process; it is possible that GLMakie
        wasn't built correctly.  To rerun the build process, run `Pkg.build("GLMakie"),
        or enter the Pkg REPL mode (`]`) and then type `build GLMakie`.
        """)
end

if WORKING_OPENGL
     # don't put this into try catch, to not mess with normal errors
    include("gl_backend.jl")
end

function activate!(use_display=true)
    b = GLBackend()
    AbstractPlotting.register_backend!(b)
    AbstractPlotting.set_glyph_resolution!(AbstractPlotting.High)
    AbstractPlotting.current_backend[] = b
    AbstractPlotting.inline!(!use_display)
end

function __init__()
    if WORKING_OPENGL
        activate!()
    else
        error("Loaded OpenGL Backend, but OpenGL isn't working")
    end
end

export set_window_config!
export assetpath, loadasset

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end


text:
cd(dirname(@__FILE__))
include("./functions/load_inputfunctions.jl")
J,δ =1000,.01
cs,al,q,w =  get_params()
include("./functions/load_stuff.jl")

B,S = generate_structures(cs,al,q,w,J,δ);
al == 0 ? solve_deterministic!(B,cs) : solve_Benders!(B,S)
println("")


text:
function get_pv_coordinates(obj::FieldPVCoordinatesProvider, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end


text:
# state transitions

transition(t::TrebuchetState, s::Array) = transition(t::TrebuchetState, s::Array, t.stage)

function transition(t::TrebuchetState, s::Array, ::Val{:Ground})
    t.aw = AnglularVelocities(s[4:6]...)
    t.a = Angles(s[1:3]...)
    if isa(t.terminated, Val{:Ground})
        t.stage = Val{:Hang}()
    elseif t.Tn < zero(typeof(t.Tn))
        @warn "Weight is too light ($(t.m.w))"
        t.stage = Val{:End}()
    end
end

function transition(t::TrebuchetState, s::Array, ::Val{:Hang})
    t.aw = AnglularVelocities(s[4:6]...)
    t.a = Angles(s[1:3]...)
    t.p = sling_end(t)
    t.v = sling_velocity(t)
    t.sol.ReleaseVelocity = t.p
    t.sol.ReleasePositon = t.v
    if isa(t.terminated, Val{:Hang})
        t.stage = Val{:Released}()
    end
end

function transition(t::TrebuchetState, s::Array, ::Val{:Released})
    t.p = Vec(s[1], s[2])
    t.v = Vec(s[3], s[4])
    if isa(t.terminated, Val{:Released})
        t.stage = Val{:End}()
    end
end


text:
#
## Boolean methods
#

function excl(b::Sidef_Types_Bool_Bool)
    b.value ? FALSE : TRUE
end

function !(b::Sidef_Types_Bool_Bool)
    b.value ? FALSE : TRUE
end

function string(b::Sidef_Types_Bool_Bool)
    b.value ? "true" : "false"
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


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

"""
    CartesianGrid(dims, origin, spacing)

A Cartesian grid with dimensions `dims`, lower left corner at `origin`
and cell spacing `spacing`. The three arguments must have the same length.

    CartesianGrid(start, finish, dims=dims)

Alternatively, construct a Cartesian grid from a `start` point (lower left)
to a `finish` point (upper right).

    CartesianGrid(start, finish, spacing)

Alternatively, construct a Cartesian grid from a `start` point to a `finish`
point using a given `spacing`.

    CartesianGrid(dims)
    CartesianGrid(dim1, dim2, ...)

Finally, a Cartesian grid can be constructed by only passing the dimensions
`dims` as a tuple, or by passing each dimension `dim1`, `dim2`, ... separately.
In this case, the origin and spacing default to (0,0,...) and (1,1,...).

## Examples

Create a 3D grid with 100x100x50 locations:

```julia
julia> CartesianGrid(100,100,50)
```

Create a 2D grid with 100x100 locations and origin at (10.,20.) units:

```julia
julia> CartesianGrid((100,100),(10.,20.),(1.,1.))
```

Create a 1D grid from -1 to 1 with 100 locations:

```julia
julia> CartesianGrid((-1.,),(1.,), dims=(100,))
```
"""
struct CartesianGrid{Dim,T} <: Mesh{Dim,T}
  dims::Dims{Dim}
  origin::Point{Dim,T}
  spacing::SVector{Dim,T}

  function CartesianGrid{Dim,T}(dims, origin, spacing) where {Dim,T}
    @assert all(dims .> 0) "dimensions must be positive"
    @assert all(spacing .> 0) "spacing must be positive"
    new(dims, origin, spacing)
  end
end

CartesianGrid(dims::Dims{Dim}, origin::Point{Dim,T},
              spacing::SVector{Dim,T}) where {Dim,T} =
  CartesianGrid{Dim,T}(dims, origin, spacing)

CartesianGrid(dims::Dims{Dim}, origin::NTuple{Dim,T},
              spacing::NTuple{Dim,T}) where {Dim,T} =
  CartesianGrid{Dim,T}(dims, Point(origin), SVector(spacing))

function CartesianGrid(start::Point{Dim,T}, finish::Point{Dim,T},
                       spacing::SVector{Dim,T}) where {Dim,T}
  dims = Tuple(ceil.(Int, (finish - start) ./ spacing))
  CartesianGrid{Dim,T}(dims, start, spacing)
end

CartesianGrid(start::NTuple{Dim,T}, finish::NTuple{Dim,T},
              spacing::NTuple{Dim,T}) where {Dim,T} =
  CartesianGrid(Point(start), Point(finish), SVector(spacing))

CartesianGrid(start::Point{Dim,T}, finish::Point{Dim,T};
              dims::Dims{Dim}=ntuple(i->100, Dim)) where {Dim,T} =
  CartesianGrid{Dim,T}(dims, start, (finish - start) ./ dims)

CartesianGrid(start::NTuple{Dim,T}, finish::NTuple{Dim,T};
              dims::Dims{Dim}=ntuple(i->100, Dim)) where {Dim,T} =
  CartesianGrid(Point(start), Point(finish); dims=dims)

CartesianGrid{T}(dims::Dims{Dim}) where {Dim,T} =
  CartesianGrid{Dim,T}(dims, ntuple(i->zero(T), Dim), ntuple(i->one(T), Dim))

CartesianGrid{T}(dims::Vararg{Int,Dim}) where {Dim,T} = CartesianGrid{T}(dims)

CartesianGrid(dims::Dims{Dim}) where {Dim} = CartesianGrid{Float64}(dims)

CartesianGrid(dims::Vararg{Int,Dim}) where {Dim} = CartesianGrid{Float64}(dims)

==(g1::CartesianGrid, g2::CartesianGrid) =
  g1.dims    == g2.dims    &&
  g1.origin  == g2.origin  &&
  g1.spacing == g2.spacing

Base.size(g::CartesianGrid) = g.dims
Base.minimum(g::CartesianGrid) = g.origin
Base.maximum(g::CartesianGrid) = g.origin + g.dims .* g.spacing
Base.extrema(g::CartesianGrid) = minimum(g), maximum(g)
spacing(g::CartesianGrid) = g.spacing

function vertices(g::CartesianGrid)
  inds = CartesianIndices(g.dims .+ 1)
  ivec(g.origin + (ind.I .- 1) .* g.spacing for ind in inds)
end

elements(g::CartesianGrid) = (g[i] for i in 1:nelements(g))

topology(g::CartesianGrid) = GridTopology(size(g))

# -----------------
# DOMAIN INTERFACE
# -----------------

function element(g::CartesianGrid{Dim}, ind::Int) where {Dim}
  I = CartesianIndices(g.dims)[ind]
  o = coordinates(g.origin)
  s = g.spacing
  i = I.I

  if Dim == 1 # segment
    p1 = (o[1] + (i[1] - 1) * s[1],)
    p2 = (o[1] + (    i[1]) * s[1],)
    Segment(p1, p2)
  elseif Dim == 2 # quadrangle
    p1 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2] - 1) * s[2])
    p2 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2] - 1) * s[2])
    p3 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2]    ) * s[2])
    p4 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2]    ) * s[2])
    Quadrangle(p1, p2, p3, p4)
  elseif Dim == 3 # hexahedron
    p1 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2] - 1) * s[2],
          o[3] + (i[3] - 1) * s[3])
    p2 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2] - 1) * s[2],
          o[3] + (i[3] - 1) * s[3])
    p3 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2]    ) * s[2],
          o[3] + (i[3] - 1) * s[3])
    p4 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2]    ) * s[2],
          o[3] + (i[3] - 1) * s[3])
    p5 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2] - 1) * s[2],
          o[3] + (i[3]    ) * s[3])
    p6 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2] - 1) * s[2],
          o[3] + (i[3]    ) * s[3])
    p7 = (o[1] + (i[1]    ) * s[1],
          o[2] + (i[2]    ) * s[2],
          o[3] + (i[3]    ) * s[3])
    p8 = (o[1] + (i[1] - 1) * s[1],
          o[2] + (i[2]    ) * s[2],
          o[3] + (i[3]    ) * s[3])
    Hexahedron(p1, p2, p3, p4, p5, p6, p7, p8)
  else
    throw(ErrorException("not implemented"))
  end
end

nelements(g::CartesianGrid) = prod(g.dims)

function centroid(g::CartesianGrid{Dim}, ind::Int) where {Dim}
  intcoords = CartesianIndices(g.dims)[ind]
  neworigin = coordinates(g.origin) .+ g.spacing ./ 2
  Point(ntuple(i -> neworigin[i] + (intcoords[i] - 1)*g.spacing[i], Dim))
end

Base.eltype(g::CartesianGrid) = typeof(g[1])

# ----------------
# OPTIONAL TRAITS
# ----------------

isgrid(::Type{<:CartesianGrid}) = true

# ----------------------------
# ADDITIONAL INDEXING METHODS
# ----------------------------

"""
    grid[istart:iend,jstart:jend,...]

Return a subgrid of the Cartesian `grid` using integer ranges
`istart:iend`, `jstart:jend`, ...
"""
Base.getindex(g::CartesianGrid{Dim}, r::Vararg{UnitRange{Int},Dim}) where {Dim} =
  getindex(g, CartesianIndex(first.(r)):CartesianIndex(last.(r)))

function Base.getindex(g::CartesianGrid{Dim}, I::CartesianIndices{Dim}) where {Dim}
  start  = coordinates(g.origin) .+ (first(I).I .- 1) .* g.spacing
  finish = coordinates(g.origin) .+ (last(I).I      ) .* g.spacing
  dims   = size(I)
  CartesianGrid(Point(start), Point(finish), dims=dims)
end

Base.view(g::CartesianGrid{Dim}, I::CartesianIndices{Dim}) where {Dim} = getindex(g, I)

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, g::CartesianGrid{Dim,T}) where {Dim,T}
  dims = join(g.dims, "×")
  print(io, "$dims CartesianGrid{$Dim,$T}")
end

function Base.show(io::IO, ::MIME"text/plain", g::CartesianGrid)
  println(io, g)
  println(io, "  minimum: ", minimum(g))
  println(io, "  maximum: ", maximum(g))
  print(  io, "  spacing: ", Tuple(spacing(g)))
end


text:
"""
 Just counting the number of edges
 in the corresponding clique-expansion
 of each hypergraph
"""
using Pkg
Pkg.activate(".")
using LTMSim
using SimpleHypergraphs

project_path = dirname(pathof(LTMSim))
data_path = joinpath(project_path, "..", "data", "hgs")

hg_files = readdir(data_path)

hgs = [prune_hypergraph!(hg_load(joinpath(data_path, hg_file))) for hg_file in hg_files]

for (index, hg) in enumerate(hgs)
    s = 0

    for he in 1:nhe(hg)
        vs = length(getvertices(hg, he))
        s += (vs * (vs-1)) / 2
    end

    println(hg_files[index], " ", size(hg), " ", s)
end


text:
module ErlangTermFormat

export encode, decode

const VERSION_MAGIC = 0x83

include("decode.jl")
include("encode.jl")

end


text:
#
# Database queries.
#

# Binding table associates identifiers with query constructors.

abstract AbstractBinding

typealias MaybeBinding Nullable{AbstractBinding}
typealias BindingTable ImmutableDict{Tuple{Symbol, Int}, AbstractBinding}

const NO_BINDING_TABLE = BindingTable()

lookup(tbl::BindingTable, ident::Symbol, arity::Int) =
    let mb = get(tbl, (ident, arity))
        isnull(mb) ? get(tbl, (ident, -1)) : mb
    end

# Database can produce global and local binding tables.

abstract AbstractDatabase

globals(db::AbstractDatabase) =
    error("not implemented")
locals(db::AbstractDatabase, dom::Domain) =
    error("not implemented")

immutable EmptyDatabase <: AbstractDatabase
end

const NO_DB = EmptyDatabase()

globals(db::EmptyDatabase) = NO_BINDING_TABLE
locals(db::EmptyDatabase, dom::Domain) = locals(db, dom.desc)

locals(db::EmptyDatabase, desc) = NO_BINDING_TABLE

function locals(db::EmptyDatabase, desc::Vector{Output})
    tbl = NO_BINDING_TABLE
    for pos in eachindex(desc)
        fld = FieldQuery(Domain(desc), pos)
        tbl = assoc(tbl, (Symbol(pos), 0), QueryBinding(fld))
        tag = decoration(desc[pos], :tag, Symbol, Symbol(""))
        if tag != Symbol("")
            tbl = assoc(tbl, (tag, 0), QueryBinding(fld))
        end
    end
    return tbl
end

# Scope is a pair of global and local binding tables.

immutable Scope
    db::AbstractDatabase
    dom::Domain
    globs::BindingTable
    locs::BindingTable
end

Scope(db::AbstractDatabase, dom::Domain) =
    Scope(db, dom, globals(db), locals(db, dom))

Scope(db::AbstractDatabase, desc::Union{Symbol, Type}) =
    Scope(db, Domain(desc))

database(scope::Scope) = scope.db
domain(scope::Scope) = scope.dom
globals(scope::Scope) = scope.globs
locals(scope::Scope) = scope.locs

lookup(scope::Scope, ident::Symbol, arity::Int) =
    let mb = lookup(scope.locs, ident, arity)
        isnull(mb) ? lookup(scope.globs, ident, arity) : mb
    end

# Signature of the query expression.

abstract AbstractSignature
abstract AbstractPrimitive <: AbstractSignature

# Binding tables for the input and the output scopes.

immutable Namespace
    db::AbstractDatabase
    ilocs::BindingTable
    olocs::BindingTable
    globs::BindingTable
end

const NO_NAMESPACE =
    Namespace(NO_DB, NO_BINDING_TABLE, NO_BINDING_TABLE, NO_BINDING_TABLE)

# Query interface.

immutable Query
    sig::AbstractSignature
    args::Vector{Query}
    ity::Input
    oty::Output
    ns::Namespace
    src::Any
end

const NO_ARGUMENTS = Query[]

Query(sig::AbstractPrimitive) =
    Query(sig, NO_ARGUMENTS, input(sig), output(sig), NO_NAMESPACE, nothing)
Query(sig::AbstractSignature, args::Vector{Query}) =
    Query(sig, args, input(sig, args), output(sig, args), NO_NAMESPACE, nothing)
Query(sig::AbstractSignature, args::Vector{Query}, ity::Union{Domain,Input}, oty::Union{Domain,Output}) =
    Query(sig, args, convert(Input, ity), convert(Output, oty), NO_NAMESPACE, nothing)
Query(sig::AbstractPrimitive, ity::Union{Domain,Input}, oty::Union{Domain,Output}) =
    Query(sig, NO_ARGUMENTS, convert(Input, ity), convert(Output, oty))

show(io::IO, q::Query) =
    print(io, "$(q.ity) -> $(q.oty)")

setsignature(q::Query, sig::AbstractSignature) =
    Query(sig, q.args, q.ity, q.oty, q.ns, q.src)
setarguments(q::Query, args::Vector{Query}) =
    Query(q.sig, args, q.ity, q.oty, q.ns, q.src)
setinput(q::Query, ity::Input) =
    Query(q.sig, q.args, ity, q.oty, q.ns, q.src)
setoutput(q::Query, oty::Output) =
    Query(q.sig, q.args, q.ity, oty, q.ns, q.src)
setnamespace(q::Query, ns::Namespace) =
    Query(q.sig, q.args, q.ity, q.oty, ns, q.src)
setnamespace(q::Query, db::AbstractDatabase) =
    let globs = globals(db),
        ilocs = locals(db, domain(q.ity)),
        olocs = locals(db, domain(q.oty)),
        ns = Namespace(db, ilocs, olocs, globs)
        setnamespace(q, ns)
    end
setsource(q::Query, src::Any) =
    Query(q.sig, q.args, q.ity, q.oty, q.ns, src)

setsignature(sig::AbstractSignature) =
    q -> setsignature(q, sig)
setarguments(args::Vector{Query}) =
    q -> setarguments(q, args)
setinput(itype::Input) =
    q -> setinput(q, itype)
setoutput(oty::Output) =
    q -> setoutput(q, oty)
setnamespace(ns::Union{Namespace, AbstractDatabase}) =
    q -> setnamespace(q, ns)
setsource(src::Any) =
    q -> setsource(q, src)

signature(q::Query) = q.sig

input(q::Query) = q.ity
output(q::Query) = q.oty

iscope(q::Query) = Scope(q.ns.db, domain(q.ity), q.ns.globs, q.ns.ilocs)
oscope(q::Query) = Scope(q.ns.db, domain(q.oty), q.ns.globs, q.ns.olocs)

isplain(q::Query) = isplain(output(q))
isoptional(q::Query) = isoptional(output(q))
isplural(q::Query) = isplural(output(q))

input(sig::AbstractSignature, args::Vector{Query}) =
    input(sig, args...)
output(sig::AbstractSignature, args::Vector{Query}) =
    output(sig, args...)

istub(q::Query) =
    let dom = domain(q.ity),
        ns = Namespace(q.ns.db, q.ns.ilocs, q.ns.ilocs, q.ns.globs)
        Query(ItSig(), NO_ARGUMENTS, Input(dom), Output(dom), ns, nothing)
    end

ostub(q::Query) =
    let dom = domain(q.oty),
        ns = Namespace(q.ns.db, q.ns.olocs, q.ns.olocs, q.ns.globs)
        Query(ItSig(), NO_ARGUMENTS, Input(dom), Output(dom), ns, nothing)
    end

decorate(q::Query, decor::Decoration) =
    let dom = domain(q.oty),
        ns = Namespace(q.ns.db, q.ns.olocs, q.ns.olocs, q.ns.globs)
        q >> Query(ItSig(), NO_ARGUMENTS, Input(dom), Output(decorate(dom, decor)), ns, nothing)
    end

# Executing queries.

ev(q::Query, ctx::InputContext=InputContext(), dom=Domain(Void), vals::AbstractVector=[nothing]) =
    ev(q, InputFlow(ctx, dom, vals))

ev(q::Query, dom, vals::AbstractVector) =
    ev(q, InputFlow(InputContext(), dom, vals))

ev{T}(q::Query, vals::AbstractVector{T}) =
    ev(q, InputFlow(InputContext(), T, vals))

ev(q::Query, iflow::InputFlow) =
    ev(q.sig, q.args, q.ity, q.oty, iflow)

ev(sig::AbstractPrimitive, ::Vector{Query}, ity::Input, oty::Output, iflow::InputFlow) =
    ev(sig, ity, oty, iflow)

ev(sig::AbstractPrimitive, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(oty, ev(sig, values(iflow)))

const OPTIMIZE_CACHE = Dict{Query, Query}()

function optimize(q::Query)
    if q in keys(OPTIMIZE_CACHE)
        return OPTIMIZE_CACHE[q]
    end
    oq = q
    OPTIMIZE_CACHE[q] = oq
    return oq
end

format(q::Query) =
    format(q, decoration(domain(output(q)), :fmt, Union{Void, Symbol, Vector{Symbol}}, nothing))

format(q::Query, ::Void) = q

format(q::Query, name::Symbol) =
    let scope = oscope(q),
        binding = get(lookup(scope, name, 0)),
        tag = decoration(output(q), :tag, Symbol, Symbol(""))
        q = q >> format(compile(binding, scope, Syntax[]))
        if tag != Symbol("")
            q = q |> decorate(:tag => tag)
        end
        q
    end

format(q::Query, names::Vector{Symbol}) =
    let scope = oscope(q),
        bindings = [get(lookup(scope, name, 0)) for name in names],
        tag = decoration(output(q), :tag, Symbol, Symbol(""))
        q = q >> RecordQuery(
                    [format(compile(binding, scope, Syntax[])) for binding in bindings])
        if tag != Symbol("")
            q = q |> decorate(:tag => tag)
        end
        q
    end

function execute{T}(q::Query, val::T=nothing; paramvals...)
    q = optimize(q)
    itype = input(q)
    otype = output(q)
    if isrelative(itype)
        frameoffs = InputFrame(OneTo(2))
    else
        frameoffs = InputFrame()
    end
    pflows = InputSlotFlow[]
    parammap = Dict{Symbol,Any}(paramvals)
    for (pname, ptype) in slots(itype)
        paramval = get(parammap, pname, nothing)
        paramcol = Column(OneTo(2), [paramval])
        paramflow = InputSlotFlow(pname, OutputFlow(ptype, paramcol))
        push!(pflows, paramflow)
    end
    iflow = InputFlow(InputContext(), T, T[val], frameoffs, pflows)
    oflow = ev(q, iflow)
    return oflow[nothing]
end

include("query/const.jl")
include("query/it.jl")
include("query/collection.jl")
include("query/mapping.jl")
include("query/compose.jl")
include("query/record.jl")
include("query/field.jl")
include("query/lift.jl")
include("query/count.jl")
include("query/exists.jl")
include("query/aggregate.jl")
include("query/filter.jl")
include("query/sort.jl")
include("query/take.jl")
include("query/connect.jl")
include("query/group.jl")
include("query/rollup.jl")
include("query/unique.jl")
include("query/given.jl")
include("query/parameter.jl")
include("query/around.jl")
include("query/frame.jl")
include("query/switch.jl")
include("query/merge.jl")
include("query/expect.jl")
include("query/sql.jl")

describe(obj) = describe(STDOUT, obj)

function describe(io::IO, q::Query)
    describe_query(io, q.sig, q.args)
end

function describe_query(io::IO, sig::AbstractSignature, args::Vector{Query})
    describe(io, sig)
    print(io, "(")
    first = true
    for arg in args
        if !first
            print(io, ", ")
        end
        first = false
        describe_query(io, arg.sig, arg.args)
    end
end

function describe_query(io::IO, sig::AbstractPrimitive, ::Vector{Query})
    describe(io, sig)
end

function describe_query(io::IO, ::ComposeSig, args::Vector{Query})
    first = true
    for arg in args
        if !first
            print(io, ".")
        end
        first = false
        describe(io, arg)
    end
end

function describe_query(io::IO, ::RecordSig, args::Vector{Query})
    print(io, "{")
    first = true
    for arg in args
        if !first
            print(io, ", ")
        end
        first = false
        describe(io, arg)
    end
    print(io, "}")
end

function describe(io::IO, sig::AbstractSignature)
    print(io, typeof(sig))
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# Unit tests for Varpro
#

using DelimitedFiles
using Test
using Random
using Varpro

include("set1.jl")
include("gaussmix.jl")

@testset "Varpro" begin

    verbose = false
    @test runone(rexp, verbose=verbose)
    @test runone(cexp, verbose=verbose)
    @test runone(example, verbose=verbose, atol=1e-4)
    @test runone(double_exponential, verbose=verbose)
    @test runone(ctoo, verbose=verbose)
    @test runone(h1_ringdown, verbose=verbose, atol=1e-2)
    @test testgauss(;verbose=verbose)

end


text:
module TestStemmer
using Test
using SnowballStemmer

algs = stemmer_types()

@test !isempty(algs)

for alg in algs
    stmr = Stemmer(alg)
    release(stmr)
end

test_cases = Dict{String, Any}(
    "english" => Dict{String, String}(
        "working" => "work",
        "worker" => "worker",
        "aβc" => "aβc",
        "a∀c" => "a∀c"
    )
)

for (alg, test_words) in test_cases
    stmr = Stemmer(alg)
    for (n,v) in test_words
        @test v == stem(stmr, n)
    end
end

end


text:
## The grid2d domain implements the environment used for testing in our reference
## CPP repository; see the examples in https://github.com/whoenig/libMultiRobotPlanning/tree/master/example
@with_kw struct Grid2DState <: MAPFState
    time::Int64
    x::Int64
    y::Int64
end
Base.isequal(s1::Grid2DState, s2::Grid2DState) = (s1.time, s1.x, s1.y) == (s2.time, s2.x, s2.y)

function equal_except_time(s1::Grid2DState, s2::Grid2DState)
    return (s1.x == s2.x) && (s1.y == s2.y)
end

@enum Action Up=1 Down=2 Left=3 Right=4 Wait=5

struct Grid2DAction <: MAPFAction
    action::Action
end

@enum ConflictType Vertex=1 Edge=2

@with_kw struct Grid2DConflict <: MAPFConflict
    time::Int64
    agent_1::Int64
    agent_2::Int64
    type::ConflictType
    x1::Int64       = -1
    y1::Int64       = -1
    x2::Int64       = -1
    y2::Int64       = -1
end

@with_kw struct VertexConstraint
    time::Int64
    x::Int64
    y::Int64
end

Base.isless(vc1::VertexConstraint, vc2::VertexConstraint) = (vc1.time, vc1.x, vc1.x) < (vc2.time, vc2.x, vc2.y)
Base.isequal(vc1::VertexConstraint, vc2::VertexConstraint) = (vc1.time, vc1.x, vc1.y) == (vc2.time, vc2.x, vc2.y)

@with_kw struct EdgeConstraint
    time::Int64
    x1::Int64
    y1::Int64
    x2::Int64
    y2::Int64
end

Base.isless(ec1::EdgeConstraint, ec2::EdgeConstraint) = (ec1.time, ec1.x1, ec1.y1, ec1.x2, ec1.y2) < (ec2.time, ec2.x1, ec2.y1, ec2.x2, ec2.y2)
Base.isequal(ec1::EdgeConstraint, ec2::EdgeConstraint) = (ec1.time, ec1.x1, ec1.y1, ec1.x2, ec1.y2) == (ec2.time, ec2.x1, ec2.y1, ec2.x2, ec2.y2)

@with_kw struct Grid2DConstraints <: MAPFConstraints
    vertex_constraints::Set{VertexConstraint}   = Set{VertexConstraint}()
    edge_constraints::Set{EdgeConstraint}       = Set{EdgeConstraint}()
end

get_empty_constraint(::Type{Grid2DConstraints}) = Grid2DConstraints()

Base.isempty(constraints::Grid2DConstraints) = (isempty(constraints.vertex_constraints) && isempty(constraints.edge_constraints))

const Grid2DLocation = NamedTuple{(:x,:y),Tuple{Int64,Int64}}


text:
# AtmosMassEnergyLoss
#
# Dump mass and energy loss

using ..Atmos
using ..BalanceLaws
using ..Mesh.Topologies
using ..Mesh.Grids
using ..MPIStateArrays

"""
    setup_atmos_mass_energy_loss(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String,
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create and return a `DiagnosticsGroup` containing the
"AtmosMassEnergyLoss" diagnostics for Atmos LES and GCM
configurations. All the diagnostics in the group will run at the
specified `interval`, be interpolated to the specified boundaries
and resolution, and written to files prefixed by `out_prefix`
using `writer`.
"""
function setup_atmos_mass_energy_loss(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String,
    writer = NetCDFWriter(),
    interpol = nothing,
)
    @assert isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosMassEnergyLoss",
        Diagnostics.atmos_mass_energy_loss_init,
        Diagnostics.atmos_mass_energy_loss_fini,
        Diagnostics.atmos_mass_energy_loss_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

# store initial mass and energy
mutable struct AtmosMassEnergyLoss
    Σρ₀::Union{Nothing, AbstractFloat}
    Σρe₀::Union{Nothing, AbstractFloat}

    AtmosMassEnergyLoss() = new(nothing, nothing)
end
const AtmosMassEnergyLoss₀ = AtmosMassEnergyLoss()

function atmos_mass_energy_loss_init(dgngrp, currtime)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    dg = Settings.dg
    bl = dg.balance_law
    Q = Settings.Q
    FT = eltype(Q)

    energy_model(bl) isa TotalEnergyModel ||
        error("Only TotalEnergyModel supported")
    ρ_idx = varsindices(vars_state(bl, Prognostic(), FT), "ρ")
    ρe_idx = varsindices(vars_state(bl, Prognostic(), FT), :(energy.ρe))
    Σρ₀ = weightedsum(Q, ρ_idx)
    Σρe₀ = weightedsum(Q, ρe_idx)

    if mpirank == 0
        AtmosMassEnergyLoss₀.Σρ₀ = Σρ₀
        AtmosMassEnergyLoss₀.Σρe₀ = Σρe₀

        # outputs are scalars
        dims = OrderedDict()

        # set up the variables we're going to be writing
        vars = OrderedDict(
            "mass_loss" => ((), FT, Variables["mass_loss"].attrib),
            "energy_loss" => ((), FT, Variables["energy_loss"].attrib),
        )

        # create the output file
        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        noov = Settings.no_overwrite
        init_data(dgngrp.writer, dfilename, noov, dims, vars)
    end

    return nothing
end

function atmos_mass_energy_loss_collect(dgngrp, currtime)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    dg = Settings.dg
    bl = dg.balance_law
    Q = Settings.Q
    FT = eltype(Q)

    energy_model(bl) isa TotalEnergyModel ||
        error("Only TotalEnergyModel supported")
    ρ_idx = varsindices(vars_state(bl, Prognostic(), FT), "ρ")
    ρe_idx = varsindices(vars_state(bl, Prognostic(), FT), :(energy.ρe))
    Σρ = weightedsum(Q, ρ_idx)
    Σρe = weightedsum(Q, ρe_idx)

    if mpirank == 0
        δρ = (Σρ - AtmosMassEnergyLoss₀.Σρ₀) / AtmosMassEnergyLoss₀.Σρ₀
        δρe = (Σρe - AtmosMassEnergyLoss₀.Σρe₀) / AtmosMassEnergyLoss₀.Σρe₀

        varvals = OrderedDict("mass_loss" => δρ, "energy_loss" => δρe)

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

end

function atmos_mass_energy_loss_fini(dgngrp, currtime) end


text:
#= module xim.jl
   started at 2020-05-09
   aim to contain all software for x-ray image treatment

   to load it use...
   using Pkg
   Pkg.add(Pkg.PackageSpec(url="https://github.com/NelsonMCLopes/xim.git"))
   using xim
   xim_basics_testsum(1,2)
   xim_haircut_testsum(1,2)

=#
module xim

include("xim_basics.jl")
include("xim_haircut.jl")
# this is a file to contain the haircut code

export xim_basics_testsum
export xim_basics_load_packages
export xim_haircut_testsum

export xim_hotHaircut
export xim_coldHaircut

end # module
#

max_stars_repo_name: NelsonMCLopes/xim

max_issues_repo_name: NelsonMCLopes/xim

max_forks_repo_name: NelsonMCLopes/xim


text:
using DynamicPolynomials

@polyvar x
@polyvar y
@polyvar z
@polyvar w

function progs(ndims, order)
    m = 1
    dims = ["x", "y", "z", "w"]

    progs = Dict{String, Vector{String}}()
    for i in 1:ndims
        d = dims[i]
        progs[d] = Vector{String}()
        push!(progs[d], "a[$m]")
        m += 1
        for i1 in 1:ndims
            push!(progs[d], "a[$m] * $(dims[i1])")
            m += 1
            for i2 in i1:ndims
                push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2])")
                m += 1
                if order > 2
                    for i3 in i2:ndims
                        push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3])")
                        m += 1
                        if order > 3
                            for i4 in i3:ndims
                                push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3]) * $(dims[i4])")
                                m += 1
                                if order > 4
                                    for i5 in i4:ndims
                                        push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3]) * $(dims[i4]) * $(dims[i5])")
                                        m += 1
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    progs, m-1
end

function exprs(ndims, order)
    ps, as = progs(ndims, order)

    es = Dict{String, Expr}()
    for d in keys(ps)
        es[d] = Meta.parse(join(ps[d], " + "))
    end

    es, as
end

function polys(es::Dict{String, Expr}, a::Vector{Float64})
    ps = Dict{String, DynamicPolynomials.Polynomial{true,Float64}}()

    for d in keys(es)
            ps[d] = eval(es[d])
    end
    ps
end

function string2args(s::String)
    [(Int(c)-77)/10.0 for c in s]
end

function write_raw(iters, w, h)
    hrange = (minimum(iters["x"]), maximum(iters["x"]))
    hscale = (w-1) / (hrange[2]-hrange[1])
    yrange = (minimum(iters["y"]), maximum(iters["y"]))
    yscale = (h-1) / (yrange[2]-yrange[1])

    println(hrange, hscale)
    println(yrange, yscale)

    pixels = zeros(Int64, w, h)

    for i in 1:length(iters["x"])
        x = floor(Int64, 1+(iters["x"][i] - hrange[1]) * hscale)
        y = floor(Int64, 1+(iters["y"][i] - yrange[1]) * yscale)
        pixels[x, y] += 1
    end

    gmax = 255 / maximum(pixels)
    f = open("/tmp/gs.data", "w+")
    for p in pixels
        write(f, floor(UInt8, p*gmax))
    end
end

es, as = exprs(2, 2)
a = string2args("EAGHNFODVNJCP"[2:end])
ps = polys(es, a)

n = 100000000
iters = Dict{String, Vector{Float64}}()
for d in keys(ps)
    iters[d] = zeros(Float64, n)
end

states = Dict{String, Float64}()
states["x"] = 0.05
states["y"] = 0.05
states["z"] = 0.05
states["w"] = 0.05

for i in 1:n
    for d in keys(iters)
        states[d] = iters[d][i] = ps[d](x=>states["x"], y=>states["y"], z=>states["z"], w=>states["w"])
    end
end

write_raw(iters, 2400, 1800)


text:
function solve()
    n = parse(Int, readline())
    as = [parse(Int, x) for x in split(readline())]

    ret = 0
    h = as[1]
    for i=2:n
        ret += max(0, h - as[i])
        h = max(h, as[i])
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia