text:
using DataStructures
using DelimitedFiles
import JSON

pkgdir = "julia/LibDLF.jl/"

# Create package directories
mkpath(pkgdir * "src")
mkpath(pkgdir * "src/lib")
mkpath(pkgdir * "test")

# Copy package license
cp("LICENSE", pkgdir * "LICENSE" ,force=true)

# Copy filter data file license
cp("../LICENSE", pkgdir * "src/lib/LICENSE",force=true)

# Get current version number in git
version = chomp(split(read(`git describe --tags`, String), "-")[1][2:end])

# Create Project.toml
iop = open(abspath(pkgdir * "Project.toml"), "w")
println(iop,"name = \"LibDLF\"")
println(iop,"uuid = \"f0c3f387-4ff6-435f-9d63-77e28b8d1347\"")
println(iop,"authors = [\"The emsig community <info@emsig.xyz> \"]")
println(iop,"version = \"$version\"")
println(iop,"\n[deps]")
println(iop,"DelimitedFiles = \"8bb1440f-4735-579b-a4ab-409b98df4dab\"")
println(iop,"\n[compat]")
println(iop,"julia = \"1\"")
println(iop,"\n[extras]")
println(iop,"Test = \"8dfed614-e22c-5e08-85e1-65c5234f0b40\"")
println(iop,"\n[targets]")
println(iop,"test = [\"Test\"]")
close(iop)

# Open runtests.jl
iotest = open(abspath(pkgdir * "test/runtests.jl"), "w")
println(iotest,"using LibDLF")
println(iotest,"using Test")

# Test settings
println(iotest,"\n# Relative tolerance for test result accuracy")
println(iotest,"test_rtol = 1e-4")

println(iotest,"\n# Test at x=1. where filters usually have small relative error")
println(iotest,"x_test = 1.")

# Add testing function pairs to runtest.jl
println(iotest,"\n# Testing function pairs (from Anderson,1975, with a=1)")
println(iotest,"rhs_j0(x)   = exp(-x^2 / 4.) / 2.")
println(iotest,"rhs_j1(x)   = exp(-x^2 / 4.) * x / 4.")
println(iotest,"lhs_j0(lam) = exp(-lam^2) * lam")
println(iotest,"lhs_j1(lam) = exp(-lam^2) * lam^2")
println(iotest,"rhs_fcos(x) = exp(-x^2 / 4.) / 2 * sqrt(π)")
println(iotest,"rhs_fsin(x) = exp(-x^2 / 4.) / 4 * sqrt(π) * x")
println(iotest,"lhs_fcos(ω) = exp(-ω^2)")
println(iotest,"lhs_fsin(ω) = exp(-ω^2) * ω")

# J0 transform test function
println(iotest,"\nfunction j0_test(base,j0)")
println(iotest,"\tlam = base/x_test")
println(iotest,"\tfht = (lhs_j0.(lam)' * j0) / x_test")
println(iotest,"\trhs =  rhs_j0(x_test)")
println(iotest,"\treturn abs(fht - rhs) ./ abs(rhs) ")
println(iotest,"end")

# J1 transform test function
println(iotest,"\nfunction j1_test(base,j1)")
println(iotest,"\tlam = base/x_test")
println(iotest,"\tfht = (lhs_j1.(lam)' * j1) / x_test")
println(iotest,"\trhs =  rhs_j1(x_test)")
println(iotest,"\treturn abs(fht - rhs) ./ abs(rhs) ")
println(iotest,"end")

# cos transform test function
println(iotest,"\nfunction fcos_test(base,fcos)")
println(iotest,"\tω = base/x_test")
println(iotest,"\tfct = (lhs_fcos.(ω)' * fcos) / x_test")
println(iotest,"\trhs =  rhs_fcos(x_test)")
println(iotest,"\treturn abs(fct - rhs) ./ abs(rhs) ")
println(iotest,"end")

# sin transform test function
println(iotest,"\nfunction fsin_test(base,fsin)")
println(iotest,"\tω = base/x_test")
println(iotest,"\tfst = (lhs_fsin.(ω)' * fsin) / x_test")
println(iotest,"\trhs =  rhs_fsin(x_test)")
println(iotest,"\treturn abs(fst - rhs) ./ abs(rhs) ")
println(iotest,"end")

println(iotest,"\n# Unit tests")
# Read in .json file listing all filters
filters = JSON.parsefile("../lib/filters.json", dicttype=DataStructures.OrderedDict)

# Create Julia module
iol = open(abspath(pkgdir * "src/LibDLF.jl"), "w")

# Module name
println(iol,"module LibDLF\n")

# Create LibDLF.jl files
for type in filters.keys

    stype = titlecase(type)

    # Make directory for this type
    mkpath(pkgdir * "src/lib/$stype")

    # Add filter type section to runtests.jl
    println(iotest, "\n@testset \"$stype\" begin")

    # Include sub module file in parent
    println(iol, "include(\"$stype.jl\")")

    # Create sub module file for filter type
    iot = open(abspath(pkgdir * "src/$stype.jl"), "w")
    println(iot, "module $stype\n")

    # Add library path variable
    println(iot, "libpath = @__DIR__")

    # Add cache
    println(iot, "\ncache = Dict() # local cache for any filters already loaded")

    # Add filter functions
    for filt in filters[type]

        # Get and write header as docstring
        iof = open(abspath("../" * filt["file"]), "r")

        # Title
        println(iot, "\n\"\"\"")
        sname = filt["name"]
        println(iot, "\t $sname()\n")
        println(iot, readline(iof)[2:end])

        # Get vals and preformat if sin & cos
        vals = replace(filt["values"], "," => ", ")
        if type == "fourier"
            vals = replace(vals, "cos" => "fcos")
            vals = replace(vals, "sin" => "fsin")
        end

        # Rest of header
        for line in eachline(iof)

            # Do not copy the title-underline; just newline
            if contains(line, "========")
                println(iot, "")

            # Empty lines: only remove comment
            elseif line == "#\n"
                println(iot, "")

            # The license is the last thing of the header
            elseif contains(line, "This file is part of libdlf")

                # Add returned vals
                println(iot, "# Returns\n")
                println(iot, "base, $vals :: Array{Float64,1}")
                println(iot, "Filter base and its values.\n")

                # Example
                println(iot, "# Example\n")
                println(iot, "```julia")
                println(iot, "base, $vals = LibDLF.$stype.$sname()")
                println(iot, "```")

                # Finish header
                println(iot, "\n\"\"\"")

                # Stop header loop
                break

            # Print entire line
            else
                println(iot, line[2:end])

            end

        end

        # Read filter data from text file
        dat = readdlm("../" * filt["file"], comments=true)

        # Save filter data to binary file
        bfile = pkgdir * "src/" * filt["file"] * ".bin"
        iob = open(bfile, "w")
        write(iob, htol.(dat)) # htol converts from host endian from to little endian
        close(iob)

        # Write function to load binary data
        bfile  =  filt["file"] * ".bin"
        println(iot, "function $sname()")
        println(iot, "\tif !haskey(cache,\"$sname\") # read and add to cache")
        println(iot, "\t\tbfile = joinpath(libpath,\"$bfile\")")
        println(iot, "\t\tbdat = Array{Float64}(undef,", size(dat, 1), ",", size(dat, 2), ")");
        println(iot, "\t\tread!(open(bfile,\"r\"),bdat)")
        println(iot, "\t\tbdat .= ltoh(bdat) # file save with little endian so convert to host's endian")
        println(iot, "\t\tcache[\"$sname\"]= tuple([bdat[:,c] for c in 1:size(bdat,2)]...)")
        println(iot, "\tend")
        println(iot, "\treturn cache[\"$sname\"]")
        println(iot, "end")

        # Close file
        close(iof)

        # Add filter tests to runtests.jl
        println(iotest, "\n\t@testset \"$sname\" begin")
        println(iotest, "\t\tbase, $vals = LibDLF.$stype.$sname()")
        for val in strip.(split(vals, ","))
            ft = val * "_test"
            println(iotest, "\t\t@test $ft(base,$val) < test_rtol")
        end
        println(iotest, "\tend")
    end

    println(iotest, "end")

    # Close filter type sub module
    println(iot, "\nend")
    close(iot)
end

# Close LibDLF module
println(iol,"\nend")
close(iol)

# Close runtest.jl
close(iotest)

# Create README
iord = open(abspath(pkgdir * "README.md"), "w")
println(iord,"\n# Julia package for LibDLF\n")
str = "[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5172893.svg)]"
str *= "(https://doi.org/10.5281/zenodo.5172893)"
println(iord,"\n$str")
str = "[![package-julia]"
str *= "(https://github.com/emsig/libdlf/actions/workflows/package-julia.yml/badge.svg)]"
str *= "(https://github.com/emsig/libdlf/actions/workflows/package-julia.yml)"
println(iord,"\n$str")
println(iord,"\nThis package is auto-generated. See
[emsig/libdlf](https://github.com/emsig/libdlf) for installation and usage
instructions.")

close(iord)

max_stars_repo_name: emsig/libdlf

max_issues_repo_name: emsig/libdlf

max_forks_repo_name: emsig/libdlf


text:
"""
    @display expr

Display `expr` using [`print_boxed`](@ref) along with the location
from where the macro was called (useful for debugging); a mixture of
`@show`, `display`, and `@debug`.

# Examples

```julia-repl
julia> @display sin.(1:10)
┌ sin.(1:10) =  10-element Array{Float64,1}:
│    0.8414709848078965
│    0.9092974268256817
│    0.1411200080598672
│   -0.7568024953079282
│   -0.9589242746631385
│   -0.27941549819892586
│    0.6569865987187891
│    0.9893582466233818
│    0.4121184852417566
│   -0.5440211108893698
└  @ REPL[16]:1
```

"""
macro display(a,io=:stdout)
    aname = "$a"
    suffix = "@ "*string(__source__.file)*":"*string(__source__.line)
    quote
        buf = IOBuffer()
        show(buf, MIME"text/plain"(), $(esc(a)))
        msglines = split(String(take!(buf)), "\n")
        length(msglines) > 1 && push!(msglines, "")
        print_boxed($(esc(io)), msglines, $aname*" =", $suffix, color=:green)
    end
end

"""
    display_matrix(io, A)

Display the matrix `A` on `io` wrapped in square brackets. Quite
simple implementation that assumes the first line contains type
information and the rest of the lines are the matrix elements.

# Examples

```jldoctest
julia> n = 10; o = ones(n);

julia> T = Tridiagonal(o[2:end], -2o, o[2:end]);

julia> display_matrix(stdout, T)
 10×10 Tridiagonal{Float64,Array{Float64,1}}:
⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎤
⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢
⎣   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦
```

"""
function display_matrix(io::IO, A::AbstractMatrix)
    buf = IOBuffer()
    (height, width) = displaysize(io)
    show(IOContext(buf, :limit => true, :displaysize => (height, width-2)),
         MIME"text/plain"(), A)

    lines = split(String(take!(buf)), "\n")
    println(io, " ", lines[1])
    body = lines[2:end]
    if length(body) == 1
        println(io, "[", body[1], "]")
    else
        println(io, "⎡", body[1], "⎤")
        for l in body[2:end-1]
            println(io, "⎢", l, "⎢")
        end
        println(io, "⎣", body[end], "⎦")
    end
end

ext: jl


text:
ReLU(x::R) where {R<:Real} = max(x, 0)
ReLU(x::C) where {C<:Complex} = max(real(x), 0) + im * max(imag(x), 0)
logabs(x) = log.(abs.(x))

"""
    Tanh(x::Complex) = tanh(real(x))+im*tanh(imag(x))
"""
Tanh(x::T) where {T<:Complex} = tanh(real(x)) + im * tanh(imag(x)) # note that this is not the analytic function tanh
Tanh(x::T) where {T<:Real} = tanh(x)
"""
    softplus(x::Real)
    softplus(x::Complex)
Not the standard softplus, but one biased so that we return zero at zero. The complex operates on each component separately instead of using complex log/exponential.
"""
softplus(x::T) where {T<:Real} = log(1 + exp(x)) - log(2)
softplus(x::T) where {T<:Complex} = log(1 + exp(real(x))) - log(2) + im * (log(1 + exp(imag(x))) - log(2))
piecewiseLinear(x::T) where {T<:Real} = (x >= 0 ? x : x / 2)
piecewiseLinear(x::T) where {T<:Complex} = ((real(x) >= 0 && imag(x >= 0)) ? x : x / 2)

"""
    spInverse(x::Union{<:Real, <:Complex})
inverse of the softplus function (note that it isn't accurate for input smaller
than -36 because of numerical issues in the softplus function itself).
"""
spInverse(x::Real) = x > 0 ? log(exp(x) - 1) : -36.737
spInverse(x::Complex) = spInverse(real(x)) + im * spInverse(imag(x))
aTanh(x::T) where {T<:Real} = abs.(x) < 1 ? atanh(real(x)) : sign(x) * 10
aTanh(x::Complex) = aTanh(real(x)) + im * aTanh(imag(x))
plInverse(x::Real) = (x > 0 ? x : 2 * x)


text:
function generatePopPlot( popReport::DataFrame, timeFactor::Float64,
    showPlot::Bool, savePlot::Bool, filename::String )

    # Generate the plot.
    plotTitle = string( "Node population plot" )
    plotData = Matrix( popReport[:, 2:end] )
    labels = hcat( string.( names( popReport )[2:end] )... )
    ymax = max( maximum( plotData ), 20 )

    plt = plot( popReport[:, :timePoint] / timeFactor, plotData, show=false,
        size=(960, 540), title=plotTitle, labels=labels, lw=2, ylim=[0, ymax] )

    # Show plot if needed.
    if showPlot
        gui( plt )
    end  # if showPlot

    # Save plot if needed.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            filename = string( "Node population plot  ", dateStr )
        end  # if filename == ""

        if !ispath( dirname( filename ) )
            mkpath( dirname( filename ) )
        end  # if !ispath( dirname( filename ) )

        extension = split( filename, "." )
        extension = length( extension ) == 1 ? "" : extension[end]
        filename = string( filename, extension ∈ extensions ? "" : ".svg" )
        savefig( plt, filename )
    end  # if savePlot

end  # generatePopPlot( popReport, timeFactor, showPlot, savePlot, filename )

function generateEvolutionPlot( popReport::DataFrame, mpSim::MPsim,
    fluxReports::Dict{String,Tuple}, nodes::NTuple{N,String},
    timeFactor::Float64, showPlot::Bool, savePlot::Bool, filename::String,
    plotFunction::Function, plotFunction2::Function ) where N

    fileroot = nothing
    extension = nothing

    # Get file name ready.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            fileroot = string( "Population evolution plot (", dateStr, ") of " )
            extension = "svg"
        else
            extension = split( filename, "." )
            fileroot = length( extension ) == 1 ? string( extension[1] ) :
                join( extension[1:(end - 1)], "." )
            extension = length( extension ) == 1 ? "" : extension[end]
            extension = extension ∈ extensions ? extension : "svg"
        end  # if filename == ""

        if !ispath( dirname( fileroot ) )
            mkpath( dirname( fileroot ) )
        end  # if !ispath( dirname( fileroot ) )
    end  # if savePlot

    timePoints = popReport[:, :timePoint] / timeFactor

    # Generate plots.
    for node in unique( nodes )
        if !haskey( fluxReports, node )
            continue
        end  # if !haskey( fluxReports, node )

        plotfilename = string( fileroot, " ", node, ".", extension )
        fluxReport = fluxReports[node]

        # Data of overview plot.
        pop = popReport[:, Symbol( node )]
        inFlux = fluxReport[1][:, 3]
        outFlux = fluxReport[2][:, 3]
        netFlux = inFlux - outFlux
        plotData = hcat( pop, inFlux, outFlux, netFlux )
        ymin = minimum( vcat( 0, netFlux ) )
        ymax = maximum( pop )
        ymax = ymax - 20 < ymin ? ymin + 20 : ymax

        # Plots.
        plotTitle = string( "Node population plot" )
        popPlot = plot( timePoints, plotData, title=plotTitle, show=false,
            labels=["population" "in flux" "out flux" "net flux"],
            color=[:black :green :red :blue], lw=2,
            ylim=[ymin, ymax] )
        inFluxPlot = plotFunction( fluxReport[1], string( node ), :in,
            timeFactor )
        outFluxPlot = plotFunction( fluxReport[2], string( node ), :out,
            timeFactor )

        if haskey( mpSim.compoundNodeList, node )
            withinFluxPlot = plotFunction( fluxReport[3], string( node ),
                :within, timeFactor )
            compPlot = plotFunction2( fluxReport[4], string( node ),
                timeFactor )

            layout = @layout [
                grid( 1, 2 ){0.5625h}
                grid( 1, 3 )
            ]
            plt = plot( popPlot, compPlot, inFluxPlot, outFluxPlot,
                withinFluxPlot, show=false, size=(1440, 960), layout=layout )
        else
            layout = @layout [
                a{0.5625h}
                grid( 1, 2 )
            ]
            plt = plot( popPlot, inFluxPlot, outFluxPlot, show=false,
                size=(960, 960), layout=layout )
        end  # if haskey( mpSim.compoundNodeList, node )

        # Show plot if needed.
        if showPlot
            gui( plt )
        end  # if showPlot

        # Save plot if needed.
        if savePlot
            savefig( plt, plotfilename )
        end  # if savePlot
    end  # for node in unique( nodes )

end  # generateEvolutionPlot( popReport, mpSim, fluxReports, nodes, timeFactor,
     #   showPlot, savePlot, fileroot, extension, plotFunction, plotFunction2 )


text:
action_types_dict = Dict(
:RIGHT_ONLY => [
    (:NOOP,),
    (:right,),
    (:right, :A),
    (:right, :B),
    (:right, :A, :B)
],

:SIMPLE_MOVEMENT => [
    (:NOOP,),
    (:right,),
    (:right, :A),
    (:right, :B),
    (:right, :A, :B),
    (:A,),
    (:left,)
],

:COMPLEX_MOVEMENT => [
    (:NOOP,),
    (:right,),
    (:right, :A),
    (:right, :B),
    (:right, :A, :B),
    (:A,),
    (:left, :A),
    (:left, :B),
    (:left, :A, :B),
    (:down,),
    (:up,)
])
button_map = Dict{Symbol, UInt8}(
    :right  => 0b10000000,
    :left   => 0b01000000,
    :down   => 0b00100000,
    :up     => 0b00010000,
    :start  => 0b00001000,
    :select => 0b00000100,
    :B      => 0b00000010,
    :A      => 0b00000001,
    :NOOP   => 0b00000000,
)


text:
x = create_var(3, 2)
r = repmat(x, 4, 5)
test_expr(r, (x,))

x = create_var(1, 1)
r = repmat(x, 1, 1)
test_expr(r, (x,))

x = create_var(3, 3)
y = x + [1 2 3; 4 5 6; 7 8 9]
r = repmat(y, 1, 3)
test_expr(r, (x,))

x = create_var(3, 3)
y = x + [1 2 3; 4 5 6; 7 8 9]
r = repmat(y, 4, 1)
test_expr(r, (x,))

info("All repmat tests passed")


text:
using Test, CacheServers

include("TestFragment.jl")

s = DefaultServer{Scalar, Variable{Float64}}()
var = Scalar(2.0)

@testset "alloc" begin
    alloc!(s, var, Variable(var.val))
    @test s.storage[objectid(var)] isa Variable{Float64}
    @test iscacheable(s, var) == true
    @test iscached(s, var) == true # always cached
end

@testset "push & pull" begin
    push!(s, Grad(3.0), var)
    @test pull(s, var, Grad) == 3.0
    @test pull(s, var, Param) == 2.0

    push!(s, Param(1.0), var)
    @test pull(s, var, Grad) == 3.0
    @test pull(s, var, Param) == 1.0
end

@testset "delete & clear" begin
    clear!(s, var)
    @test pull(s, var, Grad) == 0.0
    delete!(s, var)
    @test_throws KeyError pull(s, var, Grad)
end


text:
using nevanlinna_ac
using PyPlot
using DelimitedFiles
using Printf

fac = [0, 1]
A = x-> x*gaussian(x,fac)
x = [i for i in range(-5,5,length=1000)]
y = [A(i) for i in x ]

plot(x,y)
PyPlot.display_figs()

open("./data/ag3_bose_b20.txt","w") do file
    for i = 1:1000
        writedlm(file,[x[i] A(x[i])])
    end
end

num = 30
β = 20

B = x -> gaussian(x,fac)
ωn = [Masubara_freq(i,β,type=:b) for i=1:num]
GF = [Masubara_GF(i,B,β,type=:b) for i=1:num]

open("./data/ig3_bose_b20.txt","w") do file
    for i = 1:num
        writedlm(file,[ωn[i] real(GF[i]) imag(GF[i])])
    end
end


text:
"""
    OffsetInteger{O, T}

OffsetInteger type mainly for indexing.
* `O` - The offset relative to Julia arrays. This helps reduce copying when
communicating with 0-indexed systems such ad OpenGL.
"""
struct OffsetInteger{O, T <: Integer} <: Integer
    i::T
    OffsetInteger{O, T}(x::Integer) where {O, T <: Integer} = new{O, T}(T(x + O))
end

const ZeroIndex{T <: Integer} = OffsetInteger{-1, T}
const OneIndex{T <: Integer} = OffsetInteger{0, T}
const GLIndex = ZeroIndex{Cuint}

raw(x::OffsetInteger) = x.i
raw(x::Integer) = x
value(x::OffsetInteger{O, T}) where {O, T} = raw(x) - O
value(x::Integer) = x

function show(io::IO, oi::OffsetInteger{O, T}) where {O, T}
    print(io, "|$(raw(oi)) (indexes as $(value(oi))|")
end

Base.eltype(::Type{OffsetInteger{O, T}}) where {O, T} = T
Base.eltype(oi::OffsetInteger) = eltype(typeof(oi))

# constructors and conversion
OffsetInteger{O1, T1}(x::OffsetInteger{O2, T2}) where {O1, O2, T1 <: Integer, T2 <: Integer} = OffsetInteger{O1, T1}(T2(x))

OffsetInteger{O}(x::Integer) where {O} = OffsetInteger{O, eltype(x)}(x)
OffsetInteger{O}(x::OffsetInteger) where {O} = OffsetInteger{O, eltype(x)}(x)
# This constructor has a massive method invalidation as a consequence,
# and doesn't seem to be needed, so let's remove it!
(::Type{IT})(x::OffsetInteger) where {IT <: Integer} = IT(value(x))

Base.@pure pure_max(x1, x2) = x1 > x2 ? x1 : x2
Base.promote_rule(::Type{T1}, ::Type{OffsetInteger{O, T2}}) where {T1 <: Integer, O, T2} = T1
Base.promote_rule(::Type{OffsetInteger{O1, T1}}, ::Type{OffsetInteger{O2, T2}}) where {O1, O2, T1, T2} = OffsetInteger{pure_max(O1, O2), promote_type(T1, T2)}

#Need to convert to Int here because of: https://github.com/JuliaLang/julia/issues/35038
Base.to_index(I::OffsetInteger) = convert(Int, raw(OneIndex(I)))
Base.to_index(I::OffsetInteger{0}) = convert(Int, raw(I))

# basic operators
for op in (:(-), :abs)
    @eval Base.$op(x::T) where {T <: OffsetInteger} = T($(op)(value(x)))
end

for op in (:(+), :(-), :(*), :(/), :div)
    @eval begin
        @inline function Base.$op(x::OffsetInteger{O}, y::OffsetInteger{O}) where O
            OffsetInteger{O}($op(value(x), value(y)))
        end
    end
end

for op in (:(==), :(>=), :(<=), :(<) , :(>), :sub_with_overflow)
    @eval begin
        @inline function Base.$op(x::OffsetInteger{O}, y::OffsetInteger{O}) where O
            $op(x.i, y.i)
        end
    end
end


text:
using HMRowUnification
using Documenter

makedocs(;
    modules=[HMRowUnification],
    authors="thautwarm <twshere@outlook.com> and contributors",
    repo="https://github.com/thautwarm/HMRowUnification.jl/blob/{commit}{path}#L{line}",
    sitename="HMRowUnification.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://thautwarm.github.io/HMRowUnification.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/thautwarm/HMRowUnification.jl",
)


text:
#=

   The Euro coin problem.

   From Think Bayes, page 33ff
   """
   A statistical statement appeared in "The Guardian" on Friday January 4, 2002:
      When spun on edge 250 times, a Belgian one-euro coin
      came up heads 140 times and tails 110. 'It looks very
      suspicious to me,' said Barry Blight, a statistics lecturer
      at the London School of Economics. 'If the coin were
      unbiased, the chance of getting a result as extreme as
      that would be less than 7%.'

   But do these data give evidence that the coin is biased rather than fair?
   """

   Continues on page 41:
   """
   Exercise 4.1. Suppose that instead of observing coin tosses directly, you measure
   the outcome using an instrument that is not always correct. Specifically, suppose
   there is a probability y that an actual heads is reported as tails, or actual tails re-
   ported as heads.

   Write a class that estimates the bias of a coin given a series of outcomes and the
   value of y .
   How does the spread of the posterior distribution depend on y ?
   """

   mean(prob): 0.5213151098904925
   Distributions of variable probLt0_5 (num:0)
   0.00000 =>    7020  (0.702000)
   1.00000 =>    2980  (0.298000)

   Distributions of variable error (num:0)
   0.00000 =>    7051  (0.705100)
   1.00000 =>    2949  (0.294900)

   Cf euro_coin_problem.jl
   Cf ~/blog/euro_coin_unreliable_measurements.blog
      ~/webppl/euro_coin_unreliable_measurements.wppl
=#

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

@model function euro_coin_problem_unreliable_measurements()
    n = 250

    # Probability of throwing head
    prob ~ Beta(2,2)

    # We measure incorrect with probability 0.2
    error ~ flip(0.2)

    throwCoin ~ filldist(flip(prob),n)

    coin = tzeros(n)
    for i in 1:n
        coin[i] ~ Dirac(error == true ? 1-throwCoin[i] : throwCoin[i])
    end

    sum250 ~ Dirac(sum([coin[i] == true ? 1 : 0 for i in 1:n]))
    true ~ Dirac(sum250 == 140)
    probLt0_5 ~ Dirac(prob < 0.5)

end

model = euro_coin_problem_unreliable_measurements()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

println("mean(prob): ", mean(chns[:prob]))
show_var_dist_pct(chns, :probLt0_5)
show_var_dist_pct(chns, :error)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using BinaryBuilder

# Collection of sources required to build openspecfun
sources = [
    "https://github.com/JuliaLang/openspecfun/archive/v0.5.3.tar.gz" =>
    "1505c7a45f9f39ffe18be36f7a985cb427873948281dbcd376a11c2cd15e41e7",
]

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

# It needs to be told it's on Windows
if [[ ${target} == *mingw* ]]; then
    OS=WINNT
    libdir=$prefix/bin
elif [[ ${target} == *darwin* ]]; then
    OS=Darwin
fi

# Build it
make OS=${OS} CC="$CC" CXX="$CXX" FC="$FC" -j${nproc}

# Install it
mkdir -p $libdir $prefix/include
mv libopenspecfun*.${dlext}* $libdir/
"""

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

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libopenspecfun", :libopenspecfun)
]

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

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

max_stars_repo_name: andreasnoack/OpenspecfunBuilder

max_issues_repo_name: andreasnoack/OpenspecfunBuilder

max_forks_repo_name: andreasnoack/OpenspecfunBuilder


text:
#
# Defaults
#
CIRCLE_DEFAULT_PLOT_POINTS = 720
WIREFRAME_DEFAULT = 32
SURFACE_RESOLUTION_DEFAULT = 32
#
# Plotting Recipe – Poincaré Ball
#
@recipe function f(
    M::Hyperbolic{2},
    pts::AbstractVector{P},
    vecs::Union{AbstractVector{T},Nothing}=nothing;
    circle_points=CIRCLE_DEFAULT_PLOT_POINTS,
    geodesic_interpolation=-1,
    hyperbolic_border_color=RGBA(0.0, 0.0, 0.0, 1.0),
) where {P<:PoincareBallPoint,T<:PoincareBallTVector}
    @series begin
        φr = range(0, stop=2 * π, length=circle_points)
        x = [cos(φ) for φ in φr]
        y = [sin(φ) for φ in φr]
        seriestype := :path
        seriescolor := hyperbolic_border_color
        label := ""
        x, y
    end
    show_axis --> false
    framestyle -> :none
    axis --> false
    xlims --> (-1.01, 1.01)
    ylims --> (-1.01, 1.01)
    grid --> false
    aspect_ratio --> :equal
    tickfontcolor --> RGBA(1.0, 1.0, 1.0, 1.0)
    x = []
    y = []
    if vecs === nothing
        if geodesic_interpolation < 0
            seriestype --> :scatter
            x = [p.value[1] for p in pts]
            y = [p.value[2] for p in pts]
        else
            lpts = empty(pts)
            for i in 1:(length(pts) - 1)
                # push interims points on geodesics between two points.
                push!(
                    lpts,
                    shortest_geodesic(
                        M,
                        pts[i],
                        pts[i + 1],
                        collect(range(0, 1, length=geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                    )...,
                )
            end
            push!(lpts, last(pts)) # add last end point
            # split into x, y and plot as curve
            seriestype --> :path
            x = [p.value[1] for p in lpts]
            y = [p.value[2] for p in lpts]
        end
    else
        seriestype := :quiver
        quiver := ([v.value[1] for v in vecs], [v.value[2] for v in vecs])
        x = [p.value[1] for p in pts]
        y = [p.value[2] for p in pts]
    end
    return x, y
end
#
# Plotting Recipe – Poincaré Half plane
#
@recipe function f(
    M::Hyperbolic{2},
    pts::AbstractVector{P},
    vecs::Union{AbstractVector{T},Nothing}=nothing;
    geodesic_interpolation=-1,
) where {P<:PoincareHalfSpacePoint,T<:PoincareHalfSpaceTVector}
    aspect_ratio --> :equal
    framestyle --> :origin
    x = []
    y = []
    if vecs === nothing
        if geodesic_interpolation < 0
            seriestype --> :scatter
            x = [p.value[1] for p in pts]
            y = [p.value[2] for p in pts]
        else
            lpts = empty(pts)
            for i in 1:(length(pts) - 1)
                # push interims points on geodesics between two points.
                push!(
                    lpts,
                    shortest_geodesic(
                        M,
                        pts[i],
                        pts[i + 1],
                        collect(range(0, 1, length=geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                    )...,
                )
            end
            push!(lpts, last(pts)) # add last end point
            # split into x, y, z and plot as curve
            seriestype --> :path
            x = [p.value[1] for p in lpts]
            y = [p.value[2] for p in lpts]
        end
    else
        quiver := ([v.value[1] for v in vecs], [v.value[2] for v in vecs])
        seriestype := :quiver
        x = [p.value[1] for p in pts]
        y = [p.value[2] for p in pts]
    end
    return x, y
end
#
# Plotting Recipe – Hyperboloid
#
@recipe function f(
    M::Hyperbolic{2},
    pts::Union{AbstractVector{P},Nothing}=nothing,
    vecs::Union{AbstractVector{T},Nothing}=nothing;
    geodesic_interpolation=-1,
    wireframe=true,
    wires=WIREFRAME_DEFAULT,
    wires_x=wires,
    wires_y=wires,
    wireframe_color=RGBA(0.0, 0.0, 0.0, 1.0),
    surface=false,
    surface_resolution=SURFACE_RESOLUTION_DEFAULT,
    surface_resolution_x=surface_resolution,
    surface_resolution_y=surface_resolution,
    surface_color=RGBA(0.9, 0.9, 0.9, 0.8),
) where {P,T}
    px = isnothing(pts) ? [-2.0, 2.0] : [p[1] for p in pts]
    py = isnothing(pts) ? [-2.0, 2.0] : [p[2] for p in pts]
    # part I: wire
    if wireframe
        x = range(min(px...), max(px...), length=wires_x)
        y = range(min(py...), max(py...), length=wires_y)
        z = sqrt.(1 .+ (x .^ 2)' .+ y .^ 2)
        @series begin
            seriestype := :wireframe
            seriescolor := wireframe_color
            x, y, z
        end
    end
    # part II: solid hyperboloid
    if surface
        x = range(min(px...), max(px...), length=surface_resolution_x)
        y = range(min(py...), max(py...), length=surface_resolution_y)
        z = sqrt.(1 .+ (x .^ 2)' .+ y .^ 2)
        @series begin
            seriestype := :surface
            seriescolor := surface_color
            x, y, z
        end
    end
    show_axis --> false
    framestyle -> :none
    axis --> false
    xlims --> (min(px...), max(px...))
    ylims --> (min(py...), max(py...))
    grid --> false
    colorbar --> false
    tickfontcolor --> RGBA(1.0, 1.0, 1.0, 1.0)
    #
    # just pts given -> plot points
    x = []
    y = []
    z = []
    if pts !== nothing
        if vecs === nothing
            if geodesic_interpolation < 0
                seriestype --> :scatter
                x = [p[1] for p in pts]
                y = [p[2] for p in pts]
                z = [p[3] for p in pts]
            else
                lpts = empty(pts)
                for i in 1:(length(pts) - 1)
                    # push interims points on geodesics between two points.
                    push!(
                        lpts,
                        shortest_geodesic(
                            M,
                            pts[i],
                            pts[i + 1],
                            collect(range(0, 1, length=geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                        )...,
                    )
                end
                push!(lpts, last(pts)) # add last end point
                # split into x, y and plot as curve
                seriestype --> :path
                x = [p[1] for p in lpts]
                y = [p[2] for p in lpts]
                z = [p[3] for p in lpts]
            end
        else
            seriestype := :quiver
            quiver := ([v[1] for v in vecs], [v[2] for v in vecs], [v[3] for v in vecs])
            x = [p[1] for p in pts]
            y = [p[2] for p in pts]
            z = [p[3] for p in pts]
        end
        return x, y, z
    else
        return nothing
    end
end
#
# Plotting Recipe – Sphere
#
@recipe function f(
    M::Sphere{2,ℝ},
    pts::Union{AbstractVector{P},Nothing}=nothing,
    vecs::Union{AbstractVector{T},Nothing}=nothing;
    geodesic_interpolation=-1,
    wireframe=true,
    wires=WIREFRAME_DEFAULT,
    wires_lat=wires,
    wires_lon=wires,
    wireframe_color=RGBA(0.0, 0.0, 0.0, 1.0),
    surface=false,
    surface_resolution=SURFACE_RESOLUTION_DEFAULT,
    surface_resolution_lat=surface_resolution,
    surface_resolution_lon=surface_resolution,
    surface_color=RGBA(0.9, 0.9, 0.9, 0.8),
) where {P,T}
    # part I: wire
    if wireframe
        u = range(0, 2π, length=wires_lon + 1)
        v = range(0, π, length=wires_lat + 1)
        x = cos.(u) * sin.(v)'
        y = sin.(u) * sin.(v)'
        z = repeat(cos.(v)', outer=[wires_lon + 1, 1])
        @series begin
            seriestype := :wireframe
            seriescolor := wireframe_color
            return x, y, z
        end
    end
    # part II: solid sphere
    if surface
        u = range(0, 2π, length=surface_resolution_lon + 1)
        v = range(0, π, length=surface_resolution_lat + 1)
        x = cos.(u) * sin.(v)'
        y = sin.(u) * sin.(v)'
        z = repeat(cos.(v)', outer=[wires_lon + 1, 1])
        @series begin
            seriestype := :surface
            seriescolor := surface_color
            return x, y, z
        end
    end
    show_axis --> false
    framestyle -> :none
    axis --> false
    xlims --> (-1.05, 1.05)
    ylims --> (-1.05, 1.05)
    zlims --> (-1.05, 1.05)
    grid --> false
    colorbar --> false
    tickfontcolor --> RGBA(1.0, 1.0, 1.0, 1.0)
    x = []
    y = []
    z = []
    if pts !== nothing
        if vecs === nothing #plot points/geodesics
            if geodesic_interpolation < 0
                seriestype --> :scatter
                x = [p[1] for p in pts]
                y = [p[2] for p in pts]
                z = [p[3] for p in pts]
            else
                lpts = empty(pts)
                for i in 1:(length(pts) - 1)
                    # push interims points on geodesics between two points.
                    push!(
                        lpts,
                        shortest_geodesic(
                            M,
                            pts[i],
                            pts[i + 1],
                            collect(range(0, 1, length=geodesic_interpolation + 2))[1:(end - 1)], # omit end point
                        )...,
                    )
                end
                push!(lpts, last(pts)) # add last end point
                # split into x, y and plot as curve
                seriestype --> :path
                x = [p[1] for p in lpts]
                y = [p[2] for p in lpts]
                z = [p[3] for p in lpts]
            end
        else
            seriestype := :quiver
            quiver := ([v[1] for v in vecs], [v[2] for v in vecs], [v[3] for v in vecs])
            x = [p[1] for p in pts]
            y = [p[2] for p in pts]
            z = [p[3] for p in pts]
        end
    end
    return x, y, z
end


text:
let user_data_dir
    ENV["PATH"] = JULIA_HOME*";"*joinpath(JULIA_HOME,"..","Git","bin")*";"*ENV["PATH"]
    haskey(ENV,"JULIA_EDITOR") || (ENV["JULIA_EDITOR"] = "start")
    user_data_dir = abspath(ENV["AppData"],"julia")
    isdir(user_data_dir) || mkdir(user_data_dir)
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module MixedTests

using Interpolations, Base.Test

N = 10

for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
    A2 = rand(Float64, N, N) * 100
    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
        itp_a = @inferred(constructor(copier(A2), (BSpline(Linear()), BSpline(Quadratic(BC()))), GT()))
        itp_b = @inferred(constructor(copier(A2), (BSpline(Quadratic(BC())), BSpline(Linear())), GT()))
        @test @inferred(size(itp_a)) == size(A2)
        @test @inferred(size(itp_b)) == size(A2)

        for j = 2:N-1, i = 2:N-1
            @test_approx_eq_eps itp_a[i,j] A2[i,j] sqrt(eps(A2[i,j]))
            @test_approx_eq_eps itp_b[i,j] A2[i,j] sqrt(eps(A2[i,j]))
        end

        for i = 1:10
            dx, dy = rand(), rand()
            @test_approx_eq itp_a[2+dx,2] (1-dx)*A2[2,2]+dx*A2[3,2]
            @test_approx_eq itp_b[2,2+dy] (1-dy)*A2[2,2]+dy*A2[2,3]
        end
    end
end

# AbstractArrays
function copyshared(A)
    B = SharedArray(eltype(A), size(A))
    copy!(B, A)
end

for (constructor, copier) in ((interpolate, x->x), (interpolate!, copyshared))
    A2 = SharedArray(Float64, (N,N), init=A->rand!(A))
    for i = 1:length(A2)
        A2[i] *= 100
    end
    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
        itp_a = @inferred(constructor(copier(A2), (BSpline(Linear()), BSpline(Quadratic(BC()))), GT()))
        itp_b = @inferred(constructor(copier(A2), (BSpline(Quadratic(BC())), BSpline(Linear())), GT()))
        if constructor == interpolate!
            @test isa(itp_a.coefs, SharedArray)
            @test isa(itp_b.coefs, SharedArray)
        end
        @test @inferred(size(itp_a)) == size(A2)
        @test @inferred(size(itp_b)) == size(A2)

        for j = 2:N-1, i = 2:N-1
            @test_approx_eq_eps itp_a[i,j] A2[i,j] sqrt(eps(A2[i,j]))
            @test_approx_eq_eps itp_b[i,j] A2[i,j] sqrt(eps(A2[i,j]))
        end

        for i = 1:10
            dx, dy = rand(), rand()
            @test_approx_eq itp_a[2+dx,2] (1-dx)*A2[2,2]+dx*A2[3,2]
            @test_approx_eq itp_b[2,2+dy] (1-dy)*A2[2,2]+dy*A2[2,3]
        end
    end
end

end


text:
"""
    isinside(chart, point) -> Bool

Returns true is the given point is in the image of the given chart, false otherwise.
"""
function isinside(cell, point)

  u = carttobary(cell, point)
  T = eltype(u)

  tol = eps(T) * 1e3
  w = one(T)
  for i in 1:dimension(cell)
    0+tol < u[i] < 1-tol || return false
    w -= u[i]
  end

  0+tol < w < 1-tol || return false
  return true

end

"""
	isinclosure(simplex, point) -> Bool

Determine whether point is in the closure of simplex. False positives are possible for points just outside of the simplex.
"""
function isinclosure(cell, point)

  u = carttobary(cell, point)
  T = eltype(u)
  dim = dimension(cell)
  udim = length(point)

  tol = eps(T) * 1e3
  w = one(T)
  for i in 1:dim
    0-tol < u[i] < 1+tol || return false
    w -= u[i]
  end

  0-tol < w < 1+tol || return false

  # finally, check that the point indeed is in the plane of the cell
  A = zeros(T, udim, dim+1)
  origin = cell[dim+1]
  for i in 1:dim
      for j in 1:udim
          A[j,i] = cell[i][j] - origin[j]
      end
  end
  for j in 1:udim; A[j,end] = point[j] - origin[j]; end
  rank(A) == dim+1 && return false

  return true

end

max_stars_repo_name: shakisparki/CompScienceMeshes

max_issues_repo_name: shakisparki/CompScienceMeshes

max_forks_repo_name: shakisparki/CompScienceMeshes


text:
@testset "Checks" begin
    X = randn(5, 3)
    y = randn(5)
    @test R.check_nrows(X, y) === nothing
    y = randn(4)
    @test_throws DimensionMismatch R.check_nrows(X, y)
    @test R.check_pos(1)
    @test_throws ArgumentError R.check_pos(-1)
end

@testset "Augment" begin
    X = randn(5, 3)
    X1 = R.augment_X(X, false)
    @test X1 === X
    X1 = R.augment_X(X, true)
    @test X1 == hcat(X, ones(5, 1))
end

@testset "Apply" begin
    n, p = 5, 3
    X = randn(n, p)
    X1 = R.augment_X(X, true)
    θ = randn(p)
    θ1 = randn(p+1)
    @test R.apply_X(X, θ) == X*θ
    @test R.apply_X(X, θ1) ≈ X1*θ1
    # multiclass
    c = 3
    θa  = randn(p)
    θb  = randn(p)
    θc  = randn(p)
    θ   = vcat(θa, θb, θc)
    θa1 = randn(p+1)
    θb1 = randn(p+1)
    θc1 = randn(p+1)
    θ1  = vcat(θa1, θb1, θc1)
    @test R.apply_X(X, θ, c) ≈  X * hcat(θa, θb, θc)
    @test R.apply_X(X, θ1, c) ≈ X1 * hcat(θa1, θb1, θc1)

    θ1 = randn(p+1)
    Xθ = zeros(n)
    R.apply_X!(Xθ, X, θ1)
    @test Xθ ≈ X1 * θ1

    z = randn(n)
    Xtz = zeros(p)
    R.apply_Xt!(Xtz, X, z)
    @test Xtz ≈ X'z

    Xtz2 = zeros(p+1)
    R.apply_Xt!(Xtz2, X, z)
    @test Xtz2 ≈ X1'z

    v = view(Xtz2, 1:p)
    R.apply_Xt!(v, X, z)
    @test Xtz2 ≈ X1'z
end

@testset "Sigmoid" begin
    @test R.sigmoid(zero(Float32)) == 0.5f0
    @test R.sigmoid(zero(Float64)) == 0.5
    @test R.logsigmoid(zero(Float32)) == log(0.5f0)
    @test R.logsigmoid(zero(Float64)) == log(0.5)

    @test R.sigmoid(50) == 1.0
    @test R.sigmoid(50f0) == 1.0f0
    @test R.sigmoid(-50) == 0.0
    @test R.sigmoid(-50f0) == 0.0f0

    @test R.logsigmoid(50) == 0.0
    @test R.logsigmoid(50f0) == 0.0f0
    @test R.logsigmoid(-50) == -50.0
    @test R.logsigmoid(-50f0) == -50.0f0

    x = randn()
    @test -R.σ(-x) ≈ (R.σ(x) - 1.0)
end

@testset "Hat matrix" begin
    λ = 3
    X = randn(5, 3)
    X1 = R.augment_X(X, true)
    @test R.form_XtX(X, true, λ) ≈ X1'X1 + λ*I
end

@testset "Soft-thresh" begin
    x = randn(50)
    η = 0.5
    y = R.soft_thresh.(x, η)

    z = copy(x)
    m1 = x .> η
    m2 = abs.(x) .<= η
    m3 = x .< -η

    z[m1] .= x[m1] .- η
    z[m2] .= 0.0
    z[m3] .= x[m3] .+ η

    @test z == y
end

@testset "Clip" begin
    x = randn(10)
    x[[1, 4, 7]] .= 1e-7 .* rand(3)
    tau = 1e-5
    z = R.clip.(x, tau)
    mask = abs.(x) .>= tau
    @test z[mask] == x[mask]
    @test all(zi == tau for zi in z[.!mask])
end


text:
#using JuMP
using Ipopt
include("mds_metric.jl")

contains_inf(array) = mapreduce(x->(isinf(x)), &, array)

# A data structure representing the problem data and constraints
type MDSProblem
    num_points::Int
    dist::Array{Float64, 2}
    sphere_radius::Float64
    lads::Array{Int, 1}
    nads::Array{Int, 1}
    nucleolus_radius::Float64

    MDSProblem(n::Int, dist::Array{Float64, 2}) =
        new(n, dist, 1.0, Array(Int,0), Array(Int,0), 0)
    MDSProblem(n::Int, dist::Array{Float64, 2}, r::Float64) =
        new(n, dist, r, Array(Int,0), Array(Int,0), 0)
    # New problem with just a radius constraint
    MDSProblem(dist::Array{Float64, 2}, r::Float64) =
        new(size(dist)[1], dist, r, Array(Int,0), Array(Int,0), 0)
    MDSProblem(dist::Array{Float64, 2}, r::Float64, lads_file::String) =
        new(size(dist)[1], dist, r, load_lads(lads_file), Array(Int,0), 0)
    MDSProblem(dist::Array{Float64, 2}, r::Float64, lads_file::String, nads_file::String, nucleolus_r::Float64) =
        new(size(dist)[1], dist, r, load_lads(lads_file), load_lads(nads_file), nucleolus_r)
end

function load_data(filename::String; exponent=1/3, scale=1)
    # loads data, convert to NxN array of "wish distances"
    contact_data = readcsv(filename, Float64)
    #dist = contact_freq_to_distance_matrix(filename)
    dist = scale./(contact_data).^exponent
    return dist
end

function load_lads(filename::String)
    # should return a Vector{Float64}
    # a "lads file" is just the length followed by a list of indices.
    f = open(filename, "r")
    length = int(strip(readline(f)))
    lads = Array(Int, length)
    for (i, line) in enumerate(readlines(f))
        lads[i] = int(strip(line))
    end
    return lads
end

function remove_zeros(dist::Array{Float64, 2})
    # Removes all-zero and all-inf rows/columns from a contact map
    dist_index = 1
    while dist_index <= size(dist)[1]
        if mapreduce(x->(isinf(x) || x==0), &, dist[dist_index, 1:end-1])
            dist = [dist[1:dist_index-1, :], dist[dist_index+1:end, :]]
            dist = [dist[:, 1:dist_index-1] dist[:,dist_index+1:end]]
        else
            dist_index += 1
        end
    end
    return dist
end

function remove_infs(prob::MDSProblem)
    # Removes all rows from the distance matrix that have all infinities.
    dist_index = 1
    dist = prob.dist
    while dist_index <= size(dist)[1]
        if mapreduce(x->(isinf(x) || x==0), &, dist[dist_index, 1:end-1])
            dist = [dist[1:dist_index-1, :], dist[dist_index+1:end, :]]
            dist = [dist[:, 1:dist_index-1] dist[:,dist_index+1:end]]
            new_lads = filter(x->x!=dist_index, prob.lads)
            for (i, lad) in enumerate(new_lads)
                if lad > dist_index
                    new_lads[i] -=1
                end
            end
            prob.lads = new_lads
            new_nads = filter(x->x!=dist_index, prob.nads)
            for (i, nad) in enumerate(new_nads)
                if nad > dist_index
                    new_nads[i] -=1
                end
            end
            prob.nads = new_nads
        else
            dist_index += 1
        end
    end
    prob.dist = dist
    prob.num_points = size(dist)[1]
    return prob
end

function make_ipopt_problem(mds::MDSProblem; radius_constraint=false, lad_constraint=false, nad_constraint=false)
    # Creates an Ipopt problem for the given problem instance

    # 3 vars for eaxh point- x, y, z
    num_vars = 3*mds.num_points
    num_constraints = 0
    n, n2 = size(mds.dist)
    dist = mds.dist
    rad = mds.sphere_radius
    nuc_rad = mds.nucleolus_radius
    initial_x = [rand()*2*rad - rad for i in 1:num_vars]
    # m is the number of constraints
    m = 0
    if radius_constraint
        m = mds.num_points
    end
    #println(m)
    # constraint lower bound
    g_L = zeros(m)
    # constraint upper bound
    g_U = fill(rad^2, m)
    if lad_constraint
        for l in mds.lads
            g_L[l] = 0.9*rad^2
        end
    end
    # TODO: add constraint for nucleolus center
    # nuc_c is the index fo the x-element of the nucleolus coordinate.
    nuc_c = 3*mds.num_points + 1
    num_jac_elements=m*3
    if nad_constraint
        num_vars += 3
        m += 3 + length(mds.nads)
        append!(initial_x, [rand()*2*rad - rad for i in 1:3])
        append!(g_L, zeros(3 + length(mds.nads)))
        append!(g_U, fill((rad - nuc_rad)^2, 3))
        append!(g_U, fill(nuc_rad^2, length(mds.nads)))
        num_jac_elements += 3 + 6*length(mds.nads)
    end

    interactions = Array((Int, Int), 0)
    for i in 1:n-1
        for j in i+1:n
            if dist[i,j] != Inf
                push!(interactions, (i, j))
            end
        end
    end
    println("Number of interactions: ", length(interactions))
    println("Number of variables: ", num_vars)
    println("Number of constraints: ", m)
    println("Number of nonzero Jacobian elements: ", num_jac_elements)

    ###### Callback Functions

    function eval_f(x::Vector{Float64})
        val = 0
        for (i, j) in interactions
            i_c = 3*(i-1) + 1
            j_c = 3*(j-1) + 1
            s = (x[i_c] - x[j_c])^2 + (x[i_c+1] - x[j_c+1])^2 + (x[i_c+2] - x[j_c+2])^2
            w = 1/dist[i,j]^2
            val += w*(sqrt(s) - dist[i,j])^2
        end
        #println("x max: ", maximum(x))
        return val
    end

    function eval_g(x::Vector{Float64}, g::Vector{Float64})
        g[:] = zeros(size(g))
        if !nad_constraint
            for i in 1:m
                i_c = 3*(i-1) + 1
                g[i] = x[i_c]^2 + x[i_c+1]^2 + x[i_c+2]^2
            end
        else
            for i in 1:mds.num_points
                i_c = 3*(i-1) + 1
                g[i] = x[i_c]^2 + x[i_c+1]^2 + x[i_c+2]^2
            end
            g_index = mds.num_points+1
            g[g_index] = x[nuc_c]^2 + x[nuc_c+1]^2+x[nuc_c+2]^2
            g_index += 1
            for i in mds.nads
                i_c = 3*(i-1) + 1
                g[g_index] = (x[i_c] - x[nuc_c])^2 + (x[i_c+1] - x[nuc_c+1])^2
                    + (x[i_c+2] - x[nuc_c+2])^2
                g_index += 1
            end
        end
    end

    function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
        grad_f[:] = zeros(size(grad_f))
        for (i, j) in interactions
            # PROBLEm: we get huge gradients
            i_c = 3*(i-1) + 1
            j_c = 3*(j-1) + 1
            #s = sum((x[i_c:i_c+2] .- x[j_c:j_c+2]).^2)
            s = (x[i_c] - x[j_c])^2 + (x[i_c+1] - x[j_c+1])^2 + (x[i_c+2] - x[j_c+2])^2
            s2 = sqrt(s)
            f = s2-dist[i,j]
            if s2 > 0.0000001
                w = 1/dist[i,j]^2
                common_term = w*(1-dist[i,j]/s2)
                #if abs(common_term) > 50
                #    krintln("index: ", i, j)
                #    println("x: ", x[i_c:i_c+2], " ", x[j_c:j_c+2])
                #end
                grad_f[i_c] += common_term*2*(x[i_c] - x[j_c])
                grad_f[i_c+1] += common_term*2*(x[i_c+1] - x[j_c+1])
                grad_f[i_c+2] += common_term*2*(x[i_c+2] - x[j_c+2])
                grad_f[j_c] -= common_term*2*(x[i_c] - x[j_c])
                grad_f[j_c+1] -= common_term*2*(x[i_c+1] - x[j_c+1])
                grad_f[j_c+2] -= common_term*2*(x[i_c+2] - x[j_c+2])
            end
        end
        #println(grad_f)
        #println("grad max: ", maximum(grad_f))
        #println("grad min: ", minimum(grad_f))
        #mi = indmax(grad_f)
        #println("index of grad max: ", mi)
        #println("x max: ", maximum(x))
        #println("x min: ", minimum(x))
    end

    function eval_jac_g(
        x::Vector{Float64},         # Current solution
        mode,                       # Either :Structure or :Values
        rows::Vector{Int32},        # Sparsity structure - row indices
        cols::Vector{Int32},        # Sparsity structure - column indices
        values::Vector{Float64})
        #println("eval jac g")
        if mode == :Structure
            #println("eval jac g structure")
            if !nad_constraint
                for i in 1:m
                    i_c = 3*(i-1) + 1
                    rows[i_c] = i
                    rows[i_c+1] = i
                    rows[i_c+2] = i
                    cols[i_c] = i_c
                    cols[i_c+1] = i_c+1
                    cols[i_c+2] = i_c+2
                end
            else
                for i in 1:mds.num_points
                    i_c = 3*(i-1) + 1
                    #println(i_c)
                    rows[i_c] = i
                    rows[i_c+1] = i
                    rows[i_c+2] = i
                    cols[i_c] = i_c
                    cols[i_c+1] = i_c+1
                    cols[i_c+2] = i_c+2
                end
                constraint_index = mds.num_points + 1
                i_c = 3*(mds.num_points)+1
                rows[i_c] = constraint_index
                rows[i_c+1] = constraint_index
                rows[i_c+2] = constraint_index
                cols[i_c] = nuc_c
                cols[i_c+1] = nuc_c+1
                cols[i_c+2] = nuc_c+2
                constraint_index += 1
                i_c += 3
                for i in mds.nads
                    #println(i_c)
                    x_index = 3*(i-1)+1
                    rows[i_c] = constraint_index
                    rows[i_c+1] = constraint_index
                    rows[i_c+2] = constraint_index
                    rows[i_c+3] = constraint_index
                    rows[i_c+4] = constraint_index
                    rows[i_c+5] = constraint_index
                    cols[i_c] = x_index
                    cols[i_c+1] = x_index+1
                    cols[i_c+2] = x_index+2
                    cols[i_c+3] = nuc_c
                    cols[i_c+4] = nuc_c+1
                    cols[i_c+5] = nuc_c+2
                    i_c += 6
                    constraint_index += 1
                end
            end
        else
            if !nad_constraint
                for i in 1:m
                    i_c = 3*(i-1) + 1
                    values[i_c] = 2*x[i_c]
                    values[i_c+1] = 2*x[i_c+1]
                    values[i_c+2] = 2*x[i_c+2]
                end
            else
                for i in 1:mds.num_points
                    i_c = 3*(i-1) + 1
                    values[i_c] = 2*x[i_c]
                    values[i_c+1] = 2*x[i_c+1]
                    values[i_c+2] = 2*x[i_c+2]
                end
                #constraint_index = mds.num_points+1
                i_c = 3*(mds.num_points)+1
                values[i_c] = 2*x[nuc_c]
                values[i_c+1] = 2*x[nuc_c+1]
                values[i_c+2] = 2*x[nuc_c+1]
                i_c += 3
                for i in mds.nads
                    x_index = 3*(i-1)+1
                    values[i_c] = 2*(x[x_index]-x[nuc_c])
                    values[i_c+1] = 2*(x[x_index+1]-x[nuc_c+1])
                    values[i_c+2] = 2*(x[x_index+2]-x[nuc_c+2])
                    values[i_c+3] = -2*(x[x_index]-x[nuc_c])
                    values[i_c+4] = -2*(x[x_index+1]-x[nuc_c+1])
                    values[i_c+5] = -2*(x[x_index+2]-x[nuc_c+2])
                    i_c += 6
                end
            end
        end
        #println("finished eval jac g")
    end

    function eval_h(
        x::Vector{Float64},         # Current solution
        mode,                       # Either :Structure or :Values
        rows::Vector{Int32},        # Sparsity structure - row indices
        cols::Vector{Int32},        # Sparsity structure - column indices
        obj_factor::Float64,        # Lagrangian multiplier for objective
        lambda::Vector{Float64},    # Multipliers for each constraint
        values::Vector{Float64})    # The values of the Hessian

        if mode == :Structure
            # rows[...] = ...
            # ...
            # cols[...] = ...
        else
            # values[...] = ...
        end
    end

    function intermediate(
        alg_mod::Int,
        iter_count::Int,
        obj_value::Float64,
        inf_pr::Float64, inf_du::Float64,
        mu::Float64, d_norm::Float64,
        regularization_size::Float64,
        alpha_du::Float64, alpha_pr::Float64,
        ls_trials::Int)
        # ...
        return true  # Keep going
    end

    problem = createProblem(
        num_vars::Int,                     # Number of variables
        fill(-rad, num_vars)::Vector{Float64},       # Variable lower bounds
        fill(rad, num_vars)::Vector{Float64},       # Variable upper bounds
        m::Int,                     # Number of constraints
        g_L,       # Constraint lower bounds
        g_U,       # Constraint upper bounds
        num_jac_elements,              # Number of non-zeros in Jacobian
        0,             # Number of non-zeros in Hessian
        eval_f,                     # Callback: objective function
        eval_g,                     # Callback: constraint evaluation
        eval_grad_f,                # Callback: objective function gradient
        eval_jac_g,                 # Callback: Jacobian evaluation
        nothing)           # Callback: Hessian evaluation
    addOption(problem, "hessian_approximation", "limited-memory")
    addOption(problem, "tol", 1e-1)
    addOption(problem, "acceptable_tol", 1e-1)
    addOption(problem, "max_iter", 500)
    problem.x = initial_x
    return problem

end


text:
module cdtest_

using IUP
using IUP_IM
using IUP_CD

include("../src/libiup_h.jl")
include("../src/im_process_h.jl")
include("../src/cdiup_h.jl")

include("cdtest_h.jl")
#include("cdtest.jl")

export
  cdtest

type tCTC_l
    iup_canvas::Ptr{cdCanvas}
end

global antialias = 1

global ctgc = tCTC(
        C_NULL,
        int32(0), int32(0), 0.0, int32(0),
        C_NULL,
        int32(0),
        C_NULL,
        C_NULL,
        int32(0),
        C_NULL,
        int32(0),
        C_NULL,
        C_NULL,
        uint32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        0.0,
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
#       Array_100_Cuchar,
#        Array_100_Clong,
#        Array_4_Cint,
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        int32(0),
        C_NULL,
        C_NULL,
        C_NULL,
        int32(0),
        int32(0),
        C_NULL,
        int32(0),
#        Array_300_tPoint,
#        Array_256_Uint8,
#        Array_80_Uint8,
        int32(0),
        int32(0),
#        Array_40_Uint8,
        int32(0),
        int32(0)
#        Ptr{tList}
    )

global box_pos   = tBoxPos(int32(0), int32(0), int32(0), int32(0), int32(0), int32(0))
global arc_pos   = tArcPos(int32(0), int32(0), int32(0), int32(0), 0.0, 0.0)
global line_pos  = tLinePos(int32(0), int32(0), int32(0), int32(0))
global mark_pos  = tMarkPos(int32(0), int32(0), int32(0))
global pixel_pos = tPixelPos(int32(0), int32(0))
global fl_nx, fl_ny, fl_mx, fl_my
global line_x1, line_x2, line_y1, line_y2
global line_lastwhat = CLOSE
global poly_lastwhat = CLOSE
global poly_x1, poly_x2, poly_y1, poly_y2
global box_x1, box_x2, box_y1, box_y2
global box_lastwhat = CLOSE
global arc_xc, arc_yc, arc_y1, arc_x1
global arc_lastwhat = CLOSE

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

	# carrega o LED
	err = IupLoad("C:/programs/Gits/IUP.jl/examples/cdtest.led")
	#cdtest_loadled()

	if (err == C_NULL)	# inicializa o CDTest
		CDTestInit()
	else			# imprime uma mensagem de erro
		IupMessage("LED Error:", err)
		IupClose()
		return
	end

	IupMainLoop()
	#CDTestClose()
	IupClose()
end

#-------------------------------------------------------------------------
# Inicializa o CD Test.
#-------------------------------------------------------------------------
function CDTestInit()

	# inicializaccao dos drivers
	#DriversInit()

	# mostra o dialogo principal
	IupShow(IupGetHandle("dlgMain"))

	# inicializacao da barra de cores */
	ctgc.foreground = int32(0);		ctgc.background = int32(0)
	#ColorBarInit(IupGetHandle("dlgMain"), IupGetHandle("cnvColorBar"), convert(Ptr{Clong},ctgc.foreground), convert(Ptr{Clong},ctgc.background))

	# cria o canvas do CD associado ao canvas do IUP */
	ctgc.iup_canvas = cdCreateCanvas(cdContextIup(), IupGetHandle("cnvMain"))

w = convert(Ptr{Cint}, [int32(0)])
h = convert(Ptr{Cint}, [int32(0)])
cdCanvasGetSize(ctgc.iup_canvas, w, h, C_NULL, C_NULL);
w = unsafe_load(w)
h = unsafe_load(h)
@show(w,h)
	if (antialias != 0) cdCanvasSetAttribute(ctgc.iup_canvas, "ANTIALIAS", "0")	end
	cdActivate(ctgc.iup_canvas);

	# associa os call-backs
	setcallbacks()

	# os call-backs do canvas devem ser associados depois de sua criacao */
	#IupSetFunction("cmdRepaint", (Icallback) fRepaint);
	IupSetFunction("cmdMotionCB", cfunction(fMotionCB, Int, (Ptr{Ihandle}, Cint, Cint, Ptr{Uint8})))
	IupSetFunction("cmdButtonCB", cfunction(fButtonCB, Int, (Ptr{Ihandle}, Char, Cint, Cint, Cint, Cint)))
	#IupSetFunction("cmdResizeCB", (Icallback) fResizeCB);
	#IupSetFunction("cmdGetFocusCB", (Icallback) fGetFocusCB);

	# inicializacao do contexto */
	#ctgc.write_mode = 0;
	#ctgc.line_style = 0;
	#ctgc.line_cap = 0;
	#ctgc.line_join = 0;
	ctgc.fill_mode = IUP_CD.CD_EVENODD;
	ctgc.line_width = 1;
	ctgc.font_style = IUP_CD.CD_PLAIN;
	ctgc.font_typeface = IUP_CD.CD_SYSTEM;
	ctgc.font_size = IUP_CD.CD_STANDARD;
	ctgc.text_alignment = IUP_CD.CD_BASE_LEFT;
	#ctgc.text_orientation = 0
	ctgc.back_opacity =  IUP_CD.CD_TRANSPARENT;
	ctgc.mark_type = IUP_CD.CD_STAR;
	ctgc.poly_mode = IUP_CD.CD_OPEN_LINES;
	ctgc.interior_style = IUP_CD.CD_SOLID;
	ctgc.hatch = IUP_CD.CD_HORIZONTAL;
	#ctgc.following = 0
	ctgc.foreground = IUP_CD.CD_BLACK
	ctgc.background = IUP_CD.CD_WHITE
#	ctgc.head = C_NULL;
	#ctgc.test_image = C_NULL;
	#ctgc.sim = 0;
	ctgc.stretch_play = 1;
	ctgc.dlg_x = IUP_CENTER;
	ctgc.dlg_y = IUP_CENTER;
	#ctgc.visible = 0;

	# inicializa os vetores stipple e pattern */
	#initsamples();

	# inicializa o CDTest com a primitiva LINE */
	ctgc.dlg_cur_prim = IupGetHandle("dlgLB");
	ctgc.bt_cur_prim = IupGetHandle("btCurPrim");
	ctgc.cur_prim = LINE;

	# inicializaccao do Canvas do IUP */
	cdActivate(ctgc.iup_canvas)
	cdFont(ctgc.font_typeface,ctgc.font_style,ctgc.font_size);
	cdBackground(ctgc.background);
	cdClear()

	w = convert(Ptr{Cint}, [ctgc.w])
	h = convert(Ptr{Cint}, [ctgc.h])
	cdGetCanvasSize(w, h, convert(Ptr{Cdouble}, C_NULL), convert(Ptr{Cdouble}, C_NULL))
	ctgc.w = unsafe_load(w)
	ctgc.h = unsafe_load(h)
	ctgc.bpp = cdGetColorPlanes()

	mm = convert(Ptr{Cdouble}, [0.0])
	cdPixel2MM(1, 0, mm, convert(Ptr{Cdouble}, C_NULL))
	mm = unsafe_load(mm)
	xres = 1.0/mm;
	ctgc.res = xres;

	# inicializa o canvas off-screen de double-bufering */
	#ctgc.buffering = int32(0)
	#ctgc.buffer_canvas = C_NULL

	# inicializa o clipping
	#ctgc.clip_mode = IUP_CD.CD_CLIPOFF

	x1 = convert(Ptr{Cint}, [ctgc.clip_xmin]);	x2 = convert(Ptr{Cint}, [ctgc.clip_xmax])
	y1 = convert(Ptr{Cint}, [ctgc.clip_ymin]);	y2 = convert(Ptr{Cint}, [ctgc.clip_ymax])
	cdGetClipArea(x1, x2, y1, y2)
	ctgc.clip_xmin = unsafe_load(x1);	ctgc.clip_xmax = unsafe_load(x2)
	ctgc.clip_ymin = unsafe_load(y1);	ctgc.clip_ymax = unsafe_load(y2)
	#cdGetClipArea(&(ctgc.clip_xmin), &(ctgc.clip_xmax), &(ctgc.clip_ymin), &(ctgc.clip_ymax));

	# zera o buffer temporario de pontos */
	#ctgc.num_points = Cint(0)

	# a lista de primitivas jah estah vazia, nao hah o que apagar */
	IupSetAttribute(IupGetHandle("itEditUndo"), IUP_ACTIVE, IUP_NO);

	# atualiza o tamanho do canvas em pixels na barra de titulo */
	#s = @sprintf(ctgc.title, "CDTest 5.3 (%dx%d - %dbpp)", ctgc.w, ctgc.h, ctgc.bpp)
	#IupSetAttribute(IupGetHandle("dlgMain"), IUP_TITLE, ctgc.title);

	# inicializa a barra de status */
	#s = @sprintf(ctgc.status_line, "LEFT click and drag.");
	#set_status();

	# inicializa a posiccao do mouse */
	mouse_pos(int32(0), int32(0))

	# constroi os dialogos do CDTest sem mostra-los */
	IupMap(IupGetHandle("dlgLB"));
	IupMap(IupGetHandle("dlgAS"));
	IupMap(IupGetHandle("dlgPixel"));
	IupMap(IupGetHandle("dlgImage"));
	IupMap(IupGetHandle("dlgImageRGB"));
	IupMap(IupGetHandle("dlgMark"));
	IupMap(IupGetHandle("dlgText"));
	IupMap(IupGetHandle("dlgClip"));
	IupMap(IupGetHandle("dlgAttributes"));
	IupMap(IupGetHandle("dlgWDCanvas"));
	IupMap(IupGetHandle("dlgPICCanvas"));

	# cria o canvas WD
	use_contextplus = 1
	if (use_contextplus != 0) cdUseContextPlus(1)	end
	ctgc.wd_canvas  = cdCreateCanvas(cdContextIup(), IupGetHandle("cnvWDCanvas"));
	ctgc.pic_canvas = cdCreateCanvas(cdContextIup(), IupGetHandle("cnvPICCanvas"));
	ctgc.picture    = cdCreateCanvas(cdContextIup(), convert(Ptr{Void}, convert(Ptr{Uint8},"")))
	if (use_contextplus != 0) cdUseContextPlus(0)	end

	# CDTEST default values */
	cdActivate(ctgc.picture);
	#cdLineStyleDashes(ctgc.dashes, 4);
	#cdPattern(10, 10, ctgc.pattern);
	#cdStipple(10, 10, ctgc.stipple);
	cdInteriorStyle(0);

	cdActivate(ctgc.wd_canvas);
	#cdLineStyleDashes(ctgc.dashes, 4);
	#cdPattern(10, 10, ctgc.pattern);
	#cdStipple(10, 10, ctgc.stipple);
	cdInteriorStyle(0);

	# reativa o canvas IUP */
	cdActivate(ctgc.iup_canvas);
	#cdLineStyleDashes(ctgc.dashes, 4);
	#cdPattern(10, 10, ctgc.pattern);
	#cdStipple(10, 10, ctgc.stipple);
	cdInteriorStyle(0);
end

# ----------------------------------------------------------------------------
function CDTestClose()
	#dellist()

	#ColorBarClose()

	if (ctgc.buffer_canvas != 0) cdKillCanvas(ctgc.buffer_canvas)	end
	if (ctgc.test_image != 0) cdKillImage(ctgc.test_image)	end
	cdKillCanvas(ctgc.picture);
	cdKillCanvas(ctgc.pic_canvas);
	cdKillCanvas(ctgc.wd_canvas);
	cdKillCanvas(ctgc.iup_canvas);

	IupDestroy(IupGetHandle("dlgLB"));
	IupDestroy(IupGetHandle("dlgAS"));
	IupDestroy(IupGetHandle("dlgPixel"));
	IupDestroy(IupGetHandle("dlgImage"));
	IupDestroy(IupGetHandle("dlgImageRGB"));
	IupDestroy(IupGetHandle("dlgMark"));
	IupDestroy(IupGetHandle("dlgText"));
	IupDestroy(IupGetHandle("dlgClip"));
	IupDestroy(IupGetHandle("dlgAttributes"));
	IupDestroy(IupGetHandle("dlgWDCanvas"));
	IupDestroy(IupGetHandle("dlgPICCanvas"));

	IupDestroy(IupGetHandle("dlgMain"))

	if (use_contextplus != 0)
		cdFinishContextPlus()
	end
end

# --------------------------------------------------------------------------
function setcallbacks()
	IupSetFunction("cmdFileExit",  cfunction(fFileExit, Int, (),))
	IupSetFunction("cmdHelpAbout", cfunction(fHelpAbout,Int, (),))

	IupSetFunction("cmdLine", cfunction(fLine, Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdPoly", cfunction(fPoly, Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdRect", cfunction(fRect, Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdBox",  cfunction(fBox,  Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdArc",  cfunction(fArc,  Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdSector", cfunction(fSector,  Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdText", cfunction(fText, Int, (Ptr{Ihandle},)))
	IupSetFunction("cmdMark", cfunction(fMark, Int, (Ptr{Ihandle},)))

	IupSetFunction("cmdDraw", cfunction(fDraw, Int, (),))
end

# Tchau
function fFileExit()
	global ctgc
	IupHide(ctgc.dlg_cur_prim)
	return IUP_CLOSE
end

# Show About window
function fHelpAbout()
	IupSetAttribute(IupGetHandle("lblVersion"), IUP_TITLE, cdVersion())
	IupShow(IupGetHandle("dlgHelpAbout"))
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------*/
# Muda a primitiva corrente.                                              */
#-------------------------------------------------------------------------*/
function fLine(self::Ptr{Ihandle})
@show("MERDA_fLine")
	if (ctgc.cur_prim != LINE)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgLine");
		if (ctgc.cur_prim == POLY)	# termina o poligono em andamento */
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = LINE;
		ctgc.following = 0;
		ctgc.visible = iscurvisible();
		IupHide(ctgc.dlg_cur_prim);      # esconde o dialogo anterior */
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X);
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y);
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgLB");
		IupSetAttribute(ctgc.dlg_cur_prim, IUP_TITLE, "Line Parameters");
	end
	#sprintf(ctgc.status_line, "LEFT click and drag.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------
function fPoly(self::Ptr{Ihandle})
	if (ctgc.cur_prim != POLY)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgPoly");
		ctgc.cur_prim = POLY;
		ctgc.following = 0
		ctgc.visible = iscurvisible()
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X);
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y);
		end
    	ctgc.dlg_cur_prim = IupGetHandle("dlgPoly");
    	# zera o buffer temporario de pontos
    	ctgc.num_points = 0;
	end
	#sprintf(ctgc.status_line, "LEFT click to add vertex. RIGHT click to end.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------*/
# Muda a primitiva corrente.                                              */
#-------------------------------------------------------------------------*/
function fRect(self::Ptr{Ihandle})
	if (ctgc.cur_prim != RECT)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgRect")
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = RECT;
		ctgc.following = 0;
		ctgc.visible = iscurvisible();
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X);
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y);
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgLB");
		IupSetAttribute(ctgc.dlg_cur_prim, IUP_TITLE, "Rect Parameters");
	end
	#sprintf(ctgc.status_line, "LEFT click and drag.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------
function fBox(self::Ptr{Ihandle})
	if (ctgc.cur_prim != BOX)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgBox")
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = BOX;
		ctgc.following = 0;
		ctgc.visible = iscurvisible();
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X)
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y)
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgLB");
		IupSetAttribute(ctgc.dlg_cur_prim, IUP_TITLE, "Box Parameters");
	end
	#sprintf(ctgc.status_line, "LEFT click and drag.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------
function fArc(self::Ptr{Ihandle})
	if (ctgc.cur_prim != ARC)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgArc");
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = ARC
		ctgc.following = 0
		ctgc.visible = iscurvisible();
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X);
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y);
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgAS");
		IupSetAttribute(ctgc.dlg_cur_prim, IUP_TITLE, "Arc Parameters");
	end
	#sprintf(ctgc.status_line, "LEFT click at center and drag.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------
function fSector(self::Ptr{Ihandle})

	if (ctgc.cur_prim != SECTOR)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgSector");
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = SECTOR
		ctgc.following = 0
		ctgc.visible = iscurvisible()
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X)
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y)
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgAS");
		IupSetAttribute(ctgc.dlg_cur_prim, IUP_TITLE, "Sector Parameters");
	end
	#sprintf(ctgc.status_line, "LEFT click at center and drag.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
function fText(self::Ptr{Ihandle})

	if (ctgc.cur_prim != TEXT)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgText");
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0);
		end
		ctgc.cur_prim = TEXT;
		ctgc.following = false;
		ctgc.visible = iscurvisible();
		IupHide(ctgc.dlg_cur_prim);
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X);
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y);
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgText");
	end
	#sprintf(ctgc.status_line, "LEFT click.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
function fMark(self::Ptr{Ihandle})

	if (ctgc.cur_prim != MARK)
		IupSetAttribute(ctgc.bt_cur_prim, IUP_IMAGE, "imgMark")
		if (ctgc.cur_prim == POLY)
			fButtonCB(C_NULL, IUP_BUTTON3, 0, 0, 0, 0)
		end
		ctgc.cur_prim = MARK
		ctgc.following = false
		ctgc.visible = iscurvisible()
		IupHide(ctgc.dlg_cur_prim)
		if (ctgc.visible != 0)
			ctgc.dlg_x = IupGetInt(ctgc.dlg_cur_prim, IUP_X)
			ctgc.dlg_y = IupGetInt(ctgc.dlg_cur_prim, IUP_Y)
		end
		ctgc.dlg_cur_prim = IupGetHandle("dlgMark")
	end
	#sprintf(ctgc.status_line, "LEFT click.");
	#set_status();
	if (ctgc.visible != 0 || IupGetAttribute(self, "ISMENU") != C_NULL)
		IupShowXY(ctgc.dlg_cur_prim, ctgc.dlg_x, ctgc.dlg_y)
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
# Desenha a primitiva entrada na caixa de dialogo corrente.
#-------------------------------------------------------------------------
function fDraw()
	# atualiza os dados entrados na caixa de dialogo no contexto da primitiva corrente
	if (ctgc.cur_prim == LINE)
		line_pos.x1 = IupGetInt(IupGetHandle("txtLBX1"), IUP_VALUE)
		line_pos.x2 = IupGetInt(IupGetHandle("txtLBX2"), IUP_VALUE)
		line_pos.y1 = IupGetInt(IupGetHandle("txtLBY1"), IUP_VALUE)
		line_pos.y2 = IupGetInt(IupGetHandle("txtLBY2"), IUP_VALUE)
	elseif (ctgc.cur_prim == RECT || ctgc.cur_prim == BOX)
		box_pos.xmin = IupGetInt(IupGetHandle("txtLBX1"), IUP_VALUE)
		box_pos.xmax = IupGetInt(IupGetHandle("txtLBX2"), IUP_VALUE)
		box_pos.ymin = IupGetInt(IupGetHandle("txtLBY1"), IUP_VALUE)
		box_pos.ymax = IupGetInt(IupGetHandle("txtLBY2"), IUP_VALUE)
	elseif (ctgc.cur_prim == ARC || ctgc.cur_prim == SECTOR || ctgc.cur_prim == CHORD)
		arc_pos.xc = IupGetInt(IupGetHandle("txtASXC"), IUP_VALUE)
		arc_pos.yc = IupGetInt(IupGetHandle("txtASYC"), IUP_VALUE)
		arc_pos.w = IupGetInt(IupGetHandle("txtASW"), IUP_VALUE)
		arc_pos.h = IupGetInt(IupGetHandle("txtASH"), IUP_VALUE)
	elseif (ctgc.cur_prim == PIXEL)
		pixel_pos.x = IupGetInt(IupGetHandle("txtPixelX"), IUP_VALUE)
		pixel_pos.y = IupGetInt(IupGetHandle("txtPixelY"), IUP_VALUE)
	elseif (ctgc.cur_prim == MARK)
		mark_pos.x = IupGetInt(IupGetHandle("txtMarkX"), IUP_VALUE)
		mark_pos.y = IupGetInt(IupGetHandle("txtMarkY"), IUP_VALUE)
	end

	# efetivamente desenha a primitiva
	draw()

	return IUP_DEFAULT
end

#-------------------------------------------------------------------------*/
# Funcao associada ao botao do mouse.                                     */
#-------------------------------------------------------------------------*/
function fButtonCB(self::Ptr{Ihandle}, b::Char, e::Integer, x::Integer, y::Integer, r::Integer)

	cdActivate(ctgc.iup_canvas)
	cdWriteMode(IUP_CD.CD_NOT_XOR)
	cdForeground(convert(Clong, IUP_CD.CD_BLACK))
	cdLineStyle(IUP_CD.CD_CONTINUOUS)
	cdLineWidth(1)
	cdClip(IUP_CD.CD_CLIPOFF)

	yp = convert(Ptr{Cint}, [y]);
	cdUpdateYAxis(yp)
	y = int32(unsafe_load(yp))
	x = int32(x)				# <=========== NAO DEVIA SER Int64
	mouse_pos(x, y)

	if (b == IUP_BUTTON1)
		if (e != 0)
			if (ctgc.cur_prim == LINE)
				follow(x, y)
				line(NEWPOINT, x, y)
				line_pos.x1 = x
				line_pos.y1 = y
				ctgc.following = 1
			elseif ((ctgc.cur_prim == RECT) || (ctgc.cur_prim == BOX))
				follow(x, y)
				box(NEWPOINT, x, y)
				box_pos.x = x
				box_pos.y = y
				ctgc.following = 1
			elseif (ctgc.cur_prim == CLIP)
				box(NEWPOINT, x, y)
				follow(x, y)
				ctgc.following = 1
			elseif (ctgc.cur_prim == ARC || ctgc.cur_prim == SECTOR || ctgc.cur_prim == CHORD)
				follow(x, y);
				arc(CENTER, x, y);
				arc_pos.xc = x;
				arc_pos.yc = y;
				ctgc.following = true;
			elseif (ctgc.cur_prim == MARK || ctgc.cur_prim == PIXEL || ctgc.cur_prim == TEXT)
				follow(x, y)
				draw()
			end
		else
			if (ctgc.cur_prim == LINE)
				if (ctgc.following != 0)
					ctgc.following = 0
					line(CLOSE, x, y)
					cdClip(ctgc.clip_mode)
					draw()
				end
			elseif (ctgc.cur_prim == RECT || ctgc.cur_prim == BOX)
				if (ctgc.following != 0)
					ctgc.following = 0
					box(CLOSE, x, y)
					cdClip(ctgc.clip_mode)
					draw()
				end
			elseif (ctgc.cur_prim == CLIP)
				if (ctgc.following != 0)
					ctgc.following = 0
					box(CLOSE, x, y)
					fClipArea()
				end
			elseif (ctgc.cur_prim == ARC || ctgc.cur_prim == SECTOR || ctgc.cur_prim == CHORD)
				if (ctgc.following != 0)
					ctgc.following = false
					arc(CLOSE, x, y)
					cdClip(ctgc.clip_mode)
					draw()
					arc_pos.w = 0
					arc_pos.h = 0
				end
			elseif (ctgc.cur_prim == POLY)
				cdActivate(ctgc.iup_canvas);
				polygon(NEWPOINT, x, y);
				#newpolypoint(x, y);
				#sprintf(ctgc.status_line, "cdVertex( %d, %d)", x, y);
				#set_status();
				ctgc.following = 1;
 			end
 		end
	elseif (b == IUP_BUTTON3)
		if (e != 0)
			if (ctgc.cur_prim == IMAGE)
			elseif (ctgc.cur_prim == RGB)
			elseif (ctgc.cur_prim == CLIP)
			end
		else
			if (ctgc.cur_prim == POLY)
				cdActivate(ctgc.iup_canvas)
				ctgc.following = 0;
				polygon(CLOSE, x, y)
				#sprintf(ctgc.status_line, "cdEnd()");
				#set_status();
				fDraw()
			end
		end
	end

	return IUP_DEFAULT
end

# ------------------------------------------------------------------------
function iscurvisible()
	vis = IupGetAttribute(ctgc.dlg_cur_prim, IUP_VISIBLE)
	if (vis == C_NULL)
		return 0
	end
	return (bytestring(vis) == "YES") ? 1 : 0
end

# ------------------------------------------------------------------------
function fMotionCB(self::Ptr{Ihandle}, x::Int32, y::Int32, r::Ptr{Uint8})

	if (ctgc.iup_canvas == C_NULL)
		return IUP_DEFAULT
	end

	cdActivate(ctgc.iup_canvas)
	yp = convert(Ptr{Cint}, [y]);
	cdUpdateYAxis(yp)
	y = int32(unsafe_load(yp))
	mouse_pos(x, y)

	if (ctgc.following != 0)
		if (ctgc.cur_prim == LINE)
			line(MOVE, x, y)
			follow(x, y)
		elseif (ctgc.cur_prim == RECT || ctgc.cur_prim == BOX)
			box(MOVE, x, y)
			follow(x, y)
		elseif (ctgc.cur_prim == CLIP)
		elseif (ctgc.cur_prim == ARC || ctgc.cur_prim == SECTOR || ctgc.cur_prim == CHORD)
			arc(MOVE, x, y)
			follow(x, y)
		elseif (ctgc.cur_prim == IMAGE)
		elseif (ctgc.cur_prim == RGB)
		elseif (ctgc.cur_prim == POLY)
			polygon(MOVE, x, y)
		end
	end
	return IUP_DEFAULT
end

#-------------------------------------------------------------------------
# Atualiza a posiccao do mouse no dialogo principal.
#-------------------------------------------------------------------------
function mouse_pos(x::Int32, y::Int32)
	# salva a posiccao do mouse no contexto
	ctgc.x = x
	ctgc.y = y

	# atualiza a posiccao do mouse no feedback para o usuario
	s = @sprintf("%4d,%4d", x, y)
	IupSetAttribute(IupGetHandle("lbMousePos"), "TITLE", s);
end

#-------------------------------------------------------------------------*/
# Segue as coordenadas do mouse atualizando o TXT apropriado.             */
#-------------------------------------------------------------------------*/
function follow(x::Integer, y::Integer)
	fl_nx = @sprintf("%d", x)
	fl_ny = @sprintf("%d", y)

	if (ctgc.cur_prim == PIXEL)
		# atualiza os parametros do pixel */
		pixel_pos.x = x;
		pixel_pos.y = y;
		# atualiza a caixa de dialogo */
		fl_mx = @sprintf("%d", x);
		fl_my = @sprintf("%d", y);
		IupSetAttribute(IupGetHandle("txtPixelX"), IUP_VALUE, fl_mx);
		IupSetAttribute(IupGetHandle("txtPixelY"), IUP_VALUE, fl_my);
	elseif (ctgc.cur_prim == MARK)
		# atualiza os parametros da mark */
		mark_pos.x = x;
		mark_pos.y = y;
		# atualiza a caixa de dialogo */
		fl_mx = @sprintf("%d", x)
		fl_my = @sprintf("%d", y)
		IupSetAttribute(IupGetHandle("txtMarkX"), IUP_VALUE, fl_mx)
		IupSetAttribute(IupGetHandle("txtMarkY"), IUP_VALUE, fl_my)
	elseif (ctgc.cur_prim == RECT || ctgc.cur_prim == BOX)
		# atualiza os parametros da box
		if (ctgc.following != 0)
			if (x < box_pos.x)
				box_pos.xmin = x
				box_pos.xmax = box_pos.x
			else
				box_pos.xmax = x
				box_pos.xmin = box_pos.x
			end

			if (y < box_pos.y)
				box_pos.ymin = y
				box_pos.ymax = box_pos.y
			else
				box_pos.ymax = y
				box_pos.ymin = box_pos.y
			end
		else
			box_pos.xmax = box_pos.xmin = x
			box_pos.ymax = box_pos.ymin = y
		end
		# atualiza a caixa de dialogo
		fl_mx = @sprintf("%d", box_pos.xmin);
		fl_nx = @sprintf("%d", box_pos.xmax);
		fl_my = @sprintf("%d", box_pos.ymin);
		fl_ny = @sprintf("%d", box_pos.ymax);
		IupSetAttribute(IupGetHandle("txtLBX1"), IUP_VALUE, fl_mx);
		IupSetAttribute(IupGetHandle("txtLBX2"), IUP_VALUE, fl_nx);
		IupSetAttribute(IupGetHandle("txtLBY1"), IUP_VALUE, fl_my);
		IupSetAttribute(IupGetHandle("txtLBY2"), IUP_VALUE, fl_ny);
	elseif (ctgc.cur_prim == LINE)
		fl_mx = @sprintf("%d", x);
		fl_my = @sprintf("%d", y);
		line_pos.x2 = x;
		line_pos.y2 = y;
		if (ctgc.following != 0)
			IupSetAttribute(IupGetHandle("txtLBX2"), IUP_VALUE, fl_mx);
			IupSetAttribute(IupGetHandle("txtLBY2"), IUP_VALUE, fl_my);
		else
			line_pos.x1 = x;
			line_pos.y1 = y;
			IupSetAttribute(IupGetHandle("txtLBX1"), IUP_VALUE, fl_mx);
			IupSetAttribute(IupGetHandle("txtLBX2"), IUP_VALUE, fl_mx);
			IupSetAttribute(IupGetHandle("txtLBY1"), IUP_VALUE, fl_my);
			IupSetAttribute(IupGetHandle("txtLBY2"), IUP_VALUE, fl_my);
		end
	elseif (ctgc.cur_prim == ARC || ctgc.cur_prim == SECTOR || ctgc.cur_prim == CHORD)
		if (ctgc.following != 0)
			# atualiza os parametros do arc
			arc_pos.w = 2*abs(arc_pos.xc-x+1)
			arc_pos.h = 2*abs(arc_pos.yc-y+1)
			# atualiza a caixa de dialogo
			fl_mx = @sprintf("%d", arc_pos.w)
			fl_my = @sprintf("%d", arc_pos.h)
			IupSetAttribute(IupGetHandle("txtASW"), IUP_VALUE, fl_mx)
			IupSetAttribute(IupGetHandle("txtASH"), IUP_VALUE, fl_my)
		else
			arc_pos.xc = x;
			arc_pos.xc = y;
			fl_mx = @sprintf("%d", x)
			fl_my = @sprintf("%d", y)
			IupSetAttribute(IupGetHandle("txtASXC"), IUP_VALUE, fl_mx);
			IupSetAttribute(IupGetHandle("txtASYC"), IUP_VALUE, fl_my);
		end
	elseif (ctgc.cur_prim == TEXT)
		IupSetAttribute(IupGetHandle("txtTextX"), IUP_VALUE, fl_nx)
		IupSetAttribute(IupGetHandle("txtTextY"), IUP_VALUE, fl_ny)
    end
end

#-------------------------------------------------------------------------*/
# Desenha uma linha em rubber band.                                       */
#-------------------------------------------------------------------------*/
function line(what::Integer, x::Integer, y::Integer)
# I have to have those redeclared here otherwise it errors saying they are not declared ??????????????????
	global line_lastwhat
	global line_x1, line_x2, line_y1, line_y2

	if (what == NEWPOINT)
		line_x1 = x;                 # novo segmento comeca no... */
		line_x2 = x;                 # novo segmento comeca no... */
		line_y1 = y;                 # ...fim do primeiro */
		line_y2 = y;                 # ...fim do primeiro */
	elseif (what == MOVE)
		if (line_lastwhat == MOVE)
			cdLine(line_x1, line_y1, line_x2, line_y2);  # apaga o segmento velho */
		end
		cdLine(line_x1, line_y1, x, y);     # desenha o novo */
		line_x2 = x;                        # o novo se... */
		line_y2 = y;                        # ...torna velho */
	elseif (what == REPAINT)
		cdLine(line_x1, line_y1, line_x2, line_y2);        # recupera o segmento perdido */
		return                        # nao modifica lastwhat */
	elseif (what == CLOSE)
		cdLine(line_x1, line_y1, line_x2, line_y2);      # apaga o ultimo segmento */
	end
	line_lastwhat = what
end

#-------------------------------------------------------------------------*/
# Desenha o poligono em rubber band.                                      */
#-------------------------------------------------------------------------*/
function polygon(what::Integer, x::Integer, y::Integer)
# I have to have those redeclared here otherwise it errors saying they are not declared ??????????????????
global poly_lastwhat
global poly_x1, poly_x2, poly_y1, poly_y2
	if (what == NEWPOINT)
		if (poly_lastwhat != CLOSE)
			cdLine(poly_x1, poly_y1, poly_x2, poly_y2);      # ...apaga a anterior e... */
			cdLine(poly_x1, poly_y1, x, y);        # desenha a definitiva */
		end
		poly_x1 = x;                        # novo segmento comeca no... */
		poly_y1 = y;                        # fim do primeiro */
	elseif (what == MOVE)
		if (poly_lastwhat == MOVE)
			cdLine(poly_x1, poly_y1, poly_x2, poly_y2);      # apaga o segmento velho */
		end
		cdLine(poly_x1, poly_y1, x, y);     # desenha o novo
		poly_x2 = x;                        # o novo se...
		poly_y2 = y;                        # torna velho
	elseif (what == REPAINT)
		cdLine(poly_x1, poly_y1, poly_x2, poly_y2);        # recupera o segmento perdido */
		return;                             # nao modifica lastwhat
	elseif (what == CLOSE)
		if (poly_lastwhat != CLOSE)
			cdLine(poly_x1, poly_y1, poly_x2, poly_y2);        # apaga o ultimo segmento */
			# apaga o poligono temporario inteiro */
			for (i = 1:ctgc.num_points-1)
				cdLine(ctgc.points[i].x, ctgc.points[i].y, ctgc.points[i+1].x, ctgc.points[i+1].y);
			end
		end
	end
	poly_lastwhat = what
end

#-------------------------------------------------------------------------*/
# Desenha uma caixa em rubber band.                                       */
#-------------------------------------------------------------------------*/
function box(what::Integer, x::Integer, y::Integer)
	global box_x1, box_x2, box_y1, box_y2
	global box_lastwhat

	if (what == NEWPOINT)
		box_x1 = x;                 # novo segmento comeca no... */
		box_x2 = x;                 # novo segmento comeca no... */
		box_y1 = y;                 # ...fim do primeiro */
		box_y2 = y;                 # ...fim do primeiro */
	elseif (what == MOVE)
		if (box_lastwhat == MOVE)
			frame(box_x1, box_y1, box_x2, box_y2);       # apaga a caixa anterior */
		end
		frame(box_x1, box_y1, x, y); # desenha a nova
		box_x2 = x;                  # o novo se...
		box_y2 = y;                  # torna velho
	elseif (what == REPAINT)
		frame(box_x1, box_y1, box_x2, box_y2);         # restaura a caixa perdida
		return;                      # nao modifica lastwhat */
	elseif (what == CLOSE)
		#frame(box_x1, box_y1, box_x2, box_y2);         # apaga a caixa definitiva
	end
	box_lastwhat = what
end

#-------------------------------------------------------------------------
# Desenha uma caixa centrada, em rubber band.
#-------------------------------------------------------------------------
function arc(what::Integer, x::Integer, y::Integer)
	global arc_xc, arc_yc, arc_y1, arc_x1
	global arc_lastwhat

	if (what == CENTER)
		arc_x1 = x;                 # novo segmento comeca no... */
		arc_xc = x;                 # novo segmento comeca no... */
		arc_y1 = y;                 # ...fim do primeiro */
		arc_yc = y;                 # ...fim do primeiro */
	elseif (what == MOVE)
		if (arc_lastwhat == MOVE)
			cdArc(arc_xc, arc_yc, 2*abs(arc_xc-arc_x1+1), 2*abs(arc_yc-arc_y1+1), 0., 360.);
		end
		cdArc(arc_xc, arc_yc, 2*abs(arc_xc-x+1), 2*abs(arc_yc-y+1), 0., 360.);
		arc_x1 = x;                        # o novo se... */
		arc_y1 = y;                        # torna velho */
	elseif (what == REPAINT)
		cdArc(arc_xc, arc_yc, 2*abs(arc_xc-arc_x1+1), 2*abs(arc_yc-arc_y1+1), 0., 360.);
		return                        # nao modifica lastwhat
	elseif (what == CLOSE)
		cdArc(arc_xc, arc_yc, 2*abs(arc_xc-arc_x1+1), 2*abs(arc_yc-arc_y1+1), 0., 360.);
	end
  	arc_lastwhat = what
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Desenha uma caixa vazia
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function  frame(x1::Cint, y1::Cint, x2::Cint, y2::Cint)
	cdLine(x1, y1, x1, y2)
	cdLine(x2, y1, x2, y2)
	cdLine(x1, y1, x2, y1)
	cdLine(x1, y2, x2, y2)
end

#-------------------------------------------------------------------------*/
# Desenha a primitiva entrada na caixa de dialogo ativa.                  */
#-------------------------------------------------------------------------*/
function draw()

	IupSetAttribute(IupGetHandle("itEditUndo"), IUP_ACTIVE, IUP_YES);
	ctgc.line_width = IupGetInt(IupGetHandle("txtLineWidth"), IUP_VALUE);
	if (ctgc.line_width < 1) ctgc.line_width = 1;	end
	ctgc.font_size = IupGetInt(IupGetHandle("txtFontSize"), IUP_VALUE);

	# escolhe entre o canvas na tela e o off-screen */
	if (ctgc.buffering != 0)
		cdActivate(ctgc.iup_canvas);
		cdClip(CD_CLIPOFF);
		cdActivate(ctgc.buffer_canvas)
	else
		cdActivate(ctgc.iup_canvas)
	end

	# atualiza os atributos */
	cdWriteMode(ctgc.write_mode);
	cdLineStyle(ctgc.line_style);
	cdLineCap(ctgc.line_cap);
	cdLineJoin(ctgc.line_join);
	cdLineWidth(ctgc.line_width);
	cdForeground(ctgc.foreground);
	cdBackground(ctgc.background);
	cdHatch(ctgc.hatch);
	cdFillMode(ctgc.fill_mode);
	cdInteriorStyle(ctgc.interior_style);
	cdBackOpacity(ctgc.back_opacity);

	cdClipArea(ctgc.clip_xmin, ctgc.clip_xmax, ctgc.clip_ymin, ctgc.clip_ymax)
	cdClip(ctgc.clip_mode);

	if (ctgc.cur_prim == LINE)
		# atualiza a linha de status */
		#sprintf(ctgc.status_line, "cdLine( %d, %d, %d, %d)", line_pos.x1,
		#	line_pos.y1, line_pos.x2, line_pos.y2);
		#set_status();
		# desenha a line na tela */
		cdLine(line_pos.x1, line_pos.y1, line_pos.x2, line_pos.y2);
		# arquiva a line */
		#newline(line_pos.x1, line_pos.y1, line_pos.x2, line_pos.y2);
	elseif (ctgc.cur_prim == RECT)
		# atualiza a linha de status */
		#sprintf(ctgc.status_line,"cdRect( %d, %d, %d, %d)",box_pos.xmin,
		#	box_pos.xmax, box_pos.ymin, box_pos.ymax);
		#set_status();
		# desenha a box na tela */
		cdRect(box_pos.xmin, box_pos.xmax, box_pos.ymin, box_pos.ymax);
		# armazena a box */
		#newrect(box_pos.xmin, box_pos.xmax, box_pos.ymin, box_pos.ymax);
	elseif (ctgc.cur_prim == BOX)
		#sprintf(ctgc.status_line,"cdBox( %d, %d, %d, %d)",box_pos.xmin,
		#	box_pos.xmax, box_pos.ymin, box_pos.ymax);
		#set_status();
		# desenha a box na tela */
		cdBox(box_pos.xmin, box_pos.xmax, box_pos.ymin, box_pos.ymax);
		# armazena a box */
		#newbox(box_pos.xmin, box_pos.xmax, box_pos.ymin, box_pos.ymax);
	elseif (ctgc.cur_prim == ARC)
		arc_pos.angle1 = IupGetFloat(IupGetHandle("txtASAngle1"),IUP_VALUE);
		arc_pos.angle2 = IupGetFloat(IupGetHandle("txtASAngle2"),IUP_VALUE);
		# atualiza a linha de status */
		#sprintf(ctgc.status_line,"cdArc( %d, %d, %d, %d, %.5G, %.5G)", arc_pos.xc,
		#	arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		#set_status();
		# desenha o arc na tela */
		cdArc(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		# armazena o arc */
		#newarc(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
	elseif (ctgc.cur_prim == SECTOR)
		arc_pos.angle1 = IupGetFloat(IupGetHandle("txtASAngle1"),IUP_VALUE);
		arc_pos.angle2 = IupGetFloat(IupGetHandle("txtASAngle2"),IUP_VALUE);
		# atualiza a linha de status */
		#sprintf(ctgc.status_line,"cdSector( %d, %d, %d, %d, %.5G, %.5G)", arc_pos.xc,
		#	arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		#set_status();
		cdSector(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		#newsector(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
	elseif (ctgc.cur_prim == CHORD)
		arc_pos.angle1 = IupGetFloat(IupGetHandle("txtASAngle1"),IUP_VALUE);
		arc_pos.angle2 = IupGetFloat(IupGetHandle("txtASAngle2"),IUP_VALUE);
		#sprintf(ctgc.status_line,"cdChord( %d, %d, %d, %d, %.5G, %.5G)", arc_pos.xc,
		#	arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		#set_status();
		cdChord(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
		#newchord(arc_pos.xc, arc_pos.yc, arc_pos.w, arc_pos.h, arc_pos.angle1, arc_pos.angle2);
	elseif (ctgc.cur_prim == PIXEL)
	elseif (ctgc.cur_prim == MARK)
		mark_pos.size = IupGetInt(IupGetHandle("txtMarkSize"),IUP_VALUE);
		#sprintf(ctgc.status_line,"cdMark( %d, %d)", mark_pos.x, mark_pos.y);
		#set_status();
		# desenha a marca na tela */
		cdMarkType(ctgc.mark_type);
		cdMarkSize(mark_pos.size);
		cdMark(mark_pos.x, mark_pos.y);
		#newmark(mark_pos.x, mark_pos.y, mark_pos.size);
	elseif (ctgc.cur_prim == TEXT)
		if (IupGetAttribute(IupGetHandle("txtTextS"),IUP_VALUE) != C_NULL)
			a=IupGetInt(IupGetHandle("txtTextX"),IUP_VALUE);
			b=IupGetInt(IupGetHandle("txtTextY"),IUP_VALUE);
			#sprintf(ctgc.status_line," cdText( %d, %d, ""%.3s""...)", a, b,
			#	IupGetAttribute(IupGetHandle("txtTextS"),IUP_VALUE));
			#set_status();
			ctgc.text_orientation = IupGetInt(IupGetHandle("txtTextOrientation"),IUP_VALUE);
			cdFont(ctgc.font_typeface,ctgc.font_style,ctgc.font_size);
			cdTextAlignment(ctgc.text_alignment);
			cdTextOrientation(ctgc.text_orientation);
			cdText(a,b,IupGetAttribute(IupGetHandle("txtTextS"),IUP_VALUE));
			#newtext(a,b,IupGetAttribute(IupGetHandle("txtTextS"),IUP_VALUE));
		end
	elseif (ctgc.cur_prim == POLY)
		if (ctgc.num_points > 1)
			cdBegin(ctgc.poly_mode);
			for (a = 1:ctgc.num_points)
				cdVertex(ctgc.points[a].x,ctgc.points[a].y);
			end
			cdEnd()
			if (ctgc.poly_mode != CD_CLIP)
				#newpoly()
			end
			ctgc.num_points = 0
		end
	end

	if (ctgc.buffering != 0)
		cdFlush()
		cdActivate(ctgc.iup_canvas);
		cdClip(ctgc.clip_mode)
	end

end

end 	# module


text:
using SafeTestsets
@safetestset "My aI tests_" begin include("myAiTests.jl") end


text:
"""
    @with(
        resource₁ = source₁,
        resource₂ = source₂,
        ...
        resourceₙ = sourceₙ,
    ) do
        use(resource₁, resource₂, ..., resourceₙ)
    end

Open resources, run the do block body, and cleanup the resources.
"""
:(ContextManagers.@with)

"""
    ContextManagers.with(sources...) do resources...
        use(resources...)
    end

Open resources, run the do block body, and cleanup the resources.
"""
ContextManagers.with

"""
    ContextManagers.maybeenter(source) -> context or nothing

Start a `context` managing the resource. Or return `nothing` when `source` does
not implement the context manager interface.

Default implementation returns `nothing`; i.e., no context manager interface.
"""
ContextManagers.maybeenter

"""
    ContextManagers.value(context) -> resource

Default implementation is a pass-through (identity) function.
"""
ContextManagers.value

"""
    ContextManagers.exit(context)
    ContextManagers.exit(context, err) -> nothing or ContextManagers.Handled()

Cleanup the `context`.  Default implementation is `close`.

Roughly speaking,

```julia
@with(resource = f(args...)) do
    use(resource)
end
```

is equivalent to

```julia
context = something(ContextManagers.maybeenter(source))
try
    resource = ContextManagers.value(context)

    use(resource)

finally
    ContextManagers.exit(context)
end
```

In the two-argument version `ContextManagers.exit(context, err)` (where `err` is
`nothing` or an `Exception`), the error can be suppressed by returning
`ContextManagers.Handled()`.
"""
ContextManagers.exit


text:
"Import filter from file"
function readFilter(filename, skipstart = 2)
	saved = readdlm(filename,',', skipstart = skipstart)
	z1  = saved[1,:] + im*saved[2,:]
	a1  = saved[3,:] + im*saved[4,:]

	return z1,a1
end

"Import weight function from file"
function readWS(filename)
	saved = readdlm(filename,',')
	intv  = saved[1,:]
	wght  = saved[2,:]

	return intv,wght
end


text:
# Quadcopter
"""
# References
[1] F. A. Goodarzi, D. Lee, and T. Lee,
“Geometric stabilization of a quadrotor UAV with a payload connected by flexible cable,”
in 2014 American Control Conference, Jun. 2014, pp. 4925–4930, doi: 10.1109/ACC.2014.6859419.
# Variables
p ∈ R^3: inertial position
v ∈ R^3: inertial velocity
R ∈ so(3): inertial-to-body rotation matrix (R_I2B)
ω ∈ R^3: angular rate of body in inertial frame
u ∈ R^4: [f, M...] where f ∈ R: total thrust, M ∈ R^3: moment
"""
Base.@kwdef struct GoodarziQuadcopter <: Quadcopter
    m = 0.5  # kg
    J = 1e-2 * Diagonal([0.557, 0.557, 1.050])  # kg m²
    J_inv = inv(J)
    g = 9.81  # m/s²
    # input limits
    dim_input = 4
    u_min = zeros(dim_input)
    u_max = (m*g) * ones(dim_input)
end

function saturate(env::GoodarziQuadcopter, u)
    @unpack u_min, u_max = env
    u_saturated = clamp.(u, u_min, u_max)
end

function input_to_force_moment(env::GoodarziQuadcopter, u)
    ν = u
end