text:
## selection_tools.jl
export expand_dims, select_region_view, select_region, select_region!, flatten_trailing_dims
export slice

"""
    slice(arr, dim, index)

Return a `N` dimensional slice (where one dimensions has size 1) of the N-dimensional `arr` at the index position
`index` in the `dim` dimension of the array.
It holds `size(out)[dim] == 1`.

Examples
```jldoctest
julia> x = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> NDTools.slice(x, 1, 1)
1×3 view(::Matrix{Int64}, 1:1, 1:3) with eltype Int64:
 1  2  3
```
"""
function slice(arr::AbstractArray{T, N}, dim::Integer, index::Integer) where {T, N}
    inds = slice_indices(axes(arr), dim, index)
    return @view arr[inds...]
end

"""
    slice_indices(a, dim, index)

 # Arguments:
`a` should be the axes obtained by `axes(arr)` of an array.
`dim` is the dimension to be selected and `index` the index of it.

Examples
```jldoctest
julia> NDTools.slice_indices((1:10, 1:20, 1:12, 1:33), 1, 3)
(3:3, 1:20, 1:12, 1:33)
```
"""
function slice_indices(a::NTuple{N, T}, dim::Integer, index::Integer) where {T, N}
    inds = ntuple(i -> i == dim ? (a[i][index]:a[i][index])
                                : (first(a[i]):last(a[i])),
                  Val(N))
    return inds
end

"""
    expand_dims(x, ::Val{N})

Expands the dimensions of an array to a given number of dimensions.

Examples
The result is a 5D array with singleton dimensions at the end
```jldoctest
julia> expand_dims(ones((1,2,3)), Val(5))
1×2×3×1×1 Array{Float64, 5}:
[:, :, 1, 1, 1] =
 1.0  1.0

[:, :, 2, 1, 1] =
 1.0  1.0

[:, :, 3, 1, 1] =
 1.0  1.0
```
"""
function expand_dims(x, ::Val{N}) where N
    return reshape(x, (size(x)..., ntuple(x -> 1, Val(N - ndims(x)))...))
end

"""
    flatten_trailing_dims(arr, max_dim=Val(max_dim))

Flattens (squeezes) the trailing dims. `max_dim` denotes the last dimension to keep.
The implementation uses reshape and thus returns a modified view of the array
referring to the same data.
By default max_dim is adjusted such that a 2N array is squeezed into an N+1 array as needed for a scan.
"""
function flatten_trailing_dims(arr::AbstractArray,
                               ::Val{max_dim}) where max_dim
    reshape(arr,(size(arr)[1:max_dim-1]...,prod(size(arr)[max_dim:end])))
end

"""
    select_region_view(src; new_size=size(src), center=ft_center_diff(size(src)).+1, pad_value=zero(eltype(src)))

selects (extracts) a region of interest (ROI), defined by `new_size` and centered at `center` in the source image. Note that
the number of dimensions can be smaller in `new_size` and `center`, in which case the default values will be insterted
into the missing dimensions. `new_size` does not need to fit into the source array and missing values will be replaced with `pad_value`.

Arguments:
+ `src`. The source array to select from.
+ `new_size`. The size of the array view after the operation finished. By default the original size is assumed
+ `center`. Specifies the center of the new view in coordinates of the old view. By default an alignment of the Fourier-centers is assumed.
+ `pad_value`. Specifies the value which is inserted in case the ROI extends to outside the source area.

The returned results is a mutable view, which allows this method to also be used for writing into a ROI

See also
+ [`select_region`](@ref)

Examples
```jldoctest
julia> select_region_view(ones(3,3),new_size=(7,7),center=(1,3))
7×7 NDTools.MutablePaddedView{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function select_region_view(src::Array{T,N}; new_size=size(src), center=ft_center_diff(size(src)).+1, pad_value=zero(eltype(src))) where {T,N}
    new_size = Tuple(expand_size(new_size, size(src)))
    center = Tuple(expand_size(center, ft_center_diff(size(src)).+1))
    oldcenter = ft_center_diff(new_size).+1
    MutablePaddedView(PaddedView(pad_value, src,new_size, oldcenter .- center.+1)) :: MutablePaddedView{T, N, NTuple{N,Base.OneTo{Int64}}, OffsetArrays.OffsetArray{T, N, Array{T, N}}}
end

"""
    get_src_dst_range(src_size, dst_size, new_size, center)

A helpfer function to calculate the index ranges to copy from source size `src_size` to destination size `dst_size` with the
integer center position of the destination aligning with the position in the source as specified by center.
"""
function get_src_dst_range(src_size, dst_size, new_size, src_center, dst_ctr=dst_size .÷2 .+1)
    ROI_center = (new_size.÷2 .+1)
    src_start = src_center .- ROI_center  .+1 # start of the first pixel to copy (without clipping)
    src_end = src_start .+ new_size .- 1 # the last pixel to copy
    src_start_clip = max.(1, src_start)
    src_end_clip = min.(src_end, src_size)
    if any(src_start_clip .> src_size) || any(src_end_clip .< 1)
        return (1:0), (1:0)
    end
    extra_src_start = max.(0, src_start_clip .- src_start)
    extra_src_end = max.(0, src_end .- src_end_clip)
    copy_size = new_size .- extra_src_start .- extra_src_end

    dst_start = dst_ctr .- ROI_center .+1 .+ extra_src_start
    dst_end = dst_start .+ copy_size .- 1
    dst_end_clip = min.(dst_end, dst_size)
    dst_start_clip = max.(1, dst_start)
    if any(dst_start_clip .> dst_size) || any(dst_end_clip .< 1)
        return (1:0), (1:0)  # returns an empty range for all coordinates
    end

    extra_dst_start = max.(0, dst_start_clip .- dst_start)
    extra_dst_end = max.(0, dst_end .- dst_end_clip)
    src_start_clip = src_start_clip .+ extra_dst_start
    src_end_clip = max.(0, src_end_clip .- extra_dst_end)

    range_src = Tuple((src_start_clip[d]:src_end_clip[d]) for d in 1:length(src_start))
    range_dst = Tuple((dst_start_clip[d]:dst_end_clip[d]) for d in 1:length(dst_start))
    return range_src, range_dst
end

"""
    assign_to!(a,b)

assignes array (or value) `b` to array `a` pointwise.
Helper function to be passed as an operator to functions such as select_region!
"""
function assign_to!(a,b)
    a .= b
end

"""
    add_to!(a,b)

adds array (or value) `b` to array `a` pointwise.
Hlper function to be passed as an operator to functions such as select_region!
"""
function add_to!(a,b)
    a  .+= b
end

"""
    sub_to!(a,b)

subtractss array (or value) `b` from array `a` pointwise and assigns into a.
Hlper function to be passed as an operator to functions such as select_region!
"""
function sub_to!(a,b)
    a  .-= b
end

"""
    mul_to!(a,b)

multiplies array (or value) `b` to array `a` pointwise.
Helper function to be passed as an operator to functions such as select_region!
"""
function mul_to!(a,b)
    a .*= b
end

"""
    div_to!(a,b)

divides array `a` by array (or value) `b` to array a pointwise.
Helper function to be passed as an operator to functions such as select_region!
"""
function div_to!(a,b)
    a ./= b
end

"""
    select_region!(src, dst=nothing, new_size=size(src), center=size(src).÷2 .+1, dst_center=nothing, pad_value=zero(eltype(mat), operator!=assign_to!))

selects (extracts, pads, shifts) a region of interest (ROI), defined by `new_size` and centered with the destination center aligned at
the position `center` in the source image. Note that the number of dimensions in `new_size`,  `center` and `dst_center` can be smaller ,
in which case default values (see below) will be insterted into the missing dimensions. `new_size` does not need to fit into the source array
and missing values will be replaced with `pad_value`, if no `dst` is provided.

As opposed to `select_region()`, this version returns a copy rather than a view or, alternatively, also writes into a destination array `dst`
(`new_size` is then interpreted to refer to the maximally assigned region).
If `nothing` is provided for `dst`, a new array of size `new_size` is created.

Arguments:
+ `src`. The source array to select from.
+ `dst`. The destination array to write into, if provided. By default `dst=nothing` a new array is created. The `dst`array (or new array) is returned.
+ `new_size`. The size of the array view after the operation finished. By default the original size is assumed
+ `center`. Specifies the center of the new view in coordinates of the old view. By default an alignment of the Fourier-center (right center) is assumed.
+ `dst_center`. defines the center coordinate in the destination array which should align with the above source center. If nothing is provided, the right center pixel of the `dst` array or new array is used.
+ `pad_value`. specifies the value which is inserted in case the ROI extends to outside the source area. This is only used, if no `dst` array is provided.
+ `operator!`. allows to provide a user_defined array assignment function. The function my_op!(dst,src) should operator on array views and typically perform the assignment elementwise, overwriting the entries in dst.
                Five such functions are exported by NDTools: `assign_to!`, `add_to!`, `sub_to!`, `mul_to!`, `div_to!`, representing the operations `.=`, `.+=`, `.-=`, `.*=` and `./=` respectively.

The returned results is the destination (or newly created) array.
Note that this version is rather fast, since it consists of only a sinlge sub-array assigment on views, avoiding copy operations.

See also
+ [`select_region`](@ref)

Examples:
```jdoctest
julia> a = ones(5,6);

julia> dst=select_region(a,new_size=(10,10), dst_center=(1,1)) # pad a with zeros to a size of (10,10), but place original center at the corner
10×10 Matrix{Float64}:
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

 julia> select_region!(2 .*a,dst, dst_center=size(dst)) # write a doubled version into the bottom right corner
10×10 Matrix{Float64}:
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  2.0  2.0  2.0  2.0
 0.0  0.0  0.0  0.0  0.0  0.0  2.0  2.0  2.0  2.0
 0.0  0.0  0.0  0.0  0.0  0.0  2.0  2.0  2.0  2.0
```
"""
function select_region!(src, dst; new_size=size(dst),
                        center=size(src).÷2 .+1, dst_center=size(dst).÷ 2 .+1, operator! =assign_to!)
    new_size = Tuple(expand_size(new_size, size(dst)))
    center = Tuple(expand_size(center, size(src).÷2 .+1))
    dst_center = Tuple(expand_size(dst_center, size(dst).÷ 2 .+1))

    range_src, range_dst = get_src_dst_range(size(src),size(dst),new_size,center, dst_center)
    if !isempty(range_dst)
        v_src = @view src[range_src...]
        v_dst = @view dst[range_dst...]
        # more generic than the line below
        operator!(v_dst, v_src)
        # dst[range_dst...] .+= view(src, range_src...)
    end
    return dst
end

"""
    select_region(src; new_size=size(mat), center=ft_center_diff(size(mat)).+1, pad_value=zero(eltype(mat)))

selects (extracts) a region of interest (ROI), defined by `new_size` and centered at `center` in the source image. Note that
the number of dimensions can be smaller in `new_size` and `center`, in which case the default values will be insterted
into the missing dimensions. `new_size` does not need to fit into the source array and missing values will be replaced with `pad_value`.

Arguments:
+ `src`. The array to extract the region from.
+ `new_size`. The size of the array view after the operation finished. By default the original size is assumed
+ `center`. Specifies the center of the new view in coordinates of the old view. By default an alignment of the Fourier-centers is assumed.
+ `pad_value`. Specifies the value which is inserted in case the ROI extends to outside the source area.

The returned result is a newly allocated array of the same type as the src. This is currently faster that select_region_view().

See also
+ [`select_region_view`](@ref)

Examples
```jldoctest
julia> select_region(ones(3,3),new_size=(7,7),center=(1,3))
7×7 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> a = ones((3,3))
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> dst=select_region(a,new_size=(10,10), dst_center=(1,1)) # pad a with zeros to a size of (10,10), but place original center at the corner
10×10 Matrix{Float64}:
 1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function select_region(src::AbstractArray{T,N}; new_size=size(src), center=size(src).÷2 .+1, pad_value=zero(eltype(src)), dst_center = new_size .÷ 2 .+1) where {T,N}
    new_size = Tuple(expand_size(new_size, size(src)))
    # replace missing coordinates with the new center position
    dst_center = Tuple(expand_size(dst_center, new_size .÷ 2 .+1))

    pad_value = eltype(src)(pad_value)

    # if the new_size is fully enclosed, we can simply use similar
    # since there is no pad
    dst = let
        # check whether it is fully enclosing and if center is really in the middle
        # otherwise some indices can be outside of the original array and would need
        # to be padded
        if all(new_size .<= size(src)) && center == size(src) .÷ 2 .+ 1
            similar(src, new_size)
        else
            arr_n = similar(src, new_size)
            fill!(arr_n, pad_value)
        end
    end
    select_region!(src,dst;new_size=new_size, center=center, dst_center=dst_center)
    return dst
end


text:
# Randomized matrix initialization

function randinit(m::Int,
                  n::Int,
                  type::String="Gaussian")
    if type == "Gaussian"
        return randn(Float64, (m,n))
    elseif type == "Uniform"
        return rand(Float64, (m, n))
    elseif type == "Rademacher"
        return -1
    end
end

ext: jl


text:
k = 2 # order of interaction to investigate
r = 2 # number of gRNAs per combinatorial gRNA/Cas construct

x = 20 # number of target genes
g = 6 # number of gRNAs per target gene
n_gRNA_total = x * g  # total number of gRNAs in the experiment
ϵ_KO = 0.8  # global knockout efficiency; fraction of mutations leading to effective gene knockout

Random.seed!(1)
ρ = 2; # ratio of the frequency of the most abundant gRNA and the frequency of the least abundant gRNA
l = 50; u = ρ*l
m = (l+u)/2;
sd = (u-l)/2;
p_gRNA_freq = gRNA_frequency_distribution(m, sd, l, u, n_gRNA_total; normalize=true)  # generate gRNA frequency distribution in construct library

Random.seed!(1)
f_act = 0.9  # fraction of all gRNAs that is active
ϵ_edit_act = 0.95;  # average genome editing efficiency of active gRNAs
ϵ_edit_inact = 0.1; # average genome editing efficiency of inactive gRNAs
sd_act = 0.01  # standard deviation
p_gRNA_edit = gRNA_edit_distribution(f_act, ϵ_edit_act, ϵ_edit_inact, sd_act, n_gRNA_total; visualize=false); # generate genome editing distribution of the gRNAs

max_stars_repo_name: kirstvh/MultiplexCrisprDOE

max_issues_repo_name: kirstvh/MultiplexCrisprDOE

max_forks_repo_name: kirstvh/MultiplexCrisprDOE


text:
using Printf
using Logging
using Oceananigans

Logging.global_logger(OceananigansLogger())

function simulate_lid_driven_cavity(; Re, N, end_time)
    topology = (Flat, Bounded, Bounded)
    domain = (y=(0, 1), z=(0, 1))
    grid = RectilinearGrid(topology=topology, size=(N, N); domain...)

    v_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(1),
                                    bottom = ValueBoundaryCondition(0))

    w_bcs = FieldBoundaryConditions(north = ValueBoundaryCondition(0),
                                    south = ValueBoundaryCondition(0))

    model = NonhydrostaticModel(
                       grid = grid,
                   buoyancy = nothing,
                    tracers = nothing,
                   coriolis = nothing,
        boundary_conditions = (v=v_bcs, w=w_bcs),
                    closure = ScalarDiffusivity(ν=1/Re)
    )

    u, v, w = model.velocities
    ζ = Field(∂y(w) - ∂z(v))

    fields = (; v, w, ζ)
    global_attributes = Dict("Re" => Re)
    output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s"))
    field_output_writer =
        NetCDFOutputWriter(model, fields, filepath="lid_driven_cavity_Re$Re.nc", schedule=TimeInterval(0.1),
                           global_attributes=global_attributes, output_attributes=output_attributes)

    max_Δt = 0.25 * model.grid.Δyᵃᶜᵃ^2 * Re / 2  # Make sure not to violate diffusive CFL.
    wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=max_Δt)

    cfl = AdvectiveCFL(wizard)
    dcfl = DiffusiveCFL(wizard)

    simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=print_progress,
                            iteration_interval=20, parameters=(cfl=cfl, dcfl=dcfl))

    simulation.output_writers[:fields] = field_output_writer

    run!(simulation)

    return simulation
end

function print_progress(simulation)
    model = simulation.model
    cfl, dcfl = simulation.parameters

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / simulation.stop_time)

    # Find maximum velocities.
    vmax = maximum(abs, model.velocities.v)
    wmax = maximum(abs, model.velocities.w)

    i, t = model.clock.iteration, model.clock.time
    @info @sprintf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e",
                   progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt)

    return nothing
end

simulate_lid_driven_cavity(Re=100,   N=128, end_time=15)
simulate_lid_driven_cavity(Re=400,   N=128, end_time=20)
simulate_lid_driven_cavity(Re=1000,  N=128, end_time=25)
simulate_lid_driven_cavity(Re=3200,  N=128, end_time=50)
simulate_lid_driven_cavity(Re=5000,  N=256, end_time=50)
simulate_lid_driven_cavity(Re=7500,  N=256, end_time=75)
simulate_lid_driven_cavity(Re=10000, N=256, end_time=100)


text:
module WarpFlow # module for exporing Guassian kernel function and derivatives

export Flow, TrFlow

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

# Custom Types

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

typealias Array1 Array{Array{Float64,1},1} # for lists of vectors

typealias Array2 Array{Array{Float64,2},1} # for lists of jacobians

abstract AbstractFlow

immutable Flow <: AbstractFlow
	kappa::Array1
	eta_coeff::Array1
	phix::Array1
	Dphix::Array2
	dim::Int64
end

immutable TrFlow <: AbstractFlow
	kappa::Array1
	eta_coeff::Array1
	phix::Array1
	Dphix::Array2
	dlkappa::Array1
	dleta_coeff::Array1
	dlphix::Array1
	dlDphix::Array2
	dim::Int64
end

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

# type constructors

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

# zero Flow constructor
Flow(dim, nKap, nPhi) = Flow(array1(dim, nKap), array1(dim, nKap), array1(dim, nPhi), array2(dim, nPhi), dim)

# zero TrFlow constructor
TrFlow(dim, nKap, nPhi) = TrFlow(array1(dim, nKap), array1(dim, nKap), array1(dim, nPhi), array2(dim, nPhi), array1(dim, nKap), array1(dim, nKap), array1(dim, nPhi), array2(dim, nPhi), dim) # zero Flow constructor

# this is used to initialize the transpose flow at time 1 in get_grad
function TrFlow(yin::Flow, target::Function)
	nKap = length(yin.kappa)
	nPhi = length(yin.phix)
	dleta_coeff = array1(yin.dim, nKap)
	dlkappa     = array1(yin.dim, nKap)
	dlphix      = Array{Float64,1}[target(yin.phix)[2][i]/nPhi for i = 1:nPhi]
	dlDphix     = Array{Float64,2}[inv(yin.Dphix[i])/nPhi for i = 1:nPhi]
	TrFlow(yin.kappa, yin.eta_coeff, yin.phix, yin.Dphix, dlkappa, dleta_coeff, dlphix, dlDphix, yin.dim)
end

array1(dim, k) = Array{Float64,1}[zeros(dim) for i=1:k] # constructor for Array1

array2(dim, k) = Array{Float64,2}[zeros(dim, dim) for i=1:k] # constructor for Array2

array2eye(dim, k) = Array{Float64,2}[eye(dim) for i=1:k] # constructor for Array2

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

#  extend basic functions to Flow types

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

import Base: +, -, *, /

# defines + between Flow and TrFlow... used to update with  the return value of get_grad
function +(yin::Flow, zin::TrFlow)
	nKap = length(yin.kappa) # TODO: thro an error if dims don't match
	nPhi = length(yin.phix)
	yout = Flow(yin.dim, 0, 0) # initialize an empty Flow
	for i = 1:nKap
		push!(yout.kappa,     yin.kappa[i]     + zin.dlkappa[i])
		push!(yout.eta_coeff, yin.eta_coeff[i] + zin.dleta_coeff[i])
	end
	# phix and Dphix don't get updated from TrFlow
	for i = 1:nPhi
		push!(yout.phix,   yin.phix[i])
		push!(yout.Dphix,  yin.Dphix[i])
	end
	yout
end

# defines +, -, *, / between Flow and Flow, used in ode
for opt = (:+, :-, :*, :/)
	@eval begin
		function $opt(yin::Flow, zin::Flow)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = Flow(yin.dim, 0, 0)
			for i = 1:nKap
				push!(yout.kappa,     $opt(yin.kappa[i], zin.kappa[i]))
				push!(yout.eta_coeff, $opt(yin.eta_coeff[i], zin.eta_coeff[i]))
			end
			for i = 1:nPhi
				push!(yout.phix,  $opt(yin.phix[i],  zin.phix[i]))
				push!(yout.Dphix, $opt(yin.Dphix[i], zin.Dphix[i]))
			end
			yout
		end # end function definition
	end
end

# defines +, -, *, / between TrFlow and TrFlow, used in ode
for opt = (:+, :-, :*, :/)
	@eval begin
		function $opt(yin::TrFlow, zin::TrFlow)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = TrFlow(yin.dim, 0, 0)
			for i = 1:nKap
				push!(yout.kappa,       $opt(yin.kappa[i], zin.kappa[i]))
				push!(yout.eta_coeff,   $opt(yin.eta_coeff[i], zin.eta_coeff[i]))
				push!(yout.dlkappa,     $opt(yin.dlkappa[i], zin.dlkappa[i]))
				push!(yout.dleta_coeff, $opt(yin.dleta_coeff[i], zin.dleta_coeff[i]))
			end
			for i = 1:nPhi
				push!(yout.phix,   $opt(yin.phix[i], zin.phix[i]))
				push!(yout.Dphix,  $opt(yin.Dphix[i], zin.Dphix[i]))
				push!(yout.dlphix, $opt(yin.dlphix[i], zin.dlphix[i]))
				push!(yout.dlDphix,$opt(yin.dlDphix[i], zin.dlDphix[i]))
			end
			yout
		end # end function definition
	end
end

# defines +, -, *, / between Flow and Float64
for opt = (:+, :-, :*, :/)
	@eval begin
		function $opt(yin::Flow, a::Real)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = Flow(yin.dim, 0, 0)
			for i = 1:nKap
				push!(yout.kappa, $opt(yin.kappa[i],a))
				push!(yout.eta_coeff, $opt(yin.eta_coeff[i],a))
			end
			for i = 1:nPhi
				push!(yout.phix, $opt(yin.phix[i],a))
				push!(yout.Dphix, $opt(yin.Dphix[i],a))
			end
			yout
		end # end function definition
	end
end

# defines +, -, *, / between TrFlow and Float64
for opt = (:+, :-, :*, :/)
	@eval begin
		function $opt(yin::TrFlow, a::Real)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = TrFlow(yin.dim, 0, 0)
			for i = 1:nKap
				push!(yout.kappa, $opt(yin.kappa[i],a))
				push!(yout.eta_coeff, $opt(yin.eta_coeff[i],a))
				push!(yout.dlkappa, $opt(yin.dlkappa[i],a))
				push!(yout.dleta_coeff, $opt(yin.dleta_coeff[i],a))
			end
			for i = 1:nPhi
				push!(yout.phix, $opt(yin.phix[i],a))
				push!(yout.Dphix, $opt(yin.Dphix[i],a))
				push!(yout.dlphix, $opt(yin.dlphix[i],a))
				push!(yout.dlDphix, $opt(yin.dlDphix[i],a))
			end
			yout
		end # end function definition
	end
end

# +, * are associative for both Flow and TrFlow
for opt = (:+, :*)
	@eval $opt{T<:AbstractFlow}(a::Real, yin::T) = $opt(yin, a)
end

-{T<:AbstractFlow}(a::Real, yin::T) = -1.0 * -(yin, a)

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

# code for extending ode algorithms to Flows and TrFlows

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

function get_grad(y0::Flow, target::Function, lambda, sigma)
    dydt(t,y) = d_forward_dt(y, sigma)
    (t1,y1) = ode23_abstract_end(dydt, [0,1], y0) # Flow y0 forward to time 1

    z1 = TrFlow(y1, target) # Initalize the transpose flow at time 1

    dzdt(t,z) = d_backward_dt(z, sigma)
    (t0,z0) = ode23_abstract_end(dzdt, [0,1], z1) # Flow z1 backward from time 1 to 0

    # add the regularization
    nKap = length(y0.kappa)
    for i = 1:nKap, j = 1:nKap
        z0.dleta_coeff[i] -= lambda * y0.eta_coeff[j] * R(y0.kappa[i], y0.kappa[j], sigma)
        z0.dlkappa[i]     -= lambda * dot(y0.eta_coeff[i], y0.eta_coeff[j]) * gradR(y0.kappa[i], y0.kappa[j], sigma)
    end
    z0 # return the TrFlow
end

function d_forward_dt(yin::Flow, sigma)
    nKap = length(yin.kappa)
    nPhi = length(yin.phix)
    yout = Flow(yin.dim, nKap, nPhi) # intialize the return Flow
    # eta and kappa
    for i = 1:nKap, j = 1:nKap
        yout.eta_coeff[i] -= dot(yin.eta_coeff[i], yin.eta_coeff[j]) * gradR(yin.kappa[i], yin.kappa[j], sigma)
        yout.kappa[i] += yin.eta_coeff[j] * R(yin.kappa[i], yin.kappa[j], sigma)
    end
    # phi
    for i = 1:nPhi, j = 1:nKap
        yout.phix[i] += yin.eta_coeff[j] * R(yin.phix[i], yin.kappa[j], sigma)
        yout.Dphix[i] += yin.eta_coeff[j] * transpose(gradR(yin.phix[i], yin.kappa[j], sigma)) * yin.Dphix[i]
    end
    yout
end

function d_backward_dt(yin::TrFlow, sigma)
    nKap = length(yin.kappa)
    nPhi = length(yin.phix)
    yout = TrFlow(yin.dim, nKap, nPhi)

    # the following are re-used
    gradRkk = Array{Float64,1}[WarpFlow.gradR(yin.kappa[i], yin.kappa[j],sigma) for i = 1:nKap, j = 1:nKap]
    gradRpk = Array{Float64,1}[WarpFlow.gradR(yin.phix[i], yin.kappa[j],sigma) for i = 1:nPhi, j = 1:nKap]
    Rkk = Float64[WarpFlow.R(yin.kappa[i], yin.kappa[j],sigma) for i = 1:nKap, j = 1:nKap]
    Rpk = Float64[WarpFlow.R(yin.phix[i], yin.kappa[j],sigma) for i = 1:nPhi, j = 1:nKap]
    g1g2Rkkdl = Array{Float64,1}[WarpFlow.g1g2R(yin.kappa[i],yin.kappa[j], sigma)*yin.dleta_coeff[i] for i = 1:nKap, j = 1:nKap]# note g1g1Rkk = - g1g2Rkk
    g1g1RpkD  = Array{Float64,1}[WarpFlow.g1g1R(yin.phix[i],yin.kappa[j],sigma)*yin.Dphix[i][:,col]*transpose(yin.eta_coeff[j])*yin.dlDphix[i][:,col]  for i = 1:nPhi, j = 1:nKap, col = 1:yin.dim]

    # eta and kappa
    for i = 1:nKap, j = 1:nKap
        yout.eta_coeff[i] += dot(yin.eta_coeff[i], yin.eta_coeff[j]) * gradRkk[i,j]
        yout.kappa[i] -= yin.eta_coeff[j] * Rkk[i,j]
    end
    # phi and Dphi
    for i = 1:nPhi, j = 1:nKap
        yout.phix[i] -= yin.eta_coeff[j] * Rpk[i,j]
        yout.Dphix[i] -= yin.eta_coeff[j] * transpose(gradRpk[i,j]) * yin.Dphix[i]
    end
    # dlphix and dlDphix
    for i = 1:nPhi, j = 1:nKap
		yout.dlphix[i] += gradRpk[i,j] * transpose(yin.eta_coeff[j]) * yin.dlphix[i]
        for col = 1:yin.dim
            yout.dlphix[i] +=  g1g1RpkD[i,j,col] #g1g1R(yin.phix[i], yin.kappa[j], sigma) * yin.Dphix[i][:,col] * transpose(yin.eta_coeff[j]) * yin.dlDphix[i][:,col]
            yout.dlDphix[i][:,col] +=  gradRpk[i,j] * transpose(yin.eta_coeff[j]) * yin.dlDphix[i][:,col]
        end
    end
    # dlkappa
    for i = 1:nKap, j = 1:nKap
        yout.dlkappa[i] += dot(yin.eta_coeff[i], yin.eta_coeff[j]) .* g1g2Rkkdl[i,j] #(g1g1R(yin.kappa[i], yin.kappa[j], sigma) * yin.dleta_coeff[i])
        yout.dlkappa[i] -= dot(yin.eta_coeff[i], yin.eta_coeff[j]) .* g1g2Rkkdl[j,i] #(g1g2R(yin.kappa[j], yin.kappa[i], sigma) * yin.dleta_coeff[j])
        yout.dlkappa[i] += gradRkk[i,j] * transpose(yin.eta_coeff[j]) * yin.dlkappa[i]
        yout.dlkappa[i] -= gradRkk[j,i] * transpose(yin.eta_coeff[i]) * yin.dlkappa[j]
        yout.dleta_coeff[i] -= yin.eta_coeff[j] * transpose(gradRkk[i,j]) * yin.dleta_coeff[i]
        yout.dleta_coeff[i] -= yin.eta_coeff[j] * transpose(gradRkk[j,i]) * yin.dleta_coeff[j]
        yout.dleta_coeff[i] += Rkk[j,i] * yin.dlkappa[j]
    end
    for i = 1:nKap, j = 1:nPhi
        yout.dlkappa[i] -= gradRpk[j,i] * transpose(yin.eta_coeff[i]) * yin.dlphix[j]
        yout.dleta_coeff[i] += Rpk[j,i] * yin.dlphix[j]
        for col = 1:yin.dim
            yout.dlkappa[i] -= g1g1RpkD[j,i,col] # -g1g2R(yin.phix[j], yin.kappa[i],sigma) * yin.Dphix[j][:,col] * transpose(yin.eta_coeff[i]) * yin.dlDphix[j][:,col]
            yout.dleta_coeff[i] += dot(yin.Dphix[j][:,col], gradRpk[j,i]) * yin.dlDphix[j][:,col]
        end
    end
    yout
end

function ode23_abstract_end(F::Function, tspan::AbstractVector, y0)
    rtol = 1.e-4
    atol = 1.e-6
    threshold = atol / rtol

    t = tspan[1]
    tfinal = tspan[end]
    tdir = sign(tfinal - t)
    hmax = abs(0.1*(tfinal-t))
    y = y0
    tlen = length(t)

    # Compute initial step size.
    s1 = F(t, y)
    r = norm_inf_any(s1/max(norm_inf_any(y), threshold)) + realmin()
    h = tdir*0.8*rtol^(1/3)/r

    while t != tfinal  # The main loop.
        hmin = 16*eps()*abs(t)
        if abs(h) > hmax; h = tdir*hmax; end
        if abs(h) < hmin; h = tdir*hmin; end

        # Stretch the step if t is close to tfinal.
        if 1.1*abs(h) >= abs(tfinal - t)
            h = tfinal - t
        end

        # Attempt a step.
        s2 = F(t+h/2.0, y+(h/2.0)*s1)
        s3 = F(t+3.0*h/4.0, y + (3.0*h/4.0)*s2)
        tnew = t + h
        ynew = y + (h*2.0/9.0)*s1 + (h*3.0/9.0)*s2 + (h*4.0/9.0)*s3
        s4 = F(tnew, ynew)

        # Estimate the error.
        e = (-h*5.0/72.0)*s1 + (h*6.0/72.0)*s2 + (h*8.0/72.0)*s3 + (-h*9.0/72.0)*s4
        err = norm_inf_any(e)/max(max(norm_inf_any(y), norm_inf_any(ynew)), threshold) + realmin()
        # Accept the solution if the estimated error is less than the tolerance.
        if err <= rtol
            t = tnew
            y = ynew
            s1 = s4   # Reuse final function value to start new step
        end

        # Compute a new step size.
        h = h*min(5.0, 0.8*(rtol/err)^(1/3))

        # Exit early if step size is too small.
        if abs(h) <= hmin
            println("Step size ", h, " too small at t = ", t)
            t = tfinal
        end
    end # while (t != tfinal)
    return (t, y)
end # ode23_abstract_end

# used in ode
function norm_inf_any{D}(y::Array{Float64,D})
    norm(y,Inf)
end

# used in ode
function norm_inf_any(y::Flow)
    nKap = length(y.kappa)
    nPhi = length(y.phix)
    m=-Inf
    for i=1:nKap
        m = max(norm_inf_any(y.kappa[i]), m)
        m = max(norm_inf_any(y.eta_coeff[i]), m)
    end
    for i=1:nPhi
        m = max(norm_inf_any(y.phix[i]), m)
        m = max(norm_inf_any(y.Dphix[i]), m)
    end
    m
end

# used in ode
function norm_inf_any(y::TrFlow)
    nKap = length(y.kappa)
    nPhi = length(y.phix)
    m=-Inf
    for i=1:nKap
        m = max(norm_inf_any(y.kappa[i]), m)
        m = max(norm_inf_any(y.eta_coeff[i]), m)
        m = max(norm_inf_any(y.dlkappa[i]), m)
        m = max(norm_inf_any(y.dleta_coeff[i]), m)
    end
    for i=1:nPhi
        m = max(norm_inf_any(y.phix[i]), m)
        m = max(norm_inf_any(y.Dphix[i]), m)
        m = max(norm_inf_any(y.dlphix[i]), m)
        m = max(norm_inf_any(y.dlDphix[i]), m)
    end
    m
end

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

# Kernels for the spline flows

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

r(d::Real,sigma) = exp(-0.5 * d * d / (sigma * sigma))

function rp_div_d(d::Real,sigma)
    s2 = sigma * sigma
    -exp(-0.5 * d * d / s2) / s2
end

function rpp(d::Real,sigma)
    rd = r(d,sigma)
    s2 = sigma * sigma
    d * d * rd / (s2 * s2) - rd / s2
end

function R{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
    r(norm(x-y),sigma)
end

function gradR{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
    v=x-y
    n=norm(v)
    v * rp_div_d(n,sigma)
end

outer(u,v) = u*transpose(v)

function g1g2R(x, y, sigma)
    v = x-y
    n = norm(v)
    u = v/n
    eey = eye(length(x))
    rpd = rp_div_d(n,sigma)
    if n == 0
        G = -rpp(n,sigma) * eey
        return G
    else
        G = -rpd * eey
        G += outer(u,-u) * (rpp(n,sigma) - rpd)
        return G
    end
end

g1g1R(x, y,sigma) = -g1g2R(x, y,sigma)

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

# Target densities

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

#------------------------------------
#  Uniform targets
#--------------------------------------------
function targetUnif1d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, limit = 0.4, center = 0.0)
	density = Float64[]
	gradH   = Array{Float64,1}[]
	for k =1:length(kappa)
		kap = kappa[k][1]
		if kap <= center-limit
			dentmp = exp(-((kap  - (center - limit)).^2.0)/(2.0 * targSig * targSig))
			gradHtmp = -(kap - (center - limit)) / (targSig.^2.0)
		elseif center-limit < kap < center+limit
			dentmp = 1.0
			gradHtmp = 0.0
		else
			dentmp = exp(-((kap - (center + limit)).^2.0)/(2.0*targSig.^2.0))
			gradHtmp = -(kap - (center + limit)) / (targSig.^2.0)
		end
		push!(density, dentmp / (limit + limit + targSig*sqrt(2.0*pi))  )
		push!(gradH,[gradHtmp])
	end
	density, gradH
end
function targetUnif2d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, limit = 0.4, center = 0.0)
	N = length(kappa)
	densitx, gradHx = targetUnif1d(Array{Float64,1}[[kappa[i][1]] for i=1:N]; targSig = targSig, limit = limit, center = center)
	density, gradHy = targetUnif1d(Array{Float64,1}[[kappa[i][2]] for i=1:N]; targSig = targSig, limit = limit, center = center)
	densitx .* density, Array{Float64,1}[[gradHx[i][1], gradHy[i][1]] for i=1:N]
end
function targetUnif3d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, limit = 0.4, center = 0.0)
	N = length(kappa)
	densitx, gradHx = targetUnif1d(Array{Float64,1}[[kappa[i][1]] for i=1:N]; targSig = targSig, limit = limit, center = center)
	density, gradHy = targetUnif1d(Array{Float64,1}[[kappa[i][2]] for i=1:N]; targSig = targSig, limit = limit, center = center)
	densitz, gradHz = targetUnif1d(Array{Float64,1}[[kappa[i][2]] for i=1:N]; targSig = targSig, limit = limit, center = center)
	densitx .* density .* densitz, Array{Float64,1}[[gradHx[i][1], gradHy[i][1], gradHz[i][1]] for i=1:N]
end

#------------------------------------
#  Normal Targets
#--------------------------------------------
function targetNormal1d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, center = 0.0)
	density = Float64[]
	density = Float64[]
	gradH   = Array{Float64,1}[]
	for k = 1:length(kappa)
		push!(density, exp(- (kappa[k][1] - center).^2 / (2.0 *targSig^2) ) / (sqrt(2.0*pi) * targSig) )
		push!(gradH,[-(kappa[k][1]-center)/(targSig^2)])
	end
	density, gradH
end
function targetNormal2d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, center = [0.0,0.0])
	N = length(kappa)
	densitx, gradHx = targetNormal1d(Array{Float64,1}[[kappa[i][1]] for i=1:N], targSig = targSig, center = center[1])
	density, gradHy = targetNormal1d(Array{Float64,1}[[kappa[i][2]] for i=1:N], targSig = targSig, center = center[2])
	densitx .* density, Array{Float64,1}[[gradHx[i][1], gradHy[i][1]] for i=1:N]
end
function targetNormal3d(kappa::Array{Array{Float64,1},1}; targSig = 0.1, center = [0.0,0.0,0.0])
	N = length(kappa)
	densitx, gradHx = targetNormal1d(Array{Float64,1}[[kappa[i][1]] for i=1:N], targSig = targSig, center = center[1])
	density, gradHy = targetNormal1d(Array{Float64,1}[[kappa[i][2]] for i=1:N], targSig = targSig, center = center[2])
	densitz, gradHz = targetNormal1d(Array{Float64,1}[[kappa[i][2]] for i=1:N], targSig = targSig, center = center[3])
	densitx .* density .* densitz, Array{Float64,1}[[gradHx[i][1], gradHy[i][1], gradHz[i][1]] for i=1:N]
end

#------------------------------------
# The following is old code and needs to be ported to Julia
#--------------------------------------------
function targetGaussianMixture2d(x, p,mu1,mu2,sigma1,sigma2)
	den=  p*normpdf(x,mu1,sigma1) + (1-p)*normpdf(x,mu2,sigma2)
	dden1=(  p)*normpdf(x,mu1,sigma1).*(-(x-mu1)/sigma1^2)
	dden2=(1-p)*normpdf(x,mu2,sigma2).*(-(x-mu2)/sigma2^2)
	# dden=den.*(-(kap-mu)/TargSig^2)
	gradH=(dden1+dden2)./den
	den, gradH
end

function hfun12d(DefXnew,DefYnew)
	n=2
	sigma=.5
	gradx=(-n/2)*(1/sigma).*(2*DefXnew/sigma)./(1+(DefXnew/sigma).^2)
	grady=(-n/2)*(1/sigma).*(2*DefYnew/sigma)./(1+(DefYnew/sigma).^2)

	expH=(1/sigma^2)*(gamma(n/2)/(sqrt(pi)*gamma((n-1)/2)))^2 * (1+(DefXnew/sigma).^2).^(-n/2).*(1+(DefYnew/sigma).^2).^(-n/2)
	expH, gradx, grady
end
function hfun22d(DefXnew,DefYnew)
	n=2
	sigma=3
	sigma2=1
	gradx=(-n/2)*(1/sigma).*(2*DefXnew/sigma)./(1+(DefXnew/sigma).^2)
	grady=(-1/sigma2^2).*(DefYnew)

	expH=(1/sigma)*(gamma(n/2)/(sqrt(pi)*gamma((n-1)/2))) * (1+(DefXnew/sigma).^2).^(-n/2).*(1/(sqrt(2*pi)*sigma2)).*(exp(-DefYnew.^2/(2*sigma^2)))
	expH, gradx, grady
end

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

#  misc functions

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

function meshgrid(side_x,side_y)
    x = repmat(reshape([side_x],(1,length(side_x))) ,length(side_y),1)
    y = repmat(reshape([side_y],(length(side_y),1)) ,1,length(side_x))
    x,y
end

end # module

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
using Test

include("etl.jl")

@testset "single letter" begin
    input = Dict(1=>['A'])
    output = Dict('a'=>1)
    @test transform(input) == output
end

@testset "single score with multiple letters" begin
    input = Dict(1=>['A', 'E', 'I', 'O', 'U'])
    output = Dict('a'=>1, 'e'=>1, 'i'=>1, 'o'=>1, 'u'=>1)
    @test transform(input) == output
end

@testset "multiple scores with multiple letters" begin
    input = Dict(1=>['A', 'E'], 2=>['D', 'G'])
    output = Dict('g'=>2, 'e'=>1, 'a'=>1, 'd'=>2)
    @test transform(input) == output
end

@testset "multiple scores with differing numbers of letters" begin
    input = Dict(1=>[ 'A', 'E', 'I', 'O', 'U', 'L', 'N', 'R', 'S', 'T' ],
                 2=>[ 'D', 'G' ], 3=>[ 'B', 'C', 'M', 'P' ],
                 4=>[ 'F', 'H', 'V', 'W', 'Y' ], 5=>[ 'K' ],
                 8=>[ 'J', 'X' ], 10=>[ 'Q', 'Z' ])
    output = Dict('a'=>1, 'b'=>3,  'c'=>3, 'd'=>2, 'e'=>1,
                  'f'=>4, 'g'=>2,  'h'=>4, 'i'=>1, 'j'=>8,
                  'k'=>5, 'l'=>1,  'm'=>3, 'n'=>1, 'o'=>1,
                  'p'=>3, 'q'=>10, 'r'=>1, 's'=>1, 't'=>1,
                  'u'=>1, 'v'=>4,  'w'=>4, 'x'=>8, 'y'=>4,
                  'z'=>10)
    @test transform(input) == output
end


text:
#connect: Connect to DSN, returns Connection object, also stores Connection information in global default 'conn' object and global 'Connections' connections array
function connect(dsn::String; usr::String="", pwd::String="")
	global Connections
	global conn
	global env
	dsn_number = 0
	env == C_NULL && (env = ODBCAllocHandle(SQL_HANDLE_ENV,SQL_NULL_HANDLE) )
	dbc = ODBCAllocHandle(SQL_HANDLE_DBC,env)
	ODBCConnect!(dbc,dsn,usr,pwd)
	stmt = ODBCAllocHandle(SQL_HANDLE_STMT,dbc)
	for c in Connections
		if (c.dsn==dsn)
			dsn_number+=1
		end
	end
	conn = Connection(dsn,dsn_number+1,dbc,stmt,null_resultset)
	#print("Connection $(conn.number) to $(conn.dsn) successful.")
	push!(Connections,conn)
	return conn
end
#avancedconnect:
function advancedconnect(conn_string::String="", driver_prompt::Uint16=SQL_DRIVER_NOPROMPT)
	global Connections
	global conn
	global env
	dsn_number = 0
	env == C_NULL && (env = ODBCAllocHandle(SQL_HANDLE_ENV,SQL_NULL_HANDLE) )
	dbc = ODBCAllocHandle(SQL_HANDLE_DBC,env)
	ODBCDriverConnect!(dbc,conn_string,driver_prompt)
	stmt = ODBCAllocHandle(SQL_HANDLE_STMT,dbc)
	for c in Connections
		if (c.dsn==conn_string)
			dsn_number+=1
		end
	end
	conn = Connection(conn_string,dsn_number+1,dbc,stmt,null_resultset)
	#print("Connection $(conn.number) to $(conn.dsn) successful.")
	push!(Connections,conn)
	return conn
end
#query: Sends query string to DBMS, once executed, resultset metadata is returned, space is allocated, and results are returned
function query(querystring::String,conn::Connection=conn; output::Output=DataFrame,delim::Char=',')
	if conn == null_connection
		error("[ODBC]: A valid connection was not specified (and no valid default connection exists)")
	end
	ODBCFreeStmt!(conn.stmt_ptr)
	ODBCQueryExecute(conn.stmt_ptr,querystring)
	holder = DataFrame[]
	while true
		meta = ODBCMetadata(conn.stmt_ptr,querystring)
		if meta.rows == 0
			push!(holder, DataFrame())
		else
			columns, indicator, rowset = ODBCBindCols(conn.stmt_ptr,meta)

			if output == DataFrame
				resultset = ODBCFetchDataFrame(conn.stmt_ptr,meta,columns,rowset,indicator)
			else
				resultset = ODBCDirectToFile(conn.stmt_ptr,meta,columns,rowset,output,delim,length(holder))
			end
			push!(holder,resultset)
		end
		(@FAILED SQLMoreResults(conn.stmt_ptr)) && break
	end
	conn.resultset = length(holder) == 1 ? holder[1] : holder
	ODBCFreeStmt!(conn.stmt_ptr)
	return conn.resultset
end
# sql"..." string literal for convenience; it doesn't do anything different than query right now,
#but we could potentially do some interesting things here
macro sql_str(s)
	query(s)
end
#querymeta: Sends query string to DBMS, once executed, resultset metadata is returned
#it may seem odd to include the other arguments for querymeta, but it's so switching between query and querymeta doesn't require exluding args (convenience)
function querymeta(querystring::String,conn::Connection=conn; output::Output=DataFrame,delim::Char=',')
	if conn == null_connection
		error("[ODBC]: A valid connection was not specified (and no valid default connection exists)")
	end
	ODBCFreeStmt!(conn.stmt_ptr)
	ODBCQueryExecute(conn.stmt_ptr,querystring)
	holder = Metadata[]
	while true
		push!(holder,ODBCMetadata(conn.stmt_ptr,querystring))
		(@FAILED SQLMoreResults(conn.stmt_ptr)) && break
	end
	conn.resultset = length(holder) == 1 ? holder[1] : holder
	ODBCFreeStmt!(conn.stmt_ptr)
	return conn.resultset
end
#disconnect:
function disconnect(connection::Connection=conn)
	global conn
	global Connections
	ODBCFreeStmt!(connection.stmt_ptr)
	SQLDisconnect(connection.dbc_ptr)
	for x = 1:length(Connections)
		if connection.dsn == Connections[x].dsn && connection.number == Connections[x].number
			splice!(Connections,x)
			if is(conn,connection)
				if length(Connections) != 0
					conn = Connections[end]
				else
					conn = null_connection #Reset conn to null default connection
				end
			end
		end
	end
end
#List Installed Drivers
function listdrivers()
	global env
	env == C_NULL && (env = ODBCAllocHandle(SQL_HANDLE_ENV,SQL_NULL_HANDLE) )
	descriptions = String[]
	attributes = String[]
	driver_desc = zeros(Uint8, 256)
	desc_length = Array(Int16, 1)
	driver_attr = zeros(Uint8, 256)
	attr_length = Array(Int16, 1)
	while @SUCCEEDED SQLDrivers(env,driver_desc,desc_length,driver_attr,attr_length)
		push!(descriptions,ODBCClean(driver_desc,1,desc_length[1]))
		push!(attributes,ODBCClean(driver_attr,1,attr_length[1]))
	end
	return [descriptions attributes]
end
#List defined DSNs
function listdsns()
	global env
	env == C_NULL && (env = ODBCAllocHandle(SQL_HANDLE_ENV,SQL_NULL_HANDLE) )
	descriptions = String[]
	attributes = String[]
	dsn_desc = zeros(Uint8, 256)
	desc_length = Array(Int16, 1)
	dsn_attr = zeros(Uint8, 256)
	attr_length = Array(Int16, 1)
	while @SUCCEEDED SQLDataSources(env,dsn_desc,desc_length,dsn_attr,attr_length)
		push!(descriptions,ODBCClean(dsn_desc,1,desc_length[1]))
		push!(attributes,ODBCClean(dsn_attr,1,attr_length[1]))
	end
	return [descriptions attributes]
end


text:
using FixedPolynomials
using Test

import DynamicPolynomials
const Impl = DynamicPolynomials
include("poly_test.jl")
include("config_test.jl")
include("system_test.jl")


text:
struct LCMControlPublisher{M <:MechanismState, Tau, C}
    state::M
    controller::C
    desired_controller_period_ms::Int
    τ::Tau
    robot_info::HumanoidRobotInfo{Float64}
    lcm::LCM
    robot_command_channel::String
    atlas_command_msg::atlas_command_t
    robot_state_msg::robot_state_t
    encodebuffer::FastWriteBuffer
    decodebuffer::FastReadBuffer
end

function LCMControlPublisher(robot_info::HumanoidRobotInfo, controller::C;
        robot_state_channel::String = "EST_ROBOT_STATE",
        robot_command_channel::String = "ROBOT_COMMAND",
        desired_controller_period_ms = 1) where {C}
    mechanism = robot_info.mechanism
    state = MechanismState(mechanism)
    τ = similar(velocity(state))
    lcm = LCM()
    atlas_command_msg = atlas_command_t()
    robot_state_msg = robot_state_t()
    encodebuffer = FastWriteBuffer()
    decodebuffer = FastReadBuffer()
    publisher = LCMControlPublisher(state, controller, desired_controller_period_ms, τ, robot_info,
        lcm, robot_command_channel, atlas_command_msg, robot_state_msg, encodebuffer, decodebuffer)

    atlas_command_msg.num_joints = num_actuators(robot_info)
    resize!(atlas_command_msg)
    for (i, actuator) in enumerate(actuators(robot_info))
        atlas_command_msg.joint_names[i] = actuator.name
    end

    let publisher = publisher
        sub = subscribe(lcm, robot_state_channel, (channel, data) -> handle_robot_state_msg(publisher, data))
        set_queue_capacity(sub, 1)
    end

    publisher
end

function handle(publisher::LCMControlPublisher; async=true)
    lcm = publisher.lcm
    if async
        @async while true
            LCMCore.handle(lcm)
        end
    else
        while true
            LCMCore.lcm_handle(lcm)
        end
    end
end

function handle_robot_state_msg(publisher::LCMControlPublisher, data::AbstractVector{UInt8})
    setdata!(publisher.decodebuffer, data)
    decode!(publisher.robot_state_msg, publisher.decodebuffer)
    set!(publisher.state, publisher.robot_state_msg, publisher.robot_info)
    t = publisher.robot_state_msg.utime / 1e6
    publisher.controller(publisher.τ, t, publisher.state)
    state_des = publisher.state # TODO
    set!(publisher.atlas_command_msg, publisher.τ, t, state_des, publisher.robot_info, publisher.desired_controller_period_ms)
    encode(publisher.encodebuffer, publisher.atlas_command_msg)
    bytes = take!(publisher.encodebuffer)
    publish(publisher.lcm, publisher.robot_command_channel, bytes)
    nothing
end


text:
struct Forward <: InferenceType end
const FORWARD = Forward()

forward_sampling_results() = NonparametricSamplingResults{Forward}(FORWARD,Array{Float64, 1}(), Array{Any, 1}(), Array{Trace, 1}())

@doc raw"""
    function forward_sampling(f; params = (), num_iterations = 1)

Draws samples from the model's joint density. Equivalent to calling
`f` in a loop `num_iterations` times, but results are collected in a
`NonparametricSamplingResults` for easier postprocessing.
"""
function forward_sampling(f; params = (), num_iterations = 1)
    results = forward_sampling_results()
    for n in 1:num_iterations
        t = trace()
        r = f(t, params...)
        push!(results.return_values, r)
        push!(results.traces, t)
        push!(results.log_weights, logprob(t))
    end
    results
end

export Forward, FORWARD, forward_sampling


text:
if !Sys.iswindows()

    libdir = dirname(dlpath(libmpi))

    mpidir = get(ENV, "JULIA_MPI_PATH", joinpath(libdir, ".."))
    incdir = get(ENV, "JULIA_INCLUDE_PATH", joinpath(mpidir, "include"))
    mpicc = get(ENV, "JULIA_MPICC", joinpath(mpidir, "bin", "mpicc"))

    # `cflags` includes libraries and needs to be listed at the end of the
    # compile/link command
    if haskey(ENV, "JULIA_MPI_CFLAGS")
        cflags = Base.shell_split(ENV["JULIA_MPI_CFLAGS"])
    else
        libfilename = basename(libmpi)
        libbasename, = splitext(libfilename)
        libname = replace(libbasename, r"^lib" => s"")
        cflags = `-I$incdir -L$libdir -Wl,-rpath,$libdir -l$libname`
    end

    # We expect the subdirectory `MPIconstants-1.4.0` to contain a copy of
    # <https://github.com/eschnett/MPIconstants>
    srcdir = "MPIconstants-1.4.0"

    run(`$mpicc -o generate_compile_time_mpi_constants $(joinpath(srcdir, "generate_compile_time_mpi_constants.c")) $cflags`)
    run(`./generate_compile_time_mpi_constants`)

    run(`$mpicc -fPIC -shared -o libload_time_mpi_constants.so $(joinpath(srcdir, "load_time_mpi_constants.c")) $cflags`)

else

    cp("compile_time_mpi_constants_msmpi_$(Sys.ARCH).jl", "compile_time_mpi_constants.jl"; force=true)

end


text:
"""
      Rocket3D
"""

struct Rocket3D{I, T} <: Model{I, T}
      n::Int
      m::Int
      d::Int

      mass          # mass
      inertia       # inertia matrix
      inertia_inv   # inertia matrix inverse
      gravity       # gravity
      length        # length (com to thruster)
end

function f(model::Rocket3D, z, u, w)
      # states
      x = view(z,1:3)
      r = view(z,4:6)
      v = view(z,7:9)
      ω = view(z,10:12)

      # force in body frame
      F = view(u, 1:3)

      # torque in body frame
      τ = @SVector [model.length * u[2],
                    -model.length * u[1],
                    0.0]

      SVector{12}([v;
                   0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0*(ω' * r) * r);
                   model.gravity + (1.0 / model.mass) * MRP(r[1], r[2], r[3]) * F;
                   model.inertia_inv * (τ - cross(ω, model.inertia * ω))])
end

n, m, d = 12, 3, 0

mass = 1.0
len = 1.0
inertia = Diagonal(@SVector[1.0 / 12.0 * mass * len^2.0, 1.0 / 12.0 * mass * len^2.0, 1.0e-5])
inertia_inv = Diagonal(@SVector[1.0 / (1.0 / 12.0 * mass * len^2.0), 1.0 / (1.0 / 12.0 * mass * len^2.0), 1.0 / (1.0e-5)])

model = Rocket3D{Midpoint, FixedTime}(n, m, d,
                  mass,
                  inertia,
                  inertia_inv,
                  @SVector[0.0, 0.0, -9.81],
                  len)

function visualize!(vis, p::Rocket3D, q; Δt = 0.1, mesh = true, T_off = length(q))
	default_background!(vis)

	if mesh
		obj_rocket = joinpath(pwd(), "models/starship/Starship.obj")
		mtl_rocket = joinpath(pwd(), "models/starship/Starship.mtl")
		ctm = ModifiedMeshFileObject(obj_rocket, mtl_rocket, scale=1.0)
		setobject!(vis["rocket"]["starship"], ctm)

		settransform!(vis["rocket"]["starship"],
			compose(Translation(0.0, 0.0, -p.length),
				LinearMap(0.25 * RotY(0.0) * RotZ(0.5 * π) * RotX(0.5 * π))))

        body = Cylinder(Point3f0(0.0, 0.0, -1.25),
          Point3f0(0.0, 0.0, 0.5),
          convert(Float32, 0.125))

        setobject!(vis["rocket"]["body"], body,
          MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))
	else
		body = Cylinder(Point3f0(0.0, 0.0, -1.0 * model.length),
			Point3f0(0.0, 0.0, 1.0 * model.length),
			convert(Float32, 0.15))

		setobject!(vis["rocket"], body,
			MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
			anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
	end

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for t = 1:length(q)
	    MeshCat.atframe(anim, t) do
            if t >= T_off
				setvisible!(vis["rocket"]["body"], false)
			else
				setvisible!(vis["rocket"]["body"], true)
			end
	        settransform!(vis["rocket"],
	              compose(Translation(q[t][1:3]),
	                    LinearMap(MRP(q[t][4:6]...) * RotX(0.0))))
	    end
	end

	# settransform!(vis["/Cameras/default"], compose(Translation(0.0, 0.0, 0.0),
	# LinearMap(RotZ(pi/2))))
	MeshCat.setanimation!(vis, anim)
end

# include(joinpath(pwd(), "models/visualize.jl"))
#
# vis = Visualizer()
# render(vis)
# q0 = zeros(model.n)
#
# # simulation test
# u = [0.0; 0.0; 1.0]
# T = 10
# h = 0.1
# x_hist = [q0]
#
# for t = 1:T
# 	push!(x_hist, fd(model, x_hist[end], u, zeros(model.d), h, t))
# end
#
# visualize!(vis, model, x_hist, Δt = 0.1)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
import Base: show, *, -, +, isless, ==, convert, conj, truncate, zero, one, promote_rule, A_mul_B!, dot, transpose

export variables

immutable Symbols{T<:AbstractString}
    names :: Array{T,1}
end

immutable Poly{T<:Number}
    n        :: Int             # number of symbols in monomials
    m        :: Int             # number of monomials
    syms     :: Symbols
    c        :: Array{T,1}      # coefficients in polynomial
    alpha    :: Array{Int, 2}   # powers of symbols for all monomials
    deg      :: Int

    function Poly(syms::Symbols, alpha::Array{Int,2}, c::Array{T,1})
        m, n = size(alpha)
        if length(c) != m error("incompatible dimensions") end

        deg = size(alpha, 1) > 0 ? maximum(sum(alpha,2)) : 0
        new(n, m, syms, c, alpha, deg)
    end

    # we allow a constant polynomial without a proper symbol basis
    Poly(c::T) = new(0, 1, Symbols([""]), [c], Array(T, 1, 0), 0)
end

Poly{T<:Number}(syms::Symbols, alpha::Array{Int,2}, c::Array{T,1}) = Poly{T}(syms, alpha, c)
Poly{T<:Number}(c::T) = Poly{T}(c)

convert{S,T}(::Type{Poly{S}}, p::Poly{T}) = Poly(p.syms, p.alpha, convert(Array{S}, p.c))
convert{S,T}(::Type{Poly{S}}, c::T) = Poly{promote_type(S,T)}(p.syms, p.alpha, convert(promote_type(S,T), c))
convert{S<:Number,T<:Number}(::Type{Polyopt.Poly{S}}, v::T) = convert(Polyopt.Poly{S},convert(S,v))

zero{T<:Number}(::Type{Poly{T}}) = Poly(zero(T))
zero{T<:Number}(::Poly{T}) = Poly(zero(T))
one{T<:Number}(::Type{Poly{T}}) = Poly(one(T))
one{T<:Number}(::Poly{T}) = Poly(one(T))

promote_rule{S<:Number,T<:Number}(::Type{S}, ::Type{Poly{T}}) = Poly{promote_type(S,T)}
promote_rule{S<:Number,T<:Number}(::Type{Poly{S}}, ::Type{Poly{T}}) = Poly{promote_type(S,T)}
promote_rule(::Type{Rational}, ::Type{Int}) = Rational

dot(a::Poly, b::Poly) = a*b
dot(a::Number, b::Poly) = a*b
dot(a::Poly, b::Number) = a*b

function show{T}(io::IO, p::Poly{T})
    if p.m == 0 print(io, zero(T)) end
    for i=1:p.m
        if i==1
            print(io, p.c[i]>=0 ? "" : "-")
        else
            print(io, p.c[i]>=0 ? "+" : "-")
        end

        ci = abs(p.c[i])
        if sum(p.alpha[i,:]) == 0
            print(io, ci)
        else
            if ci != one(ci) print(io, "$(ci)*") end
            for k=1:p.n
                if p.alpha[i,k] > 0
                    print(io, p.syms.names[k])
                    if p.alpha[i,k] > 1
                        print(io, "^$(p.alpha[i,k])")
                    end
                    if k < p.n && sum(p.alpha[i,k+1:end]) > 0
                        print(io, "*")
                    end
                end
            end
        end
    end
end

-{T<:Number}(p::Poly{T}) = Poly{T}(p.syms, p.alpha, -p.c)
+{S<:Number,T<:Number}(p::Poly{S}, a::T) = p + Poly{T}(p.syms, zeros(Int64, 1, p.n), T[a])
+{S<:Number,T<:Number}(a::S, p::Poly{T}) = Poly{S}(p.syms, zeros(Int64, 1, p.n), S[a]) + p
-{S<:Number,T<:Number}(p::Poly{T}, a::S) = p - Poly{S}(p.syms, zeros(Int64, 1, p.n), S[a])
-{S<:Number,T<:Number}(a::S, p::Poly{T}) = Poly{S}(p.syms, zeros(Int64, 1, p.n), S[a]) - p

*{S<:Number,T<:Number}(p::Poly{S}, a::T) = simplify(Poly{promote_type(S,T)}(p.syms, p.alpha, p.c*a))
*{S<:Number,T<:Number}(a::S, p::Poly{T}) = *(p::Poly{T}, a::S)

function add{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    U = Array{promote_type(S, T),1}
    p1, p2 = promote_poly(p1, p2)
    Poly{promote_type(S,T)}(p1.syms, vcat(p1.alpha, p2.alpha), vcat(convert(U,p1.c), convert(U,p2.c)))
end

function sub{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    U = Array{promote_type(S, T),1}
    p1, p2 = promote_poly(p1, p2)
    Poly{promote_type(S,T)}(p1.syms, vcat(p1.alpha, p2.alpha), vcat(convert(U,p1.c), -convert(U,p2.c)))
end

+{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T}) = simplify( add(p1, p2) )
-{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T}) = simplify( sub(p1, p2) )

function *{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    p1, p2 = promote_poly(p1, p2)
    if p2.m == 1
        Poly(p1.syms, p1.alpha .+ p2.alpha, p1.c*p2.c[1])
    elseif p1.m == 1
        Poly(p1.syms, p2.alpha .+ p1.alpha, p2.c*p1.c[1])
    else
        if p1.m < p2.m
            #r = Poly(p1.syms, p2.alpha .+ p1.alpha[1:1,:], p2.c*p1.c[1])
            r = Poly(p1.syms, p2.alpha .+ view(p1.alpha,1:1,:), p2.c*p1.c[1])
            for k=2:p1.m
                #r = add(r, Poly(p1.syms, p2.alpha .+ p1.alpha[k:k,:], p2.c*p1.c[k]))
                r = add(r, Poly(p1.syms, p2.alpha .+ view(p1.alpha,k:k,:), p2.c*p1.c[k]))
            end
        else
            #r = Poly(p1.syms, p1.alpha .+ p2.alpha[1:1,:], p1.c*p2.c[1])
            r = Poly(p1.syms, p1.alpha .+ view(p2.alpha,1:1,:), p1.c*p2.c[1])
            for k=2:p2.m
                #r = add(r, Poly(p1.syms, p1.alpha .+ p2.alpha[k:k,:], p1.c*p2.c[k]))
                r = add(r, Poly(p1.syms, p1.alpha .+ view(p2.alpha,k:k,:), p1.c*p2.c[k]))
            end
        end
        simplify(r)
    end
end

function ^{T<:Number}(p::Poly{T}, a::Int)
    if a<0 error("power must be nonnegative") end
    if a==0 return Poly{T}(p.syms, zeros(Int, 1, p.n), [one(T)]) end
    if a==1 return Poly{T}(p.syms, p.alpha, p.c) end
    if p.m == 1 return Poly{T}(p.syms, p.alpha*a, p.c.^a) end

    powers = Array(Poly{T}, Int(1+ceil(log2(a))))

    @inbounds begin
    powers[1] = p
    k = 1
    while 2^k<=a
        powers[k+1] = powers[k]*powers[k]
        k += 1
    end

    b = 2^(k-1)
    r = powers[k]
    while ( b != a)
        k -= 1
        if (b + 2^(k-1) <= a)
            b += 2^(k-1)
            r = r*powers[k]
        end
    end

    end
    r
end

typealias MatOrVec{T} Union{Array{T,1},Array{T,2}}

*{T<:Number,S<:Number}(a::T, v::MatOrVec{Poly{S}}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(v::MatOrVec{Poly{T}}, a::S) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(a::Poly{T}, v::MatOrVec{Poly{S}}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(v::MatOrVec{Poly{S}}, a::Poly{T}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))

function =={S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    r = p1-p2
    r.m == 0 || r.n == 0
end

=={S<:Number,T<:Number}(p::Poly{S}, c::T) = p == Poly(c)
=={S<:Number,T<:Number}(c::S, p::Poly{T}) = p == Poly(c)

conj{T<:Number}(p::Poly{T}) = Poly{T}(p.syms, p.alpha, conj(p.c))
convert{T<:Number}(::Type{Poly{T}}, a::T) = Poly{T}(convert(T,a))
#isconst{T<:Number}(p::Poly{T}) = (p.m == 0 || p.n == 0)
isconst{T<:Number}(p::Poly{T}) = p.n == 0

function A_mul_B!{S<:Number,T<:Number,U<:Number,V<:Number}(alpha::Poly{S}, A::SparseMatrixCSC{T,Int}, x::Array{Poly{U},1}, beta::Poly{V}, y::Array{Poly{V},1})

    @inbounds begin
    for i=1:length(y)
        if beta == zero(y[i])
            y[i] = zero(y[i])
        else
            y[i] *= beta
        end
    end

    for j=1:size(A,2)
        for k=A.colptr[j]:A.colptr[j+1]-1
            y[A.rowval[k]] += alpha*A.nzval[k]*x[j]
        end
    end
    end
    y
end

transpose{T<:Number}(p::Poly{T}) = p

# promote polynomial to share the same Symbol basis
function promote_poly{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    if (p1.syms == p2.syms) || (isconst(p1) && isconst(p2)) return (p1, p2) end

    if isconst(p1)  && ~isconst(p2)
        return (Poly{S}(p2.syms, zeros(Int, p1.m, p2.n), p1.c), p2)
    elseif ~isconst(p1) && isconst(p2)
        return (p1, Poly{T}(p1.syms, zeros(Int, p2.m, p1.n), p2.c))
    else
        error("incompatible monomial bases")
    end
end

function evalpoly{T<:Number,S<:Number}(p::Poly{T}, x::AbstractArray{S})
    val = zero(S)
    for j=1:p.m
        val += p.c[j]*prod(vec(x) .^ vec(p.alpha[j,:]))
    end
    val
end

# combine identical terms and remove zero terms
function simplify{T<:Number}(p::Poly{T})

    # handle simple cases explicitly to speed function up
    if p.m == 1 && p.c[1] == zero(T)
        return Poly{T}(p.syms, Array(Int, 0, p.n), Array(T, 0))
    end

    if p.m <= 1
        return p
    end

    #y = ordermap(p)
    #perm = sortperm(y)

    c = 1:p.n
    rows = [ view(p.alpha,i,c) for i=1:p.m ]
    perm = sortperm(rows, order=Base.Order.Lexicographic)
    y = p.alpha*rand(p.n)

    first = Array(Int, p.m)
    c = Array(T, p.m)

    l = 1
    @inbounds begin
    first[l] = 1
    c[l] = p.c[perm[1]]

    for k=2:p.m
        if y[perm[k]] == y[perm[first[l]]]
            c[l] += p.c[perm[k]]
        else
            l += 1
            first[l] = k
            c[l] = p.c[perm[k]]
        end
    end

    I = view(c,1:l) .!= 0
    y = Poly{T}(p.syms, p.alpha[perm[view(first,1:l)[I]],:], view(c,1:l)[I])
    end
    y
end

function truncate{T<:Number}(p::Poly{T}, threshold=1e-10)
    lgt = 0
    labeled = falses(p.m)

    for k=1:p.m
        if abs(p.c[k]) < threshold
            labeled[k] = true
        else
            lgt += 1
        end
    end

    # the unlabeled terms form the reduced polynomial
    alpha = Array(Int, lgt, p.n)
    c = Array(T, lgt)
    lgt = 1
    for k=1:p.m
        if ~labeled[k]
            alpha[lgt,:] = p.alpha[k,:]
            c[lgt] = p.c[k]
            lgt += 1
        end
    end
    perm = sortperm([ vec(alpha[i,1:p.n]) for i=1:size(alpha,1)], lt=lexless)
    Poly(p.syms, alpha[perm,:], c[perm])
end

# ordering for polynomials
function isless{T<:Number}(p1::Poly{T}, p2::Poly{T})
    lexless(p1.alpha, p2.alpha)
end

# function ordermap{T<:Number}(p::Poly{T})
#
#     @inbounds begin
#
#     r = zeros(Float64, p.n)
#     r[p.n] = 1
#     t = p.deg
#
#     for l=p.n-1:-1:1
#         r[l] = t+1
#         t += p.deg*r[l]
#     end
#
#     end
#
#     p.alpha*r
# end

# function order{T<:Number}(p::Poly{T})
#     perm = sortperm(ordermap(p))
#     Poly{T}(p.syms, p.alpha[perm,:], p.c[perm])
# end

function order{T<:Number}(p::Poly{T})
    c = 1:size(p.alpha,2)
    rows = [ view(p.alpha,i,c) for i=1:size(p.alpha,1) ]
    perm = sortperm(rows, order=Base.Order.Lexicographic)
    Poly{T}(p.syms, p.alpha[perm,:], p.c[perm])
end

variables(syms::Symbols) = [Poly{Int}(syms, [zeros(Int,1,k-1) 1 zeros(Int,1,length(syms.names)-k)], [1]) for k=1:length(syms.names)]
variables{T<:AbstractString}(syms::Vector{T}) = variables(Symbols(syms))
variables{T<:AbstractString}(sym::T, n::Int) = variables([string(sym,i) for i=1:n])


text:
"""
The Kernels module defines a collection of kernels to avoid namespace collisions:

- `biweight`
- `cosine`
- `epanechnikov`
- `triangular`
- `tricube`
- `triweight`
- `uniform`
- `gaussian`
- `logistic`
"""
module Kernels
inrange(u::Float64) = abs(u) <= 1.0

biweight(u::Float64)        = inrange(u) ? (1.0 - u ^ 2) ^ 2 : 0.0
cosine(u::Float64)          = inrange(u) ? cos(0.5 * π * u) : 0.0
epanechnikov(u::Float64)    = inrange(u) ? 1.0 - u ^ 2 : 0.0
triangular(u::Float64)      = inrange(u) ? 1.0 - abs(u) : 0.0
tricube(u::Float64)         = inrange(u) ? (1.0 - abs(u) ^ 3) ^ 3 : 0.0
triweight(u::Float64)       = inrange(u) ? (1.0 - u ^ 2) ^ 3 : 0.0
uniform(u::Float64)         = inrange(u) ? 0.5 : 0.0

gaussian(u::Float64) = exp(-0.5 * u ^ 2)
logistic(u::Float64) = 1.0 / (exp(u) + 2.0 + exp(-u))
end


text:
module BidirectionalEncoder

using Flux
using Requires
using Requires: @init
using BSON

using ..Basic
using ..Pretrain: isbson, iszip, istfbson, zipname, zipfile, findfile
export Bert, load_bert_pretrain, bert_pretrain_task, masklmloss

include("./bert.jl")
include("./tfckpt2bson.jl")
include("./load_pretrain.jl")
include("./tokenizer.jl")
include("./wordpiece.jl")
include("./pretrain_helper.jl")

end


text:
"""
    vcat!!(s::SparseMatrixCSC, J,V)

*WARNING*: `s` is invalid after this operation!!! Its symbol should be reassigned
to the output of `vcat!!`.
Append new rows to `s` with 1 value `V[r]` at column `J[r]` each.
Because `SparseMatrixCSC` is immutable, return a new, valid
`SparseMatrixCSC` that points to the original's data structures.
For large matrices, this is much faster than `vcat`.
"""
function vcat!!(s::SparseMatrixCSC, J,V)
  length(V)==length(J) || error("[vcat!!] J,V must have the same length")
  k= length(V); k_s= nnz(s)
  m= s.m+k
  resize!(s.nzval, k_s+k)
  resize!(s.rowval, k_s+k)

  # Rearrange values-to-add by ascending column
  col_order= sortperm(J)
  J= J[col_order]; V= V[col_order]
  # Calculate how many steps forward each column start moves
  colptr_cat= copy(s.colptr)
  for c in J
     @inbounds colptr_cat[c+1:end].+= 1
  end

  ## Fill in the new values at the correct places
  # NOTE start from the end, because that's where the empty places are
  for c= s.n:-1:1
    colrange= s.colptr[c] : s.colptr[c+1]-1
    colrange_cat= colptr_cat[c] : colptr_cat[c+1]-1
    # 1: Transport previous values to new places`
    @inbounds s.rowval[colrange_cat[1:length(colrange)]].= s.rowval[colrange]
    @inbounds s.nzval[colrange_cat[1:length(colrange)]].= s.nzval[colrange]

    # 2: add new values
    if length(colrange_cat) > length(colrange)
      # OPTIMIZE J.==c is a lot of comparisons!
       @inbounds s.rowval[colrange_cat[length(colrange)+1:end]].=
          col_order[J .== c] .+ s.m
       @inbounds s.nzval[colrange_cat[length(colrange)+1:end]].= V[J .== c]
    end
  end

  # Construct the new SparseMatrixCSC
  SparseMatrixCSC(m,s.n,colptr_cat,s.rowval,s.nzval)
end
"""
    hcat!!(s::SparseMatrixCSC, I,V)

*WARNING*: `s` is invalid after this operation!!! Its symbol should be reassigned
to the output of `hcat!!`.
Append new columns to `s` with 1 value `V[c]` at row `I[c]` each.
Because `SparseMatrixCSC` is immutable, return a new, valid
`SparseMatrixCSC` that points to the original's data structures.
For large matrices, this is *much* faster than `hcat`.
"""
function hcat!!(s::SparseMatrixCSC, I,V)
  length(V)==length(I) || error("[hcat!!] I,V must have the same length")
  k= length(V); k_s= nnz(s)
  n= s.n+k
  resize!(s.nzval, k_s+k)
  resize!(s.rowval, k_s+k)
  resize!(s.colptr, n+1)
  s.nzval[k_s+1:end].= V
  s.rowval[k_s+1:end].= I
  s.colptr[s.n+2:end].= s.colptr[s.n+1] .+ (1:k)
  # Construct the new SparseMatrixCSC
  SparseMatrixCSC(s.m,n,s.colptr,s.rowval,s.nzval)
end
# Change just the number of columns, without adding any new values
function hcat!!(s::SparseMatrixCSC,k)
  n= s.n+k
  resize!(s.colptr,n+1)
  s.colptr[s.n+2:end].= s.colptr[s.n+1]
  SparseMatrixCSC(s.m,n,s.colptr,s.rowval,s.nzval)
end


text:
# conversion examples

# lossy

#=
julia> convert(Float64, 1//3)
0.3333333333333333

julia> convert(Rational, convert(Float64, 1//3))
6004799503160661//18014398509481984
=#

#= REPL

------------------------------------------------------------------------------
Use case #1

julia> x = rand(3)
3-element Array{Float64,1}:
 0.7472603457673705
 0.08245417518187859
 0.06299555866248618

julia> x[1] = 1
1

julia> x
3-element Array{Float64,1}:
 1.0
 0.08245417518187859
 0.06299555866248618

------------------------------------------------------------------------------
Use case #2

julia> mutable struct Foo
           x::Float64
       end

julia> foo = Foo(1.0)
Foo(1.0)

julia> foo.x = 2
2

julia> foo
Foo(2.0)

------------------------------------------------------------------------------
Use case #3

julia> struct Foo
           x::Float64
           Foo(v) = new(v)
       end

julia> Foo(1)
Foo(1.0)

------------------------------------------------------------------------------
Use case #4

julia> function foo()
           local x::Float64
           x = 1
           println(x, " has type of ", typeof(x))
       end

julia> foo()
1.0 has type of Float64

------------------------------------------------------------------------------
Use case #5

julia> function foo()::Float64
           return 1
       end

julia> foo()
1.0

------------------------------------------------------------------------------
Use case #6

julia> ccall((:exp, "libc"), Float64, (Float64,), 2)
7.38905609893065

------------------------------------------------------------------------------

julia> subtypetree(AbstractFloat)
AbstractFloat
    BigFloat
    Float16
    Float32
    Float64

julia> twice(x::AbstractFloat) = 2x;

julia> twice(1.0)
2.0

julia> BigFloat("1.5e1234")
1.500000000000000000000000000000000000000000000000000000000000000000000000000007e+1234

------------------------------------------------------------------------------

julia> twice(2)
ERROR: MethodError: no method matching twice(::Int64)

------------------------------------------------------------------------------

julia> twice(x) = 2x
twice (generic function with 1 method)

julia> twice(1.0)
2.0

julia> twice(1)
2

------------------------------------------------------------------------------

julia> twice(x::Number) = 2x
twice (generic function with 1 method)

julia> twice(2.0)
4.0

julia> twice(2)
4

julia> twice(2//3)
4//3
=#


text:
struct CompiledOctTreeNode{T<:AbstractPoint3D} <: SpatialTreeNode
    point::T
    l::Float64
    next::Int64
end

CompiledOctTreeNode(point::T, l::Float64, next::Int64) where T <: AbstractPoint3D = CompiledOctTreeNode{T}(point, l, next)

CompiledOctTreeNode(n::OctTreeNode{T}) where T <: AbstractPoint3D = CompiledOctTreeNode{T}(n.point, isleaf(n) ? -1.0 : 2.0*n.r, -1)

withnext(cn::CompiledOctTreeNode{T}, next::Int64) where T<:AbstractPoint3D = CompiledOctTreeNode{T}(cn.point, cn.l, next)

mutable struct CompiledOctTree{T<:AbstractPoint3D} <: SpatialTree
    nodes::SharedArray{CompiledOctTreeNode{T}, 1}
    number_of_nodes_used::Int64
    faststack::Array{Int64, 1}
end

CompiledOctTree(n::Int64, ::Type{T}) where T <: AbstractPoint3D =
    CompiledOctTree{T}(SharedArray{CompiledOctTreeNode{T}}(2*n), 0, Array{Int64}(undef,100000))

function stop_cond(q::OctTreeNode{T}, ct::CompiledOctTree{T}) where T <: AbstractPoint3D
    isemptyleaf(q) && return true # nothing to do
    ct.number_of_nodes_used += 1
    ct.nodes[ct.number_of_nodes_used] = CompiledOctTreeNode(q)
    q.id = ct.number_of_nodes_used
    return false
end

function compile!(ct::CompiledOctTree{T}, t::OctTree{T}) where T<:AbstractPoint3D
    ct.number_of_nodes_used = 0
    map(t, ct)
    childs = Array{OctTreeNode{T}}(undef, 8)

    # fix neighbours
    @inbounds for i in 1:t.number_of_nodes_used
        q=t.nodes[i]
        !q.is_divided && continue
        childs[1] = q.lxlylz
        childs[2] = q.lxlyhz
        childs[3] = q.lxhylz
        childs[4] = q.lxhyhz
        childs[5] = q.hxlylz
        childs[6] = q.hxlyhz
        childs[7] = q.hxhylz
        childs[8] = q.hxhyhz
        for a in 1:7
            qa = childs[a]
            qa.id <= 0 && continue
            # searching qa neighbor...
            qb = qa
            for b in (a+1):8
                if childs[b].id > 0
                    qb = childs[b]
                    break
                end
            end
            (qb === qa) && continue # no neighbor found lets leave next=-1
            # qa and qb are neighbors -- fix that in the compiled tree
            ct.nodes[qb.id] = withnext(ct.nodes[qb.id], qa.id)
        end
    end
    ct
end

function map(t::CompiledOctTree{T}, cond_data) where T<:AbstractPoint3D
    curr_stack_ix = 1
    t.faststack[1] = 1
    @inbounds while curr_stack_ix > 0
        node_ix = t.faststack[curr_stack_ix]
        curr_stack_ix -= 1
        if node_ix>0 && node_ix<=t.number_of_nodes_used && !stop_cond(t.nodes[node_ix], cond_data)
            curr_stack_ix += 1
            next_ix = node_ix+1
            t.faststack[curr_stack_ix] = next_ix
            while next_ix > 0
                curr_stack_ix += 1
                next_ix = t.nodes[next_ix].next
                t.faststack[curr_stack_ix] = next_ix
            end
        end
    end
end


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Script for executing Gurobi on selected datasets on the DRO problem with Wasserstein metric
# based ambiguity sets. Command line usage:
# julia run_solver.jl <dataset> <method>
# method: -1=automatic, 0=primal simplex, 1=dual simplex, 2=barrier, 3=concurrent, 4=deterministic concurrent, 5=deterministic concurrent simplex.

using Gurobi
using JuMP
using LinearAlgebra
using SparseArrays

BLAS.set_num_threads(1)

include("../src/problems/standardLP.jl")
include("../src/algorithms/utils/exitcriterion.jl")
include("../src/algorithms/utils/results.jl")
include("../src/algorithms/utils/helper.jl")
include("../src/algorithms/clvr_lazy.jl")
include("../src/problems/dro/utils/libsvm_parser.jl")
include("../src/problems/dro/wasserstein.jl")

DATASET_INFO = Dict([
    ("a1a", (123, 1605)),
    ("a9a", (123, 32561)),
    ("gisette", (5000, 6000)),
    ("news20", (1355191, 19996)),
    ("rcv1", (47236, 20242)),
])

# ARGS
dataset = ARGS[1]
gurobimethod = parse(Int, ARGS[2])

# Dataset parameters
if !haskey(DATASET_INFO, dataset)
    throw(ArgumentError("Invalid dataset name supplied."))
end
dim_dataset, num_dataset = DATASET_INFO[dataset]
filepath = "./data/$(dataset).txt"

# Problem instance parameters
κ = 0.1
ρ = 10.

# Problem instance instantiation
yX_T = read_libsvm_into_yXT_sparse(filepath, dim_dataset, num_dataset)
A_T, b, c = droreformuation_wmetric_hinge_standardformnormalized(yX_T, κ, ρ)
problem = StandardLinearProgram(A_T, b, c)

# Gurobi
println("Setting up LP")
starttime = time()
vector_model = Model(Gurobi.Optimizer)
set_optimizer_attribute(vector_model, "Threads", 1)
set_optimizer_attribute(vector_model, "Method", gurobimethod)
@variable(vector_model, x[1:size(A_T)[1]] >= 0)
@constraint(vector_model, equality, A_T' * x .== b)
@objective(vector_model, Min, c' * x)
endtime = time()
println("=====> Setting up: $(endtime - starttime)")

println("Solving LP")
starttime = time()
optimize!(vector_model)
endtime = time()
println("=====> Solve time: $(endtime - starttime)")

println(objective_value(vector_model))
println(termination_status(vector_model))
println(primal_status(vector_model))

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
const TA{F, N} = Tracker.TrackedArray{F, N, Array{F, N}}
const TR{F} = Tracker.TrackedReal{F}

mutable struct State{A1, A2, A3}
  delta0::A1 # L0
  delta1::A1 # L1
  sig2::A1 # I
  W::A2 # I x K
  eta0::A3 # I x J x K
  eta1::A3 # I x J x K
  v::A1 # K
  H::A2 # J x K
  alpha::A1 # 1 (F won't work, TrackedReals don't work as expected)
  eps::A1 # I
  y_m::TA{Float64, 2} # I x J
  y_log_s::TA{Float64, 2} # I x J

  State(A::Type) = new{A{1}, A{2}, A{3}}()
end

const StateMP = State{ADVI.MPA{Float64, 1}, ADVI.MPA{Float64, 2}, ADVI.MPA{Float64, 3}}

function State(c::Constants)
  s = State(ADVI.MPA{Float64})
  s.delta0 = ADVI.ModelParam(c.L[0], "positive", m=ones(c.L[0]), s=ones(c.L[0]))
  s.delta1 = ADVI.ModelParam(c.L[1], "positive", m=ones(c.L[1]), s=ones(c.L[1]))
  s.eps = ADVI.ModelParam(c.I, "unit", m=fill(log(1e-6), c.I), s=fill(1e-7, c.I))
  s.sig2 = ADVI.ModelParam(c.I, "positive", m=fill(log(.1), c.I), s=fill(.1, c.I))
  s.W = ADVI.ModelParam((c.I, c.K - 1), "simplex")
  s.eta0 = ADVI.ModelParam((c.I, c.J, c.L[0] - 1), "simplex")
  s.eta1 = ADVI.ModelParam((c.I, c.J, c.L[1] - 1), "simplex")
  s.v = ADVI.ModelParam(c.K, "unit")
  s.H = ADVI.ModelParam((c.J, c.K), "unit")
  s.alpha = VB.ADVI.ModelParam("positive")
  s.y_m = param(fill(-3.0, c.I, c.J))
  s.y_log_s = param(fill(log(.1), c.I, c.J))

  return s
end

function rsample(s::StateMP; AT::Type=TA{Float64})
  real = State(AT)
  tran = State(AT)

  for key in fieldnames(State)
    if !(key in (:y_m, :y_log_s))
      f = getfield(s, key)
      if typeof(f) <: Array
        println("In State.jl: This message should not be printing!")
        # TODO: optimize
        rs = []
        ts = []
        for each_f in f
          println(f)
          r = ADVI.rsample(f)
          t = ADVI.transform(f, r)
          append!(rs, r)
          append!(ts, t)
        end
        setfield!(real, key, rs)
        setfield!(tran, key, ts)
      else
        r = ADVI.rsample(f)
        t = ADVI.transform(f, r)
        setfield!(real, key, r)
        setfield!(tran, key, t)
      end
    end
  end

  return real, tran
end

function rsample(s::StateMP, y::Vector{M}, c::Constants; AT::Type=TA{Float64}) where M
  # Sample real and transformed parameters
  real, tran = rsample(s, AT=AT)

  # Draw y and compute log q(y|m)
  log_qy = 0.0
  yout = [begin
    vae = VAE(s.y_m[i:i, :], s.y_log_s[i:i, :])
    yi, log_qyi = vae(y[i], c.N[i])
    @assert size(yi) == size(y[i])
    log_qy += log_qyi
    yi
  end for i in 1:c.I]

  return real, tran, yout, log_qy
end