text:
using NPZ
using Distributions
using LinearAlgebra

"""
    generate_gaussian_data(N::Int64, D::Int64, K::Int64,MixtureVar::Number)

Generate `N` observations, generated from `K` `D` dimensions Gaussians, with the Gaussian means sampled from a `Normal` distribution with mean `0` and `MixtureVar` variance.

Returns `(Samples, Labels, Clusters_means, Clusters_cov)`

# Example
```julia
julia> x,y,clusters = generate_gaussian_data(10000,2,6,100.0)
[3644, 2880, 119, 154, 33, 3170]
...
```
"""
function generate_gaussian_data(N::Int64, D::Int64, K::Int64, MixtureVar::Int64)
    x = randn(Float32, D, N)
    tpi = rand(Dirichlet(ones(K)*0.01))
    tzn = rand(Multinomial(N, tpi))
    tz = zeros(Float32, N)

    tmean = zeros(Float32, D, K)
    tcov = zeros(Float32, D, D, K)

    ind = 1
    println(tzn)
    for i = 1:length(tzn)
        indices = ind:ind+tzn[i]-1
        tz[indices] .= i
        tmean[:, i] .=
            rand(MvNormal(zeros(Float32, D), MixtureVar * Matrix{Float32}(I, D, D)))
        tcov[:, :, i] .= rand(InverseWishart(D + 2, Matrix{Float32}(I, D, D)))
        d = MvNormal(tmean[:, i], tcov[:, :, i])
        for j in indices
            x[:, j] = rand(d)
        end
        ind += tzn[i]
    end
    x, tz, tmean, tcov
end

"""
     generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)

Generate `N` observations, generated from `K` `D` features Multinomial vectors, with `trials` draws from each vector.

Returns `(Samples, Labels, Vectors)`

# Example
```julia
julia> generate_mnmm_data(10000, 10, 5, 100)
...
```
"""
function generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)
    clusters = zeros(Float64, D, K)
    x = zeros(Float32, D, N)
    labels = rand(1:K, (N,))
    for i = 1:K
        alphas = rand(1:20, (D,))
        alphas[i] = rand(30:100)
        clusters[:, i] = rand(Dirichlet(alphas))
    end
    for i = 1:N
        x[:, i] = rand(Multinomial(trials, clusters[:, labels[i]]))
    end
    return x, labels, clusters
end


text:
using PyPlot, SegmentIntersections

function plot_random_lines()
    random_coordinates = 100 .* rand(100, 4);
    segments = Segment{Float64}[]
    for datap in eachrow(random_coordinates)
        push!(segments, Segment(datap...))
    end

    intersections_brute = find_intersections_brute(segments);

    intersections = find_intersections(segments);

    fig, ax = plt.subplots(1,2, sharex=true, sharey=true, figsize=(10,4))
    for segment in segments
        for axis in ax
            axis.plot(
                [segment.p.x, segment.q.x],
                [segment.p.y, segment.q.y],
                color = "C0",
                linewidth = 1,
            )
        end
    end
    ax[1].set_title("Brute force algorithm")
    ax[2].set_title("Bentley-Ottmann algorithm")
    xs = [point.x for point in intersections];
    ys = [point.y for point in intersections];
    ax[1].scatter(xs, ys, color = "red", s = 30)
    xs = [point.x for point in intersections_brute];
    ys = [point.y for point in intersections_brute];
    ax[2].scatter(xs, ys, color = "red", s = 30)
    plt.subplots_adjust(wspace=0, hspace=0)
    fig.savefig("random_lines.png", dpi=300, bbox_inches="tight")
end

function plot_complete_graph(n=16, r=100)
    theta_values = range(0, 2π, length=n+1)[1:end-1]
    coordinates = [[r*cos(theta), r*sin(theta)] for theta in theta_values]
    segments = Segment{Float64}[]
    for i in 1:n
        for j in (i+1):n
            push!(segments, Segment(Point(coordinates[i]...), Point(coordinates[j]...)))
        end
    end

    intersections_brute = find_intersections_brute(segments);

    intersections_bo = find_intersections(segments, tol=1e-9);

    fig, ax = plt.subplots(1,2, sharex=true, sharey=true, figsize=(10,4))
    for segment in segments
        for axis in ax
            axis.plot(
                [segment.p.x, segment.q.x],
                [segment.p.y, segment.q.y],
                color = "C0",
                linewidth = 1,
            )
        end
    end
    ax[1].set_title("Brute force algorithm")
    xs = [point.x for point in intersections_brute];
    ys = [point.y for point in intersections_brute];
    ax[1].scatter(xs, ys, color = "red", s = 30)
    xs = [point.x for point in intersections_bo];
    ys = [point.y for point in intersections_bo];
    ax[2].scatter(xs, ys, color = "red", s = 30)
    ax[2].set_title("Bentley-Ottmann algorithm")
    plt.subplots_adjust(wspace=0, hspace=0)
    fig.savefig("complete_graph.png", dpi=300, bbox_inches="tight")
end

plot_random_lines()
plot_complete_graph(14)

ext: jl


text:
module ThreadsAPIBenchmarksAnalyzers

include("analyze_task.jl")
include("analyze_notify.jl")

end # module


text:
# this function escapes the ampersand, which is probably the most common use case
function escape_ampersand(s::String)
    repl_dict = Dict("&" => "\\&")
    for (old, new) in repl_dict
        s = replace.(s, Ref(old => new))
    end
    return s
end


text:
using ADCME
using LinearAlgebra
using AdFem
using DelimitedFiles

A = readdlm("fenics/A.txt")
mesh = Mesh(8,8,1/8)
ρ = ones(get_ngauss(mesh))
B = compute_fem_mass_matrix1(ρ, mesh)
sess = Session(); init(sess)
B0 = run(sess, B)
B0 = Array(B0)
@show norm(A - B0)

A = readdlm("fenics/A2.txt")
mesh = Mesh(8, 8, 1.0 / 8,degree=2)
ρ = ones(get_ngauss(mesh))
B = compute_fem_mass_matrix1(ρ, mesh)
sess = Session(); init(sess)
B0 = run(sess, B)
B0 = Array(B0)

E = Int64.(readdlm("fenics/edges.txt"))
Edof = get_edge_dof(E, mesh)
DOF = [1:mesh.nnode; Edof .+ mesh.nnode]
B0 = B0[DOF, DOF]

@show norm(A - B0)


text:
"""
        set_kit(proj_name, proj_path)

    Create a folder named `proj_name` in `proj_path`, as well as a Pluto notebook and an internal file structure inside `proj_path/proj_name`.

    File structure and use of the notebook are detailed [here](https://github.com/ludmillafigueiredo/computational_notebooks).

    See also: [`Pluto`](https://juliapackages.com/p/pluto)

    # Example
    ```julia-repl
     julia> include("set_kit.jl")
     julia> set_kit("example_project", pwd())
    ```
    """
function set_kit(proj_name::String = "project", proj_path::String = pwd())
    ## Check project path and whether it is empty
    println(string("Will set project at ", joinpath(proj_path, proj_name)))
    println("Is that correct? Answer yes or no:")
    user_ans = readline()
    if user_ans == "no"
        println("Check help('set_kit') to set up your project's location.")
    else
        if isdir(joinpath(proj_path, proj_name))
            println("There already is a project with this name in this path. Move it or change it's name.")
        else
            ## Create project in path
            cd(proj_path)
            ## Create file structure
            ### first level
            mkdir(proj_name)
            ### main folders of a project
            map(mkdir, map(x -> joinpath(proj_name, x),
                           ["results","text", "submission"]))
            ### main folders of results
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"results"), x),
                           ["data", "scripts"]))
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"results", "data"), x),
                           ["raw", "processed"]))
            ### main folders of text
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"text"), x),
                           ["figures","tables", "supplementary"]))

            ### main folder
            touch(joinpath(proj_path, proj_name, "README.txt"))
            open(joinpath(proj_path, proj_name, "README.txt"), "w") do io
                write(io, "This folder contains the set up for a reproducible workflow as described by https://github.com/ludmillafigueiredo/computational_notebooks.git\n")
                write(io, "\n")
                write(io, "The file structure is organized as such:\n")
                write(io, "\n")
                write(io, "projet_name\n")
                write(io, "|-- README.txt\n")
                write(io, "|-- main.Rmd\n")
                write(io, "|-- results\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- data\n")
                write(io, "|   |   |-- raw\n")
                write(io, "|   |   |-- processed\n")
                write(io, "|   |-- scripts\n")
                write(io, "|-- text\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- main.doc\n")
                write(io, "|   |-- figures\n")
                write(io, "|   |-- tables\n")
                write(io, "|   |-- supplementary\n")
                write(io, "|-- submission\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- journal1\n")
                write(io, "|       |-- first\n")
                write(io, "|   |-- journal2\n")
                write(io, "|       |-- first\n")
                write(io, "|       |-- revisions\n")
            end

            ### results folder
            touch(joinpath(proj_path, proj_name, "results/README.txt"))
            open(joinpath(proj_path, proj_name, "results/README.txt"), "w") do io
                write(io, "This folder contains all files of results or their processing, organized in the following subfolders:")
                write(io, "\n")
                write(io, "`data/raw`: your raw data files. These should not be protected against any change after the first storage")
                write(io, "\n")
                write(io, "`data/processed`: these are files generated by processing the raw data and these are the ones used in the analysis. It can be a copy of the raw data if that is already ready to use. Any processing of the raw data to generate the files here should be documented in the notebook.")
                write(io, "\n")
                write(io, "`scripts`: all code used to process the data, and which, for some reason or another is not included in the notebook because they too cumbersome or not of upmost relevancy for comprehension.")
            end

            ### text folder
            touch(joinpath(proj_path, proj_name, "text/README.txt"))
            open(joinpath(proj_path, proj_name, "text/README.txt"), "w") do io
                write(io, "This folder contains the main text of the manuscript, folders containing the figures and tables (unformatted) to be included in it, as well as a folder with the supplementary material.")
            end

            ### submission folder
            touch(joinpath(proj_path, proj_name, "submission/README.txt"))
            open(joinpath(proj_path, proj_name, "submission/README.txt"), "w") do io
                write(io, "This folder contains the files specific to journal submissions, e.g. cover letters, submitted versions.")
            end

            ## Create the minimal notebook
            touch(joinpath(proj_path, proj_name, string(proj_name, ".jl")))
            open(joinpath(proj_path, proj_name, string(proj_name, ".jl")), "w") do io
                write(io, "### A Pluto.jl notebook ###")
                write(io, "\n\n")
                write(io, "using Markdown")
                write(io, "\n")
                write(io, "using InteractiveUtils")
                write(io, "\n\n")
                write(io, "# ╔═╡ 307335ba-cf46-11eb-28b0-f199cf048ae6")
                write(io, "\n")
                write(io, "md\" ## Your title here\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ 37c7bffc-cf46-11eb-19b8-4175e38e818b")
                write(io, "\n")
                write(io, "## Paths for easy call from the Pluto notebook")
                write(io, "\n")
                write(io, "begin")
                write(io, "\n")
                write(io, "\tdata_dir = joinpath(\"results\", \"clean_data\") ## do NOT play with stuff in raw_data. That is your back-up")
                write(io, "\n")
                write(io, "\tscripts_dir = joinpath(\"results\", \"scripts\")")
                write(io, "\n")
                write(io, "\tsuppl_dir = joinpath(\"results\", \"supplementary\")")
                write(io, "\n")
                write(io, "\tsemiprods_dir = joinpath(\"results\", \"semi_products\")")
                write(io, "\n")
                write(io, "\tfigures_dir = joinpath(\"text\", \"figures\")")
                write(io, "\n")
                write(io, "\ttables_dir = joinpath(\"text\", \"tables\")")
                write(io, "\n")
                write(io, "end")
                write(io, "\n\n\n")
                write(io, "# ╔═╡ Cell order:")
                write(io, "\n")
                write(io, "# ╠═307335ba-cf46-11eb-28b0-f199cf048ae6")
                write(io, "\n")
                write(io, "# ╠═37c7bffc-cf46-11eb-19b8-4175e38e818b")

                println("Your project is ready to go")
            end
        end
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function temp_pkg_dir(fn::Function)
    local env_dir
    local old_load_path
    local old_depot_path
    local old_home_project
    local old_active_project
    try
        old_load_path = copy(LOAD_PATH)
        old_depot_path = copy(DEPOT_PATH)
        old_home_project = Base.HOME_PROJECT[]
        old_active_project = Base.ACTIVE_PROJECT[]
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        Base.HOME_PROJECT[] = nothing
        Base.ACTIVE_PROJECT[] = nothing
        withenv("JULIA_PROJECT" => nothing, "JULIA_LOAD_PATH" => nothing) do
            mktempdir() do env_dir
                mktempdir() do depot_dir
                    push!(LOAD_PATH, "@", "@v#.#", "@stdlib")
                    push!(DEPOT_PATH, depot_dir)
                    fn(env_dir)
                end
            end
        end
    finally
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        append!(LOAD_PATH, old_load_path)
        append!(DEPOT_PATH, old_depot_path)
        Base.HOME_PROJECT[] = old_home_project
        Base.ACTIVE_PROJECT[] = old_active_project
    end
end

function cd_tempdir(f)
    mktempdir() do tmp
        cd(tmp) do
            f(tmp)
        end
    end
end

isinstalled(pkg) = Base.locate_package(Base.PkgId(pkg.uuid, pkg.name)) !== nothing

function write_build(path, content)
    build_filename = joinpath(path, "deps", "build.jl")
    mkpath(dirname(build_filename))
    write(build_filename, content)
end

function with_current_env(f)
    Pkg.activate(".")
    try
        f()
    finally
        Pkg.activate()
    end
end

function with_temp_env(f, env_name::AbstractString="Dummy")
    env_path = joinpath(mktempdir(), env_name)
    Pkg.generate(env_path)
    Pkg.activate(env_path)
    try
        applicable(f, env_path) ? f(env_path) : f()
    finally
        Pkg.activate()
    end
end

function with_pkg_env(fn::Function, path::AbstractString="."; change_dir=false)
    Pkg.activate(path)
    try
        if change_dir
            cd(fn, path)
        else
            fn()
        end
    finally
        Pkg.activate()
    end
end

import LibGit2
using UUIDs
const TEST_SIG = LibGit2.Signature("TEST", "TEST@TEST.COM", round(time()), 0)
const TEST_PKG = (name = "Example", uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a"))


text:
using Test

using DesignAutomation

@testset "bracket" begin
    fn(x) = 2 - 4 * x + exp(x)
    res = [round(i, digits=3) for i in bracket(fn, .5)]
    @test res == [.5, 1.309, 2.618]
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
dense_conv(n, k, k0) = Chain(
    BatchNorm(n * k + k0, relu),
    # The 1x1 convolutions produce 4K feature maps
    Conv( (1,1), n * k + k0 => 4 * k),
    BatchNorm(4 * k, relu),
    Conv( (3,3), 4 * k => k; pad = 1)
)

struct DenseBlock
    layers::Vector
    transition
    apply_transition::Bool
end

Flux.@functor DenseBlock

function (D::DenseBlock)(x)
    for f in D.layers
        y = f(x)
        x = cat(x, y; dims = 3)
    end
    # Transition Layer
    ret = D.apply_transition ? D.transition(x) : x
    return ret
end

function DenseBlock(k::Int, k0, N::Int; apply_transition = true)
    layers = [dense_conv(i-1, k, k0) for i in 1:N]

    transition = Chain(
        BatchNorm( N * k + k0, relu),
        Conv( (1,1), N * k + k0 => k0),
        x -> meanpool(x, (2,2); pad = 0, stride = 2),
    )
    return DenseBlock(layers, transition, apply_transition)
end

function DenseNet(k)
    counts = (6, 12, 64, 48)
    return Chain(
        Conv( (7,7), 3 => 2 * k, relu; pad = 3, stride = 2),
        x -> maxpool(x, (3,3); pad = 1, stride = 2),
        DenseBlock(k, 2*k, counts[1]),
        DenseBlock(k, 2*k, counts[2]),
        DenseBlock(k, 2*k, counts[3]),
        DenseBlock(k, 2*k, counts[4]; apply_transition = false),
        x -> meanpool(x, (7,7); pad = 0),
        x -> reshape(x, :, size(x,4)),
        Dense(counts[4] * k + 2 * k, 1000),
        x -> log.(max.(x, Float32(1e-9))),
        softmax
    )
end

function densenet_training(batchsize; growth = 32)
    x = rand(Float32, 224, 224, 3, batchsize)
    x = (x .- mean(x)) ./ std(x)

    y = rand(Float32, 1000, batchsize)
    random_labels!(y)

    f = ForwardLoss(DenseNet(growth), Flux.crossentropy)
    return Actualizer(f, x, y; optimizer = nGraph.SGD(Float32(0.001)))
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
# QP example
#
#    minimize 2 x^2 + y^2 + xy + x + y
#
#       s.t.  x, y >= 0
#             x + y = 1
#
#    solution: (0.25, 0.75), objv = 1.875
#

using CPLEX, Base.Test

@testset "QP 01" begin
    env = CPLEX.Env()
    CPLEX.set_param!(env, "CPX_PARAM_SCRIND", 0)

    model = CPLEX.Model(env, "qp_02")

    CPLEX.add_vars!(model, [1., 1.], 0., Inf)

    CPLEX.add_qpterms!(model, [1, 1, 2], [1, 2, 2], [2., 1., 1.])
    CPLEX.add_constr!(model, [1., 1.], '=', 1.)

    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test abs(sol[1]) < 1e-4
    @test isapprox(sol[2], 1, rtol=1e-4)
    @test isapprox(CPLEX.get_objval(model), 1.5, rtol=1e-8)
end


text:
#=
String classification, normalization and width functions

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

# Recommended by deprecate
text_width(str::Str) = mapreduce(text_width, +, str; init=0)
Base.Unicode.normalize(str::Str, opt::Symbol) = normalize(str, opt)

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

text_width(str::Str{Union{ASCIICSE,Latin_CSEs}}) = length(str)

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

function is_latin(str::MaybeSub{String})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + siz
        while pnt < fin
            cu = get_codeunit(pnt)
            # cu must be 1) 0-0x7f, or 2) 0xc2 or 0xc3 followed by 0x80-0xbf
            (cu < 0x7f ||
             ((cu - 0xc2) < 0x02 &&
              (pnt += 1) < fin && is_valid_continuation(get_codeunit(pnt)))) ||
              return false
            pnt += 1
        end
        true
    end
end

@inline function check_3byte(cu, pnt)
    b2 = get_codeunit(pnt-1)
    b3 = get_codeunit(pnt)
    is_valid_continuation(b2) && is_valid_continuation(b3) &&
        !is_surrogate_codeunit(_mskup32(cu, 0xf, 12) | _mskup32(b2, 0x3f, 6) | (b3 & 0x3f))
end

function is_bmp(str::MaybeSub{String})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + siz
        while pnt < fin
            cu = get_codeunit(pnt)
            # cu must be 1) 0-0x7f, or 2) 0xc2 or 0xc3 followed by 0x80-0xbf
            # c2-df -> de,df
            (cu < 0x7f ||
             ((cu - 0xc2) < 0x1e && (pnt += 1) < fin && checkcont(pnt)) ||
             ((cu - 0xe0) < 0x0f && (pnt += 2) < fin && check_3byte(cu, pnt))) ||
             return false
             pnt += 1
        end
        true
    end
end

function is_latin(str::AbstractString)
    @inbounds for ch in str
        is_latin(ch) || return false
    end
    true
end

function is_bmp(str::AbstractString)
    @inbounds for ch in str
        is_bmp(ch) || return false
    end
    true
end

function is_unicode(str::AbstractString)
    @inbounds for ch in str
        is_unicode(ch) || return false
    end
    true
end


text:
using BinDeps
@BinDeps.setup

visa = library_dependency("visa", aliases = ["visa64","VISA","/Library/Frameworks/VISA.framework/VISA"])
visa_path_found = BinDeps._find_library(visa)

if ~isempty(visa_path_found)
    @BinDeps.install Dict(:visa=>:libvisa)
else
    @warn "No VISA libraries found, please check VISA installation for visa64 library!"
end


text:
abstract type AbstractCoefficientCache{T} end

function Base.checkbounds(::Type{Bool}, u::AbstractVector, cache::AbstractCoefficientCache)
    length(u) > length(cache.inv_left_weights)+length(cache.inv_right_weights)
end

"""
    VarCoefDerivativeCoefficients

The coefficients of a variable coefficient derivative operator on a nonperiodic grid.
"""
struct VarCoefDerivativeCoefficients{T,CoefficientCache<:AbstractCoefficientCache{T},
                                     LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients} <: AbstractDerivativeCoefficients{T}
    # coefficients defining the operator and its action
    coefficient_cache::CoefficientCache
    left_weights::SVector{LeftWidth, T}
    right_weights::SVector{RightWidth, T}
    mode::ExecutionMode
    # corresponding orders etc.
    derivative_order::Int
    accuracy_order  ::Int
    source_of_coefficients::SourceOfCoefficients

    function VarCoefDerivativeCoefficients(coefficient_cache::CoefficientCache,
                                           left_weights::SVector{LeftWidth, T},
                                           right_weights::SVector{RightWidth, T},
                                           mode::ExecutionMode,
                                           derivative_order::Int, accuracy_order::Int,
                                           source_of_coefficients::SourceOfCoefficients) where {T,CoefficientCache<:AbstractCoefficientCache{T},
                                                                                LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients}
        new{T,CoefficientCache,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients}(
            coefficient_cache, left_weights, right_weights, mode, derivative_order, accuracy_order, source_of_coefficients)
    end
end

function Base.show(io::IO, coefficients::VarCoefDerivativeCoefficients)
    if derivative_order(coefficients) == 1
        print(io, "Coefficients of the variable-coefficient first-derivative operator")
    elseif  derivative_order(coefficients) == 2
        print(io, "Coefficients of the variable-coefficient second-derivative operator")
    elseif  derivative_order(coefficients) == 3
        print(io, "Coefficients of the variable-coefficient third-derivative operator")
    else
        print(io, "Coefficients of the variable-coefficient ", derivative_order(coefficients),
              "-derivative operator")
    end
    print(io, " of order ", accuracy_order(coefficients), " of ")
    print(io, source_of_coefficients(coefficients))
end

# Compute `α*D*u + β*dest` using the coefficients `b` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                u::AbstractVector, b::AbstractVector, α, β)
    @unpack coefficient_cache, mode = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(dest) == length(b) DimensionMismatch
        @argcheck checkbounds(Bool, dest, coefficient_cache) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, coefficient_cache, u, b, α, β)
    convolve_interior_coefficients!(dest, coefficient_cache, u, b, α, β, mode)
end

# Compute `α*D*u` using the coefficients `b` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                u::AbstractVector, b::AbstractVector, α)
    @unpack coefficient_cache, mode = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(dest) == length(b) DimensionMismatch
        @argcheck checkbounds(Bool, dest, coefficient_cache) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, coefficient_cache, u, b, α)
    convolve_interior_coefficients!(dest, coefficient_cache, u, b, α, mode)
end

@inline function left_length(cache::AbstractCoefficientCache)
    length(cache.inv_left_weights)
end

@inline function right_length(cache::AbstractCoefficientCache)
    length(cache.inv_right_weights)
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, mode)
    for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, mode::ThreadedMode)
    Threads.@threads for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, β, mode)
    for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval + β*dest[i]
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, β, mode::ThreadedMode)
    Threads.@threads for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval + β*dest[i]
    end end
end

abstract type AbstractVariableCoefficientNonperiodicDerivativeOperator{T} <: AbstractNonperiodicDerivativeOperator{T} end
abstract type AbstractVariableCoefficientPeriodicDerivativeOperator{T} <: AbstractPeriodicDerivativeOperator{T} end

@inline source_of_coefficients(D::Union{AbstractVariableCoefficientNonperiodicDerivativeOperator,AbstractVariableCoefficientNonperiodicDerivativeOperator}) = source_of_coefficients(D.coefficients)

@inline function lower_bandwidth(D::AbstractVariableCoefficientNonperiodicDerivativeOperator)
    lower_bandwidth(D.coefficients.coefficient_cache)
end

@inline function upper_bandwidth(D::AbstractVariableCoefficientNonperiodicDerivativeOperator)
    upper_bandwidth((D.coefficients.coefficient_cache))
end

"""
    DissipationOperator

A dissipation operator on a nonperiodic finite difference grid.
See [`dissipation_operator`](@ref).
"""
struct DissipationOperator{T,Coefficients<:VarCoefDerivativeCoefficients{T},Grid} <: AbstractVariableCoefficientNonperiodicDerivativeOperator{T}
    factor::T
    coefficients::Coefficients
    grid::Grid
    b::Vector{T}

    function DissipationOperator(factor::T,
                                 coefficients::Coefficients,
                                 grid::Grid, b::Vector{T}) where {T,Coefficients<:VarCoefDerivativeCoefficients{T},Grid}
        @argcheck checkbounds(Bool, grid, coefficients.coefficient_cache) DimensionMismatch
        @argcheck length(grid) == length(b)
        factor < 0 && @warn("Negative dissipation strength shouldn't be used.")

        new{T,Coefficients,Grid}(factor,coefficients, grid, b)
    end
end

function Base.show(io::IO, D::DissipationOperator{T}) where {T}
    if  derivative_order(D) == 2
        print(io, "SBP second-derivative dissipation operator")
    else
        print(io, "SBP ", derivative_order(D),
              "-derivative dissipation operator")
    end
    print(io, " of order ", accuracy_order(D))
    if get(io, :compact, false) == false
        print(io, " on a grid in [", first(grid(D)), ", ", last(grid(D)),
                "] using ", length(grid(D)), " nodes \n")
        print(io, "and coefficients")
    end
    print(io, " of ", source_of_coefficients(D))
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                       u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α, β)
end

# Compute `α*D*u` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                       u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α)
end

"""
    dissipation_operator(source_of_coefficients, order, xmin, xmax, N,
                         left_weights, right_weights, mode=FastMode())

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a weighted `order`-th derivative on a grid between `xmin` and
`xmax` with `N` grid points up to order of accuracy 2 with coefficients given
by `source_of_coefficients`.
The norm matrix is given by `left_weights` and `right_weights`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.
"""
function dissipation_operator(source_of_coefficients, order, xmin, xmax, N,
                              left_weights, right_weights,
                              strength=one(xmin+xmax),
                              mode=FastMode(),
                              parallel=nothing)
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :dissipation_operator)
        mode = _parallel_to_mode(parallel)
    end
    grid = construct_grid(source_of_coefficients, order, xmin, xmax, N)
    coefficients, b = dissipation_coefficients(source_of_coefficients, order, grid, left_weights, right_weights, mode)
    DissipationOperator(strength, coefficients, grid, b)
end

"""
    dissipation_operator([source_of_coefficients=MattssonSvärdNordström2004()],
                         D::DerivativeOperator{T};
                         strength=one(T),
                         order::Int=accuracy_order(D),
                         mode=D.coefficients.mode)

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a weighted `order`-th derivative adapted to the derivative
operator `D` with coefficients given in `source_of_coefficients`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.
"""
function dissipation_operator(source_of_coefficients, D::DerivativeOperator{T};
                              strength=one(T),
                              order::Int=accuracy_order(D),
                              mode=D.coefficients.mode,
                              parallel=nothing) where {T}
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :dissipation_operator)
        mode = _parallel_to_mode(parallel)
    end
    x = grid(D)
    dissipation_operator(source_of_coefficients, order, first(x), last(x), length(x), D.coefficients.left_weights, D.coefficients.right_weights, strength, mode)
end

function dissipation_operator(D::DerivativeOperator; kwargs...)
    source = MattssonSvärdNordström2004()
    dissipation_operator(source, D; kwargs...)
end


text:
# Model
include_model("hopper")

# Free-time model
model_ft = free_time_model(model)

# Horizon
T = 49
Tm = 25

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= 250.0 #model_ft.uU
_uu[end] = 5.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -250.0 #model_ft.uL
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.5 , 0.0, 0.5]
q_right = [-0.25, 0.5 + 0.25, pi / 2.0, 0.25]
q_top = [-0.5, 0.5 + 0.5, pi, 0.25]
q_left = [-0.75, 0.5 + 0.25, 3.0 * pi / 2.0, 0.25]
qT = [-1.0, 0.5,  2.0 * pi, 0.5]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; q1],
		xT = [qT; qT])
xl[Tm][nq .+ (1:nq)] = q_top
xu[Tm][nq .+ (1:nq)] = q_top

q_ref = [linear_interpolation(q1, q_right, 13)...,
         linear_interpolation(q_right, q_top, 13)[2:end]...,
         linear_interpolation(q_top, q_left, 13)[2:end]...,
         linear_interpolation(q_left, qT, 13)[2:end]...]

x_ref = configuration_to_state(q_ref)

# Objective
include_objective(["velocity"])
obj_tracking = quadratic_time_tracking_objective(
    [100.0 * Diagonal([1.0; 10.0; 0.1; 0.1; 1.0; 10.0; 0.1; 0.1]) for t = 1:T],
    [Diagonal([1.0, 1.0,
		zeros(model_ft.m - model_ft.nu)...]) for t = 1:T-1],
    [x_ref[t] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T-1],
    1.0)

obj_contact_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)
obj_velocity = velocity_objective([Diagonal(ones(model_ft.nq)) for t = 1:T],
	model_ft.nq, idx_angle = (3:3))
obj = MultiObjective([obj_tracking, obj_contact_penalty, obj_velocity])

# Constraints
include_constraints(["contact", "free_time"])
con_contact = contact_constraints(model_ft, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_free_time, con_contact])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = deepcopy(x_ref) # linear interpolation on state
u0 = [[1.0e-3 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
# include_snopt()

@time z̄ , info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

x̄, ū = unpack(z̄, prob)

@show ū[4][end]
@show check_slack(z̄, prob)

using Plots
tf, t, h̄ = get_time(ū)
plot(t[1:end-1], hcat(ū...)[1:2,:]', linetype=:steppost,
	xlabel="time (s)", ylabel = "control",
	label = ["angle" "length"],
	width = 2.0, legend = :top)
plot(t[1:end-1], h̄, linetype=:steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model_ft,
	state_to_configuration([[x̄[1] for i = 1:50]...,x̄..., [x̄[end] for i = 1:50]...]),
	Δt = h̄[1],
	scenario = :flip)


text:
## Computes a direct plug-in selector of the
## bandwidth for local linear regression as
## described in the 1996 J. Amer. Statist. Assoc.
## paper by Ruppert, Sheather and Wand.

function dpill(x::Vector{Float64}, y::Vector{Float64};
               blockmax::Int = 5,
               divisor::Int = 20,
               trim::Float64 = 0.01,
               proptrun::Float64 = 0.05,
               gridsize::Int = 401,
               range_x::Vector{Float64} = Float64[],
               truncate = true)
    ## Trim the 100(trim)% of the data from each end (in the x-direction).
    sp = sortperm(x)
    x = x[sp] #don't sort in place so original x and y are not muted
    y = y[sp]
    indlow = ifloor(trim*length(x)) + 1
    indupp = length(x) - ifloor(trim*length(x))

    x = x[indlow:indupp]
    y = y[indlow:indupp]

    ## Rename common parameters
    n = length(x)
    M = gridsize
    if range_x == Float64[]
        range_x = [minimum(x), maximum(x)]
    end

    a = range_x[1]
    b = range_x[2]

    ## Bin the data

    gpoints = linspace(a, b, M)
    xcounts, ycounts = rlbin(x, y, gpoints, truncate)

    ## Choose the value of N using Mallow's C_p
    Nmax = max(min(ifloor(n/divisor), blockmax), 1)
    Nval = cpblock(x, y, Nmax, 4)

    ## Estimate sig^2, theta_22 and theta_24 using quartic fits
    ## on "Nval" blocks.

    (sigsqQ, th22Q, th24Q) = blkest(x, y, Nval, 4)

    ## Estimate theta_22 using a local cubic fit
    ## with a "rule-of-thumb" bandwidth: "gamseh"

    gamseh = (sigsqQ*(b-a)/(abs(th24Q)*n))

    if th24Q < 0.0
        gamseh = (3.0 * gamseh / (8.0 * sqrt(pi)))^(1.0/7.0)
    end
    if th24Q > 0.0
        gamseh = (15.0 * gamseh / (16.0 * sqrt(pi)))^(1.0/7.0)
    end

    mddest = locpoly(xcounts, ycounts, gamseh, drv=2, range_x=range_x, binned=true)[2]

    llow = ifloor(proptrun*M) + 1
    lupp = M - ifloor(proptrun*M)
    th22kn = sum((mddest[llow:lupp].^2) .* xcounts[llow:lupp])/n

    ## Estimate sigma^2 using a local linear fit
    ## with a "direct plug-in" bandwidth: "lamseh"
    C3K = (1.0/2.0) + 2.0*sqrt(2.0) - (4.0/3.0)*sqrt(3.0)
    C3K = (4.0*C3K/(sqrt(2.0*pi)))^(1.0/9.0)
    lamseh = C3K*(((sigsqQ^2)*(b-a)/((th22kn*n)^2))^(1.0/9.0))

    ## Now compute a local linear kernel estimate of
    ## the variance.
    mest = locpoly(xcounts, ycounts, lamseh, range_x=range_x, binned=true)[2]
    Sdg = sdiag(xcounts, lamseh, range_x=range_x, binned=true)[2]
    SSTdg = sstdiag(xcounts, lamseh, range_x=range_x, binned=true)[2]
    sigsqn = sum(y.^2) .- 2*sum(mest.*ycounts) + sum((mest.^2).*xcounts)
    sigsqd = n - 2*sum(Sdg.*xcounts) + sum(SSTdg.*xcounts)
    sigsqkn = sigsqn/sigsqd

    ## Combine to obtain final answer.
    (sigsqkn*(b-a)/(2.0*sqrt(pi)*th22kn*n))^(1.0/5.0)
end


text:
using Merlin

function maxrows(m::Matrix)
  _, inds = findmax(m, 1)
  map!(i -> ind2sub(size(m), i)[1], inds)
  vec(inds)
end

function train(path)
  # data
  worddict = read_wordlist("$(path)/words.lst")
  chardict, tagdict = Dict(), Dict()
  traindata = read_conll("$(path)/wsj_00-18.conll", true, worddict, chardict, tagdict)
  println("#word: $(length(worddict)), #char: $(length(chardict)), #tag: $(length(tagdict))")
  traindata = traindata[1:5000]
  testdata = read_conll("$(path)/wsj_22-24.conll", false, worddict, chardict, tagdict)

  # model
  model = POSModel(path)
  opt = SGD(0.0075)

  for iter = 1:10
    println("iter: $(iter)")
    golds, preds = Int[], Int[]
    opt.learnrate = 0.0075 / iter
    loss = 0.0

    for i in randperm(length(traindata))
      tokens = traindata[i]
      append!(golds, map(t -> t.tagid, tokens))

      # forward & prediction
      out = forward(model, tokens)
      maxidx = maxrows(out.value)
      append!(preds, maxidx)

      # loss function
      px = zeros(Float32, size(out.value,1), length(tokens))
      for j = 1:length(tokens)
        px[tokens[j].tagid, j] = 1.0
      end
      out = CrossEntropy(px)(out)
      loss += sum(out.value)

      # backward & update
      backward!(out)
      update!(model, opt)
    end
    println("loss: $(loss)")
    acc = accuracy(golds, preds)
    println("train acc: $(acc)")
    decode(model, testdata)
    println("")
  end
  println("finish")
end

function decode(m::POSModel, data::Vector{Vector{Token}})
  golds, preds = Int[], Int[]
  for i = 1:length(data)
    tokens = data[i]
    append!(golds, map(t -> t.tagid, tokens))
    out = forward(m, tokens)
    append!(preds, maxrows(out.value))
  end
  acc = accuracy(golds, preds)
  println("test acc: $(acc)")
end

max_stars_repo_name: hshindo/POSTagging-benchmark

max_issues_repo_name: hshindo/POSTagging-benchmark

max_forks_repo_name: hshindo/POSTagging-benchmark


text:
@testset "Lasso small (v. split, $T)" for T in [Float32, Float64, ComplexF32, ComplexF64]
    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra
    using AbstractOperators: MatrixOp

    A1 = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
    ]
    A2 = T[
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    A = vcat(MatrixOp(A1), MatrixOp(A2))
    b1 = T[1.0, 2.0]
    b2 = T[3.0, 4.0]

    m1, n = size(A1)
    m2, _ = size(A2)
    m = m1 + m2

    R = real(T)

    lam = R(0.1) * norm([A1; A2]' * [b1; b2], Inf)
    @test typeof(lam) == R

    f = SeparableSum(Translate(SqrNormL2(R(1)), -b1), Translate(SqrNormL2(R(1)), -b2))
    f2 = Sum(LeastSquares(A1, b1), LeastSquares(A2, b2))
    g = NormL1(lam)

    Lf = opnorm([A1; A2])^2

    x_star = T[-3.877278911564627e-01, 0, 0, 2.174149659863943e-02, 6.168435374149660e-01]

    TOL = R(1e-4)

    @testset "ForwardBackward" begin

        ## Nonfast/Nonadaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.ForwardBackward(tol = TOL)
        x, it = solver(x0, f = f2, g = g, Lf = Lf)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 150
        @test x0 == x0_backup

        # Nonfast/Adaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.ForwardBackward(tol = TOL, adaptive = true)
        x, it = solver(x0, f = f2, g = g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 300
        @test x0 == x0_backup

        # Fast/Nonadaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.FastForwardBackward(tol = TOL)
        x, it = solver(x0, f = f2, g = g, Lf = Lf)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 100
        @test x0 == x0_backup

        # Fast/Adaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver =
            ProximalAlgorithms.FastForwardBackward(tol = TOL, adaptive = true)
        x, it = solver(x0, f = f2, g = g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 200
        @test x0 == x0_backup
    end

    @testset "ZeroFPR" begin

        # ZeroFPR/Nonadaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.ZeroFPR(tol = TOL)
        x, it = solver(x0, f = f, A = A, g = g, Lf = opnorm([A1; A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20
        @test x0 == x0_backup

        # ZeroFPR/Adaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.ZeroFPR(adaptive = true, tol = TOL)
        x, it = solver(x0, f = f, A = A, g = g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20
        @test x0 == x0_backup

    end

    @testset "PANOC" begin

        # PANOC/Nonadaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.PANOC(tol = TOL)
        x, it = solver(x0, f = f, A = A, g = g, Lf = opnorm([A1; A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20
        @test x0 == x0_backup

        ## PANOC/Adaptive

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.PANOC(adaptive = true, tol = TOL)
        x, it = solver(x0, f = f, A = A, g = g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20
        @test x0 == x0_backup

    end

end


text:
@testset "QP Example from Maros" begin
  @testset "parseqps($T, filename)" for T in subtypes(AbstractFloat)
    # MPS description
    mps   = parseqps(T, joinpath(dirname(@__FILE__), "qp-example.qps"))

    @test name(mps) == "QP example"

    @test variables(mps) == [:c1, :c2]

    Q, q₁, q₂ = objective(mps)
    @test Q ≈ [8 2; 2 10] && eltype(Q) == T
    @test q₁ ≈ [1.5; -2] && eltype(q₁) == T
    @test q₂ ≈ 4 && eltype(q₂) == T

    A, b = equalities(mps)
    @test isempty(A) == true && eltype(A) == T
    @test isempty(b) == true && eltype(b) == T

    C, c₁, c₂ = inequalities(mps)
    @test C ≈ [2 1; -1 2] && eltype(C) == T
    @test c₁ ≈ [2; -Inf] && eltype(c₁) == T
    @test c₂ ≈ [Inf, 6] && eltype(c₂) == T

    lb, ub = bounds(mps)
    @test lb ≈ zeros(2) && eltype(lb) == T
    @test ub ≈ [20, Inf] && eltype(ub) == T
    # Canonical description
    canon = canonical(mps)

    @test name(canon) == name(mps)

    @test variables(canon) == variables(mps)

    QQ, qq1, qq2 = objective(canon)
    @test QQ ≈ Q && qq1 ≈ q₁ && qq2 ≈ q₂

    AA, bb = equalities(canon)
    @test AA ≈ A && bb ≈ b

    C̃, c̃ = inequalities(canon)
    @test C̃ ≈ [-2 -1; -1 2] && eltype(C̃) == T
    @test c̃ ≈ [-2; 6] && eltype(c̃) == T

    llb, uub = bounds(canon)
    @test llb ≈ lb && uub ≈ ub
  end

  @testset "parse_sparseqps($T, filename)" for T in subtypes(AbstractFloat)
    # MPS description
    mps   = parse_sparseqps(T, joinpath(dirname(@__FILE__), "qp-example.qps"))

    @test name(mps) == "QP example"

    @test variables(mps) == [:c1, :c2]

    Q, q₁, q₂ = objective(mps)
    @test typeof(Q) == SparseMatrixCSC{T,Int} && Matrix(Q) ≈ [8 2; 2 10]
    @test q₁ ≈ [1.5; -2] && eltype(q₁) == T
    @test q₂ ≈ 4 && eltype(q₂) == T

    A, b = equalities(mps)
    @test typeof(A) == SparseMatrixCSC{T,Int} && isempty(A) == true
    @test isempty(b) == true && eltype(b) == T

    C, c₁, c₂ = inequalities(mps)
    @test typeof(C) == SparseMatrixCSC{T,Int} && Matrix(C) ≈ [2 1; -1 2]
    @test c₁ ≈ [2; -Inf] && eltype(c₁) == T
    @test c₂ ≈ [Inf, 6] && eltype(c₂) == T

    lb, ub = bounds(mps)
    @test lb ≈ zeros(2) && eltype(lb) == T
    @test ub ≈ [20, Inf] && eltype(ub) == T
    # Canonical description
    canon = canonical(mps)

    @test name(canon) == name(mps)

    @test variables(canon) == variables(mps)

    QQ, qq1, qq2 = objective(canon)
    @test QQ ≈ Q && qq1 ≈ q₁ && qq2 ≈ q₂

    AA, bb = equalities(canon)
    @test AA ≈ A && bb ≈ b

    C̃, c̃ = inequalities(canon)
    @test typeof(C̃) == SparseMatrixCSC{T,Int} && Matrix(C̃) ≈ [-2 -1; -1 2]
    @test c̃ ≈ [-2; 6] && eltype(c̃) == T

    llb, uub = bounds(canon)
    @test llb ≈ lb && uub ≈ ub
  end
end


text:
#
# Protein - Water
#

using PDBTools
include("../../src/ComplexMixtures.jl")

dir="/home/leandro/Documents/Emilia/corrected_concentration/50"
pdb="$dir/SC9_50%gly_corrected.pdb"
traj="$dir/01/production.dcd"

atoms = PDBTools.readPDB(pdb)
solute_indexes = PDBTools.select(atoms,"protein")
solute = ComplexMixtures.Selection( solute_indexes, nmols=1 )
solvent_indexes = PDBTools.select(atoms,"resname TIP3")
solvent = ComplexMixtures.Selection( solvent_indexes, natomspermol=3 )
trajectory = ComplexMixtures.Trajectory(traj,solute,solvent)
options = ComplexMixtures.Options(binstep=0.2,n_random_samples=10)
@time R = ComplexMixtures.mddf(trajectory,options)


text:
support_scitype(::Type) = Unknown
support_scitype(d) = support_scitype(typeof(d))
support_scitype(::Type{<:UnivariateFinite}) = Finite

# NOTE: see also interface/univariate_finite.jl

"""
    levels(d::UnivariateFinite)

A list of the raw levels in the common pool of classes used to
construct `d`, equal to `get.(classes(d))`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    levels(d) # Array{String, 1}["maybe", "no", "yes"]
"""
levels(d::UnivariateFinite)  = d.decoder.pool.levels

# get the internal integer representations of the support
raw_support(d::UnivariateFinite) =
    sort!(collect(keys(d.prob_given_class)))

"""
    Distributions.support(d::UnivariateFinite)

Those classes associated with non-zero probabilities.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    support(d) # CategoricalArray{String,1,UInt32}["maybe", "no"]
"""
Distributions.support(d::UnivariateFinite) = map(d.decoder, raw_support(d))

function Base.show(stream::IO, d::UnivariateFinite)
    raw = sort!(raw_support(d)) # reflects order of pool at
                                # instantiation of d
    x1 = d.decoder(first(raw))
    p1 = d.prob_given_class[first(raw)]
    str = "UnivariateFinite($x1=>$(round(p1, sigdigits=3))"
    pairs = (d.decoder(r)=>d.prob_given_class[r] for r in raw[2:end])
    for pair in pairs
        str *= ", $(pair[1])=>$(round(pair[2], sigdigits=3))"
    end
    str *= ")"
    print(stream, str)
end

"""
    isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)

Returns `true` if and only if `Set(classes(d1) == Set(classes(d2))` and the
corresponding probabilities are approximately equal. The key-word arguments
`kwargs` are passed through to each call of `isapprox` on probability pairs.
Returns `false` otherwise.
"""
function Base.isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)
    classes1 = classes(d1)
    classes2 = classes(d2)
    for c in classes1
        c in classes2 || return false
        isapprox(pdf(d1, c), pdf(d2, c); kwargs...) ||
            return false # pdf defined below
    end
    return true
end

function average(dvec::AbstractVector{UnivariateFinite{L,U,T}};
                 weights=nothing) where {L,U,T}

    n = length(dvec)

    Dist.@check_args(UnivariateFinite, weights == nothing || n==length(weights))

    # check all distributions have consistent pool:
    first_index = first(dvec).decoder.pool.index
    for d in dvec
        d.decoder.pool.index == first_index ||
            error("Averaging UnivariateFinite distributions with incompatible"*
                  " pools. ")
    end

    # get all refs:
    refs = Tuple(reduce(union, [keys(d.prob_given_class) for d in dvec]))

    # initialize the prob dictionary for the distribution sum:
    prob_given_class = LittleDict{U,T}(refs, zeros(T, length(refs)))

    # make vector of all the distributions dicts padded to have same common keys:
    prob_given_class_vec = map(dvec) do d
        merge(prob_given_class, d.prob_given_class)
    end

    # sum up:
    if weights == nothing
        scale = 1/n
        for x in refs
            for k in 1:n
                prob_given_class[x] += scale*prob_given_class_vec[k][x]
            end
        end
    else
        scale = 1/sum(weights)
        for x in refs
            for k in 1:n
                prob_given_class[x] +=
                    weights[k]*prob_given_class_vec[k][x]*scale
            end
        end
    end

    return UnivariateFinite(first(dvec).decoder, prob_given_class)
end

function _pdf(d::UnivariateFinite{L,U,T}, ref) where {L,U,T}
    return get(d.prob_given_class, ref, zero(T))
end

"""
    Distributions.pdf(d::UnivariateFinite, x)

Probability of `d` at `x`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    pdf(d, "yes")     # 0.3
    pdf(d, v[1])      # 0.3
    pdf(d, "no")      # 0.0
    pdf(d, "house")   # throws error

Other similar methods are available too:

    mode(d)    # CategoricalString{UInt32} "maybe"
    rand(d, 5) # CategoricalArray{String,1,UInt32}["maybe", "no", "maybe", "maybe", "no"] or similar
    d = fit(UnivariateFinite, v)
    pdf(d, "maybe") # 0.25

One can also do weighted fits:

    w = [1, 4, 5, 1] # some weights
    d = fit(UnivariateFinite, v, w)
    pdf(d, "maybe") ≈ 4/11 # true

*Warning:* The `pdf` function will give wrong answers if the order of
 levels of any categorical element passed to the UnivariateFinite
 constructor is changed.

See also `classes`, `support`.
"""
Distributions.pdf(d::UnivariateFinite, x::CategoricalElement) = _pdf(d, int(x))

function Distributions.pdf(d::UnivariateFinite, level)
    p = d.decoder.pool
    if level in p.levels
        ref = p.order[p.invindex[level]]
        return _pdf(d, ref)
    else
        throw(ArgumentError(""))
    end
end

function Distributions.mode(d::UnivariateFinite)
    dic = d.prob_given_class
    p = values(dic)
    max_prob = maximum(p)
    m = first(first(dic)) # mode, just some ref for now
    for (x, prob) in dic
        if prob == max_prob
            m = x
            break
        end
    end
    return d.decoder(m)
end

"""
    _cumulative(d::UnivariateFinite)

Return the cumulative probability vector `[0, ..., 1]` for the
distribution `d`, using whatever ordering is used in the dictionary
`d.prob_given_class`. Used only for to implement random sampling from
`d`.

"""
function _cumulative(d::UnivariateFinite{L,U,T}) where {L,U,T<:Real}
    p = collect(values(d.prob_given_class))
    K = length(p)
    p_cumulative = Array{T}(undef, K + 1)
    p_cumulative[1] = zero(T)
    p_cumulative[K + 1] = one(T)
    for i in 2:K
        p_cumulative[i] = p_cumulative[i-1] + p[i-1]
    end
    return p_cumulative
end

"""
_rand(p_cumulative)

Randomly sample the distribution with discrete support `1:n` which has
cumulative probability vector `p_cumulative=[0, ..., 1]` (of length
`n+1`). Does not check the first and last elements of `p_cumulative`
but does not use them either.

"""
function _rand(p_cumulative)
    real_sample = rand()
    K = length(p_cumulative)
    index = K
    for i in 2:K
        if real_sample < p_cumulative[i]
            index = i - 1
            break
        end
    end
    return index
end

function Base.rand(d::UnivariateFinite)
    p_cumulative = _cumulative(d)
    classes = d.decoder.pool.valindex[collect(keys(d.prob_given_class))]
    return classes[_rand(p_cumulative)]
end

function Base.rand(d::UnivariateFinite, n::Int)
    p_cumulative = _cumulative(d)
    classes = d.decoder.pool.valindex[collect(keys(d.prob_given_class))]
    return [classes[_rand(p_cumulative)] for i in 1:n]
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{L}) where L
    L <: CategoricalElement ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "CategoricalValue or CategoricalString type. ")
    y = skipmissing(v) |> collect
    isempty(y) && error("No non-missing data to fit. ")
    N = length(y)
    count_given_class = Dist.countmap(y)
    classes = Tuple(keys(count_given_class))
    probs = values(count_given_class)./N
    prob_given_class = LittleDict(classes, probs)
    return UnivariateFinite(prob_given_class)
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{L},
                           weights::Nothing) where L
    return Distributions.fit(d, v)
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{L},
                           weights::AbstractVector{<:Real}) where L
    L <: CategoricalElement ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "CategoricalValue or CategoricalString type. ")
    y = broadcast(identity, skipmissing(v))
    isempty(y) && error("No non-missing data to fit. ")
    classes_seen = filter(in(unique(y)), classes(y[1]))

    # instantiate and initialize prob dictionary:
    prob_given_class = LittleDict{L,Float64}()
    for c in classes_seen
        prob_given_class[c] = 0
    end

    # compute unnormalized  probablilities:
    for i in eachindex(y)
        prob_given_class[y[i]] += weights[i]
    end

    # normalize the probabilities:
    S = sum(values(prob_given_class))
    for c in keys(prob_given_class)
        prob_given_class[c] /=S
    end

    return UnivariateFinite(prob_given_class)
end

## NORMAL DISTRIBUTION ARITHMETIC

# *(lambda::Number, d::Dist.Normal) = Dist.Normal(lambda*d.μ, lambda*d.σ)
# function +(ds::Dist.Normal...)
#     μ = sum([d.μ for d in ds])
#     σ = sum([d.σ^2 for d in ds]) |> sqrt
#     return Dist.Normal(μ, σ)
# end


text:
module ANSEL

using Unicode

# For reference, see https://en.wikipedia.org/wiki/ANSEL

ansel_map = [0x0000 0x0001 0x0002 0x0003 0x0004 0x0005 0x0006 0x0007 0x0008 0x0009 0x000A 0x000B 0x000C 0x000D 0x000E 0x000F
             0x0010 0x0011 0x0012 0x0013 0x0014 0x0015 0x0016 0x0017 0x0018 0x0019 0x001A 0x001B 0x001C 0x001D 0x001E 0x001F
             0x0020 0x0021 0x0022 0x0023 0x0024 0x0025 0x0026 0x0027 0x0028 0x0029 0x002A 0x002B 0x002C 0x002D 0x002E 0x002F
             0x0030 0x0031 0x0032 0x0033 0x0034 0x0035 0x0036 0x0037 0x0038 0x0039 0x003A 0x003B 0x003C 0x003D 0x003E 0x003F
             0x0040 0x0041 0x0042 0x0043 0x0044 0x0045 0x0046 0x0047 0x0048 0x0049 0x004A 0x004B 0x004C 0x004D 0x004E 0x004F
             0x0050 0x0051 0x0052 0x0053 0x0054 0x0055 0x0056 0x0057 0x0058 0x0059 0x005A 0x005B 0x005C 0x005D 0x005E 0x005F
             0x0060 0x0061 0x0062 0x0063 0x0064 0x0065 0x0066 0x0067 0x0068 0x0069 0x006A 0x006B 0x006C 0x006D 0x006E 0x006F
             0x0070 0x0071 0x0072 0x0073 0x0074 0x0075 0x0076 0x0077 0x0078 0x0079 0x007A 0x007B 0x007C 0x007D 0x007E 0x007F
             0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000
             0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000
             0x0000 0x0141 0x00D8 0x0110 0x00DE 0x00C6 0x0152 0x02B9 0x00B7 0x266D 0x00AE 0x00B1 0x01A0 0x01AF 0x02BC 0x0000
             0x02BB 0x0142 0x00F8 0x0111 0x00FE 0x00E6 0x0153 0x02BA 0x0131 0x00A3 0x00F0 0x0000 0x01A1 0x01B0 0x0000 0x0000
             0x00B0 0x2113 0x2117 0x00A9 0x266F 0x00BF 0x00A1 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000
             0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000
             0x0309 0x0300 0x0301 0x0302 0x0303 0x0304 0x0306 0x0307 0x0308 0x030C 0x030A 0xFE20 0xFE21 0x0315 0x030B 0x0310
             0x0327 0x0328 0x0323 0x0324 0x0325 0x0333 0x0332 0x0326 0x031C 0x032E 0xFE22 0xFE23 0x0000 0x0000 0x0313 0x0000]

gedcom_map = Dict(0xbe => 0x25a1,
                  0xbf => 0x25a0,
                  0xcd => 0x0065,
                  0xce => 0x006f,
                  0xcf => 0x00df,
                  0xfc => 0x0338)

function ansel2unicode(b::UInt8)
    i = Integer(b >> 4) + 1
    j = Integer(b & 0xf) + 1
    ansel_map[i,j]
end

gedcom2unicode(b::UInt8) = get(gedcom_map, b, ansel2unicode(b))

function anselgedcom2unicode(data::Vector{UInt8}, t::Function)
    utf8data = t.(data)
    # In the ANSEL specification, diacritic marks precede the
    # decorated letters, whereas in Unicode they follow the decorated
    # letters instead.
    for i ∈ findall(b -> b & 0xf0 ≥ 0xe0 , data)
        utf8data[[i,i+1]] = utf8data[[i+1,i]]
    end
    Unicode.normalize(transcode(String, utf8data), :NFC)
end

ansel2unicode(data) = anselgedcom2unicode(data, ansel2unicode)
gedcom2unicode(data) = anselgedcom2unicode(data, gedcom2unicode)

export ansel2unicode, gedcom2unicode

end # module


text:
using SHA  # security instincts say do not write bare passwords to a shared file even in toy code :)

mutable struct Personnel
    fullname::String
    office::String
    extension::String
    homephone::String
    email::String
    Personnel(ful,off,ext,hom,ema) = new(ful,off,ext,hom,ema)
end

mutable struct Passwd
     account::String
     password::String
     uid::Int32
     gid::Int32
     personal::Personnel
     directory::String
     shell::String
     Passwd(acc,pas,uid,gid,per,dir,she) =  new(acc,pas,uid,gid,per,dir,she)
end

function writepasswd(filename, passrecords)
    if(passrecords isa Array) == false
        passrecords = [passrecords]
    end
    fh = open(filename, "a") # should throw an exception if cannot open in a locked or exclusive mode for append
    for pas in passrecords
        record = join([pas.account, bytes2hex(sha256(pas.password)), pas.uid, pas.gid,
                 join([pas.personal.fullname, pas.personal.office, pas.personal.extension,
                 pas.personal.homephone, pas.personal.email], ','),
                 pas.directory, pas.shell], ':')
        write(fh, record, "\n")
    end
    close(fh)
end

const jsmith = Passwd("jsmith","x",1001, 1000, Personnel("Joe Smith", "Room 1007", "(234)555-8917", "(234)555-0077", "jsmith@rosettacode.org"), "/home/jsmith", "/bin/bash")
const jdoe = Passwd("jdoe","x",1002, 1000, Personnel("Jane Doe", "Room 1004", "(234)555-8914", "(234)555-0044", "jdoe@rosettacode.org"), "/home/jdoe", "/bin/bash")
const xyz = Passwd("xyz","x",1003, 1000, Personnel("X Yz", "Room 1003", "(234)555-8913", "(234)555-0033", "xyz@rosettacode.org"), "/home/xyz", "/bin/bash")

const pfile = "pfile.csv"
writepasswd(pfile, [jsmith, jdoe])
println("Before last record added, file is:\n$(readstring(pfile))")
writepasswd(pfile, xyz)
println("After last record added, file is:\n$(readstring(pfile))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JuMP
using CPLEX
include("parser.jl")
include("relaxed.jl")
include("base.jl")

export DW

function DW(fileName :: String, X_i::Array{Float64,2},  X_fi::Array{Float64,3}, X_ikf::Array{Float64,4}, E::Array{Float64,5})

	#Data acquisition, model definition
	cout_ouverture, Fct_commod, func_cost, func_capacity, nb_nodes, nb_arcs, nb_commodities, latency, node_capacity, commodity, nb_func, exclusion = read_instance(fileName)

	#state sub-problem definition
	m0 = Model(CPLEX.Optimizer)
	set_silent(m0)
	@variable(m0, x_i[1:nb_nodes] >= 0)
	@variable(m0, x_fi[1:nb_func,1:nb_nodes]>=0)
	@constraint(m0,[i in 1:nb_nodes], x_i[i] <= 1) # Relaxed binary constraint
	@constraint(m0,[i in 1:nb_nodes], sum(x_fi[f,i] for f in 1:nb_func) <= node_capacity[i], base_name = "cap" ) # (4) capacité de noeud

	#flow sub problem definition

	m_k=Dict{Int,JuMP.Model}()
	m_k[0]=m0

	for k in 1:nb_commodities

		#model definition
		mk = Model(CPLEX.Optimizer)
		set_silent(mk)
		@variable(mk, x_ikf[1:nb_nodes,1:nb_func]>=0, base_name="x_ikf_"*string(k))
		@variable(mk, e[1:nb_nodes,1:nb_nodes,1:nb_func+1] >= 0,base_name="e_"*string(k))

		for i in 1:nb_nodes
			@constraint(mk,[f in 1:nb_func], x_ikf[i,f] <= 1) # Relaxed binary constraint
			for j in 1:nb_nodes
				@constraint(mk,[f in 1:nb_func+1], e[i,j,f] <= 1) # Relaxed binary constraint
			end
		end

		#Constraints of exclusion and possible arc
		for i in 1:nb_nodes
			if exclusion[k,1]>0
				@constraint(mk,sum(x_ikf[i,w] for w in exclusion[k,:]) <= 1, base_name = "exclusion_"*string(k)*"_"*string(i)) #exclusion
			end
			for j in 1:nb_nodes
				if latency[i,j] == 0
					@constraint(mk,[f in 1:nb_func+1],e[i,j,f]==0, base_name = "no_arc_"*string(i)*"_"*string(j)) #absence d'arcs
				end
			end
		end

		#Contrainte de flots successifs
		size_fk=length( findall( y -> y > 0, Fct_commod[k,:]))
		tab_fk=sortperm( Fct_commod[k,:])[1:size_fk]
		@constraint(mk,sum(e[commodity[k,1],j,tab_fk[1]] - e[j,commodity[k,1],tab_fk[1]] for j in 1:nb_nodes)-x_ikf[commodity[k,1],tab_fk[1]]==1, base_name = "init_flot_"*string(k)) #initialisation du flot à la source
		@constraint(mk,[i in 1:nb_nodes ; i!=commodity[k,1]],sum(e[i,j,tab_fk[1]] - e[j,i,tab_fk[1]]  for j in 1:nb_nodes)+x_ikf[i,tab_fk[1]]==0, base_name = "1er_flot_"*string(k)) #flot pour traiter la première fonction de puis la source
		for t in 2:size_fk
			@constraint(mk,[i in 1:nb_nodes],sum(e[i,j,tab_fk[t]] - e[j,i,tab_fk[t]] for j in 1:nb_nodes)-x_ikf[i,tab_fk[t-1]]+x_ikf[i,tab_fk[t]]==0, base_name = "mid_flow_"*string(t)*"_"*string(k)) #flot d'un traitement au suivant
		end
		@constraint(mk,[i in 1:nb_nodes ; i!=commodity[k,2]],sum(e[i,j,nb_func+1] - e[j,i,nb_func+1] for j in 1:nb_nodes)-x_ikf[i,tab_fk[end]]==0, base_name = "last_flow_"*string(k)) #flot du dernier traitement vers le puit
		@constraint(mk,sum(e[j,commodity[k,2],end] for j in 1:nb_nodes)+x_ikf[commodity[k,2],tab_fk[end]]==1, base_name = "end_flow_"*string(k)) #fin de flot sur le puit

		for f in 1:nb_func
			@constraint(mk,sum(x_ikf[i,f] for i in 1:nb_nodes)<=1, base_name = "concatenation_"*string(k)*"_"*string(f)) #concatenation des problemes de flot
		end

		#Contraintes de latence
		@constraint(mk, sum( sum( sum(e[i,j,f]*latency[i,j] for f in 1:nb_func+1) for j in 1:nb_nodes) for i in 1:nb_nodes) <= commodity[k,4], base_name = "latence_"*string(k))

		m_k[k]=mk
	end

	#Column generation
	not_opti=true
	n=0
	while not_opti

		#Probleme maitre
		m = Model(CPLEX.Optimizer)
		set_silent(m)
		n=size(X_ikf,1)
		@variable(m, lbd[1:n] >= 0)
		@constraint(m,sum(lbd[w] for w in 1:n)>=1, base_name="convex") #Convexity constraint
		for i in 1:nb_nodes
			@constraint(m,[f in 1:nb_func],sum(lbd[w]*X_fi[w,f,i] for w in 1:n)*func_capacity[f] >=  sum( sum(lbd[w]*X_ikf[w,i,k,f] for w in 1:n) * commodity[k,3] for k in 1:nb_commodities), base_name = "n_"*string(i)) # (3) Nombre de fonctions f à placer en i
			for k in 1:nb_commodities
				@constraint(m,[f in 1:nb_func],sum(lbd[w]*X_i[w,i] for w in 1:n) >= sum(lbd[w]*X_ikf[w,i,k,f] for w in 1:n), base_name = "ouverture_"*string(i)*"_"*string(k)) # (2) ouverture du noeud en i
			end
		end
		@objective(m, Min, sum(sum(lbd[w]*X_i[w,i] for w in 1:n)*cout_ouverture + sum(sum(lbd[w]*X_fi[w,f,i] for w in 1:n)*func_cost[f,i] for f in 1:nb_func) for i in 1:nb_nodes) )
		optimize!(m)

		#Get dual values and solutions
		#println(JuMP.value.(lbd))
		eta=dual(constraint_by_name(m, "convex"))
		mu_2=Array{Float64,3}(zeros(nb_nodes,nb_commodities,nb_func))
		mu_3=Array{Float64,2}(zeros(nb_nodes,nb_func))
		for i in 1:nb_nodes
			for f in 1:nb_func
				mu_3[i,f]=dual(constraint_by_name(m,"n_"*string(i)*"["*string(f)*"]"))
				for k in 1:nb_commodities
					mu_2[i,k,f]=dual(constraint_by_name(m,"ouverture_"*string(i)*"_"*string(k)*"["*string(f)*"]"))
				end
			end
		end

		#setup for next column to add on first index
		min_cr=-eta #Minimum ecart complementaire
		X_i=vcat(X_i,X_i[1:1,:])
		X_fi=vcat(X_fi,X_fi[1:1,:,:])
		X_ikf=vcat(X_ikf,X_ikf[1:1,:,:,:])
		E=vcat(E,E[1:1,:,:,:,:])

		#state sub-problem resolution
		@objective(m_k[0], Min, sum(x_i[i]*(cout_ouverture-sum(sum(mu_2[i,k,f] for f in 1:nb_func) for k in 1:nb_commodities))+sum(x_fi[f,i]*(func_cost[f,i]-mu_3[i,f]*func_capacity[f]) for f in 1:nb_func)  for i in 1:nb_nodes) )
		optimize!(m_k[0])
		min_cr=min_cr+objective_value(m_k[0])
		X_i[end,:]=JuMP.value.(x_i)
		X_fi[end,:,:]=JuMP.value.(x_fi)

		#flow sub-problem resolution
		for k in 1:nb_commodities
			x_ikf=object_dictionary(m_k[k])[:x_ikf]
			e=object_dictionary(m_k[k])[:e]
			@objective(m_k[k], Min, - sum(sum( (-mu_3[i,f]*commodity[k,3] - mu_2[i,k,f])*x_ikf[i,f] for i in 1:nb_nodes) for f in 1:nb_func) )
			optimize!(m_k[k])
			min_cr = min_cr + objective_value(m_k[k])
			X_ikf[end,:,k,:]=JuMP.value.(x_ikf)
			E[end,:,:,k,:]=JuMP.value.(e)
		end

		#Check optimality
		if min_cr>=-1e-7
			not_opti=false
		end
	end

	#Last solve of model to get results
	m = Model(CPLEX.Optimizer)
	n=size(X_ikf,1)
	@variable(m, lbd[1:n] >= 0)
	@constraint(m, lbd[1] == 0)

	@constraint(m,sum(lbd[w] for w in 1:n)==1, base_name="convex")  #Convexity constraint
	for i in 1:nb_nodes
		@constraint(m,[f in 1:nb_func],sum(lbd[w]*X_fi[w,f,i] for w in 1:n)*func_capacity[f] >=  sum( sum(lbd[w]*X_ikf[w,i,k,f] for w in 1:n) * commodity[k,3] for k in 1:nb_commodities), base_name = "n_"*string(i)) # (3) Nombre de fonctions f à placer en i
		for k in 1:nb_commodities
			@constraint(m,[f in 1:nb_func],sum(lbd[w]*X_i[w,i] for w in 1:n) >= sum(lbd[w]*X_ikf[w,i,k,f] for w in 1:n), base_name = "ouverture_"*string(i)*"_"*string(k)) # (2) ouverture du noeud en i
		end
	end
	@objective(m, Min, sum(sum(lbd[w]*X_i[w,i] for w in 1:n)*cout_ouverture + sum(sum(lbd[w]*X_fi[w,f,i] for w in 1:n)*func_cost[f,i] for f in 1:nb_func) for i in 1:nb_nodes) )
	optimize!(m)
	println(JuMP.value.(lbd))
	mu_2_sol=Array{Float64,3}(zeros(nb_nodes,nb_commodities,nb_func))
	mu_3_sol=Array{Float64,2}(zeros(nb_nodes,nb_func))
	for i in 1:nb_nodes
		for f in 1:nb_func
			mu_3_sol[i,f]=dual(constraint_by_name(m,"n_"*string(i)*"["*string(f)*"]"))
			for k in 1:nb_commodities
				mu_2_sol[i,k,f]=dual(constraint_by_name(m,"ouverture_"*string(i)*"_"*string(k)*"["*string(f)*"]"))
			end
		end
	end
	println(solution_summary(m))

	#with dual
	m = Model(CPLEX.Optimizer)
	set_silent(m)
	@variable(m,mu_2[1:nb_nodes,1:nb_commodities,1:nb_func]>=0)
	@variable(m,mu_3[1:nb_nodes,1:nb_func]>=0)
	#@constraint(m, [ i in 1:nb_nodes*nb_func], mu_3[i] == mu_3_sol[i])
	#@constraint(m, [ i in 1:nb_nodes*nb_func*nb_commodities], mu_2[i] == mu_2_sol[i])
	@constraint(m, [ i in 1:nb_nodes*nb_func], mu_3[i] <= 30000)
	@constraint(m, [ i in 1:nb_nodes*nb_func*nb_commodities], mu_2[i] <= 30000)

	@variable(m,eta)
	for idx in 2:size(X_ikf,1)
		@constraint(m, sum(cout_ouverture*X_i[idx,i]+sum(func_cost[f,i]*X_fi[idx,f,i] for f in 1:nb_func) for i in 1:nb_nodes) - sum( X_i[idx,i]*sum(sum(mu_2[i,k,f] for k in 1:nb_commodities) for f in 1:nb_func) + sum(func_capacity[f]*mu_3[i,f]*X_fi[idx,f,i] for f in 1:nb_func) for i in 1:nb_nodes) - sum(sum(sum( (-commodity[k,3]*mu_3[i,f]-mu_2[i,k,f])*X_ikf[idx,i,k,f] for f in 1:nb_func) for k in 1:nb_commodities) for i in 1:nb_nodes) - eta >=0, base_name="con_"*string(idx))
		#@constraint(m, sum(X_i[idx,i]*(cout_ouverture-sum(sum(mu_2[i,k,f] for f in 1:nb_func) for k in 1:nb_commodities))+sum(X_fi[idx,f,i]*(func_cost[f,i]-mu_3[i,f]*func_capacity[f]) for f in 1:nb_func)  for i in 1:nb_nodes) - sum( sum(sum( (-mu_3[i,f]*commodity[k,3] - mu_2[i,k,f])*X_ikf[idx,i,k,f] for i in 1:nb_nodes) for f in 1:nb_func) for k in 1:nb_commodities) - eta >=0, base_name="con_"*string(idx))
	end
	@objective(m, Max, eta )
	optimize!(m)
	for idx in 2:size(X_ikf,1)
		print(dual(constraint_by_name(m,"con_"*string(idx))))
		print(" ")
	end
	println(" ")

	#Results
	println(solution_summary(m))
	lbd_sol=JuMP.value.(lbd)
	x_i=sum(lbd_sol[w].*X_i[w,:] for w in 1:n)
	x_fi=sum(lbd_sol[w].*X_fi[w,:,:] for w in 1:n)
	x_ikf=sum(lbd_sol[w].*X_ikf[w,:,:,:] for w in 1:n)
	e=sum(lbd_sol[w].*E[w,:,:,:,:] for w in 1:n)

	return x_i, x_fi, x_ikf, e
end

io, x_i, x_fi, x_ikf, e = Static("di-yuan/di-yuan_1/",true)
X_i=Array{Float64,2}(zeros(1,size(x_i,1)))
X_i[1,:]=x_i
X_fi=Array{Float64,3}(zeros(1,size(x_fi,1),size(x_fi,2)))
X_fi[1,:,:]=x_fi
X_ikf=Array{Float64,4}(zeros(1,size(x_ikf,1),size(x_ikf,2),size(x_ikf,3)))
X_ikf[1,:,:,:]=x_ikf
E=Array{Float64,5}(zeros(1,size(e,1),size(e,2),size(e,3),size(e,4)))
E[1,:,:,:,:]=e

x_i, x_fi, x_ikf, e=DW("di-yuan/di-yuan_1/", X_i,X_fi,X_ikf,E)
println("fin")
#x_i, x_fi, x_ikf, e=Relaxed("di-yuan/di-yuan_1/")

max_stars_repo_name: VSpitzer/RORT

max_issues_repo_name: VSpitzer/RORT

max_forks_repo_name: VSpitzer/RORT


text:
module BenchSum

using BenchmarkTools
using Transducers

suite = BenchmarkGroup()

function slow_computation(n)
    x = 0.4
    for _ in 1:n
        x = 3 / π * sin(x * π)
    end
    return x
end

xf = Map(slow_computation)

n = 2^14
maxworkload = 100
inputtable = [
    ("uniform", fill(maxworkload ÷ 2, n)),
    ("random", rand(1:maxworkload, n)),
    #=
    ("increasing", ceil.(Int, range(1, maxworkload, length = n))),
    ("decreasing", ceil.(Int, range(maxworkload, 1, length = n))),
    (
        "mountain",
        vcat(
            ceil.(Int, range(1, maxworkload, length = n ÷ 2)),
            ceil.(Int, range(maxworkload, 1, length = n ÷ 2)),
        ),
    ),
    =#
    (
        "valley",
        vcat(
            ceil.(Int, range(maxworkload, 1, length = n ÷ 2)),
            ceil.(Int, range(1, maxworkload, length = n ÷ 2)),
        ),
    ),
]

for (label, xs) in inputtable
    s1 = suite[label] = BenchmarkGroup()
    s1["foldl"] = @benchmarkable foldl(+, $xf, $xs; simd = true)
    s2 = s1["reduce"] = BenchmarkGroup()
    for basesize in 2 .^ (7:9)
        s2["basesize=$basesize"] = @benchmarkable reduce(
            +,
            $xf,
            $xs;
            basesize = $basesize,
            simd = true,
        )
    end
end

end  # module
BenchSum.suite


text:
# ---
# title: 1013. Partition Array Into Three Parts With Equal Sum
# id: problem1013
# author: Indigo
# date: 2021-01-30
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/partition-array-into-three-parts-with-equal-sum/description/>
# hidden: true
# ---
#
# Given an array `A` of integers, return `true` if and only if we can partition
# the array into three **non-empty** parts with equal sums.
#
# Formally, we can partition the array if we can find indexes `i+1 < j` with
# `(A[0] + A[1] + ... + A[i] == A[i+1] + A[i+2] + ... + A[j-1] == A[j] + A[j-1]
# + ... + A[A.length - 1])`
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [0,2,1,-6,6,-7,9,1,2,0,1]
#     Output: true
#     Explanation: 0 + 2 + 1 = -6 + 6 - 7 + 9 + 1 = 2 + 0 + 1
#
#
# **Example 2:**
#
#
#
#     Input: A = [0,2,1,-6,6,7,9,-1,2,0,1]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: A = [3,3,6,5,-2,2,5,1,-9,4]
#     Output: true
#     Explanation: 3 + 3 = 6 = 5 - 2 + 2 + 5 + 1 - 9 + 4
#
#
#
#
# **Constraints:**
#
#   * `3 <= A.length <= 50000`
#   * `-10^4 <= A[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

function can_three_parts_equal_sum(arr::Vector{Int})
    function can_three_parts_equal_sum(arr::AbstractVector{Int}, n::Int)
        targ = sum(arr)
        if targ % n != 0 || length(arr) == 0
            return false
        elseif n == 1
            return true
        end
        targ ÷= n
        acc = 0
        for i in 1:length(arr)
            acc += arr[i]
            targ == acc && return can_three_parts_equal_sum(@view(arr[(i + 1):end]), n - 1)
        end
        return false
    end
    return can_three_parts_equal_sum(arr, 3)
end
## @lc code=end


text:
# Test chordal extension
workspace()
include("qdldl.jl")
include("/Users/Micha/Dropbox/Research/OSSDP/Code/src/Graph.jl")

using QDLDL,GraphModule, Base.Test

rng = Random.MersenneTwister(131123)

numProblems = 10

# g = Graph(A,false)
# println("Original Graph is chordal?: $(GraphModule.isPerfectOrdering(g))")
# # A = A  +(2*5+1)*speye(5);
# F = QDLDL.qdldl(A,logical=true,perm=nothing)
# L = F.L

# C = full(L+L') + eye(5)

# g_c = Graph(sparse(C))
# println("Extended Graph is chordal?: $(GraphModule.isPerfectOrdering(g_c))")

@testset "QDLDL routine for chordal embedding" begin

  for iii=1:numProblems

    dim = rand(rng,50:1:200);
    density=rand(rng,0.1:0.1:0.4)
    A = sprand(rng,dim,dim,density)

    A = 0.5*(A + A')
    A = A + (2*dim+1)*speye(dim);

    g = Graph(A,false)

    if (GraphModule.isPerfectOrdering(g) == true)
      continue
    else
      F = QDLDL.qdldl(A,logical=true,perm=nothing)
      L = F.L
      C = L+L' + speye(dim)
      g_c = Graph(C)
    end

    # check results
    @test GraphModule.isPerfectOrdering(g_c)

  end
end


text:
import Pkg

@inline function _get_project_recursive_dependencies()
    result = Vector{String}(undef, 0)
    for pkginfo in values(Pkg.dependencies())
        push!(result, pkginfo.name)
    end
    append!(result, collect(values(Pkg.Types.stdlibs())))
    unique!(result)
    return result
end

@inline function _add_all_project_recursive_dependencies()
    pkg_names_to_add = Vector{String}(undef, 0)
    for pkginfo in values(Pkg.dependencies())
        push!(pkg_names_to_add, pkginfo.name)
    end
    append!(pkg_names_to_add, collect(values(Pkg.Types.stdlibs())))
    unique!(pkg_names_to_add)
    Pkg.add(pkg_names_to_add)
    return nothing
end


text:
using DynamicHMCModels

ProjDir = @__DIR__

df = DataFrame(CSV.read(joinpath(ProjDir, "..", "..", "data",  "Kline.csv"), delim=';'))

# New col logpop, set log() for population data
df[!, :logpop] = map((x) -> log(x), df[:, :population]);
df[!, :society] = 1:10;

Base.@kwdef mutable struct KlineModel{Ty <: AbstractVector,
  Tx <: AbstractMatrix, Ts <: AbstractVector}
    "Observations (total_tools)."
    y::Ty
    "Covariates (logpop)"
    x::Tx
    "Society"
    s::Ts
    "Number of observations (10)"
    N::Int
    "Number of societies (also 10)"
    N_societies::Int
end

function make_transformation(model::KlineModel)
    as( (β = as(Array, size(model.x, 2)), α = as(Array, model.N_societies), σ = asℝ₊) )
end

# Instantiate the model with data and inits.

N = size(df, 1)
N_societies = length(unique(df[:, :society]))
x = hcat(ones(Int64, N), df[:, :logpop]);
s = df[:, :society]
y = df[:, :total_tools]
model = KlineModel(; y=y, x=x, s=s, N=N, N_societies=N_societies)

# Make the type callable with the parameters *as a single argument*.

function (model::KlineModel)(θ)
    @unpack y, x, s, N, N_societies = model   # data
    @unpack β, α, σ = θ  # parameters
    ll = 0.0
    ll += logpdf(Cauchy(0, 1), σ)
    ll += sum(logpdf.(Normal(0, σ), α)) # α[1:10]
    ll += logpdf.(Normal(0, 10), β[1]) # a
    ll += logpdf.(Normal(0, 1), β[2]) # a
    ll += sum(
      [loglikelihood(Poisson(exp(α[s[i]] + dot(x[i, :], β))), [y[i]]) for i in 1:N]
    )
    ll
end

println()
θ = (β = [1.0, 0.25], α = rand(Normal(0, 1), N_societies), σ = 0.2)
model(θ) |> display
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display

println()
DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

# Set varable names

parameter_names = ["a", "bp", "sigma_society"]
pooled_parameter_names = ["a_society[$i]" for i in 1:10]

# Create a3d

a3d = Array{Float64, 3}(undef, 1000, 13, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[i].β)
    a3d[i, 3, j] = values(posterior[i].σ)
    a3d[i, 4:13, j] = values(posterior[i].α)
  end
end

chns = MCMCChains.Chains(a3d,
  vcat(parameter_names, pooled_parameter_names),
  Dict(
    :parameters => parameter_names,
    :pooled => pooled_parameter_names
  )
);

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
                            Mean                SD               Naive SE             MCSE            ESS
            a          1.076167468  0.7704872560 0.01218247319 0.0210530022 1000.000000
           bp         0.263056273  0.0823415805 0.00130193470 0.0022645077 1000.000000
  a_society.1   -0.191723568  0.2421382537 0.00382854195 0.0060563054 1000.000000
  a_society.2    0.054569029  0.2278506876 0.00360263570 0.0051693148 1000.000000
  a_society.3   -0.035935050  0.1926364647 0.00304584994 0.0039948433 1000.000000
  a_society.4    0.334355037  0.1929971201 0.00305155241 0.0063871707  913.029080
  a_society.5    0.049747513  0.1801287716 0.00284808595 0.0043631095 1000.000000
  a_society.6   -0.311903245  0.2096126337 0.00331426674 0.0053000536 1000.000000
  a_society.7    0.148637507  0.1744680594 0.00275858223 0.0047660246 1000.000000
  a_society.8   -0.164567976  0.1821341074 0.00287979309 0.0034297298 1000.000000
  a_society.9    0.277066965  0.1758237250 0.00278001719 0.0055844175  991.286501
 a_society.10   -0.094149204  0.2846206232 0.00450024719 0.0080735022 1000.000000
sigma_society    0.310352849  0.1374834682 0.00217380450 0.0057325226  575.187461
";

# Describe the chain

describe(chns) |> display
println()

# Describe the chain

describe(chns, sections=[:pooled])


text:
using BinaryBuilder

# This package builds debian glibc 2.33 for use as a replacement glibc in
# debian bookworm in order to address the lse initialization order issue.
# It is hopefully temporary until the upstream bug is fixed.

include("../common.jl")

name = "DebianGlibc"
version = v"2.33"

sources = glibc_sources(version)
push!(sources,
    ArchiveSource("http://deb.debian.org/debian/pool/main/g/glibc/glibc_2.33-7.debian.tar.xz",
    "c4dd741a08918861796cf62a8d3209a56d41d807344a7db481b73a4d17f726bc"),)
script = glibc_script()
platforms =  filter(p -> libc(p) == "glibc" && arch(p) == "aarch64",
    supported_platforms())
products = Product[] # glibc_products()
dependencies = glibc_dependencies()

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"11", lock_microarchitecture=false)


text:
module mirone_toy_

export
	mirone_toy

using IUP
using IUP_IM
using IUP_CD
using GMT
using Compat

const NEWPOINT = 0
const MOVE = 1
const CLOSE = 2
const CENTER = 3
const REPAINT = 4

type Handles
	figure1::Ptr{Ihandle}
	iup_canvas::Ptr{Ihandle}
	cd_canvas::Ptr{cdCanvas}	# cdCanvas is a composite type
end

type LineSegPix
	x1::Cint
	x2::Cint
	y1::Cint
	y2::Cint
	active::Bool
end
#function LineSegPix(x1, x2, y1, y2, active)
#	x1 = 0; x2 = 0; y1 = 0; y2 = 0; active = false
#end

include("../src/handlegraphics.jl")

global poly_lastwhat = CLOSE

# ----------------------------------------------------------------------------------
function mirone_toy()
	IupOpen()	#Initializes IUP

	imlabCreateButtonImages()

 	mainMenu = CreateMainMenu()

	# Associates handle "mainMenu" with mainMenu
	IupSetHandle("mainMenu", mainMenu)

	btNew = IupButton("", "imImgWinNew" )
	IupSetAttribute(btNew,"TIP","Creates a new image.");
	IupSetAttribute(btNew,"IMAGE","imImgWinNewButton");

	btOpen = IupButton("", "imImgWinOpen" )
	IupSetAttribute(btOpen,"TIP","Loads an image file from disk.")
	IupSetAttribute(btOpen,"IMAGE","imImgWinOpenButton")

	btLine = IupButton("", "imImgWinLine" )
	IupSetAttribute(btLine,"TIP","Draw line.")
	IupSetAttribute(btLine,"IMAGE","imImgWinLineButton")

	btCloseAll = IupButton("", "imImgWinCloseAll" );
	IupSetAttribute(btCloseAll,"TIP","Close all windows.");
	IupSetAttribute(btCloseAll,"IMAGE","imImgWinCloseAllButton");
	IupSetAttribute(btCloseAll, "ACTIVE", "NO");

	hbToolBar = IupHbox(
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),	# Blank space before first icon
		btNew,
		btOpen,
		btLine,
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),
		#IupSetAttributes(IupFill(),"EXPAND=YES"),
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO")
	)
	IupSetAttribute(hbToolBar, "btCloseAll", btCloseAll);
	IupSetAttribute(hbToolBar, "FLAT", "YES");

	cnvStatus = IupCanvas();
	IupSetAttribute(cnvStatus, "EXPAND", "HORIZONTAL");
	IupSetAttribute(cnvStatus, "RASTERSIZE", "1x25");
	IupSetAttribute(cnvStatus, "BORDER", "YES");

	iup_canvas = blabla()

	vbMainDesktop = IupVbox(
		IupSetAttributes(IupFill(),"SIZE=2, EXPAND=NO"),
		hbToolBar,
		IupSetAttributes(IupFill(),"SIZE=4, EXPAND=NO"),
		IupHbox(iup_canvas),
		IupSetAttributes(IupFill(),"SIZE=4, EXPAND=NO"),
		cnvStatus,
		IupSetAttributes(IupFill(),"SIZE=2, EXPAND=NO")
	)

	# Creates a dialog
	main_window = IupDialog(vbMainDesktop)

	handles = Handles(C_NULL, C_NULL, C_NULL)
	handles.figure1 = main_window
	handles.iup_canvas = iup_canvas
	guidata(main_window, handles)

	IupSetAttribute(main_window, "MENU","mainMenu")
	#IupSetAttribute(main_window, "CLOSE_CB","imImgWinMainExit");
	#IupSetAttribute(main_window, "DROPFILES_CB","imImgWinDropFile");
	#IupSetAttribute(main_window, "GETFOCUS_CB", "imImgWinGetFocus");

	IupSetAttribute(main_window,"MainToolBar", hbToolBar);

	#IupSetAttribute(main_window,"PLACEMENT", "MAXIMIZED");
	IupSetAttribute(main_window,"BRINGFRONT", "YES");
	IupSetAttribute(main_window, "K_cN","imImgWinWindowNext");
	IupSetAttribute(main_window, "K_cO","imImgWinOpen");

	# Sets title and size of the dialog and associates a menu to it
	IupSetAttributes(main_window, "TITLE=\"Mirone toy\", RASTERSIZE = 730x91, MENU = mainMenu");

	# associa os callbacks
	setcallbacks(handles)

	IupShowXY(main_window, IUP_CENTER, IUP_TOP) # Displays the dialog (at screen center)
	IupMainLoop()                               # Initializes IUP main loop
	IupClose()                                  # And close it when ready

end

# ---------------------------------------------------------------------------------------------
function blabla()
	iup_canvas = IupCanvas()
	IupSetAttribute(iup_canvas, "BORDER", "NO");
	IupSetCallback(iup_canvas, "ACTION", cfunction(cbCanvasRepaint, Int, (Ptr{Ihandle},)))
	IupSetCallback(iup_canvas, "MAP_CB", cfunction(cbCanvasMap, Int, (Ptr{Ihandle},)))
	return iup_canvas
end

# --------------------------------------------------------------------------------
function cbCanvasMap(iup_canvas::Ptr{Ihandle})
	handles = guidata(iup_canvas)
	cd_canvas = cdCreateCanvas(cdContextIup(), iup_canvas)
	#cd_canvas = cdCreateCanvas(cdContextDBuffer(), iup_canvas)
	IupSetAttribute(IupGetDialog(iup_canvas), "cdCanvas", cd_canvas)
	IupSetAttribute(iup_canvas, "cdCanvas", cd_canvas)			# Redundant with handles
	handles.cd_canvas = cd_canvas
	guidata(handles.figure1, handles)
	return IUP_DEFAULT
end

# ---------------------------------------------------------------------------------------------
function setcallbacks(handles::Handles)
	IupSetFunction("imImgWinOpen",  cfunction(cmOpen, Int, (Ptr{Ihandle},)))
	IupSetFunction("imImgWinLine",  cfunction(cmdDrawPoly, Int, (Ptr{Ihandle},)))

	#IupSetFunction("cmdRepaint", (Icallback) fRepaint);
	WindowButtonDownFcn(handles.iup_canvas, fButtonCB)
	WindowButtonMotionFcn(handles.iup_canvas, fMotionCB)
	#IupSetFunction("cmdResizeCB", (Icallback) fResizeCB);
end

# ------------------------------------------------------------------------
function fMotionCB(hand::Ptr{Ihandle}, x::Int32, y::Int32, r::Ptr{Uint8})
	handles = guidata(hand)
	yp = convert(Ptr{Cint}, pointer([y]))
	cdCanvasUpdateYAxis(handles.cd_canvas, yp)
	y = unsafe_load(yp)
	polygon(handles, MOVE, x, y)
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
# Function associated to the mouse buttons
#-------------------------------------------------------------------------
function fButtonCB(hand::Ptr{Ihandle}, b::Char, e::Integer, x::Integer, y::Integer, r::Ptr{Uint8})
	handles = guidata(hand)

	cdCanvasActivate(handles.cd_canvas)
	#cdCanvasWriteMode(handles.cd_canvas, IUP_CD.CD_NOT_XOR)
	cdCanvasForeground(handles.cd_canvas, convert(Clong, IUP_CD.CD_BLACK))
	cdCanvasLineStyle(handles.cd_canvas, IUP_CD.CD_CONTINUOUS)
	cdCanvasLineWidth(handles.cd_canvas, 1)
	cdCanvasClip(handles.cd_canvas, IUP_CD.CD_CLIPOFF)

	yp = convert(Ptr{Cint}, pointer([y]))
	cdCanvasUpdateYAxis(handles.cd_canvas, yp)
	y = unsafe_load(yp)
	x = int32(x)				# <=========== NAO DEVIA SER Int64

	set_current_point(handles.iup_canvas, [x,y])

	if (b == IUP_BUTTON1)		# Left button
		if (e != 0)				# button was pressed
			line_seg = getappdata(handles.iup_canvas, "lineSeg")
			if ((try isempty(line_seg) end)) return IUP_DEFAULT end
@show(line_seg)
			line_seg.x1 = x
			line_seg.y1 = y
			line_seg.active = true
			setappdata(handles.iup_canvas, "lineSeg", line_seg)
		end
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
# Desenha o poligono em rubber band.
#-------------------------------------------------------------------------
function polygon(handles::Handles, what::Integer, x::Integer, y::Integer)
# I have to have those redeclared here otherwise it errors saying they are not declared ??????????????????
	global poly_lastwhat

	line_seg = getappdata(handles.iup_canvas, "lineSeg")
#@show(line_seg)
	if ((try isempty(line_seg) end) || !line_seg.active) return end
#@show(what)

	if (what == NEWPOINT)
		if (poly_lastwhat != CLOSE)
			cdCanvasLine(line_seg.x1, line_seg.y1, line_seg.x2, line_seg.y2);      # ...apaga a anterior e...
			cdCanvasLine(line_seg.x1, line_seg.y1, x, y);        # desenha a definitiva
		end
		line_seg.x1 = x;                        # novo segmento comeca no...
		line_seg.y1 = y;                        # fim do primeiro
	elseif (what == MOVE)
		if (poly_lastwhat == MOVE)
			cdCanvasLine(handles.cd_canvas, line_seg.x1, line_seg.y1, line_seg.x2, line_seg.y2);      # apaga o segmento velho
		end
		cdCanvasLine(handles.cd_canvas, line_seg.x1, line_seg.y1, x, y);     # desenha o novo
		line_seg.x2 = x;                        # e o novo se...
		line_seg.y2 = y;                        # torna velho
	elseif (what == REPAINT)
		cdCanvasLine(handles.cd_canvas, line_seg.x1, line_seg.y1, line_seg.x2, line_seg.y2);        # recupera o segmento perdido
		return                             # nao modifica lastwhat
	elseif (what == CLOSE)
		if (poly_lastwhat != CLOSE)
			cdLine(handles.cd_canvas, line_seg.x1, line_seg.y1, line_seg.x2, line_seg.y2);        # apaga o ultimo segmento
			# apaga o poligono temporario inteiro
			for (i = 1:ctgc.num_points-1)
				cdCanvasLine(handles.cd_canvas, ctgc.points[i].x, ctgc.points[i].y, ctgc.points[i+1].x, ctgc.points[i+1].y);
			end
		end
	end
	poly_lastwhat = what
	setappdata(handles.iup_canvas, "lineSeg", line_seg)
end

# ---------------------------------------------------------------------------------------------
function cmdDrawPoly(iup_canvas::Ptr{Ihandle})
	handles = guidata(iup_canvas)
	line_seg = LineSegPix(0, 0, 0, 0, true)

	setappdata(handles.iup_canvas, "lineSeg", line_seg)

	cdCanvasBackground(handles.cd_canvas, IUP_CD.CD_WHITE);
	cdCanvasClear(handles.cd_canvas);
	cdCanvasLineWidth(handles.cd_canvas, 3);
	cdCanvasLineStyle(handles.cd_canvas, IUP_CD.CD_CONTINUOUS);
	cdCanvasForeground(handles.cd_canvas, cdEncodeAlpha(IUP_CD.CD_DARK_MAGENTA, char(128)))
	cdCanvasWriteMode(handles.cd_canvas, IUP_CD.CD_NOT_XOR)
	cdCanvasActivate(handles.cd_canvas)
	return IUP_DEFAULT
end

# ---------------------------------------------------------------------------------------------
function cmOpen(iup_canvas::Ptr{Ihandle})
	# Retrieve a file name
	FileName = "*.*" * repeat(" ",253)
	if (IupGetFile(FileName) != 0)
		return IUP_DEFAULT
	end

	ShowImage(FileName, IupGetDialog(iup_canvas))
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function ShowImage(FileName::String, iup_dialog::Ptr{Ihandle})
	image = IupGetAttribute(iup_dialog, "imImage")		# typeof(image) => Ptr{Uint8}
	image = convert(Ptr{imImage}, image)		# If I use Ptr{imImage} it Booms???
	if (image != C_NULL)
		imImageDestroy(image)
	end
	IupSetAttribute(iup_dialog, "imImage")

#
	error = pointer([0])
	image = imFileImageLoadBitmap(FileName, 0, error)
	error = unsafe_load(error)
	if (error != 0)	PrintError(error)	end
	if (image == C_NULL) return end
#
	#image = load_with_gmt(FileName)

	img = unsafe_load(image)		# Need to get access to the composite type, not to its pointer
@show(pointer_to_array(img.palette,256*3))
	handles = guidata(iup_dialog)
	IupSetAttribute(handles.figure1, "CLIENTSIZE", @sprintf("%dx%d", max(img.width,730), img.height+72))
	IupRefresh(handles.figure1)
@show(bytestring(IupGetAttribute(handles.figure1, "RASTERSIZE")))

	IupSetAttribute(iup_dialog, "imImage", image);
	cbCanvasRepaint(iup_dialog)   # we can do this because canvas inherit attributes from the dialog */
	IupShow(iup_dialog)
end

# --------------------------------------------------------------------------------
function cbCanvasRepaint(iup_canvas::Ptr{Ihandle})
	cd_canvas = convert(Ptr{cdCanvas}, IupGetAttribute(iup_canvas, "cdCanvas"))
	if (cd_canvas == C_NULL)
		return IUP_DEFAULT
	end

	cdCanvasActivate(cd_canvas);
	cdCanvasClear(cd_canvas);

	image = convert(Ptr{imImage}, IupGetAttribute(iup_canvas, "imImage"))
	if (image == C_NULL)
		return IUP_DEFAULT
	end

	#handles = guidata(iup_canvas)
	img = unsafe_load(image)		# Need to get access to the composite type, not to its pointer
	x_off = div(730 - img.width, 2)
	imcdCanvasPutImage(cd_canvas, img, x_off, 0, img.width, img.height, 0, 0, 0, 0)

	IupRefresh(iup_canvas)
	cdCanvasFlush(cd_canvas)
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function load_with_gmt(FileName::String)
	I = gmt("read -Ti " * FileName)
#	I = gmt("read C://progs_cygw//GMTdev//gmt5//branches//5.2.0//test//grdimage//gdal//needle.jpg -Ti")
	n = I.n_columns * I.n_rows
	if (I.LayerCount >= 3)
		data = pointer([pointer(I.image[:,:,1]), pointer(I.image[:,:,2]), pointer(I.image[:,:,3])])
		cmap = C_NULL
	else
		data = pointer([pointer(I.image)])
	end
	palette_count = size(I.colormap, 1)
	color_space = IM_RGB
	if (size(I.colormap,1) > 1)
		cmap = I.colormap
		color_space = IM_MAP
	end
	im = imImage(I.n_columns, I.n_rows, color_space, 0, 0, I.LayerCount, I.n_columns, n,
	             n*I.LayerCount, n, data, pointer(cmap), palette_count, C_NULL)
	image = pointer([im])
	return gc_ref(image)
end

# ----------------------------------------------------------------------------------------------
function imlabCreateButtonImages()
	new_bits = Uint8[
		2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
		,3,2,2,1,3,2,2,1,2,2,2,2,2,2,2,2
		,1,3,2,1,4,2,1,3,2,2,2,2,2,2,2,2
		,2,1,3,1,3,1,3,0,0,0,0,0,0,0,0,2
		,2,3,1,3,4,1,1,1,1,4,4,4,4,4,0,2
		,1,1,1,4,3,3,4,4,4,4,4,4,4,4,0,2
		,2,2,4,1,4,4,1,4,4,4,4,4,4,4,0,2
		,1,2,2,1,3,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,1,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,2
	]

	new_colors = [
		"0 0 0",
		"132 132 132",
		"BGCOLOR",
		"255 255 0",
		"255 255 255",
		""
	]

	open_bits = Uint8[
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,
		1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,0,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
		1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,0,
		1,0,1,3,1,0,0,0,0,0,0,0,1,1,1,1,
		1,0,3,1,3,1,3,1,3,1,3,0,1,1,1,1,
		1,0,1,3,1,3,1,3,1,3,1,0,1,1,1,1,
		1,0,3,1,3,1,0,0,0,0,0,0,0,0,0,0,
		1,0,1,3,1,0,2,2,2,2,2,2,2,2,2,0,
		1,0,3,1,0,2,2,2,2,2,2,2,2,2,0,1,
		1,0,1,0,2,2,2,2,2,2,2,2,2,0,1,1,
		1,0,0,2,2,2,2,2,2,2,2,2,0,1,1,1,
		1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
	]

	open_colors = [
		"0 0 0",
		"BGCOLOR",
		"128 128 0",
		"255 255 0",
		""
	]

	line_bits = Uint8[
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,
		1,1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,
		1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,
		1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,
		1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,
		1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,
		1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,
		1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,
		1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,
		1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,
		1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,
		1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,
		1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
	]

	CreateButtonImage(16,16, new_bits, new_colors, "imImgWinNewButton");
	CreateButtonImage(16,16, open_bits, open_colors, "imImgWinOpenButton");
	CreateButtonImage(16,16, line_bits, open_colors, "imImgWinLineButton");
end

# ----------------------------------------------------------------------------------------------
function CreateButtonImage(w::Int, h::Int, bits, colors, name::String)
	global total_images, images

	color_str = ["0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15"]
	i = 1

	iup_image = IupImage(w, h, bits)
	IupSetHandle(name, iup_image)

	while (colors[i] != "")
		if (i > 16)
			aux_color_str = @sprintf("%d", i)
			IupStoreAttribute(iup_image, aux_color_str, colors[i]);
		else
			IupStoreAttribute(iup_image, color_str[i], colors[i]);
		end
		i += 1
	end

end

# ---------------------------------------------------------------------------------------------
function CreateMainMenu()

	# Creates items of menu File
	item_new   = IupItem ("New","");
	item_open  = IupItem ("Open Grid/Image", "");
	item_close = IupItem ("Close", "");
	item_exit  = IupItem ("Exit", "item_exit_act");

	# Creates items of menu Image
	item_copy  = IupItem ("Copy", "");
	item_paste = IupItem ("Paste", "")

	# Creates items for menu triangle
	item_scalenus = IupItem ("Scalenus","")
	item_isoceles = IupItem ("Isoceles", "")
	item_equilateral = IupItem ("Equilateral", "")

	# Create menu triangle
	menu_triangle = IupMenu(item_equilateral, item_isoceles, item_scalenus)

	# Creates submenu triangle
	submenu_triangle = IupSubmenu("Triangle", menu_triangle)

	# Creates items for menu create
	item_line = IupItem ("Line", "")
	item_circle = IupItem ("Circle", "")

	# Creates menu create
	menu_create = IupMenu(item_line, item_circle, submenu_triangle)

	# Creates submenu create
	submenu_create = IupSubmenu ("Create", menu_create)

	# Creates items of menu help
	item_help = IupItem ("Help", "item_help_act")

	# Creates menus
	menu_file  = IupMenu (item_new, item_open, item_close, IupSeparator(), item_exit)
	menu_image = IupMenu (item_copy, item_paste, IupSeparator(), submenu_create)
	menu_tools = IupMenu (IupItem("Tools",""))
	menu_help  = IupMenu (item_help);

	# Creates submenus
	submenu_file  = IupSubmenu ("File", menu_file)
	submenu_image = IupSubmenu ("Image", menu_image)
	submenu_tools = IupSubmenu ("Tools", menu_tools)
	submenu_help  = IupSubmenu ("Help", menu_help)

	# The menu Draw and its descendents
	menu_draw = IupMenu(IupItem("Draw",""))
	submenu_draw = IupSubmenu ("Draw", menu_draw)

	# The menu Geography and its descendents
	menu_geog = IupMenu(IupItem("Plot-coastline",""))
	submenu_geog = IupSubmenu ("Geography", menu_geog)

	# The menu Plates and its descendents
	menu_plates = IupMenu(IupItem("Plate calculator",""))
	submenu_plates = IupSubmenu ("Plates", menu_plates)

	# The menu Mag/Grav and its descendents
	menu_mag = IupMenu(IupItem("IGRF calculator",""))
	submenu_mag = IupSubmenu ("Mag/Grav", menu_mag)

	# The menu Seismology and its descendents
	menu_seismic = IupMenu(IupItem("Seismicity",""))
	submenu_seismic = IupSubmenu ("Seismology", menu_seismic)

	# The menu Tsunamis and its descendents
	menu_tsu = IupMenu(IupItem("Tsunami Travel Time",""))
	submenu_tsu = IupSubmenu("Tsunamis", menu_tsu)

	# The menu GMT and its descendents
	menu_gmt = IupMenu(IupItem("grdfilter",""))
	submenu_gmt = IupSubmenu("GMT", menu_gmt)

	# The menu Grid Tools and its descendents
	menu_grdtool = IupMenu(IupItem("grid calculator",""))
	submenu_grdtool = IupSubmenu("Grid Tools", menu_grdtool)

	# The menu Projections and its descendents
	menu_proj = IupMenu(IupItem("None",""))
	submenu_proj = IupSubmenu("Projections", menu_proj)

	# Creates main menu with file menu
	mainMenu = IupMenu(submenu_file, submenu_image, submenu_tools, submenu_draw, submenu_geog, submenu_plates,
	                   submenu_mag, submenu_seismic, submenu_tsu, submenu_gmt, submenu_grdtool, submenu_proj,
	                   submenu_help)

	return mainMenu
end

end 	# module