text:
using BayesianSumProductNetworks
using StatsBase
using StatsFuns
using AxisArrays
using LinearAlgebra
using Random
using DelimitedFiles
using ArgParse
using ProgressMeter
using Logging
using SparseArrays
using Dates
using FileIO
using JLD2
using Serialization
using Dates

import BayesianSumProductNetworks.resample!

s = ArgParseSettings()
@add_arg_table s begin
    "--alpha"
        help = "Concentration parameter of sum nodes"
        default = 1.0
        arg_type = Float64
    "--beta"
        help = "Concentration parameter of product nodes"
        default = 10.0
        arg_type = Float64
    "--gamma"
        help = "Concentration parameter of leaf mixtures"
        default = 0.1
        arg_type = Float64
    "--Ksum"
        help = "Number of children under a sum node"
        default = 4
        arg_type = Int
    "--Kprod"
        help = "Number of children under a product node"
        default = 2
        arg_type = Int
    "--Kleaf"
        help = "Number of children under a product node"
        default = 10
        arg_type = Int
    "--J"
        help = "Number of partitions."
        default = 2
        arg_type = Int
    "--L"
        help = "Number of layers."
        default = 2
        arg_type = Int
    "--numsamples"
        help = "Write numsamples many samples to disk (default 100%)."
        default = 0
        arg_type = Int
    "--seed"
        help = "Random seed = seed * Int(first(dataset)) * chain."
        default = 42
        arg_type = Int
    "--chain"
        help = "Chain"
        default = 1
        arg_type = Int
    "dataset"
        help = "Dataset name, assuming naming convention of 20 terrible datasets."
        required = true
        arg_type = String
    "datafolder"
        help = "Path to folder containing datasets."
        required = true
        arg_type = String
    "outputfolder"
        help = "Path to folder in which output files should be stored."
        required = true
        arg_type = String
    "burnin"
        required = true
        arg_type = Int
    "iterations"
        help = "Number of iterations."
        required = true
        arg_type = Int
end

args = parse_args(s)

# Load the dataset
dataset = args["dataset"]
datafolder = args["datafolder"]

# results dir
resultsdir = args["outputfolder"]
if !isdir(resultsdir)
    mkpath(resultsdir)
end

donefile = joinpath(resultsdir, "done")
if isfile(donefile)
    @info "MCMC sampling already done, exiting."
    exit(-2)
end

# Output files.
burninfile = joinpath(resultsdir, "bspn_burnin.jld2")
modelfile = joinpath(resultsdir, "bspn_model.jld2")
llhfile = joinpath(resultsdir, "llh.jld2")
tracefile = joinpath(resultsdir, "trace.jld2")
nparamfile = joinpath(resultsdir, "nparam.jld2")
iterfile = joinpath(resultsdir, "iteration")

# check if results are locked and exit if this is the case
lockfile = joinpath(resultsdir, "filelock.lck")

# Check if we should continue the sampling...
if (isfile(modelfile) || isfile(burninfile)) && isfile(iterfile)
    d = now() - Dates.unix2datetime(stat(iterfile).mtime)
    if round(d, Minute) > Minute(60)
        # remove lock file
        if isfile(lockfile)
            rm(lockfile)
        end
    end
end

# Check if we should continue the sampling...
if isfile(lockfile)
    @info "MCMC sampling is locked, exiting."
    exit(-2)
else

    t = rand(1:15) + rand(1:15)
    sleep(t)

    if isfile(lockfile)
        @info "MCMC sampling is locked, exiting."
        exit(-2)
    end

    touch(lockfile)
end

# Write out the configuration of the run.
save(joinpath(resultsdir, "configuration.jld2"), Dict("configuration" => args))

# Start script
open(joinpath(resultsdir, "log.txt"), "w+") do io
    logger = SimpleLogger(io)
    with_logger(logger) do
        if Threads.nthreads() == 1
            @info "Run export JULIA_NUM_THREADS=nproc before running julia to use multithreading."
        else
            @info "Using $(Threads.nthreads()) threads."

            # Test if all threads are available.
            thrsa = zeros(Threads.nthreads())
            Threads.@threads for t in 1:Threads.nthreads()
                thrsa[t] = Threads.threadid()
            end

            avThreads = length(unique(thrsa[thrsa .> 0]))
            @assert avThreads == Threads.nthreads() "Only $(avThreads) threads available."
        end

        @info "Loading dataset: $dataset stored under: $datafolder"
        x_train = vcat(
            readdlm(joinpath(datafolder, dataset*".ts.data"), ','),
            readdlm(joinpath(datafolder, dataset*".valid.data"), ',')
        )

        # Set missing values to NaN
        x_train[x_train .== "missing"] .= NaN
        x_train = Float64.(x_train)

        x_test = readdlm(joinpath(datafolder, dataset*".test.data"), ',', Bool)

        # Set the random seed.
        rnd_seed = args["seed"]+args["chain"]
        rng = Random.seed!(rnd_seed)

        # Get data dimensions.
        (N_train, D) = size(x_train)
        (N_test, _) = size(x_test)

        @info "Sampling dataset of (train) size: N=$N_train, D=$D and (test) size: N=$N_test"

        # Define the base distributions
        priors = map(d -> Beta(1/2, 1/2), 1:D)
        likelihoods = map(d -> Bernoulli(rand(priors[d])), 1:D)
        sstats = map(d -> BetaSufficientStats(), 1:D)

        # Construct a second network which we use for sampling.
        @info "Constructing SPN for sampling."
        spn = templateRegion(args["alpha"],
                             args["beta"],
                             args["gamma"],
                             priors,
                             likelihoods,
                             sstats,
                             N_train, D,
                             args["Ksum"], args["Kprod"],
                             args["J"], args["Kleaf"],
                             1, args["L"], root = true)

        # Set the number of iterations.
        iterations = args["iterations"]
        burnin = args["burnin"]

        # ######### #
        # MCMC code #
        # ######### #
        ks = args["numsamples"] > 0 ? args["numsamples"] : iterations
        ks = min(iterations, ks)
        thinning = Int(iterations / ks)

        @info "Starting MCMC for $iterations iterations."

        startIter = 1
        storePredictions = true
        storeTrace = false

        # Initialize SPN.
        if isfile(modelfile)
            # Resume.
            spn = open(modelfile, "r") do fin
                deserialize(fin)
            end
            r_ = read(iterfile, Int)
            startIter = 1 + r_
            @info "Continue MCMC."
        elseif isfile(burninfile)
            spn = open(burninfile, "r") do fin
                deserialize(fin)
            end
            @info "Continue burnin or start using previously converged chain, not storing trace."
            r_ = read(iterfile, Int)
            startIter = 1 + r_
        else
            @warn "Starting from random position, not storing predictions!"
            randomscopes!(spn, x_train)
        end

        nodes = topoligicalorder(spn)

        for iteration in startIter:(burnin+iterations)

            storeTrace = iteration <= burnin
            storePredictions = iteration > burnin

            dateprint = Dates.format(now(), "yyyy-mm-ddTHH:MM:SS")

            @info "[$(dateprint)] iteration = $(iteration) "

            ancestralsampling!(spn, x_train)
            gibbssamplescopes!(spn, x_train)

            # Store the scopes, weights and parameters on the disk.
            if ((iteration-burnin) % thinning) == 0

                # Serialize model.
                fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                open(fname, "w") do fout
                    serialize(fout, spn)
                end
                if iteration > burnin
                    mv(fname, modelfile, force=true)
                else
                    mv(fname, burninfile, force=true)
                end

                if storeTrace
                    # Store parameters of root.
                    Nparam = length(spn.logweights[:])
                    θ = isfile(tracefile) ? load(tracefile, "weights") : zeros(0,Nparam)
                    fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                    save(fname, "weights", vcat(θ, spn.logweights[:]'))
                    mv(fname, tracefile, force=true)
                end

                if storePredictions
                    # Evaluate network.
                    lp_test = logpdf(spn, x_test)
                    lp_train = logpdf(spn, x_train)

                    b_ = isfile(llhfile)
                    llh_test = b_ ? load(llhfile, "test") : zeros(0,N_test)
                    llh_train = b_ ? load(llhfile, "train") : zeros(0,N_train)
                    fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                    save(fname, "train", vcat(llh_train, lp_train'), "test", vcat(llh_test, lp_test'))
                    mv(fname, llhfile, force=true)

                    nweights = isfile(nparamfile) ? load(nparamfile, "weights") : Vector{Int}()
                    nleaves = isfile(nparamfile) ? load(nparamfile, "leaves") : Vector{Int}()
                    n_ = 0
                    for n in filter(n -> n isa RegionGraphNode, nodes)
                        n_ += sum(n.active)
                    end
                    push!(nweights, n_)

                    n_ = 0
                    for n in filter(n -> n isa FactorizedDistributionGraphNode, nodes)
                        n_ += sum(any(n.obsVecs, dims=1)) * length(scope(n))
                    end
                    push!(nleaves, n_)

                    fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                    save(fname, "weights", nweights, "leaves", nleaves)
                    mv(fname, nparamfile, force=true)
                end

                write(iterfile, iteration)
            end
        end

        @info "Finished MCMC"

        # Save file indicating we are done.
        touch(donefile)
        rm(lockfile)

        tmp_files = filter(f -> isfile(joinpath(resultsdir,f)) && startswith(f, "tmp_"), readdir(resultsdir))
        for f in tmp_files
            rm(joinpath(resultsdir,f))
        end
    end
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
module spectralanalysis

export eigenstates, eigenenergies, simdiag

using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse

const nonhermitian_warning = "The given operator is not hermitian. If this is due to a numerical error make the operator hermitian first by calculating (x+dagger(x))/2 first."

"""
    eigenstates(op::Operator[, n::Int; warning=true])

Calculate the lowest n eigenvalues and their corresponding eigenstates.

This is just a thin wrapper around julia's `eig` and `eigs` functions. Which
of them is used depends on the type of the given operator. If more control
about the way the calculation is done is needed, use the functions directly.
More details can be found at
[http://docs.julialang.org/en/stable/stdlib/linalg/].

If the given operator is non-hermitian a warning is given. This behavior
can be turned off using the keyword `warning=false`.
"""
function eigenstates(op::DenseOperator, n::Int=length(basis(op)); warning=true)
    b = basis(op)
    if ishermitian(op)
        D, V = eig(Hermitian(op.data), 1:n)
        states = [Ket(b, V[:, k]) for k=1:length(D)]
        return D, states
    else
        warning && warn(nonhermitian_warning)
        D, V = eig(op.data)
        states = [Ket(b, V[:, k]) for k=1:length(D)]
        perm = sortperm(D, by=real)
        permute!(D, perm)
        permute!(states, perm)
        return D[1:n], states[1:n]
    end
end

"""
For sparse operators by default it only returns the 6 lowest eigenvalues.
"""
function eigenstates(op::SparseOperator, n::Int=length(basis(op)); warning=true)
    b = basis(op)
    if ishermitian(op)
        data = Hermitian(op.data)
    else
        warning && warn(nonhermitian_warning)
        data = op.data
    end
    D, V = eigs(data; nev=n, which=:SR)
    states = [Ket(b, V[:, k]) for k=1:length(D)]
    D, states
end

"""
    eigenenergies(op::Operator[, n::Int; warning=true])

Calculate the lowest n eigenvalues.

This is just a thin wrapper around julia's `eigvals`. If more control
about the way the calculation is done is needed, use the function directly.
More details can be found at
[http://docs.julialang.org/en/stable/stdlib/linalg/].

If the given operator is non-hermitian a warning is given. This behavior
can be turned off using the keyword `warning=false`.
"""
function eigenenergies(op::DenseOperator, n::Int=length(basis(op)); warning=true)
    b = basis(op)
    if ishermitian(op)
        D = eigvals(Hermitian(op.data), 1:n)
        return D
    else
        warning && warn(nonhermitian_warning)
        D = eigvals(op.data)
        sort!(D, by=real)
        return D[1:n]
    end
end

"""
For sparse operators by default it only returns the 6 lowest eigenvalues.
"""
eigenenergies(op::SparseOperator, n::Int=6; warning=true) = eigenstates(op, n; warning=warning)[1]

arithmetic_unary_error = operators.arithmetic_unary_error
eigenstates(op::Operator, n::Int=0) = arithmetic_unary_error("eigenstates", op)
eigenenergies(op::Operator, n::Int=0) = arithmetic_unary_error("eigenenergies", op)

"""
    simdiag(ops; atol, rtol)

Simultaneously diagonalize commuting Hermitian operators specified in `ops`.

This is done by diagonalizing the sum of the operators. The eigenvalues are
computed by ``a = ⟨ψ|A|ψ⟩`` and it is checked whether the eigenvectors fulfill
the equation ``A|ψ⟩ = a|ψ⟩``.

# Arguments
* `ops`: Vector of sparse or dense operators.
* `atol=1e-14`: kwarg of Base.isapprox specifying the tolerance of the
        approximate check
* `rtol=1e-14`: kwarg of Base.isapprox specifying the tolerance of the
        approximate check

# Returns
* `evals_sorted`: Vector containing all vectors of the eigenvalues sorted
        by the eigenvalues of the first operator.
* `v`: Common eigenvectors.
"""
function simdiag(ops::Vector{T}; atol::Real=1e-14, rtol::Real=1e-14) where T<:DenseOperator
    # Check input
    for A=ops
        if !ishermitian(A)
            error("Non-hermitian operator given!")
        end
    end

    d, v = eig(sum(ops).data)

    evals = [Vector{Complex128}(length(d)) for i=1:length(ops)]
    for i=1:length(ops), j=1:length(d)
        vec = ops[i].data*v[:, j]
        evals[i][j] = (v[:, j]'*vec)[1]
        if !isapprox(vec, evals[i][j]*v[:, j]; atol=atol, rtol=rtol)
            error("Simultaneous diagonalization failed!")
        end
    end

    index = sortperm(real(evals[1][:]))
    evals_sorted = [real(evals[i][index]) for i=1:length(ops)]
    return evals_sorted, v[:, index]
end

end # module

ext: jl


text:
using HDF5

include("../utils.jl")
include("../initializations.jl")
include("../codebook_update.jl")
include("../encodings/encode_chain.jl")

# Train a chain quantizer with viterbi decoding
function train_chainq{T <: AbstractFloat}(
  X::Matrix{T},             # d-by-n matrix of data points to train on.
  m::Integer,               # number of codebooks
  h::Integer,               # number of entries per codebook
  R::Matrix{T},             # Init rotation matrix
  B::Matrix{Int16},         # Init codes
  C::Vector{Matrix{T}},     # Init codebooks
  niter::Integer,           # number of optimization iterations
  V::Bool=false)            # whether to print progress

  d, n = size( X )
  obj  = zeros(Float32, niter+1)

  CB = zeros( Float32, size(X) )
  RX = R' * X

  # Initialize C
  C = update_codebooks_chain( RX, B, h, V )
  @printf("%3d %e\n", -2, qerror( RX, B, C ))

  # Initialize B
  B   = encoding_viterbi( RX, C, V )
  @printf("%3d %e\n", -1, qerror( RX, B, C ))

  for iter = 0:niter
    obj[iter+1] = qerror( RX, B, C  )
    @printf("%3d %e\n", iter, obj[iter+1])

    # update CB
    CB[:] = 0;
    for i = 1:m; CB += C[i][:, vec(B[i,:]) ]; end

    # update R
    U, S, VV = svd(X * CB', thin=true)
    R = U * VV';

    # update R*X
    RX = R' * X

    # Update the codebooks #
    C = update_codebooks_chain( RX, B, h, V )

    # Update the codes with lattice search
    B = encoding_viterbi( RX, C, V );

  end

  return C, B, R, obj
end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
module TestAddToField

using Test, AddToField

import Base: setproperty!

function addone(x)
	x[] += 1
	return x[]
end

struct PropertyDict
	D::Dict{Symbol, Int}
end

setproperty!(PD::PropertyDict, name::Symbol, value) = setindex!(PD.D, value, name)

@testset "nt" begin
	s = "e1"

	r = Ref(1)

	res = @addnt begin
		a = 1
		@add a
		@add b = 2
		@add :c1 c = 3
		@add "d1" d = 4
		@add "$s" e = 5
		f, g, h = 6, 7, 8

		@add f, g, h

		@add i = addone(r)
		@add "j" addone(r)
	end

	correct= (
		a = 1,
		b = 2,
		c1 = 3,
		d1 = 4,
		e1 = 5,
		f = 6,
		g = 7,
		h = 8,
		i = 2,
		j = 3
	)

	@test res == correct
end

@testset "dict" begin
	s = "e1"

	r = Ref(1)

	res = @adddict begin
		a = 1
		@add a
		@add b = 2
		@add :c1 c = 3
		@add "d1" d = 4
		@add "$s" e = 5
		f, g, h = 6, 7, 8

		@add f, g, h

		@add i = addone(r)
		@add "j" addone(r)
	end

	correct= Dict(
		:a => 1,
		:b => 2,
		:c1 => 3,
		:d1 => 4,
		:e1 => 5,
		:f => 6,
		:g => 7,
		:h => 8,
		:i => 2,
		:j => 3
	)

	@test res == correct
end

@testset "two args" begin
	D = Dict{Symbol, Int}()

	s = "e1"

	r = Ref(1)

	res = @addto! D begin
		a = 1
		@add a
		@add b = 2
		@add :c1 c = 3
		@add "d1" d = 4
		@add "$s" e = 5
		f, g, h = 6, 7, 8

		@add f, g, h

		@add i = addone(r)
		@add "j" addone(r)
	end

	correct = Dict(
		:a => 1,
		:b => 2,
		:c1 => 3,
		:d1 => 4,
		:e1 => 5,
		:f => 6,
		:g => 7,
		:h => 8,
		:i => 2,
		:j => 3
	)

	@test res == correct

	PD = PropertyDict(Dict{Symbol, Int}())

	r = Ref(1)

	res = @addto! PD begin
		a = 1
		@add a
		@add b = 2
		@add :c1 c = 3
		@add "d1" d = 4
		@add "$s" e = 5
		f, g, h = 6, 7, 8

		@add f, g, h

		@add i = addone(r)
		@add "j" addone(r)
	end

	@test res.D == correct
end

@testset "let" begin
	x = 1

	@test x == 1
	res = @addnt let
		local x = 2
		@add y = 5
	end

	@test res == (; y = 5)
	@test x == 1

	x = 1

	@test x == 1
	res = @adddict let
		local x = 2
		@add y = 5
	end

	@test res == Dict(:y => 5)
	@test x == 1

	res = @addnt let y = 1
		@add x = 1
	end;
	@test res == (;x = 1)

	res = @addnt let y = 1, z = 2
		@add x = 1
	end;
	@test res == (;x = 1)

	res = @adddict let y = 1
		@add x = 1
	end;
	@test res == Dict(:x => 1)

	res = @adddict let y = 1, z = 2
		@add x = 1
	end;
	@test res == Dict(:x => 1)

	res = Dict{Symbol, Int}()

	@addto! res let
		local x = 2
		@add y = 5
	end

	@test res == Dict([:y => 5])
	@test x == 1

	@addto! res let z = 1
		local x = 2
		@add y = 5
	end

	@test res == Dict([:y => 5])
	@test x == 1

	@addto! res let z = 1, z1 = 2
		local x = 2
		@add y = 5
	end

	@test res == Dict([:y => 5])
	@test x == 1
end

@testset "new dict" begin
	d = @addto! Dict() begin
		x = 1
		@add y = 2
	end

	@test d == Dict{Any, Any}(:y => 2)

	d = @addto! Dict() let
		x = 10
		@add y = 20
	end

	@test d == Dict{Any, Any}(:y => 20)

	d = @addto! Dict() let z = 50
		x = 100
		@add y = 200
	end

	@test d == Dict{Any, Any}(:y => 200)

	d = @addto! Dict() let z = 50, zz = 51
		x = 100
		@add y = 200
	end

	@test d == Dict{Any, Any}(:y => 200)
end

end # module


text:
function perform_thresholding(f, M, t)
    """
        Only 3 types of thresholding currently implemented
    """
    if t == "largest"
        a = sort(abs(f)[:])[end:-1:1] #sort a 1D copy of F in descending order
        T = a[M]
        y = f.*(abs(f) .> T)
    elseif t == "soft"
        s = abs(f) - M
        s = (s + abs(s))./2
        y = (f./abs(f)).*s
    elseif t == "hard"
        y = f.*(abs(f) .> M)
    end
    return y
end


text:
# Simplifications go here

# Basic maps

composedmap1(m1::IdentityMap, m2) = m2
composedmap1(m1::IdentityMap{T}, m2::Map{T}) where {T} = m2
composedmap1(m1::IdentityMap{T}, m2::Map{S}) where {S,T} = convert(Map{T}, m2)

composedmap2(m1, m2::IdentityMap) = m1
composedmap2(m1::Map{T}, m2::IdentityMap{T}) where {T} = m1
composedmap2(m1::Map{S}, m2::IdentityMap{T}) where {S,T} = convert(Map{T}, m1)

composedmap2(m1, m2::ConstantMap) = m2
composedmap2(m1::Map{T}, m2::ConstantMap) where {T} = ConstantMap{T}(constant(m2))
composedmap1(m1::ConstantMap{T}, m2) where {T} = ConstantMap{T}(m2(constant(m1)))

composedmap2(m1, m2::ZeroMap) = m2
composedmap2(m1::Map{T}, m2::ZeroMap{S,U}) where {S,T,U} = ZeroMap{T,U}()
composedmap1(m1::ConstantMap{T}, m2::ZeroMap{S,U}) where {S,T,U} = ZeroMap{T,U}()

multiply_map1(m1::ZeroMap, m2) = m1
multiply_map2(m1, m2::ZeroMap) = m2
multiply_map2(m1::ConstantMap{T}, m2::ConstantMap{S}) where {S,T} =
    ConstantMap{promote_type(S,T)}(constant(m1)*constant(m2))

sum_map1(m1::ZeroMap, m2) = m2
sum_map2(m1, m2::ZeroMap) = m1
sum_map2(m1::ConstantMap{T}, m2::ConstantMap{S}) where {S,T} =
    ConstantMap{promote_type(S,T)}(constant(m1)+constant(m2))

## Affine maps

composedmap(m1::AbstractAffineMap, m2::AbstractAffineMap) = affine_composition(m1, m2)

"""
Compute the affine map that represents map2 after map1, that is:
`y = a2*(a1*x+b1)+b2 = a2*a1*x + a2*b1 + b2`.
"""
affine_composition(map1::AbstractAffineMap, map2::AbstractAffineMap) =
    AffineMap(matrix(map2) * matrix(map1), matrix(map2)*vector(map1) + vector(map2))

affine_composition(map1::AffineMap, map2::AffineMap) =
    AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), unsafe_matrix(map2)*unsafe_vector(map1) + unsafe_vector(map2))

affine_composition(map1::LinearMap, map2::LinearMap) =
    LinearMap(unsafe_matrix(map2) * unsafe_matrix(map1))

affine_composition(map1::LinearMap, map2::AffineMap) =
    AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), unsafe_vector(map2))

affine_composition(map1::AffineMap, map2::LinearMap) =
    AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), unsafe_matrix(map2)*unsafe_vector(map1))

affine_composition(map1::Translation, map2::Translation) =
    Translation(unsafe_vector(map2) + unsafe_vector(map1))

affine_composition(map1::Translation, map2::LinearMap) =
    AffineMap(unsafe_matrix(map2), unsafe_matrix(map2)*unsafe_vector(map1))

affine_composition(map1::LinearMap, map2::Translation) =
    AffineMap(unsafe_matrix(map1), unsafe_vector(map2))

# The sum of two affine maps is again an affine map
mapsum(map1::AbstractAffineMap, map2::AbstractAffineMap) =
    AffineMap(matrix(map1)+matrix(map2), vector(map1)+vector(map2))

==(m1::ProductMap, m2::IdentityMap) = all(map(isidentity, components(m1)))
==(m1::IdentityMap, m2::ProductMap) = m2 == m1


text:
while true
	if read(ExplorerHat.input[1]) == 1
		on(ExplorerHat.motor[1], 50)
		sleep(0.5)
		off(ExplorerHat.motor[1])
	end
        sleep(0.2)
end


text:
### GENERALIZED COLUMN OPS

# We want to support arbitrary dimensional arrays.  When data comes in
# N dimensions, we assume it is an array of N-1 dimensional instances
# and the last dimension gives us the instance count.  We will refer
# to the first N-1 dimensions as generalized "columns" of the
# data. These columns are indexed by the last index of an array,
# i.e. column i corresponds to b[:,:,...,i].

# Here are some convenience functions for generalized columns:
# We consider a 1-D array a single column:

csize(a)=(ndims(a)==1 ? size(a) : size(a)[1:end-1])
csize(a,n)=tuple(csize(a)..., n) # size if you had n columns
clength(a)=(ndims(a)==1 ? length(a) : stride(a,ndims(a)))
ccount(a)=(ndims(a)==1 ? 1 : size(a,ndims(a)))
csub(a,i)=(ndims(a)==1 ? error() : sub(a, ntuple(i->(:), ndims(a)-1)..., i))
cget(a,i)=(ndims(a)==1 ? error() : getindex(a, ntuple(i->(:), ndims(a)-1)..., i))
size2(y)=(nd=ndims(y); (nd==1 ? (length(y),1) : (stride(y, nd), size(y, nd)))) # size as a matrix
size2(y,i)=size2(y)[i]

function minibatch(x, y, batchsize)
    data = Any[]
    for i=1:batchsize:ccount(x)
        j=min(i+batchsize-1,ccount(x))
        push!(data, (cget(x,i:j), cget(y,i:j)))
    end
    return data
end

# CSLICE!  Returns a slice of array b, with columns specified in range
# r, using the storage in KUarray a.  The element types need to match,
# but the size of a does not need to match, it is adjusted as
# necessary.  This is used in train and predict to get data from a raw
# array into a KUarray for minibatching.

# cslice!{A,B,T}(a::KUdense{A,T}, b::KUdense{B,T}, r::UnitRange)=cslice!(a,b.arr,r)

# function cslice!{A,T}(a::KUdense{A,T}, b::BaseArray{T}, r::UnitRange)
#     n  = clength(b) * length(r)
#     length(a.ptr) >= n || resize!(a.ptr, int(resizefactor(KUdense)*n+1))
#     b1 = 1 + clength(b) * (first(r) - 1)
#     copysync!(a.ptr, 1, b, b1, n)
#     a.arr = arr(a.ptr, csize(b, length(r)))
#     return a
# end

# # For non-contiguous columns:
# function cslice!{A,T}(a::KUdense{A,T}, b::BaseArray{T}, cols)
#     ncols = length(cols)
#     clen = clength(b)
#     n = clen * ncols
#     length(a.ptr) >= n || resize!(a.ptr, int(resizefactor(KUdense)*n+1))
#     alen = 0
#     for i=1:ncols
#         bidx = (cols[i]-1)*clen + 1
#         copysync!(a.ptr, alen+1, b, bidx, clen)
#         alen += clen
#     end
#     a.arr = arr(a.ptr, csize(b, ncols))
#     return a
# end

function cslice!{T}(a::BaseArray{T}, b::BaseArray{T}, cols)
    ncols = length(cols)
    clen = clength(b)
    n = clen * ncols
    size(a) == csize(b, ncols) || error("Size mismatch")
    alen = 0
    for i=1:ncols
        bidx = (cols[i]-1)*clen + 1
        copysync!(a, alen+1, b, bidx, clen) # t:134
        alen += clen
    end
    return a
end

function cslice!{T}(a::SparseMatrixCSC{T}, b::SparseMatrixCSC{T}, cols)
    bptr = b.colptr
    nz = 0; for i in cols; nz += bptr[i+1]-bptr[i]; end
    a.m = b.m
    a.n = length(cols)
    resize!(a.nzval, nz)
    resize!(a.rowval, nz)
    resize!(a.colptr, a.n+1)
    a.colptr[1] = a1 = aj = 1
    for bj in cols                 # copy column b[:,bj] to a[:,aj]
        b1 = bptr[bj]
        nz = bptr[bj+1]-b1
        copysync!(a.nzval, a1, b.nzval, b1, nz) # t:217
        copysync!(a.rowval, a1, b.rowval, b1, nz) # t:191
        a1 += nz
        a.colptr[aj+=1] = a1
    end
    @assert aj == a.n+1
    return a
end

# TODO: write the non-contiguous sparse version
# function cslice!{A,B,T}(a::KUsparse{A,T}, b::KUsparse{B,T}, r::UnitRange)
#     bptr = cpucopy(b.colptr)
#     nz = 0; for i in r; nz += bptr[i+1]-bptr[i]; end
#     a.m = b.m
#     a.n = length(r)
#     resize!(a.nzval, nz)
#     resize!(a.rowval, nz)
#     aptr = Array(Int32, a.n+1)
#     aptr[1] = a1 = aj = 1
#     for bj in r                 # copy column b[:,bj] to a[:,aj]
#         b1 = bptr[bj]
#         nz = bptr[bj+1]-b1
#         copysync!(a.nzval, a1, b.nzval, b1, nz)
#         copysync!(a.rowval, a1, b.rowval, b1, nz)
#         a1 += nz
#         aptr[aj+=1] = a1
#     end
#     @assert aj == a.n+1
#     a.colptr = convert(A, aptr)
#     return a
# end

# cslice!{A,T}(a::KUsparse{A,T}, b::SparseMatrixCSC{T}, r::UnitRange)=cslice!(a, convert(KUsparse, b), r)

# CCOPY! Copy n columns from src starting at column si, into dst
# starting at column di.  Used by predict to construct output.
# Don't need the sparse version, output always dense.

# ccopy!{A,T,N}(dst::BaseArray{T,N}, di, src::KUdense{A,T,N}, si=1, n=ccount(src)-si+1)=(ccopy!(dst,di,src.arr,si,n); dst)
# ccopy!{A,B,T,N}(dst::KUdense{A,T,N}, di, src::KUdense{B,T,N}, si=1, n=ccount(src)-si+1)=(ccopy!(dst.arr,di,src.arr,si,n); dst)

function ccopy!{T,N}(dst::BaseArray{T,N}, di, src::BaseArray{T,N}, si=1, n=ccount(src)-si+1)
    @assert csize(dst)==csize(src)
    clen = clength(src)
    d1 = 1 + clen * (di - 1)
    s1 = 1 + clen * (si - 1)
    copysync!(dst, d1, src, s1, clen * n)
    return dst
end

# CADD! Add n columns from src starting at column si, into dst
# starting at column di.  Used by uniq!  Don't need sparse version,
# weights always dense.

using Base.LinAlg: axpy!

# cadd!{A,T,N}(dst::BaseArray{T,N}, di, src::KUdense{A,T,N}, si=1, n=ccount(src)-si+1)=(cadd!(dst,di,src.arr,si,n); dst)
# cadd!{A,B,T,N}(dst::KUdense{A,T,N}, di, src::KUdense{B,T,N}, si=1, n=ccount(src)-si+1)=(cadd!(dst.arr,di,src.arr,si,n); dst)

function cadd!{T,N}(dst::BaseArray{T,N}, di, src::BaseArray{T,N}, si=1, n=ccount(src)-si+1)
    @assert csize(dst)==csize(src)
    @assert ccount(dst) >= di+n-1
    @assert ccount(src) >= si+n-1
    clen = clength(src)
    d1 = 1 + clen * (di - 1)
    s1 = 1 + clen * (si - 1)
    n1 = clen * n
    axpy!(n1, one(T), pointer(src, s1), 1, pointer(dst, d1), 1)
    return dst
end

# CCAT! generalizes append! to multi-dimensional arrays.  Adds the
# ability to specify particular columns to append.  Used in
# kperceptron to add support vectors.

# ccat!{A,B,T,N}(a::KUdense{A,T,N}, b::KUdense{B,T,N}, cols=(1:ccount(b)))=ccat!(a,b.arr,cols)

# function ccat!{A,T,N}(a::KUdense{A,T,N}, b::BaseArray{T,N}, cols=(1:ccount(b)))
#     @assert csize(a)==csize(b)
#     alen = length(a)
#     clen = clength(a)
#     ncols = length(cols)
#     n = alen + ncols * clen
#     length(a.ptr) >= n || resize!(a.ptr, round(Int,resizefactor(KUdense)*n+1))
#     for i=1:ncols
#         bidx = (cols[i]-1)*clen + 1
#         copysync!(a.ptr, alen+1, b, bidx, clen)
#         alen += clen
#     end
#     a.arr = arr(a.ptr, csize(a, ccount(a) + ncols))
#     return a
# end

# ccat!{A,T}(a::KUsparse{A,T}, b::SparseMatrixCSC{T}, cols=(1:ccount(b)))=ccat!(a,convert(KUsparse,b),cols)

# function ccat!{A,B,T}(a::KUsparse{A,T}, b::KUsparse{B,T}, cols=(1:ccount(b)))
#     # a: m, n, colptr, rowval, nzval
#     # colptr[i]: starting index (in rowval,nzval) of column i
#     # colptr[n+1]: nz+1
#     @assert size(a,1) == size(b,1)
#     aptr = cpucopy(a.colptr)
#     bptr = cpucopy(b.colptr)
#     na = aptr[a.n+1]-1          # count new nonzero entries in a
#     ncols = length(cols)
#     for i in cols; na += bptr[i+1]-bptr[i]; end
#     resize!(a.nzval, na)
#     resize!(a.rowval, na)
#     na = aptr[a.n+1]-1          # restart the count
#     for i=1:ncols
#         bj=cols[i]              # bj'th column of b
#         aj=a.n+i                # will become aj'th column of a
#         nz=bptr[bj+1]-bptr[bj]  # with nz nonzero values
#         @assert length(aptr) == aj
#         push!(aptr, aptr[aj]+nz) # aptr[aj+1] = aptr[aj]+nz
#         copysync!(a.nzval,na+1,b.nzval,bptr[bj],nz)
#         copysync!(a.rowval,na+1,b.rowval,bptr[bj],nz)
#         na = na+nz
#     end
#     @assert length(aptr) == a.n + ncols + 1
#     resize!(a.colptr, a.n + ncols + 1)
#     copysync!(a.colptr, a.n+2, aptr, a.n+2, ncols)
#     a.n += ncols
#     return a
# end

### UNIQ! leaves unique columns in its first argument and sums to
### corresponding columns in the remaining arguments.  Used by
### kperceptron in merging identical support vectors.

# function uniq!{A<:Array}(s::KUdense{A}, ww::KUdense...)
#     oldn = ccount(s)                                            # number of original support vectors
#     for w in ww; @assert ccount(w) == oldn; end
#     ds = Dict{Any,Int}()                                        # support vector => new index
#     newn = 0                                                    # number of new support vectors
#     for oldj=1:oldn
#         newj = get!(ds, _colkey(s,oldj), newn+1)
#         if newj <= newn                                         # s[:,oldj] already in s[:,newj]
#             @assert newj <= newn == length(ds) < oldj
#             for w in ww; cadd!(w,newj,w,oldj,1); end
#         else                                                    # s[:,oldj] to be copied to s[:,newj]
#             @assert newj == newn+1 == length(ds) <= oldj
#             newn += 1
#             if newj != oldj
#                 ccopysync!(s,newj,s,oldj,1)
#                 for w in ww; ccopy!(w,newj,w,oldj,1); end
#             end
#         end
#     end
#     @assert newn == length(ds)
#     resize!(s, csize(s, newn))
#     for w in ww; resize!(w, csize(w, newn)); end
#     return tuple(s, ww...)
# end

# function uniq!{A<:Array}(s::KUsparse{A}, ww::KUdense...)
#     oldn = ccount(s)                                            # number of original support vectors
#     for w in ww; @assert ccount(w) == oldn; end
#     ds = Dict{Any,Int}()                                        # support vector => new index
#     @assert s.colptr[1]==1
#     ncol = 0
#     nnz = 0
#     for oldj=1:oldn
#         newj = get!(ds, _colkey(s,oldj), ncol+1)
#         if newj <= ncol                                          # s[:,oldj] already in s[:,newj]
#             @assert newj <= length(ds) == ncol < oldj
#             for w in ww; cadd!(w,newj,w,oldj,1); end
#         else                                                    # s[:,oldj] to be copied to s[:,newj]
#             @assert newj == ncol+1 == length(ds) <= oldj
#             from = s.colptr[oldj]
#             nval = s.colptr[oldj+1] - from
#             to = nnz+1
#             ncol += 1
#             nnz += nval
#             if newj != oldj
#                 copysync!(s.rowval, to, s.rowval, from, nval)
#                 copysync!(s.nzval, to, s.nzval, from, nval)
#                 s.colptr[ncol+1] = nnz+1
#                 for w in ww; ccopy!(w,newj,w,oldj,1); end
#             else
#                 @assert to == from
#                 @assert s.colptr[ncol+1] == nnz+1
#             end
#         end
#     end
#     @assert length(ds) == ncol
#     s.n = ncol
#     resize!(s.colptr, ncol+1)
#     resize!(s.rowval, nnz)
#     resize!(s.nzval,  nnz)
#     for w in ww; resize!(w, csize(w, s.n)); end
#     return tuple(s, ww...)
# end

# _colkey{A<:Array}(s::KUdense{A},j)=sub(s.arr, ntuple(i->(i==ndims(s) ? (j:j) : Colon()), ndims(s))...)

# function _colkey{A<:Array}(s::KUsparse{A},j)
#     a=s.colptr[j]
#     b=s.colptr[j+1]-1
#     r=sub(s.rowval, a:b)
#     v=sub(s.nzval, a:b)
#     (r,v)
# end

# Getting columns one at a time is expensive, just copy the whole array
# CudaArray does not support sub, even if it did we would not be able to hash it
# getcol{T}(s::KUdense{CudaArray,T}, j)=(n=clength(s);copysync!(Array(T,csize(s,1)), 1, s.arr, (j-1)*n+1, n))

# we need to look at the columns, might as well copy

# function uniq!{A<:CudaArray}(s::KUdense{A}, ww::KUdense...)
#     ss = cpucopy(s)
#     uniq!(ss, ww...)
#     cslice!(s, ss, 1:ccount(ss))
#     return tuple(s, ww...)
# end

# function uniq!{A<:CudaArray}(s::KUsparse{A}, ww::KUdense...)
#     ss = cpucopy(s)
#     uniq!(ss, ww...)
#     cslice!(s, ss, 1:ccount(ss))
#     return tuple(s, ww...)
# end

# TODO: fix array types

    # ds = Dict{Any,Int}()        # dictionary of support vectors
    # ns = 0                      # number of support vectors
    # s0 = spzeros(eltype(s), Int32, size(s,1), ns) # new sv matrix
    # for j=1:size(s,2)
    #     jj = get!(ds, s[:,j], ns+1)
    #     if jj <= ns             # s[:,j] already in s0[:,jj]
    #         @assert ns == length(ds) < j
    #         u[:,jj] += u[:,j]
    #         v[:,jj] += v[:,j]
    #     else                    # s[:,j] to be added to s0
    #         @assert jj == ns+1 == length(ds) <= j
    #         ns = ns+1
    #         hcat!(s0, s, [j], 1)
    #         if jj != j
    #             u[:,jj] = u[:,j]
    #             v[:,jj] = v[:,j]
    #         end
    #     end
    # end
    # @assert ns == length(ds) == size(s0,2)
    # u = size!(u, (size(u,1),ns))
    # v = size!(v, (size(v,1),ns))
    # for f in names(s); s.(f) = s0.(f); end
    # return (s,u,v)

# function uniq!(s::SparseMatrixCSC, u::AbstractArray, v::AbstractArray)
#     ds = Dict{Any,Int}()        # dictionary of support vectors
#     ns = 0                      # number of support vectors
#     s0 = spzeros(eltype(s), Int32, size(s,1), ns) # new sv matrix
#     for j=1:size(s,2)
#         jj = get!(ds, s[:,j], ns+1)
#         if jj <= ns             # s[:,j] already in s0[:,jj]
#             @assert ns == length(ds) < j
#             u[:,jj] += u[:,j]
#             v[:,jj] += v[:,j]
#         else                    # s[:,j] to be added to s0
#             @assert jj == ns+1 == length(ds) <= j
#             ns = ns+1
#             hcat!(s0, s, [j], 1)
#             if jj != j
#                 u[:,jj] = u[:,j]
#                 v[:,jj] = v[:,j]
#             end
#         end
#     end
#     @assert ns == length(ds) == size(s0,2)
#     u = size!(u, (size(u,1),ns))
#     v = size!(v, (size(v,1),ns))
#     for f in names(s); s.(f) = s0.(f); end
#     return (s,u,v)
# end

# function uniq!(ss::CudaSparseMatrixCSC, uu::AbstractCudaArray, vv::AbstractCudaArray)
#     (s,u,v)=map(cpucopy,(ss,uu,vv))
#     (s,u,v)=uniq!(s,u,v)
#     n = size(s,2)
#     uu = size!(uu, (size(u,1),n))
#     vv = size!(vv, (size(v,1),n))
#     copysync!(uu, 1, u, 1, size(u,1)*n)
#     copysync!(vv, 1, v, 1, size(v,1)*n)
#     (ss.m, ss.n, ss.colptr, ss.rowval, ss.nzval) = (s.m, s.n, gpucopy(s.colptr), gpucopy(s.rowval), gpucopy(s.nzval))
#     return (ss,uu,vv)
# end

# function cslice!{A,T,I}(a::KUsparse{A,T,I}, b::SparseMatrixCSC{T,I}, r::UnitRange)
#     nz = 0; for i in r; nz += b.colptr[i+1]-b.colptr[i]; end
#     a.m = b.m
#     a.n = length(r)
#     resize!(a.nzval, nz)
#     resize!(a.rowval, nz)
#     aptr = Array(I, a.n+1)
#     aptr[1] = a1 = aj = 1
#     for bj in r                 # copy column b[:,bj] to a[:,aj]
#         b1 = b.colptr[bj]
#         nz = b.colptr[bj+1]-b1
#         copysync!(a.nzval.arr, a1, b.nzval, b1, nz)
#         copysync!(a.rowval.arr, a1, b.rowval, b1, nz)
#         a1 += nz
#         aptr[aj+=1] = a1
#     end
#     @assert aj == a.n+1
#     copysync!(a.colptr, aptr)
#     return a
# end


text:
using Compat

# Cairo backend for compose
import Cairo

using Cairo: CairoContext, CairoSurface, CairoARGBSurface,
             CairoEPSSurface, CairoPDFSurface, CairoSVGSurface,
             CairoImageSurface

@compat abstract type ImageBackend end
@compat abstract type PNGBackend <: ImageBackend end

@compat abstract type VectorImageBackend <: ImageBackend end
@compat abstract type SVGBackend <: VectorImageBackend end
@compat abstract type PDFBackend <: VectorImageBackend end
@compat abstract type PSBackend  <: VectorImageBackend end
@compat abstract type CairoBackend <: VectorImageBackend end

type ImagePropertyState
    stroke::RGBA{Float64}
    fill::RGBA{Float64}
    stroke_dash::Array{Float64,1}
    stroke_linecap::LineCap
    stroke_linejoin::LineJoin
    visible::Bool
    linewidth::AbsoluteLength
    fontsize::AbsoluteLength
    font::AbstractString
    clip::Nullable{ClipPrimitive}
end

type ImagePropertyFrame
    # Vector properties in this frame.
    vector_properties::Dict{Type, Property}

    # True if this property frame has scalar properties. Scalar properties are
    # emitted as a group (<g> tag) that must be closed when the frame is popped.
    has_scalar_properties::Bool

    function ImagePropertyFrame()
        return new(Dict{Type, Property}(), false)
    end
end

type Image{B <: ImageBackend} <: Backend
    out::IO
    surface::CairoSurface
    ctx::CairoContext

    width::Float64
    height::Float64

    # Current state
    stroke::RGBA{Float64}
    fill::RGBA{Float64}
    stroke_dash::Array{Float64,1}
    stroke_linecap::LineCap
    stroke_linejoin::LineJoin
    visible::Bool
    linewidth::AbsoluteLength
    fontsize::AbsoluteLength
    font::AbstractString
    clip::Nullable{ClipPrimitive}

    # Keep track of property
    state_stack::Vector{ImagePropertyState}
    property_stack::Vector{ImagePropertyFrame}
    vector_properties::Dict{Type, Nullable{Property}}

    # Close the surface when finished
    owns_surface::Bool

    # Backend is responsible for opening/closing the file
    ownedfile::Bool

    # Filename when ownedfile is true
    filename::@compat(Union{AbstractString, (@compat Void)})

    # True when finish has been called and no more drawing should occur
    finished::Bool

    # Emit on finish
    emit_on_finish::Bool

    # Points (or pixels for PNG) per mm
    ppmm::Float64

    # For use with the t/T and s/S commands in SVG-style paths
    last_ctrl1_point::Nullable{AbsoluteVec2}
    last_ctrl2_point::Nullable{AbsoluteVec2}

    @compat function (::Type{Image{B}}){B}(surface::CairoSurface, ctx::CairoContext, out::IO)
        img = new{B}()
        img.out = out
        img.width = 0
        img.height = 0
        img.surface = surface
        img.ctx = ctx

        img.stroke = default_stroke_color == nothing ?
                        RGBA{Float64}(0, 0, 0, 0) : convert(RGBA{Float64}, default_stroke_color)
        img.fill   = default_fill_color == nothing ?
                        RGBA{Float64}(0, 0, 0, 0) : convert(RGBA{Float64}, default_fill_color)
        img.stroke_dash = []
        img.stroke_linecap = LineCapButt()
        img.stroke_linejoin = LineJoinMiter()
        img.visible = true
        img.linewidth = default_line_width
        img.fontsize = default_font_size
        img.font = default_font_family
        img.clip = Nullable{ClipPrimitive}()

        img.state_stack = Array{ImagePropertyState}(0)
        img.property_stack = Array{ImagePropertyFrame}(0)
        img.vector_properties = Dict{Type, Nullable{Property}}()
        img.owns_surface = false
        img.ownedfile = false
        img.filename = nothing
        img.finished = false
        img.emit_on_finish = false
        img.ppmm = 72 / 25.4
        img.last_ctrl1_point = Nullable{AbsoluteVec2}()
        img.last_ctrl2_point = Nullable{AbsoluteVec2}()
        img
    end

    @compat function (::Type{Image{B}}){B}(surface::CairoSurface, ctx::CairoContext)
        Image{B}(surface, ctx, IOBuffer())
    end

    @compat (::Type{Image{B}}){B}(surface::CairoSurface) = Image{B}(surface, CairoContext(surface))

    @compat function (::Type{Image{B}}){B}(out::IO,
                              width::MeasureOrNumber,
                              height::MeasureOrNumber,
                              emit_on_finish::Bool=true;
                              dpi = (B == PNGBackend ? 96 : 72))

        width = size_measure(width)
        height = size_measure(height)

        if !isa(width, AbsoluteLength) || !isa(height, AbsoluteLength)
            error("Image size must be specificed in absolute units.")
        end

        ppmm = dpi / 25.4
        width = width.value * ppmm
        height = height.value * ppmm

        surface = newsurface(B, out, width, height)
        img = Image{B}(surface, CairoContext(surface), out)
        img.width = width
        img.height = height
        img.owns_surface = true
        img.emit_on_finish = emit_on_finish
        img.ppmm = ppmm
        img
    end

    @compat function (::Type{Image{B}}){B}(filename::AbstractString,
                              width::MeasureOrNumber,
                              height::MeasureOrNumber;
                              dpi = (B == PNGBackend ? 96 : 72))
        img = Image{B}(open(filename, "w"), width, height, dpi = dpi)
        img.ownedfile = true
        img.filename = filename
        img
    end

    @compat function (::Type{Image{B}}){B}(width::MeasureOrNumber,
                              height::MeasureOrNumber,
                              emit_on_finish::Bool=true;
                              dpi = (B == PNGBackend ? 96 : 72))
        img = Image{B}(IOBuffer(), width, height, emit_on_finish, dpi = dpi)
        img
    end
end

const PNG = Image{PNGBackend}
const PDF = Image{PDFBackend}
const PS =  Image{PSBackend}
const CAIROSURFACE =  Image{CairoBackend}

function canbatch(img::Image)
    for vp in values(img.vector_properties)
        if !isnull(vp)
            return false
        end
    end
    return true
end

# convert compose absolute units (millimeters) to the absolute units used by the
# cairo surface (pixels for PNG, points for all others)

function absolute_native_units{B}(img::Image{B}, u::Float64)
    img.ppmm * u
end

surface(img::Image) = img.surface

function Measures.width(img::Image)
    (Cairo.width(img.surface) / img.ppmm) * mm
end

function Measures.height(img::Image)
    (Cairo.height(img.surface) / img.ppmm) * mm
end

function iswithjs(img::Image)
    return false
end

function iswithousjs(img::Image)
    return true
end

finish{B <: ImageBackend}(::Type{B}, img::Image) = nothing

function finish(::Type{PNGBackend}, img::Image)
    Cairo.write_to_png(img.surface, img.out)
end

function finish{B <: ImageBackend}(img::Image{B})
    if img.finished
        return
    end

    if(img.owns_surface)
        Cairo.destroy(img.ctx)
    end
    finish(B, img)
    if(img.owns_surface)
        Cairo.destroy(img.surface)
    end

    img.finished = true

    if img.emit_on_finish && typeof(img.out) == IOBuffer
        display(img)
    end

    if method_exists(flush, (typeof(img.out),))
        flush(img.out)
    end

    if img.ownedfile
        close(img.out)
    end
end

function newsurface{B}(::Type{B}, out, width, height)
        local surface::CairoSurface
        if B == SVGBackend
            surface = CairoSVGSurface(out, width, height)
        elseif B == PNGBackend
            surface = CairoARGBSurface(round(Integer, width), round(Integer, height))
        elseif B == PDFBackend
            surface = CairoPDFSurface(out, width, height)
        elseif B == PSBackend
            surface = CairoEPSSurface(out, width, height)
        elseif B == CairoBackend
            surface = out
        else
            error("Unkown Cairo backend.")
        end

        if Cairo.status(surface) != Cairo.STATUS_SUCCESS
            error("Unable to create cairo surface.")
        end

        surface
end

function reset{B}(img::Image{B})
    if !img.owns_surface
        error("Backend can't be reused since an external cairo surface is being used.")
    end

    if img.ownedfile
        img.out = open(img.filename, "w")
    else
        try
            seekstart(img.out)
        catch
            error("Backend can't be reused, since the output stream is not seekable.")
        end
    end

    img.surface = newsurface(B, img.out, img.width, img.height)
    img.ctx = CairoContext(img.surface)
    img.finished = false
end

function isfinished(img::Image)
    img.finished
end

function root_box(img::Image)
    BoundingBox(width(img), height(img))
end

@compat function show(io::IO, ::MIME"image/png", img::PNG)
    write(io, String(take!(img.out)))
end

@compat function show(io::IO, ::MIME"application/pdf", img::PDF)
    write(io, String(take!(img.out)))
end

@compat function show(io::IO, ::MIME"application/postscript", img::PS)
    write(io, String(take!(img.out)))
end

# Applying Properties
# -------------------

function push_property_frame(img::Image, properties::Vector{Property})
    if isempty(properties)
        return
    end

    frame = ImagePropertyFrame()
    applied_properties = Set{Type}()
    scalar_properties = Array{Property}(0)
    for property in properties
        if isscalar(property) && !(typeof(property) in applied_properties)
            push!(scalar_properties, property)
            push!(applied_properties, typeof(property))
            frame.has_scalar_properties = true
        elseif !isscalar(property)
            frame.vector_properties[typeof(property)] = property
            img.vector_properties[typeof(property)] = property
        end
    end
    push!(img.property_stack, frame)
    if isempty(scalar_properties)
        return
    end

    save_property_state(img)
    for property in scalar_properties
        apply_property(img, property.primitives[1])
    end
end

function pop_property_frame(img::Image)
    @assert !isempty(img.property_stack)
    frame = pop!(img.property_stack)

    if frame.has_scalar_properties
        restore_property_state(img)
    end

    for (propertytype, property) in frame.vector_properties
        img.vector_properties[propertytype] = Nullable{Property}()
        for i in length(img.property_stack):-1:1
            if haskey(img.property_stack[i].vector_properties, propertytype)
                img.vector_properties[propertytype] =
                    img.property_stack[i].vector_properties[propertytype]
            end
        end
    end
end

function save_property_state(img::Image)
    push!(img.state_stack,
        ImagePropertyState(
            img.stroke,
            img.fill,
            img.stroke_dash,
            img.stroke_linecap,
            img.stroke_linejoin,
            img.visible,
            img.linewidth,
            img.fontsize,
            img.font,
            img.clip))
    Cairo.save(img.ctx)
end

function restore_property_state(img::Image)
    state = pop!(img.state_stack)
    img.stroke = state.stroke
    img.fill = state.fill
    img.stroke_dash = state.stroke_dash
    img.stroke_linecap = state.stroke_linecap
    img.stroke_linejoin = state.stroke_linejoin
    img.visible = state.visible
    img.linewidth = state.linewidth
    img.fontsize = state.fontsize
    img.font = state.font
    img.clip = state.clip
    Cairo.restore(img.ctx)
end

# Return true if the vector properties need to be pushed and popped, rather
# than simply applied.
function vector_properties_require_push_pop(img::Image)
    for (propertytype, property) in img.vector_properties
        propertytype
        if in(propertytype, [Property{FontPrimitive},
                             Property{FontSizePrimitive},
                             Property{ClipPrimitive}])
            return true
        end
    end
    return false
end

function push_vector_properties(img::Image, idx::Int)
    save_property_state(img)
    for (propertytype, property) in img.vector_properties
        if isnull(property)
            continue
        end
        primitives = get(property).primitives
        if idx > length(primitives)
            error("Vector form and vector property differ in length. Can't distribute.")
        end
        apply_property(img, primitives[idx])
    end
end

function pop_vector_properties(img::Image)
    restore_property_state(img)
end

function apply_property(img::Image, p::StrokePrimitive)
    img.stroke = p.color
end

function apply_property(img::Image, p::FillPrimitive)
    img.fill = p.color
end

function apply_property(img::Image, p::FillOpacityPrimitive)
    img.fill = RGBA{Float64}(color(img.fill), p.value)
end

function apply_property(img::Image, p::StrokeOpacityPrimitive)
    img.stroke = RGBA{Float64}(color(img.stroke), p.value)
end

function apply_property(img::Image, p::StrokeDashPrimitive)
    img.stroke_dash = map(v -> absolute_native_units(img, v.value), p.value)
end

function apply_property(img::Image, p::StrokeLineCapPrimitive)
    img.stroke_linecap = p.value
end

function apply_property(img::Image, p::StrokeLineJoinPrimitive)
    img.stroke_linejoin = p.value
end

function apply_property(img::Image, p::VisiblePrimitive)
    img.visible = p.value
end

function apply_property(img::Image, property::LineWidthPrimitive)
    img.linewidth = property.value
    Cairo.set_line_width(
        img.ctx,
        absolute_native_units(img, property.value.value))
end

function apply_property(img::Image, property::FontPrimitive)
    img.font = property.family

    font_desc = ccall((:pango_layout_get_font_description, Cairo._jl_libpango),
                      Ptr{Void}, (Ptr{Void},), img.ctx.layout)

    if font_desc == C_NULL
        size = absolute_native_units(img, default_font_size.value)
    else
        size = ccall((:pango_font_description_get_size, Cairo._jl_libpango),
                     Cint, (Ptr{Void},), font_desc)
    end

    Cairo.set_font_face(img.ctx,
        @sprintf("%s %0.2fpx", property.family, size / PANGO_SCALE))
end

function apply_property(img::Image, property::FontSizePrimitive)
    img.fontsize = property.value

    font_desc = ccall((:pango_layout_get_font_description, Cairo._jl_libpango),
                      Ptr{Void}, (Ptr{Void},), img.ctx.layout)

    if font_desc == C_NULL
        family = "sans"
    else
        family = ccall((:pango_font_description_get_family, Cairo._jl_libpango),
                       Ptr{UInt8}, (Ptr{Void},), font_desc)
        family = unsafe_string(family)
    end

    Cairo.set_font_face(img.ctx,
        @sprintf("%s %.2fpx",
            family, absolute_native_units(img, property.value.value)))
end

function apply_property(img::Image, property::ClipPrimitive)
    if isempty(property.points); return; end
    move_to(img, property.points[1])
    for point in property.points[2:end]
        line_to(img, point)
    end
    close_path(img)
    Cairo.clip(img.ctx)
    img.clip = property
end

# No-op SVG+JS only properties
function apply_property(img::Image, property::JSIncludePrimitive)
end

function apply_property(img::Image, property::JSCallPrimitive)
end

function apply_property(img::Image, property::SVGIDPrimitive)
end

function apply_property(img::Image, property::SVGClassPrimitive)
end

function apply_property(img::Image, property::SVGAttributePrimitive)
end

# Cairo Wrappers
# --------------

function current_point(img::Image)
    x = Array{Float64}(1)
    y = Array{Float64}(1)
    ccall((:cairo_get_current_point, Cairo._jl_libcairo), Void,
          (Ptr{Void}, Ptr{Float64}, Ptr{Float64}), img.ctx.ptr, x, y)
    return ((x[1] / img.ppmm)*mm, (x[2] / img.ppmm)*mm)
end

function move_to(img::Image, point::AbsoluteVec2)
    Cairo.move_to(
        img.ctx,
        absolute_native_units(img, point[1].value),
        absolute_native_units(img, point[2].value))
end

function rel_move_to(img::Image, point::AbsoluteVec2)
    Cairo.rel_move_to(
        img.ctx,
        absolute_native_units(img, point[1].value),
        absolute_native_units(img, point[2].value))
end

function line_to(img::Image, point::AbsoluteVec2)
    Cairo.line_to(
        img.ctx,
        absolute_native_units(img, point[1].value),
        absolute_native_units(img, point[2].value))
end

function rel_line_to(img::Image, point::AbsoluteVec2)
    Cairo.rel_line_to(
        img.ctx,
        absolute_native_units(img, point[1].value),
        absolute_native_units(img, point[2].value))
end

function rel_curve_to(img::Image, ctrl1::AbsoluteVec2, ctrl2::AbsoluteVec2, to::AbsoluteVec2)
    Cairo.rel_curve_to(
        img.ctx,
        absolute_native_units(img, ctrl1[1].value),
        absolute_native_units(img, ctrl1[2].value),
        absolute_native_units(img, ctrl2[1].value),
        absolute_native_units(img, ctrl2[2].value),
        absolute_native_units(img, to[1].value),
        absolute_native_units(img, to[2].value),)
end

function rectangle(img::Image, corner::AbsoluteVec2,
                   width::AbsoluteLength, height::AbsoluteLength)
    Cairo.rectangle(img.ctx,
                    absolute_native_units(img, corner[1].value),
                    absolute_native_units(img, corner[2].value),
                    absolute_native_units(img, width.value),
                    absolute_native_units(img, height.value))
end

function circle(img::Image, center::AbsoluteVec2, radius::AbsoluteLength)
    Cairo.circle(img.ctx,
                 absolute_native_units(img, center[1].value),
                 absolute_native_units(img, center[2].value),
                 absolute_native_units(img, radius.value))
end

function curve_to(img::Image, ctrl1::AbsoluteVec2, ctrl2::AbsoluteVec2, anchor::AbsoluteVec2)
    Cairo.curve_to(
        img.ctx,
        absolute_native_units(img, ctrl1[1].value),
        absolute_native_units(img, ctrl1[2].value),
        absolute_native_units(img, ctrl2[1].value),
        absolute_native_units(img, ctrl2[2].value),
        absolute_native_units(img, anchor[1].value),
        absolute_native_units(img, anchor[2].value))
end

function close_path(img::Image)
    Cairo.close_path(img.ctx)
end

function arc(img::Image, x::Float64, y::Float64, radius::Float64,
                angle1::Float64, angle2::Float64)
    Cairo.arc(
        img.ctx,
        absolute_native_units(img, x),
        absolute_native_units(img, y),
        absolute_native_units(img, radius),
        angle1, angle2)
end

function arc_negative(img::Image, x::Float64, y::Float64, radius::Float64,
                      angle1::Float64, angle2::Float64)
    Cairo.arc_negative(
        img.ctx,
        absolute_native_units(img, x),
        absolute_native_units(img, y),
        absolute_native_units(img, radius),
        angle1, angle2)
end

function translate(img::Image, tx::Float64, ty::Float64)
    Cairo.translate(
        img.ctx,
        absolute_native_units(img, tx),
        absolute_native_units(img, ty))
end

function scale(img::Image, sx::Float64, sy::Float64)
    Cairo.scale(img.ctx, sx, sy)
end

function rotate(img::Image, theta::Float64)
    Cairo.rotate(img.ctx, theta)
end

function rotate(img::Image, theta::Float64, x::Float64, y::Float64)
    ct = cos(theta)
    st = sin(theta)

    x′ = x - (ct * x - st * y)
    y′ = y - (st * x + ct * y)

    translate(img, x′, y′)
    rotate(img, theta)
end

# Convert native linecap/linejoin enums to the Cairo values.
cairo_linecap(::LineCapButt) = Cairo.CAIRO_LINE_CAP_BUTT
cairo_linecap(::LineCapRound) = Cairo.CAIRO_LINE_CAP_ROUND
cairo_linecap(::LineCapSquare) = Cairo.CAIRO_LINE_CAP_SQUARE
cairo_linejoin(::LineJoinMiter) = Cairo.CAIRO_LINE_JOIN_MITER
cairo_linejoin(::LineJoinBevel) = Cairo.CAIRO_LINE_JOIN_BEVEL
cairo_linejoin(::LineJoinRound) = Cairo.CAIRO_LINE_JOIN_ROUND

function fillstroke(img::Image, strokeonly::Bool=false)
    if !img.visible
        return
    end

    if img.fill.alpha > 0.0 && !strokeonly
        Cairo.set_source_rgba(img.ctx, img.fill.r, img.fill.g, img.fill.b,
                              img.fill.alpha)

        if img.stroke.alpha > 0.0
            Cairo.fill_preserve(img.ctx)
        else
            Cairo.fill(img.ctx)
        end
    end

    if img.stroke.alpha > 0.0
        Cairo.set_source_rgba(img.ctx, img.stroke.r, img.stroke.g,
                                     img.stroke.b, img.stroke.alpha)
        Cairo.set_dash(img.ctx, img.stroke_dash)
        Cairo.set_line_cap(img.ctx, cairo_linecap(img.stroke_linecap))
        Cairo.set_line_join(img.ctx, cairo_linejoin(img.stroke_linejoin))

        Cairo.stroke(img.ctx)
    end

    # if the path wasn't stroked or filled, we should still clear it
    Cairo.new_path(img.ctx)
end

# Form Drawing
# ------------

function draw(img::Image, form::Form)
    if vector_properties_require_push_pop(img)
        for (idx, primitive) in enumerate(form.primitives)
            push_vector_properties(img, idx)
            draw(img, primitive)
            pop_vector_properties(img)
        end
    else
        for (idx, primitive) in enumerate(form.primitives)
            for (propertytype, property) in img.vector_properties
                if isnull(property)
                    continue
                end
                primitives = get(property).primitives
                if idx > length(primitives)
                    error("Vector form and vector property differ in length. Can't distribute.")
                end
                apply_property(img, primitives[idx])
            end
            draw(img, primitive)
        end
    end
end

function draw(img::Image, prim::RectanglePrimitive)
    rectangle(img, prim.corner, prim.width, prim.height)
    fillstroke(img)
end

function draw(img::Image, prim::PolygonPrimitive)
    if length(prim.points) <= 1; return; end
    prev_ok = false
    for (i,p) in enumerate(prim.points)
        ok = isfinite(p[1].value) && isfinite(p[2].value)
        if ok && prev_ok
            line_to(img, p)
        elseif !ok && prev_ok
            close_path(img)
            fillstroke(img)
        else
            move_to(img, p)
        end
        prev_ok = ok
    end
    close_path(img)
    fillstroke(img)
end

function draw(img::Image, prim::Compose.ComplexPolygonPrimitive)
    if isempty(prim.rings); return; end

    for ring in prim.rings
        move_to(img, ring[1])
        for point in ring[2:end]
            line_to(img, point)
        end
        close_path(img)
    end
    fillstroke(img)
end

function draw(img::Image, prim::CirclePrimitive)
    circle(img, prim.center, prim.radius)
    fillstroke(img)
end

function draw(img::Image, prim::EllipsePrimitive)
    cx = prim.center[1].value
    cy = prim.center[2].value
    rx = sqrt((prim.x_point[1].value - cx)^2 +
              (prim.x_point[2].value - cy)^2)
    ry = sqrt((prim.y_point[1].value - cx)^2 +
              (prim.y_point[2].value - cy)^2)
    theta = atan2(prim.x_point[2].value - cy,
                  prim.x_point[1].value - cx)

    if !all(isfinite([cx, cy, rx, ry, theta]))
        return
    end

    save_property_state(img)
    translate(img, cx, cy)
    rotate(img, theta)
    translate(img, -rx, -ry)
    scale(img, 2rx, 2ry)
    arc(img, 0.5, 0.5, 0.5, 0.0, 2pi)
    restore_property_state(img)
    fillstroke(img)
end

function draw(img::Image, prim::LinePrimitive)
    if length(prim.points) <= 1; return; end
    prev_ok = false
    for (i,p) in enumerate(prim.points)
        ok = isfinite(p[1].value) && isfinite(p[2].value)
        if ok && prev_ok
            line_to(img, p)
        else
            move_to(img, p)
        end
        prev_ok = ok
    end
    fillstroke(img, true)
end

function get_layout_size(img::Image)
    width, height = Cairo.get_layout_size(img.ctx)
    return ((width / img.ppmm) * mm, (height / img.ppmm) * mm)
end

function draw(img::Image, prim::TextPrimitive)
    if !img.visible || (img.fill.alpha == 0.0 && img.stroke.alpha == 0.0)
        return
    end

    pos = prim.position
    Cairo.set_text(img.ctx, replace(prim.value, "&", "&"), true)
    width, height = get_layout_size(img)
    pos = (pos[1], pos[2] - height)

    if prim.halign != hleft || prim.valign != vbottom
        if prim.halign == hcenter
            pos = (pos[1] - width/2, pos[2])
        elseif prim.halign == hright
            pos = (pos[1] - width, pos[2])
        end

        if prim.valign == vcenter
            pos = (pos[1], pos[2] + height/2)
        elseif prim.valign == vtop
            pos = (pos[1], pos[2] + height)
        end
    end

    Cairo.set_source_rgba(img.ctx, img.fill.r, img.fill.g, img.fill.b,
                          img.fill.alpha)

    if prim.rot.theta != 0.0
        save_property_state(img)
        rotate(img, prim.rot.theta,
               prim.rot.offset[1].value, prim.rot.offset[2].value)
    end

    move_to(img, pos)
    Cairo.show_layout(img.ctx)

    if prim.rot.theta != 0.0
        restore_property_state(img)
    end
end

function draw(img::Image, prim::CurvePrimitive)
    move_to(img, prim.anchor0)
    curve_to(img, prim.ctrl0, prim.ctrl1, prim.anchor1)
    fillstroke(img, true)
end

function draw(img::Image, prim::BitmapPrimitive)
    error("Embedding bitmaps in Cairo backends (i.e. PNG, PDF, PS) is not supported.")
end

function draw(img::Image, prim::PathPrimitive)
    for op in prim.ops
        draw_path_op(img, op)
    end
    fillstroke(img)
end

function draw_path_op(img::Image, op::MoveAbsPathOp)
    move_to(img, op.to)
end

function draw_path_op(img::Image, op::MoveRelPathOp)
    rel_move_to(img, op.to)
end

function draw_path_op(img::Image, op::ClosePathOp)
    close_path(img)
end

function draw_path_op(img::Image, op::LineAbsPathOp)
    line_to(img, op.to)
end

function draw_path_op(img::Image, op::LineRelPathOp)
    rel_line_to(img, op.to)
end

function draw_path_op(img::Image, op::HorLineAbsPathOp)
    pos = current_point(img)
    line_to(img, (op.x, pos.y))
end

function draw_path_op(img::Image, op::HorLineRelPathOp)
    rel_line_to(img, (op.Δx, 0.0mm))
end

function draw_path_op(img::Image, op::VertLineAbsPathOp)
    pos = current_point(img)
    line_to(img, (pos.x, op.y))
end

function draw_path_op(img::Image, op::VertLineRelPathOp)
    rel_line_to(img, (0.0mm, op.Δy))
end

function draw_path_op(img::Image, op::CubicCurveAbsPathOp)
    curve_to(img, op.ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = op.ctrl2
end

function draw_path_op(img::Image, op::CubicCurveRelPathOp)
    xy = current_point(img)
    rel_curve_to(img, op.ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = (op.ctrl2[1] + xy[1], op.ctrl2[2] + xy[2])
end

function draw_path_op(img::Image, op::CubicCurveShortAbsPathOp)
    xy = current_point(img)
    ctrl1 = img.last_ctrl2_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (2*xy[1] - ctrl1[1], 2*xy[2] - ctrl1[2])
    end
    curve_to(img, ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = op.ctrl2
end

function draw_path_op(img::Image, op::CubicCurveShortRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value

    ctrl1 = img.last_ctrl2_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=(2*x1 - ctrl1[1].value) - x1),
                 Measure(abs=(2*y1 - ctrl1[2].value) - y1))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    rel_curve_to(img, ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point =
        (Measure(abs=op.ctrl2[1].value + xy.x.abs),
         Measure(abs=op.ctrl2[2].value + xy.y.abs))
end

function draw_path_op(img::Image, op::QuadCurveAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    cx, cy = op.ctrl1[1].value, op.ctrl1[2].value
    curve_to(img,
             (Measure(abs=(x1 + 2*cx)/3),
              Measure(abs=(y1 + 2*cy)/3)),
             (Measure(abs=(x2 + 2*cx)/3),
              Measure(abs=(y2 + 2*cy)/3)),
             op.to)
    img.last_ctrl1_point = op.ctrl1
end

function draw_path_op(img::Image, op::QuadCurveRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    cx, cy = op.ctrl1[1].value, op.ctrl1[2].value
    rel_curve_to(img,
                 (Measure(abs=(x1 + 2*cx)/3),
                  Measure(abs=(y1 + 2*cy)/3)),
                 (Measure(abs=(x2 + 2*cx)/3),
                  Measure(abs=(y2 + 2*cy)/3)),
             op.to)
    img.last_ctrl1_point =
        (Measure(abs=op.ctrl1[1].value + xy.x.abs),
         Measure(abs=op.ctrl1[2].value + xy.y.abs))
end

function draw_path_op(img::Image, op::QuadCurveShortAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value

    ctrl1 = img.last_ctrl1_point
    if isnull(img.last_ctrl1_point)
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=2*x1 - ctrl1[1].value),
                 Measure(abs=2*y1 - ctrl1[2].value))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    curve_to(img,
             (Measure(abs=(x1 + 2*cx)/3),
              Measure(abs=(y1 + 2*cy)/3)),
             (Measure(abs=(x2 + 2*cx)/3),
              Measure(abs=(y2 + 2*cy)/3)),
             (Measure(abs=x2), Measure(abs=y2)))
    img.last_ctrl1_point = ctrl1
end

function draw_path_op(img::Image, op::QuadCurveShortRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = x1 + op.to[1].value, y1 + op.to[2].value

    ctrl1 = img.last_ctrl1_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=(2*x1 - ctrl1[1].value) - x1),
                 Measure(abs=(2*y1 - ctrl1[2].value) - y1))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    rel_curve_to(img,
                 (Measure(abs=(x1 + 2*cx)/3),
                  Measure(abs=(y1 + 2*cy)/3)),
                 (Measure(abs=(x2 + 2*cx)/3),
                  Measure(abs=(y2 + 2*cy)/3)),
                 (Measure(abs=x2), Measure(abs=y2)))
    img.last_ctrl1_point =
        (Measure(abs=op.ctrl1[1].value + x1),
         Measure(abs=op.ctrl1[2].value + y1))
end

function draw_path_op(img::Image, op::ArcAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    rx, ry = op.rx.abs, op.ry.abs
    φ = deg2rad(op.rotation)
    draw_endpoint_arc(img, rx, ry, φ, op.largearc, op.sweep, x1, y1, x2, y2)
end

function draw_path_op(img::Image, op::ArcRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = x1 + op.to[1].value, y1 + op.to[2].value
    rx, ry = op.rx.abs, op.ry.abs
    φ = deg2rad(op.rotation)
    draw_endpoint_arc(img, rx, ry, φ, op.largearc, op.sweep, x1, y1, x2, y2)
end

# Draw an SVG style elliptical arc
function draw_endpoint_arc(img::Image, rx::Float64, ry::Float64, φ::Float64,
                           largearc::Bool, sweep::Bool,
                           x1::Float64, y1::Float64,
                           x2::Float64, y2::Float64)
    function uvangle(ux, uy, vx, vy)
        t = (ux * vx + uy * vy) / (sqrt(ux^2 + uy^2) * sqrt(vx^2 + vy^2))
        t = max(min(t, 1.0), -1.0)
        return (ux * vy - uy * vx < 0.0 ? -1 : 1.0) * acos(t)
    end

    # From: http://www.w3.org/TR/SVG/implnote.html#ArcConversionEndpointToCenter
    xm, ym = (x1 - x2)/2, (y1 - y2)/2
    x1p =  cos(φ) * xm + sin(φ) * ym
    y1p = -sin(φ) * xm + cos(φ) * ym

    u = (rx^2 * ry^2 - rx^2 * y1p^2 - ry^2 * x1p^2) / (rx^2 * y1p^2 + ry^2 * x1p^2)
    u = u >= 0.0 ? sqrt(u) : 0.0

    cxp =  u * (rx * y1p) / ry
    cyp = -u * (ry * x1p) / rx
    if sweep == largearc
        cxp = -cxp
        cyp = -cyp
    end
    cx = (x1 + x2)/2 + cos(φ) * cxp - sin(φ) * cyp
    cy = (y1 + y2)/2 + sin(φ) * cxp + cos(φ) * cyp

    θ1 = uvangle(1.0, 0.0, (x1p - cxp) / rx, (y1p - cyp) / ry)
    Δθ = uvangle((x1p - cxp) / rx, (y1p - cyp) / ry,
                 (-x1p - cxp) / rx, (-y1p - cyp) / ry) % (2.0*π)
    if Δθ > 0.0 && !sweep
        Δθ -= 2*π
    elseif Δθ < 0.0 && sweep
        Δθ += 2*π
    end

    Cairo.save(img.ctx)
    Cairo.translate(img.ctx,
                    absolute_native_units(img, cx),
                    absolute_native_units(img, cy))
    Cairo.rotate(img.ctx, φ)
    Cairo.scale(img.ctx, rx, ry)
    if sweep
        arc(img, 0.0, 0.0, 1.0, θ1, θ1 + Δθ)
    else
        arc_negative(img, 0.0, 0.0, 1.0, θ1, θ1 + Δθ)
    end
    Cairo.restore(img.ctx)
end

function draw(img::Image, batch::FormBatch)
    bounds = boundingbox(batch.primitive, img.linewidth, img.font, img.fontsize)
    width = bounds.a[1].value * img.ppmm
    height = bounds.a[2].value * img.ppmm

    xt = -bounds.x0[1].value * img.ppmm
    yt = -bounds.x0[2].value * img.ppmm

    Cairo.save(img.ctx)
    Cairo.reset_clip(img.ctx)
    Cairo.translate(img.ctx, xt, yt)
    Cairo.push_group(img.ctx)
    draw(img, batch.primitive)
    pattern = Cairo.pop_group(img.ctx)
    surface = Cairo.pattern_get_surface(pattern)
    Cairo.restore(img.ctx)

    # reapply the clipping region we just reset
    if !isnull(img.clip)
        apply_property(img, get(img.clip))
    end

    Cairo.set_antialias(img.ctx, Cairo.ANTIALIAS_NONE)
    for offset in batch.offsets
        x = offset[1].value * img.ppmm - xt
        y = offset[2].value * img.ppmm - yt
        Cairo.set_source_surface(img.ctx, surface, x, y)
        Cairo.rectangle(img.ctx, x, y, width, height)
        Cairo.fill(img.ctx)
    end
    Cairo.set_antialias(img.ctx, Cairo.ANTIALIAS_DEFAULT)
end


text:
module ImageQualityIndexes

using MappedArrays, OffsetArrays
using ImageCore, ColorVectorSpace
using ImageCore: NumberLike, Pixel, GenericImage, GenericGrayImage
using ImageDistances, ImageFiltering
using Statistics: mean, std
using Base.Iterators: repeated

include("generic.jl")
include("psnr.jl")
include("ssim.jl")
include("colorfulness.jl")

export
    # generic
    assess,

    # Peak signal-to-noise ratio
    PSNR,
    psnr,

    # Structral Similarity
    SSIM,
    ssim,

    # Colorfulness
    HASLER_AND_SUSSTRUNK_M3,
    hasler_and_susstrunk_m3,
    colorfulness

end # module


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

module FirecrackerApis

using Random
using Dates
using HTTP
using Swagger
import Swagger: set_field!, get_field, isset_field, validate_field, SwaggerApi, SwaggerModel
import Base: convert

include("modelincludes.jl")

include("api_DefaultApi.jl")

# export models
export convert, BootSource
export convert, CpuTemplate
export convert, Drive
export convert, Error
export convert, InstanceActionInfo
export convert, InstanceInfo
export convert, Logger
export convert, MachineConfiguration
export convert, NetworkInterface
export convert, PartialDrive
export convert, RateLimiter
export convert, TokenBucket

# export operations
export convert, DefaultApi

export check_required, set_field!, get_field, isset_field, convert

end


text:
# Multicollinearity Test for regression models:
# see associated pdf and data file
# overview of notes:
# if doing regression, the coefficient of each predictor (indepedent) variable tells you how much on average the response changes
# when the variable changes by a unit. however, when variables are corellated (they are not really independent) it makes inference harder

using RCall
using DataFrames
using CSV
using GLM
using Statistics

# read the first 4 columns (not the columns with S)
df = CSV.read("multicollinearity_regression_models_data.csv" , header=1, select=(i, name) -> i < 5 && return true) |> DataFrame#, select=(i, name) -> i < 5 && return true)
insertcols!(df, :pfat_weight => df.pfat .* df.weight) # create interaction term manually

rmodel = lm(@formula(neck ~ pfat + weight + activity + pfat_weight), df)
# lets calculate the Variance Inflation numbers
# using r2 of the model, calculate VIF

rmodel1 = lm(@formula(pfat ~  weight + activity + pfat_weight), df)
rmodel2 = lm(@formula(weight ~ pfat + activity + pfat_weight), df)
rmodel3 = lm(@formula(activity ~ pfat + weight + pfat_weight), df)
rmodel4 = lm(@formula(pfat_weight ~ pfat + weight + activity), df)
rmodels = (rmodel1, rmodel2, rmodel3, rmodel4)
vifs = @. 1 / (1 - r2(rmodels))

# node that for activity, VIF is 1.05 so this variable has no mulitcollinearity.
# for the other ones, there seems to be definitely multicollinearity
# first, let's remove the structural multicollinearity from the interaction term.
# to do tht, we will first mean center the data.

# there is a bug in Query that messes up the dataframe types and GLM can't read it anymore
# https://discourse.julialang.org/t/query-jl-mutate-command-does-not-preserve-column-types/41231/4
# aa = mean(df.pfat)::Float64
# ab = mean(df.weight)::Float64
# ac = mean(df.activity)::Float64
# ad = mean(df.pfat_weight)::Float64
# df_s = df |> @mutate(pfat_s .= _.pfat .- aa,
#               weight_s .= _.weight .- ab,
#               activity_s .= _.activity .- ac,
#               pfat_weight_s .= _.pfat_weight .- ad) |> DataFrame

# equivalent linq
# @from i in df begin
# @select {i..., pfat_s = mean(df.pfat)}
# @collect DataFrame
# end

# lets use standard method to add columns
df[!, :neck_s] .= mean(df.neck) .- df.neck
df[!, :pfat_s] .= mean(df.pfat) .- df.pfat
df[!, :weight_s] .= mean(df.weight) .- df.weight
df[!, :activity_s] .= mean(df.activity) .- df.activity
df[!, :pfat_weight_s] .=  df.pfat_s .* df.weight_s

smodel = lm(@formula(neck ~ pfat_s + weight_s + activity_s + pfat_weight_s), df)
# lets calculate the Variance Inflation numbers
# using r2 of the model, calculate VIF
smodel1 = lm(@formula(pfat_s ~  weight_s + activity_s + pfat_weight_s), df)
smodel2 = lm(@formula(weight_s ~ pfat_s + activity_s + pfat_weight_s), df)
smodel3 = lm(@formula(activity_s ~ pfat_s + weight_s + pfat_weight_s), df)
smodel4 = lm(@formula(pfat_weight_s ~ pfat_s + weight_s + activity_s), df)
smodels = (smodel1, smodel2, smodel3, smodel4)
vifs = @. 1 / (1 - r2(smodels))
# (3.323870441091722, 4.745648167219831, 1.0530047498114123, 1.9910631706191224)

# you see now how small the VIFs have gotten just by removing the multicollinearity

max_stars_repo_name: affans/notebooks

max_issues_repo_name: affans/notebooks

max_forks_repo_name: affans/notebooks


text:
struct GadgetPhysicalUnits{T}

    x_cgs::Quantity{T}        # position in cm
    x_physical::Quantity{T}        # position in kpc

    v_cgs::Quantity{T}     # velocity in cm/s
    v_physical::Quantity{T}     # velocity in km/s

    m_cgs::Quantity{T}     # mass in g
    m_msun::Quantity{T}     # mass in Msun
    m_physical::Quantity{T}     # mass in 10^10 Msun

    t_s::Quantity{T}          # time in sec
    t_Myr::Quantity{T}        # time in Myr

    E_cgs::Quantity{T}        # energy in erg
    E_eV::Quantity{T}         # energy in eV

    B_cgs::Quantity{T}     # magnetic field in Gauss

    rho_physical::Quantity{T}   # density in 10^10 Msun/kpc^3
    rho_cgs::Quantity{T}   # density in g/cm^3
    rho_ncm3::Quantity{T}    # density in mp/cm^3

    T_K::Quantity{T}           # temperature in K
    T_eV::Quantity{T}            # temperature in eV

    P_th_cgs::Quantity{T}      # thermal pressure in Ba
    P_CR_cgs::Quantity{T}      # cosmic ray pressure in Ba

end

"""
    GadgetPhysicalUnits(l_unit::T=3.085678e21, m_unit::T=1.989e43, v_unit::T=1.e5;
                        a_scale::T=1.0, hpar::T=1.0,
                        γ_th::T=5.0/3.0, γ_CR::T=4.0/3.0, xH::T=0.76) where T

Creates a datatype GadgetPhysicalUnits which holds the conversion factors between comoving code units and physical units.
Stores the unit information which can be converted with Unitful or UnitfulAstro.

# Keyword Arguments
- `a_scale::T = 1.0`:  Cosmological scale factor of the simulation. Can be passed with the header `h` as `h.time`.
- `hpar::T = 1.0`:     Hubble constant as 'little h'. Can be passed with header `h` as `h.h0`.
- `γ_th::T = 5.0/3.0`: Adiabatic index of gas.
- `γ_CR::T = 4.0/3.0`: Adiabatic index of cosmic ray component.
- `xH::T = 0.76`:      Hydrogen fraction of the simulation, if run without chemical model.

# Fields
| Name                         | Meaning                        |
|:---------------------------- |:-----------------------------  |
| `x_cgs::Quantity{T}`         | position in cm                 |
| `x_kpc::Quantity{T}`         | position in kpc                |
| `v_cgs::Quantity{T}`         | velocity in cm/s               |
| `v_kms::Quantity{T}`         | velocity in km/s               |
| `m_cgs::Quantity{T}`         | mass in g                      |
| `m_msun::Quantity{T}`        | mass in Msun                   |
| `t_s::Quantity{T}`           | time in sec                    |
| `t_Myr::Quantity{T}`         | time in Myr                    |
| `E_cgs::Quantity{T}`         | energy in erg                  |
| `E_eV::Quantity{T}`          | energy in eV                   |
| `B_cgs::Quantity{T}`         | magnetic field in Gauss        |
| `rho_cgs::Quantity{T}`       | density in ``g/cm^3``          |
| `rho_ncm3::Quantity{T}`      | density in ``n_p/cm^3``        |
| `T_K::Quantity{T}`           | temperature in K               |
| `T_eV::Quantity{T}`          | temperature in eV              |
| `P_th_cgs::Quantity{T}`      | thermal pressure in Ba         |
| `P_CR_cgs::Quantity{T}`      | cosmic ray pressure in Ba      |

"""
function GadgetPhysicalUnits(l_unit::T=3.085678e21, m_unit::T=1.989e43, v_unit::T=1.e5;
                            a_scale::T=1.0, hpar::T=1.0,
                            γ_th::T=5.0/3.0, γ_CR::T=4.0/3.0, xH::T=0.76) where T

    # Gadget units are given in cgs
    l_unit *= one(T)*u"cm"
    m_unit *= one(T)*u"g"
    v_unit *= one(T)*u"cm/s"

    # convert comoving output to physical units
    x_cgs      = l_unit * a_scale / hpar
    x_physical = x_cgs |> u"kpc"

    v_cgs      = v_unit * sqrt(a_scale)
    v_physical = v_cgs |> u"km/s"

    m_cgs   = m_unit / hpar
    m_msun  = m_cgs |> u"Msun"
    m_physical = m_cgs |> u"Msun"

    t_unit  = l_unit / v_unit
    t_s     = t_unit * sqrt(a_scale) / hpar  # in sec
    t_Myr   = t_s |> u"Myr"

    E_cgs = m_cgs * v_cgs^2 |> u"erg"
    E_eV = E_cgs |> u"eV"

    B_cgs = 1.0u"Gs"    # gadget outputs in cgs

    rho_cgs = m_unit/l_unit^3 * hpar^2 / a_scale^3
    rho_physical = rho_cgs |> u"Msun/kpc^3"

    yhelium = ( 1.0 - xH ) / ( 4.0 * xH )
    mean_mol_weight = (1.0 + 4.0 * yhelium) / (1.0 + 3.0 * yhelium + 1.0)

    n2ne =  (  xH + 0.5 * ( 1.0 - xH ) ) /
            ( 2.0 * xH + 0.75 * ( 1.0 - xH ) )       # conversion n_pat -> n_electrons
    umu  =  4.0 / ( 5.0 * xH + 3.0 )                 # mean molucular weight in hydr. mass

    rho_ncm3 = rho_cgs * n2ne/( umu * 1.0u"mp" ) |> u"cm^-3"
    #rho_ncm3 = rho_cgs |> u"n_p"

    yhelium = ( 1.0 - xH ) / ( 4.0 * xH )
    mean_mol_weight = (1.0 + 4.0 * yhelium) / (1.0 + 3.0 * yhelium + 1.0)

    T_cgs = (γ_th - 1.0) * v_unit^2 * 1.0u"mp" * mean_mol_weight / 1.0u"k" |> u"K"
    T_eV  = T_cgs * 1.0u"k" |> u"eV"

    P_th_cgs = a_scale^(-3) * E_cgs / l_unit^3 * hpar^2  |> u"erg/cm^3"
    P_CR_cgs = a_scale^(-4) * E_cgs / l_unit^3 * hpar^2  |> u"erg/cm^3"

    GadgetPhysicalUnits{T}(x_cgs, x_physical,
        v_cgs, v_physical,
        m_cgs, m_msun, m_physical,
        t_s, t_Myr,
        E_cgs, E_eV,
        B_cgs,
        rho_physical, rho_cgs, rho_ncm3,
        T_cgs, T_eV,
        P_th_cgs,
        P_CR_cgs
        )

end

"""
    GadgetPhysicalUnits( DT::DataType,
                            l_unit::Real=3.085678e21, m_unit::Real=1.989e43, v_unit::Real=1.e5;
                            a_scale::Real=1.0, hpar::Real=1.0,
                            γ_th::Real=5.0/3.0, γ_CR::Real=4.0/3.0, xH::Real=0.76)

Set up a unit struct with a given `DataType`.
"""
GadgetPhysicalUnits( DT::DataType,
                        l_unit::Real=3.085678e21, m_unit::Real=1.989e43, v_unit::Real=1.e5;
                        a_scale::Real=1.0, hpar::Real=1.0,
                        γ_th::Real=5.0/3.0, γ_CR::Real=4.0/3.0, xH::Real=0.76) =
    GadgetPhysicalUnits(DT(l_unit), DT(m_unit), DT(v_unit),
                        a_scale=DT(a_scale), hpar=DT(hpar),
                        γ_th=DT(γ_th), γ_CR=DT(γ_CR), xH=DT(xH))


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://nbviewer.jupyter.org/gist/antimon2/dd23b1e62d5fcd08d91fd420d7ba996d/Mandelbrot.mt.cu.jl.ipynb

# %%
using BenchmarkTools
using Plots
gr(fmt=:png)
using CUDA

plotmandelbrot(m) = heatmap(m; c=reverse(cgrad(:jet1)),
    size=(300, 300), colorbar=false, ticks=false, frame=false)

function mandelbrot(c; maxiters=2^10, threshold_abs2=2^10)
    z = zero(c)
    for i in 1:maxiters
        z = z * z + c
        abs2(z) ≥ threshold_abs2 && return i
    end
    maxiters + 1
end

n = 2^8
x = range(-0.714689, -0.714679; length=n)
y = range( 0.299872,  0.299882; length=n)
c = complex.(x', y)

@show typeof(c)
@time m = mandelbrot.(c)
@time m = mandelbrot.(c)
@show typeof(m)
plotmandelbrot(m)

# %%
c_cuda = CuMatrix{ComplexF64}(c)

@show typeof(c_cuda)
@time m_cuda = collect(mandelbrot.(c_cuda))
@time m_cuda = collect(mandelbrot.(c_cuda))
@show typeof(mandelbrot.(c_cuda))
@show typeof(collect(mandelbrot.(c_cuda)))
plotmandelbrot(m_cuda)

# %%
x32 = range(-0.714689f0, -0.714679f0; length=n)
y32 = range( 0.299872f0,  0.299882f0; length=n)
c32 = complex.(x32', y32)

@show typeof(c32)
@time m32 = mandelbrot.(c32)
@time m32 = mandelbrot.(c32)
@show typeof(m32)
plotmandelbrot(m32)

# %%
c32_cuda = cu(c32)

@show typeof(c32_cuda)
@time m32_cuda = collect(mandelbrot.(c32_cuda))
@time m32_cuda = collect(mandelbrot.(c32_cuda))
@show typeof(mandelbrot.(c32_cuda))
@show typeof(collect(mandelbrot.(c32_cuda)))
plotmandelbrot(m32_cuda)

# %%
function mandelbrot_threads(c)
    m = similar(c, typeof(mandelbrot(c[end])))
    Threads.@threads for i in keys(c)
        m[i] = mandelbrot(c[i])
    end
    m
end

@show Threads.nthreads()
@time m_th = mandelbrot_threads(c)
@time m_th = mandelbrot_threads(c)
plotmandelbrot(m_th)

# %%
@benchmark mandelbrot.($c) # CPU Float64

# %%
@benchmark collect(mandelbrot.($c_cuda)) # GPU Float64 (GPU → GPU → CPU)

# %%
@benchmark collect(mandelbrot.(CuMatrix{ComplexF64}($c))) # GPU Float64 (CPU → GPU → GPU → CPU)

# %%
@benchmark mandelbrot.($c32) # CPU Float32

# %%
@benchmark collect(mandelbrot.($c32_cuda)) # GPU Float32 (GPU → GPU → CPU)

# %%
@benchmark collect(mandelbrot.(cu($c32))) # GPU Float32 (CPU → GPU → GPU → CPU)

# %%
@show Threads.nthreads()
@benchmark mandelbrot_threads($c) # CPU Float64 multi-threaded version with nthreads = 12

# %%
N = 2^10
X = range(-0.714689, -0.714679; length=N)
Y = range( 0.299872,  0.299882; length=N)
C = complex.(X', Y)
X32 = range(-0.714689f0, -0.714679f0; length=N)
Y32 = range( 0.299872f0,  0.299882f0; length=N)
C32 = complex.(X32', Y32)

C_cuda = CuMatrix{ComplexF64}(C)
C32_cuda = cu(C32);

# %%
@benchmark mandelbrot.($C) # CPU Float64

# %%
@benchmark collect(mandelbrot.($C_cuda)) # GPU Float64 (GPU → GPU → CPU)

# %%
@benchmark collect(mandelbrot.(CuMatrix{ComplexF64}($C))) # GPU Float64 (CPU → GPU → GPU → CPU)

# %%
@benchmark mandelbrot.($C32) # CPU Float32

# %%
@benchmark collect(mandelbrot.($C32_cuda)) # GPU Float32 (GPU → GPU → CPU)

# %%
@benchmark collect(mandelbrot.(cu($C32))) # GPU Float32 (CPU → GPU → GPU → CPU)

# %%

max_stars_repo_path: 0014/mandelbrot by cuda.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/mandelbrot by cuda.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/mandelbrot by cuda.jl

max_forks_repo_name: genkuroki/public


text:
import Base: promote_rule

promote_rule(::Type{SparseMatrixCSC{T1,S1}},::Type{SparseMatrixCSC{T2,S2}}) where {T1,T2,S1,S2} = SparseMatrixCSC{promote_type(T1,T2),promote_type(S1,S2)}
promote_rule(::Type{SparseMatrixCSC{T1,S1}},::Type{<:AbstractMatrix{T2}}) where {T1,S1,T2} = Matrix{promote_type(T1,T2)}
promote_rule(::Type{<:AbstractMatrix{T1}},::Type{SparseMatrixCSC{T2,S2}}) where {T1,T2,S2} = Matrix{promote_type(T1,T2)}

promote_rule(::Type{SparseVector{T1,S1}},::Type{SparseVector{T2,S2}}) where {T1,T2,S1,S2} = SparseVector{promote_type(T1,T2),promote_type(S1,S2)}


text:
"""
    looks_like_exchange_reaction(rxn_id::String;
        exclude_biomass = false,
        biomass_strings = _constants.biomass_strings,
        exchange_prefixes = _constants.exchange_prefixes,
    )

A predicate that matches reaction identifiers that look like
exchange or biomass reactions, given the usual naming schemes in common model
repositories. Exchange reactions are identified based on matching prefixes in
the set `exchange_prefixes` and biomass reactions are identified by looking for
occurences of `biomass_strings` in the reaction id.

Also see [`find_exchange_reactions`](@ref).

# Example
```
findall(looks_like_exchange_reaction, reactions(model)) # returns indices
filter(looks_like_exchange_reaction, reactions(model)) # returns Strings

# to use the optional arguments you need to expand the function's arguments
# using an anonymous function
findall(x -> looks_like_exchange_reaction(x; exclude_biomass=true), reactions(model)) # returns indices
filter(x -> looks_like_exchange_reaction(x; exclude_biomass=true), reactions(model)) # returns Strings
```
"""
function looks_like_exchange_reaction(
    rxn_id::String;
    exclude_biomass = false,
    biomass_strings = _constants.biomass_strings,
    exchange_prefixes = _constants.exchange_prefixes,
)::Bool
    any(startswith(rxn_id, x) for x in exchange_prefixes) &&
        !(exclude_biomass && any(occursin(x, rxn_id) for x in biomass_strings))
end

"""
    find_exchange_reactions(m::MetabolicModel; kwargs...)

Shortcut for finding exchange reaction indexes in a model; arguments are
forwarded to [`looks_like_exchange_reaction`](@ref).
"""
find_exchange_reactions(m::MetabolicModel; kwargs...) =
    findall(id -> looks_like_exchange_reaction(id; kwargs...), reactions(m))

"""
    find_exchange_reaction_ids(m::MetabolicModel; kwargs...)

Shortcut for finding exchange reaction identifiers in a model; arguments are
forwarded to [`looks_like_exchange_reaction`](@ref).
"""
find_exchange_reaction_ids(m::MetabolicModel; kwargs...) =
    filter(id -> looks_like_exchange_reaction(id, kwargs...), reactions(m))

"""
    looks_like_biomass_reaction(rxn_id::String;
        exclude_exchanges = false,
        exchange_prefixes = _constants.exchange_prefixes,
        biomass_strings = _constants.biomass_strings,
    )::Bool

A predicate that matches reaction identifiers that look like biomass reactions.
Biomass reactions are identified by looking for occurences of `biomass_strings`
in the reaction id. If `exclude_exchanges` is set, the strings that look like
exchanges (from [`looks_like_exchange_reaction`](@ref)) will not match.

# Example
```
filter(looks_like_biomass_reaction, reactions(model)) # returns strings
findall(looks_like_biomass_reaction, reactions(model)) # returns indices
```
"""
function looks_like_biomass_reaction(
    rxn_id::String;
    exclude_exchanges = false,
    exchange_prefixes = _constants.exchange_prefixes,
    biomass_strings = _constants.biomass_strings,
)::Bool
    any(occursin(x, rxn_id) for x in biomass_strings) &&
        !(exclude_exchanges && any(startswith(rxn_id, x) for x in exchange_prefixes))
end

"""
    find_biomass_reactions(m::MetabolicModel; kwargs...)

Shortcut for finding biomass reaction indexes in a model; arguments are
forwarded to [`looks_like_biomass_reaction`](@ref).
"""
find_biomass_reactions(m::MetabolicModel; kwargs...) =
    findall(id -> looks_like_biomass_reaction(id; kwargs...), reactions(m))

"""
    find_biomass_reaction_ids(m::MetabolicModel; kwargs...)

Shortcut for finding biomass reaction identifiers in a model; arguments are
forwarded to [`looks_like_biomass_reaction`](@ref).
"""
find_biomass_reaction_ids(m::MetabolicModel; kwargs...) =
    filter(id -> looks_like_biomass_reaction(id, kwargs...), reactions(m))

"""
    looks_like_extracellular_metabolite(rxn_id::String;
        extracellular_suffixes = _constants.extracellular_suffixes,
        )::Bool

A predicate that matches metabolite identifiers that look like they are extracellular
metabolites. Extracellular metabolites are identified by `extracellular_suffixes` at the end of the
metabolite id.

# Example
```
filter(looks_like_extracellular_metabolite, metabolites(model)) # returns strings
findall(looks_like_extracellular_metabolite, metabolites(model)) # returns indices
```
"""
function looks_like_extracellular_metabolite(
    met_id::String;
    extracellular_suffixes = _constants.extracellular_suffixes,
)::Bool
    any(endswith(met_id, x) for x in extracellular_suffixes)
end

"""
    find_extracellular_metabolites(m::MetabolicModel; kwargs...)

Shortcut for finding extracellular metabolite indexes in a model; arguments are
forwarded to [`looks_like_extracellular_metabolite`](@ref).
"""
find_extracellular_metabolites(m::MetabolicModel; kwargs...) =
    findall(id -> looks_like_extracellular_metabolite(id; kwargs...), metabolites(m))
"""
    find_extracellular_metabolite_ids(m::MetabolicModel; kwargs...)

Shortcut for finding extracellular metabolite identifiers in a model; arguments are
forwarded to [`looks_like_extracellular_metabolite`](@ref).
"""
find_extracellular_metabolite_ids(m::MetabolicModel; kwargs...) =
    findall(id -> looks_like_extracellular_metabolite(id; kwargs...), metabolites(m))


text:
using ImageCore, ImageFiltering
using ImageTransformations, CoordinateTransformations, Interpolations
using MLDataPattern
using OffsetArrays, StaticArrays, IdentityRanges, MappedArrays, ComputationalResources
using ReferenceTests, Test, TestImages, ImageDistances, Statistics

if isdefined(OffsetArrays, :IdOffsetRange)
    OffsetRange = OffsetArrays.IdOffsetRange
else
    OffsetRange = OffsetArrays.IdentityUnitRange
end

# check for ambiguities
refambs = detect_ambiguities(ImageTransformations, Base, Core)
using Augmentor
ambs = detect_ambiguities(Augmentor, ImageTransformations, Base, Core)
# The 3 is from plain_axes with a Tuple{} (so its spurious)
#@test Set(setdiff(ambs, refambs)) == Set{Tuple{Method,Method}}()
@test length(setdiff(ambs, refambs)) == 3

str_show(obj) = @io2str show(::IO, obj)
str_showcompact(obj) = @io2str show(IOContext(::IO, :compact=>true), obj)
str_showconst(obj) = @io2str Augmentor.showconstruction(::IO, obj)

camera = testimage("cameraman")
cameras = similar(camera, size(camera)..., 2)
copyto!(view(cameras,:,:,1), camera)
copyto!(view(cameras,:,:,2), camera)
square = Gray{N0f8}[0.1 0.2 0.3; 0.4 0.5 0.6; 0.7 0.6 0.9]
square2 = rand(Gray{N0f8}, 4, 4)
rect = Gray{N0f8}[0.1 0.2 0.3; 0.4 0.5 0.6]
checkers = Gray{N0f8}[1 0 1 0 1; 0 1 0 1 0; 1 0 1 0 1]
rgb_rect = rand(RGB{N0f8}, 2, 3)

tests = [
    "tst_compat.jl",
    "tst_utils.jl",
    "operations/tst_channels.jl",
    "operations/tst_dims.jl",
    "operations/tst_convert.jl",
    "operations/tst_mapfun.jl",
    "operations/tst_noop.jl",
    "operations/tst_cache.jl",
    "operations/tst_rotation.jl",
    "operations/tst_shear.jl",
    "operations/tst_flip.jl",
    "operations/tst_crop.jl",
    "operations/tst_resize.jl",
    "operations/tst_scale.jl",
    "operations/tst_zoom.jl",
    "operations/tst_distortions.jl",
    "operations/tst_either.jl",
    "tst_operations.jl",
    "tst_pipeline.jl",
    "tst_augment.jl",
    "tst_augmentbatch.jl",
    "tst_distortedview.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
module DecimalExpansion

export decimal_expansion, to_rational

# DECIMAL EXPANSION
decimal_expansion(R::Rational) = decimal_expansion(R.num, R.den)
function decimal_expansion(numerator::Int, denominator::Int)::String
    denominator == 1 && return "$numerator"

    ordered_steps= Vector{Int}()
    unique_steps = Set{Int}()    # optimization

    decimals = Vector{Int}()

    integer_part::Int = div(numerator, denominator)
    remainder::Int = numerator - denominator*integer_part

    # iterative decimal expansion
    while true
        remainder *= 10
        push!(ordered_steps, remainder)
        length(unique_steps) == push!(unique_steps, remainder) |> length &&
            break

        digit = div(remainder, denominator)
        push!(decimals, digit)

        remainder -= digit*denominator
    end

    # separate recurring / non-recurring parts
    i = findfirst(x->x==remainder, ordered_steps)
    recurring     = decimals[i:end] |> join
    non_recurring = decimals[1:i-1] |> join

    output = integer_part |> string
    output*= "."
    output*= non_recurring
    if recurring != "0"
        if length(recurring) == 1
            output *= "$recurring"^3 * "..."
        else
            output *= "[$recurring]..."
        end
    end

    output
end

# NUMBER TO RATIONAL
function to_rational(x::T) where T<:Number
  int_type = sizeof(x) <= 8 ? Int : BigInt

  La, Lb = zero(int_type), one(int_type)
  Ua, Ub = one(int_type), zero(int_type)

  Ma, Mb = one(int_type), one(int_type)

  while true

    M = Ma/Mb
    if M < x
      La, Lb = Ma, Mb
    elseif M > x
      Ua, Ub = Ma, Mb
    else
      break
    end

    Ma, Mb = La+Ua, Lb+Ub
  end

  Ma//Mb
end

end


text:
# This file is a part of SimilaritySearch.jl

export BeamSearchSpace

@with_kw struct BeamSearchSpace <: AbstractSolutionSpace
    bsize = 8:8:64
    Δ = [0.8, 0.9, 1.0, 1.1]                  # this really depends on the dataset, be careful
    bsize_scale = (s=1.5, p1=0.8, p2=0.8, lower=2, upper=512)  # all these are reasonably values
    Δ_scale = (s=1.07, p1=0.8, p2=0.8, lower=0.6, upper=2.0)  # that should work in most datasets
end

Base.hash(c::BeamSearch) = hash((c.bsize, c.Δ, c.maxvisits))
Base.isequal(a::BeamSearch, b::BeamSearch) = a.bsize == b.bsize && a.Δ == b.Δ && a.maxvisits == b.maxvisits
Base.eltype(::BeamSearchSpace) = BeamSearch
Base.rand(space::BeamSearchSpace) = BeamSearch(bsize=rand(space.bsize), Δ=rand(space.Δ))

function combine(a::BeamSearch, b::BeamSearch)
    bsize = ceil(Int, (a.bsize + b.bsize) / 2)
    Δ = (a.Δ + b.Δ) / 2
    BeamSearch(; bsize, Δ)
end

function mutate(space::BeamSearchSpace, c::BeamSearch, iter)
    bsize = SearchModels.scale(c.bsize; space.bsize_scale...)
    Δ = SearchModels.scale(c.Δ; space.Δ_scale...)
    BeamSearch(; bsize, Δ)
end

@with_kw mutable struct OptimizeParameters <: Callback
    kind::ErrorFunction = ParetoRecall()
    initialpopulation = 16
    params = SearchParams(maxpopulation=16, bsize=4, mutbsize=16, crossbsize=8, tol=-1.0, maxiters=16)
    ksearch::Int32 = 10
    numqueries::Int32 = 64
    space::BeamSearchSpace = BeamSearchSpace()
end

optimization_space(index::SearchGraph) = BeamSearchSpace()

function setconfig!(bs::BeamSearch, index::SearchGraph, perf)
    index.search_algo.bsize = bs.bsize
    index.search_algo.Δ = bs.Δ
    index.search_algo.maxvisits = ceil(Int, perf.visited[end])
end

function runconfig(bs::BeamSearch, index::SearchGraph, q, res::KnnResult, pools)
    search(bs, index, q, res, index.hints, pools; maxvisits = 2index.search_algo.maxvisits)
    # search(bs, index, q, res, index.hints, pools)
end

"""
    execute_callback(opt::OptimizeParameters, index::SearchGraph)

SearchGraph's callback for adjunting search parameters
"""
function execute_callback(opt::OptimizeParameters, index::SearchGraph)
    optimize!(index, opt)
end

function optimize!(
    index::SearchGraph,
    opt::OptimizeParameters;
    queries=nothing,
    numqueries=opt.numqueries,
    ksearch=opt.ksearch,
    initialpopulation=opt.initialpopulation,
    params=opt.params,
    verbose=index.verbose,
)
    optimize!(index, opt.kind, opt.space; queries, ksearch, numqueries, initialpopulation, verbose, params)
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Define beam
    λ = LinRange(200,1000,600) # wavelength range [nm]
    θ = [0.] # angle of incidence [degrees]
    beam = PlaneWave(λ, θ)
    # Define layers
    l0 = LayerTMMO(RIdb.air(beam.λ))
    l1 = LayerTMMO(RIdb.dummy(beam.λ,1.45,0.0); type=:OT, d=1/4.)
    l2 = LayerTMMO(RIdb.dummy(beam.λ,3.45,0.0); type=:OT, d=1/4.)
    layers = [l0, l1, l2, l1, l2, l1, l2, l1, l2, l0]
    # Reference wavelenth [nm]
    λ0 = 700.
    # Call main script
    sol = tmm_optics(beam, layers; λ0=λ0, emfflag=true, h=10, pbgflag=true)
    return sol
end

sol = main()

# Plot spectra
plot(Spectrum1D(),
    sol.beam.λ, [sol.Spectra.Rp, sol.Spectra.Tp, 1.0.-(sol.Spectra.Rp.+sol.Spectra.Tp)],
    label=["Reflectance" "Transmittance" "Absorbance"],
    line=([:solid :dash :dashdot]),
    ylims=(0.0,1.0),
    xlims=(sol.beam.λ[1], sol.beam.λ[end]),
);
gui()

# Plot the EMF pattern for normal incidence
plot(EMF2D(),
    sol.beam.λ, sol.Misc.ℓ, log10.(sol.Field.emfp[:,1,:]),
    title=("Log of EMF intensity"),
);
gui()

# plot the refractive index profile
plot(RIprofile(), sol; plotemf=true)
gui()

# plot the photonic dispersion with custom function
plot(PBGDispersion1D(), sol.Bloch)
gui()

plot(PBGDispersion1Dalt(), sol.Bloch)
gui()

plot(PBGDispersion1Dimre(), sol.Bloch)
gui()


text:
module AbstractLattices

export ∧, ∨, dist, wedge, vee

function wedge end
function vee end

const ∧ = wedge
const ∨ = vee

wedge(x) = x
vee(x) = x

function dist end

end