text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# Example of 2d nonlinear elasticity with non-homogeneous boundary contitions

# Geometry of nodes
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])

block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true

element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
add_element!(block, element)

# boundary conditions
bc = Problem(Dirichlet, "bc", 2, "displacement")
bel1 = Element(Seg2, (1, 2))
bel2 = Element(Seg2, (3, 4))
bel3 = Element(Seg2, (4, 1))
update!((bel1, bel2, bel3), "geometry", X)
update!(bel1, "displacement 2", 0.0)
update!(bel2, "displacement 2", 0.5)
update!(bel3, "displacement 1", 0.0)
add_elements!(bc, bel1, bel2, bel3)

analysis = Analysis(Nonlinear)
add_problems!(analysis, block, bc)
run!(analysis)

# results are verified using Code Aster
eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
u3_expected =  [-2.36237356855269E-01, 5.00000000000000E-01]

u3 = block("displacement", 0.0)[3]
@test isapprox(u3, u3_expected, atol=1.0e-5)


text:
abstract type AbstractLeft end
abstract type AbstractRight end
abstract type AbstractUnit <: AbstractRight; end

struct Px <:   AbstractUnit;   val::Real; end
struct Em <:   AbstractUnit;   val::Real; end
struct Cm <:   AbstractUnit;   val::Real; end
struct In <:   AbstractUnit;   val::Real; end
struct Deg <:  AbstractUnit;   val::Real; end
struct Mm <:   AbstractUnit;   val::Real; end
struct Pt <:   AbstractUnit;   val::Real; end
struct Pc <:   AbstractUnit;   val::Real; end
struct Ex <:   AbstractUnit;   val::Real; end
struct Ch <:   AbstractUnit;   val::Real; end
struct Rem <:  AbstractUnit;   val::Real; end
struct Vw <:   AbstractUnit;   val::Real; end
struct Vh <:   AbstractUnit;   val::Real; end
struct Vmin <: AbstractUnit;   val::Real; end
struct Vmax <: AbstractUnit;   val::Real; end
struct Pcnt <: AbstractUnit;   val::Real; end

struct Hex <:   AbstractRight; val::Array{Int64}; end
struct Int_  <: AbstractRight; val::Int64;       end # These do seem redundant but actually help.
struct Float <: AbstractRight; val::Float64;     end
# ==============================================================================
struct Class <: AbstractLeft;     val::String; end
struct Tag <: AbstractLeft;       val::String; end
struct Ident <: AbstractLeft;       val::String; end
struct Attr <: AbstractLeft;      val::String; end

#"digit", "float", "value", "reference", "hex", "string", "unit"
struct Digit <: AbstractRight;     val::String; end
struct Value <: AbstractRight;     val::String; end
struct Reference <: AbstractRight; val::String; end
struct Str <: AbstractRight;       val::String; end
struct Unit <: AbstractRight;       val::String; end

struct Alpha <: AbstractRight;     val::String; end

struct Open;      val::String; end
struct Close;     val::String; end
struct End;       val::String; end

struct Hash;      val::String; end
struct Dot;       val::String; end
struct Dotdot;    val::String; end
struct Dollar;    val::String; end
struct At;        val::String; end
struct Space;     val::String; end

# ==============================================================================
units = [ "px", "em", "cm", "in", "deg", "mm", "pt", "pc",
               "ex", "ch", "rem", "vw", "vh", "vmin", "vmax"] # % & #

ext: jl


text:
using CUSOLVER
using Base.Test

@testset "dense" begin
    include("dense.jl")
end
@testset "sparse" begin
    include("sparse.jl")
end


text:
using PlotsLH
using Plots, Test

include("test_setup.jl")

@testset "All" begin
    include("xy_test.jl")
    include("xyz_test.jl")
    include("regression_test.jl")
end

# ----------------


text:
export TrainingSummary

type TrainingSummary <: Coffee
  statistic_names :: Vector{Any}

  #Default Constructor
  TrainingSummary(statistic_names...) = begin
      if !isempty(statistic_names)
          new(collect(statistic_names))
      else
          new([:iter, :obj_val])
      end
  end
end

function enjoy(lounge::CoffeeLounge, coffee::TrainingSummary, ::Net, state::SolverState)
  summaries = AbstractString[]

  for statistic_name in coffee.statistic_names
    statval = get_statistic(state, statistic_name)
    push!(summaries, format_statistic(state, statistic_name, statval))
    if statistic_name != :iter
        update_statistics(lounge, string(statistic_name), statval)
    end
  end

  @info(" TRAIN ", join(summaries, ' '))

end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

const _M = Grids._M

"""
    knl_apply_filter!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{horizontal},
                      ::Val{vertical}, Q, ::Val{states}, filtermatrix,
                      elems) where {dim, N, nstate, states, horizontal, vertical}

Computational kernel: Applies the `filtermatrix` to the `states` of `Q`.

The arguments `horizontal` and `vertical` are used to control if the filter is
applied in the horizontal and vertical reference directions, respectively.
"""
function knl_apply_filter!(::Val{dim}, ::Val{N}, ::Val{nstate},
                           ::Val{horizontal}, ::Val{vertical}, Q,
                           ::Val{states}, filtermatrix,
                           elems) where {dim, N, nstate, horizontal, vertical,
                                         states}
  FT = eltype(Q)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  filterinξ1 = horizontal
  filterinξ2 = dim == 2 ? vertical : horizontal
  filterinξ3 = dim == 2 ? false : vertical

  # Return if we are not filtering in any direction
  if !(filterinξ1 || filterinξ2 || filterinξ3)
    return
  end

  nfilterstates = length(states)

  s_filter = @shmem FT (Nq, Nq)
  s_Q = @shmem FT (Nq, Nq, Nqk, nfilterstates)
  l_Qfiltered = @scratch FT (nfilterstates, Nq, Nq, Nqk) 3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_filter[i, j] = filtermatrix[i, j]
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          @unroll for fs = 1:nfilterstates
            l_Qfiltered[fs, i, j, k] = zero(FT)
          end

          ijk = i + Nq * ((j-1) + Nq * (k-1))

          @unroll for fs = 1:nfilterstates
            s_Q[i, j, k, fs] = Q[ijk, states[fs], e]
          end
        end
      end
    end

    if filterinξ1
      @synchronize
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            @unroll for n = 1:Nq
              @unroll for fs = 1:nfilterstates
                l_Qfiltered[fs, i, j, k] += s_filter[i, n] * s_Q[n, j, k, fs]
              end
            end
          end
        end
      end

      if filterinξ2 || filterinξ3
        @loop for k in (1:Nqk; threadIdx().z)
          @loop for j in (1:Nq; threadIdx().y)
            @loop for i in (1:Nq; threadIdx().x)
              @unroll for fs = 1:nfilterstates
                s_Q[i, j, k, fs] = l_Qfiltered[fs, i, j, k]
                l_Qfiltered[fs, i, j, k] = zero(FT)
              end
            end
          end
        end
      end
    end

    if filterinξ2
      @synchronize
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            @unroll for n = 1:Nq
              @unroll for fs = 1:nfilterstates
                l_Qfiltered[fs, i, j, k] += s_filter[j, n] * s_Q[i, n, k, fs]
              end
            end
          end
        end
      end

      if filterinξ3
        @loop for k in (1:Nqk; threadIdx().z)
          @loop for j in (1:Nq; threadIdx().y)
            @loop for i in (1:Nq; threadIdx().x)
              @unroll for fs = 1:nfilterstates
                s_Q[i, j, k, fs] = l_Qfiltered[fs, i, j, k]
                (l_Qfiltered[fs, i, j, k] = zero(FT))
              end
            end
          end
        end
      end
    end

    if filterinξ3
      @synchronize
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            @unroll for n = 1:Nq
              @unroll for fs = 1:nfilterstates
                l_Qfiltered[fs, i, j, k] += s_filter[k, n] * s_Q[i, j, n, fs]
              end
            end
          end
        end
      end
    end

    # Store result
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for fs = 1:nfilterstates
            Q[ijk, states[fs], e] = l_Qfiltered[fs, i, j, k]
          end
        end
      end
    end

    @synchronize
  end
  nothing
end

function knl_apply_TMAR_filter!(::Val{nreduce}, ::Val{dim}, ::Val{N}, Q,
                                ::Val{filterstates}, vgeo,
                                elems) where {nreduce, dim, N, filterstates}
  FT = eltype(Q)

  Nq = N + 1
  Nqj = dim == 2 ? 1 : Nq

  nfilterstates = length(filterstates)

  # note that k is the second not 4th index (since this is scratch memory and
  # k needs to be persistent across threads)
  l_Q = @scratch FT (nfilterstates, Nq, Nq, Nqj) 2

  l_MJ = @scratch FT (Nq, Nq, Nqj) 2

  s_MJQ = @shmem FT (Nq * Nqj, nfilterstates)
  s_MJQclipped = @shmem FT (Nq * Nqj, nfilterstates)

  nelemperblock = 1

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        # loop up the pencil and load Q and MJ
        @unroll for k in 1:Nq
          ijk = i + Nq * ((j-1) + Nqj * (k-1))

          @unroll for sf = 1:nfilterstates
            s = filterstates[sf]
            l_Q[sf, k, i, j] = Q[ijk, s, e]
          end

          l_MJ[k, i, j] = vgeo[ijk, _M, e]
        end
      end
    end

    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for sf = 1:nfilterstates
          MJQ, MJQclipped = zero(FT), zero(FT)

          @unroll for k in 1:Nq
            MJ = l_MJ[k, i, j]
            Qs = l_Q[sf, k, i, j]
            Qsclipped = Qs ≥ 0 ? Qs : zero(Qs)

            MJQ += MJ * Qs
            MJQclipped += MJ * Qsclipped
          end

          ij = i + Nq * (j-1)

          s_MJQ[ij, sf] = MJQ
          s_MJQclipped[ij, sf] = MJQclipped
        end
      end
    end

    @synchronize

    @unroll for n = 11:-1:1
      if nreduce ≥ 2^n
        @loop for j in (1:Nqj; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            ij = i + Nq * (j-1)
            ijshift = ij + 2^(n-1)
            if ij ≤ 2^(n-1) && ijshift ≤ Nq * Nqj
              @unroll for sf = 1:nfilterstates
                s_MJQ[ij, sf] += s_MJQ[ijshift, sf]
                s_MJQclipped[ij, sf] += s_MJQclipped[ijshift, sf]
              end
            end
          end
        end
        @synchronize
      end
    end

    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for sf = 1:nfilterstates
          qs_average = s_MJQ[1, sf]
          qs_clipped_average = s_MJQclipped[1, sf]

          r = qs_average > 0 ? qs_average / qs_clipped_average : zero(FT)

          s = filterstates[sf]
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))

            Qs = l_Q[sf, k, i, j]
            Q[ijk, s, e] = Qs ≥ 0 ? r*Qs : zero(Qs)
          end
        end
      end
    end
  end

  nothing
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/05.%20Conditionals.ipynb

#=
with short-circuit evaluation
We've already seen expressions with the syntax

a & b
for two expressions or values a and b. Julia will evaluate this expression eagerly so that

false & (println("hi"); true)
prints "hi" to stdout before returning false.

On the other hand, when we replace & with &&, as in

a && b
we get short-circuit evaluation. b is only evaluated if a is true, which can help us out if evaluating b is expensive. For example,

false && (println("hi"); true)
returns false without printing "hi".

This means we can use

a && b
to conditionally evaluate b if a is true!
=#

false & (println("hi &"); true)
false && (println("hi &&"); true)

#Similarly, check out the difference between
true | (println("hi |"); true)
true || (println("hi ||"); true)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using BernsteinExpansions, Test
using DynamicPolynomials # optional polynomial interface
using StaticArrays

using BernsteinExpansions: _univariate! #, _multivariate!

# initialize some variables
@polyvar x y

@testset "Univariate monomial, floating point" begin
    # exported function
    @test univariate(x^3, 3, 1..2) == [1.0, 2.0, 4.0, 8.0]
    @test univariate(x, 5, 0..1) == [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]

    # internals
    coeffs = zeros(MVector{4, Float64})
    @test _univariate!(coeffs, 3, 3, 1.0, 2.0) == @MVector([1.0, 2.0, 4.0, 8.0])
    coeffs = zeros(MVector{6, Float64})
    @test _univariate!(coeffs, 1, 5, 0.0, 1.0) == @MVector([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
end

@testset "Univariate monomial, rational" begin
    # exported function
    @test univariate(x^3, 3, Interval(1//1, 2//1)) == [1//1, 2//1, 4//1, 8//1]
    @test univariate(x, 5, Interval(0//1, 1//1)) == [0//1, 1//5, 2//5, 3//5, 4//5, 1//1]

    # internals
    coeffs = zeros(MVector{4, Rational{Int}})
    @test _univariate!(coeffs, 3, 3, 1//1, 2//1) == [1//1, 2//1, 4//1, 8//1]
    coeffs = zeros(MVector{6, Rational{Int}})
    @test _univariate!(coeffs, 1, 5, 0//1, 1//1) == [0//1, 1//5, 2//5, 3//5, 4//5, 1//1]
end

# TODO: add tests for multivariate (ref Smith's thesis)
#

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
"""
```
imgb = binarize(Polysegment(), img)
binarize!(Polysegment(), img)
```

Uses the *polynomial segmentation* technique to group the image pixels
into two categories (foreground and background).

# Details

The approach involves constructing a univariate second-degree polynomial
such that the two roots of the polynomial represent the graylevels
of two cluster centers (i.e the foreground and background). Pixels are then
assigned to the foreground or background depending on which cluster center
is closest.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractGray` image which is to be binarized into background (`Gray(0)`)
and foreground (`Gray(1)`) pixels.

## Reference

[1] R. E. Vidal, "Generalized Principal Component Analysis (GPCA): An Algebraic Geometric Approach to Subspace Clustering and Motion Segmentation." Order No. 3121739, University of California, Berkeley, Ann Arbor, 2003.
"""
function binarize(algorithm::Polysegment,  img::AbstractArray{T,2}) where T <: AbstractGray
  binarize!(algorithm, copy(img))
end

"""
```
t = binarize!(Polysegment(), img)
```
Same as [`binarize`](@ref binarize(::Polysegment, ::AbstractArray{<:AbstractGray,2})) except that it modifies the image that was passed as an argument.
"""
function binarize!(algorithm::Polysegment,  img::AbstractArray{T,2}) where T <: AbstractGray
  # Construct data matrix for second-degree polynomial (Equation 2.3) in [1].
  x = vec(img)
  Lₙ = hcat(ones(length(x)), x, x.^2)
  F = svd(Lₙ)
  c = F.Vt[end,:]
  p = Poly(vec(c))
  μ₁, μ₂ = roots(p)
  # Binarize the image.
  map!(img,img) do val
    (val-μ₁)^2 < (val-μ₂)^2 ? Gray(0) : Gray(1)
  end
end


text:
"""
    ESw2(P::Array)

Compute the effect size `w`  for a two-way probability table
corresponding to the alternative hypothesis in the chi-squared test
of association in two-way contingency tables
"""
function ESw2(P::Array)
    pcol = sum(P,1)
    prow = sum(P,2)
    P0 = prow * pcol
    return sqrt(sum((P-P0).^2./P0))
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
# modified from Base (removed initial conversion)
function power_by_squaring(x, p::Integer)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        isone(x) && return copy(x)
        isone(-x) && return iseven(p) ? one(x) : copy(x)
        throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

max_stars_repo_name: dpsanders/LazyTaylorSeries

max_issues_repo_name: dpsanders/LazyTaylorSeries

max_forks_repo_name: dpsanders/LazyTaylorSeries


text: Dict{String,Any}("newline" => Dict{String,Any}("value" => "crlf","type" => "string"),"os" => Dict{String,Any}("value" => "DOS","type" => "string"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
export Filter_freq_band, createCirclesMask

function Filter_freq_band(Q,m,t, parm)
    n      = size(Q,1);
    Q      = [Q;zeros(Float32,size(Q))];
    d      = t[2]-t[1];
    dt     = collect(t[1]:d:(t[1]+2*n*d));
    Q_f    = zeros(Float32,size(Q));

    for i=1:size(Q,2);
	sf            = rfft(Q[:,i]);
        index         = length(sf);
        if m[1]*t[end]>0
          test        = collect(1:index);
          mask        = convert(Array{Float32,1},zeros(length(test)));
          #mask        = 1 ./ (1+exp(-parm*((test-round(m[1])))))+1 ./ (1+exp(+parm*((test-round(m[2])))));
	for j=1:length(test)
		mask[j] = 1/(1+exp(-parm*((test[j]-round(m[1])))))+1/(1+exp(+parm*((test[j]-round(m[2])))));
        end
        else
          test        = collect(1:index);
          #mask        = 1 ./ (1+exp(-parm*((test-round(m[2])))));
          mask        = convert(Array{Float32,1},zeros(length(test)));
	for j=1:length(test)
	  mask[j]     = 1/(1+exp(-parm*((test[j]-round(m[2])))));
        end
        end
        #mask                          = mask -1;
	for j=1:length(mask)
          mask[j] = mask[j]-1;
        end
        sf	                      = mask .* sf;
	q_f		              = irfft(sf,size(Q,1));
	#q_f[convert(Int64,floor(length(q_f)/2)):end] = 0;
	a = convert(Int64,floor(length(q_f)/2));
	for j = a:length(q_f)
	  q_f[j] = 0;
	end
        Q_f[:,i]                      = q_f;
   end
   Q_f = Q_f[1:n,:];
   return Q_f
end

function createCirclesMask(image,centers,radii)
# Brett Shoelson, PhD
# 9/22/2014
# Comments, suggestions welcome: brett.shoelson@mathworks.com
# Copyright 2014 The MathWorks, Inc.
#IMAGE specified
xDim,yDim   = size(image);
xc          = centers[:,1];
yc          = centers[:,2];
x           = 1:xDim;
y           = 1:yDim;
#yy          = repmat(x[:], 1, length(y));
#xx          = repmat(y[:]', length(x), 1);
yy          = repeat(x[:], 1, length(y));
xx          = repeat(y[:]', length(x), 1);
mask        = falses(xDim,yDim);
mask1       = falses(xDim,yDim);
#for ii = 1:length(radii)
#	mask = mask | convert(Array{Int,2},floor(hypot.(xx - xc[ii], yy - yc[ii]))).<=radii[ii];
#end
#added for julia 1.1
for ii = 1:length(radii)
        testx = convert(Array{Float64,2},xx);
        testxc = convert(Array{Float64,2},xc[ii]*ones(size(testx)));
        testy = convert(Array{Float64,2},yy);
        testyc = convert(Array{Float64,2},yc[ii]*ones(size(testy)));
        test11=hypot.(testx - testxc, testy - testyc);
        test12 = zeros(size(test11));
        for i=1:size(test11,1)
       		for j=1:size(test11,2)
       		test12[i,j] = floor(test11[i,j])
       		end
        end
        test13=convert(Array{Int,2},test12).<=radii[ii]
        for i=1:size(mask1,1)
       		for j=1:size(mask1,2)
       		mask[i,j] = mask1[i,j] | test13[i,j];
       		end
        end
end
return mask
end


text:
# custom callbacks
SetCustomMouseButtonCallback(x::Function) = global g_CustomCallbackMousebutton = x
SetCustomScrollCallback(x::Function) = global g_CustomCallbackScroll = x
SetCustomKeyCallback(x::Function) = global g_CustomCallbackKey = x
SetCustomCharCallback(x::Function) = global g_CustomCallbackChar = x

function ImGui_ImplGlfw_MouseButtonCallback(window::GLFW.Window, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint)
    global g_MouseJustPressed
    global g_CustomCallbackMousebutton
    g_CustomCallbackMousebutton != C_NULL && g_CustomCallbackMousebutton(window, button, action, mods)
    b = Cint(button)
    if action == GLFW.PRESS && b ≥ 0 && b < length(g_MouseJustPressed)
        g_MouseJustPressed[b+1] = true
    end
end

function ImGui_ImplGlfw_ScrollCallback(window::GLFW.Window, xoffset, yoffset)
    global g_CustomCallbackScroll
    g_CustomCallbackScroll != C_NULL && g_CustomCallbackScroll(window, xoffset, yoffset)
    io = GetIO()
    io.MouseWheelH += Cfloat(xoffset)
    io.MouseWheel += Cfloat(yoffset)
end

function ImGui_ImplGlfw_KeyCallback(window::GLFW.Window, key, scancode, action, mods)
    global g_CustomCallbackKey
    g_CustomCallbackKey != C_NULL && g_CustomCallbackKey(window, key, scancode, action, mods)

    io = GetIO()
    action == GLFW.PRESS && Set_KeysDown(io, key, true)
    action == GLFW.RELEASE && Set_KeysDown(io, key, false)
    # modifiers are not reliable across systems
    io.KeyCtrl = Get_KeysDown(io, GLFW.KEY_LEFT_CONTROL) || Get_KeysDown(io, GLFW.KEY_RIGHT_CONTROL)
    io.KeyShift = Get_KeysDown(io, GLFW.KEY_LEFT_SHIFT) || Get_KeysDown(io, GLFW.KEY_RIGHT_SHIFT)
    io.KeyAlt = Get_KeysDown(io, GLFW.KEY_LEFT_ALT) || Get_KeysDown(io, GLFW.KEY_RIGHT_ALT)
    io.KeySuper = Get_KeysDown(io, GLFW.KEY_LEFT_SUPER) || Get_KeysDown(io, GLFW.KEY_RIGHT_SUPER)
end

function ImGui_ImplGlfw_CharCallback(window::GLFW.Window, c)
    global g_CustomCallbackChar
    g_CustomCallbackChar != C_NULL && g_CustomCallbackChar(window, c)
    io = GetIO()
    0 < Cuint(c) < 0x10000 && AddInputCharacter(io, c)
end


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Heisenberg group" begin
    G = HeisenbergGroup(1)
    @test repr(G) == "HeisenbergGroup(1)"

    pts = [
        [1.0 2.0 3.0; 0.0 1.0 -1.0; 0.0 0.0 1.0],
        [1.0 4.0 -3.0; 0.0 1.0 3.0; 0.0 0.0 1.0],
        [1.0 -2.0 1.0; 0.0 1.0 1.1; 0.0 0.0 1.0],
    ]
    Xpts = [
        [0.0 2.0 3.0; 0.0 0.0 -1.0; 0.0 0.0 0.0],
        [0.0 4.0 -3.0; 0.0 0.0 3.0; 0.0 0.0 0.0],
        [0.0 -2.0 1.0; 0.0 0.0 1.1; 0.0 0.0 0.0],
    ]

    ge = allocate(pts[1])
    identity_element!(G, ge)
    @test isapprox(ge, I; atol=1e-10)

    @test check_point(G, [0.0 2.0 3.0; 0.0 1.0 -1.0; 0.0 0.0 1.0]) isa DomainError
    @test check_point(G, [1.0 2.0 3.0; 0.0 1.0 -1.0; 1.0 0.0 1.0]) isa DomainError
    @test check_point(G, [1.0 2.0 3.0; 0.0 2.0 -1.0; 0.0 0.0 1.0]) isa DomainError
    @test check_point(G, [1.0 2.0 3.0; 0.0 1.0 -1.0; 0.0 0.0 2.0]) isa DomainError
    @test check_point(G, [1.0 2.0 3.0; 0.0 1.0 -1.0; 0.0 1.0 1.0]) isa DomainError
    @test check_point(G, [1.0 2.0 3.0; 1.0 1.0 -1.0; 0.0 0.0 1.0]) isa DomainError

    @test check_vector(G, pts[1], [1.0 2.0 3.0; 0.0 0.0 -1.0; 0.0 0.0 0.0]) isa DomainError
    @test check_vector(G, pts[1], [0.0 2.0 3.0; 1.0 0.0 -1.0; 0.0 0.0 0.0]) isa DomainError
    @test check_vector(G, pts[1], [0.0 2.0 3.0; 0.0 0.0 -1.0; 0.0 0.0 2.0]) isa DomainError

    @test compose(G, pts[1], pts[2]) ≈ pts[1] * pts[2]
    test_group(
        G,
        pts,
        Xpts,
        Xpts;
        test_diff=true,
        test_invariance=true,
        test_lie_bracket=true,
        test_adjoint_action=true,
    )

    test_manifold(
        G,
        pts;
        basis_types_to_from=(DefaultOrthonormalBasis(),),
        parallel_transport=true,
        test_injectivity_radius=true,
        test_project_point=true,
        test_project_tangent=true,
        test_musical_isomorphisms=false,
    )
end


text:
"""
Loads all the tiles as per config file and generate a tileid to pass map
"""
function map_TILEID2PASS(config,group)
    #load all the files load tiles from all the pass
    if(group=="NS" || group=="")
       sky=nothing
    else
       sky=group
    end

    tiles_dic=load_DESI_tiles(config["TILES"]["tile_file"],config["TILES"]["PROGRAM"],
                  nothing,sky)
    tiles_pass=Dict{Integer,Integer}()

    for ii in 1:size(tiles_dic["TILEID"],1)
        tiles_pass[tiles_dic["TILEID"][ii]]=tiles_dic["PASS"][ii]
    end

    return tiles_pass
end

"""For each zone and each and given tracer writes the fits file including FBA information
config: input config dictionary loaded from config file
tracer: One of the tracer in the targets categories
group: The overall group we are working with this should be consistent among all the steps
mypart: which partition of the zones needs to be executed
npart_zone: How many partition we want for zones, ech partition we work in parallell cpu
columns: The list of columns tow write apart from FBA related output read from tile files
nbit_pack: How many bits to pack together, the code works for 64 or 32 but FITSIO library can-not write
UInt64 and hence only 32 will work at the moment.
"""
function PostFBA_JLD2FITS(config,tracer,group,mypart_zone,npart_zone,;
        columns=["RA","DEC","FITS_index_julia"],nbit_pack=32)
    #= for the group in config file
    scan alls the zones and write them to fits including the FBA info from tile file
    =#

    @assert ((nbit_pack==32) || (nbit_pack==64)) "Invalid bitpacking size $(nbit_pack), must be 32 or 64"

    jldfile=config["target"][tracer]["JLDfile"]
    num_fba=config["NumFBARealization"]

    #file=jldopen(jldfile,"r")
    #@show keys(file[group]["zone_77"])
    #close(file)
    #return
    #load the zones in group
    jld_order=traverse_group(jldfile,group)

    nzone=size(jld_order["zone_order"],1)
    izone_min,izone_max=partition_indices(nzone,npart_zone,mypart_zone)

    #load the listing to identify the zones which has tile
    listing=load(jldfile,"$(group)/listing")

    #maximum number of objects among the zones of interest
    max_nobj=maximum(jld_order["nobj"][izone_min:izone_max])

    #declare array to store info to be written
    tracer_targets=Dict{String,Array}("RA"=>zeros(Float32,max_nobj),"DEC"=>zeros(Float32,max_nobj))
    for col in columns
        if(col in ["FITS_index_julia","TILEID"])
            tracer_targets[col]=zeros(Int64,max_nobj)
        else
            tracer_targets[col]=zeros(Float32,max_nobj)
        end
    end

    #FBA related quantities
    tracer_targets["FBA_PASS_R1"]=zeros(Int8,max_nobj) .-1
    tracer_targets["FBA_LOCATION_R1"]=zeros(Int16,max_nobj)

    #additionally we need to load the group index to match the object in FBA
    tracer_targets["index_group"]=zeros(Int64,max_nobj)
    #To keep a count of how many times and object is assigned for R1, THis is for sanity checks
    tracer_targets["Count_R1"]=zeros(Int32,max_nobj)

    #number of bitweights are needed given fba and packing rate
    nbitweight=Int(floor(num_fba/nbit_pack))
    if(nbitweight*nbit_pack<num_fba)
        println("PostP Warning: Want to pack $(num_fba) bits in $(nbit_pack), will use only first $(nbitweight*nbit_pack)")
        println("PostP Warning: Ignoring the last $(num_fba-(nbitweight*nbit_pack)) realizations")
    end

    for ii in 1:nbitweight
        if(nbit_pack==32)
            tracer_targets["FBA_BITS$(ii)"]=zeros(UInt32,max_nobj)
        else
            tracer_targets["FBA_BITS$(ii)"]=zeros(UInt64,max_nobj)
        end
    end

    #generate and tileif 2d passs maps
    tiles_pass=map_TILEID2PASS(config,group)

    #@show izone_min,izone_max, max_nobj

    for iz in izone_min:izone_max
        fill!(tracer_targets["FBA_PASS_R1"],-1)
        fill!(tracer_targets["FBA_LOCATION_R1"],-1)
        for ii in 1:nbitweight
            fill!(tracer_targets["FBA_BITS$(ii)"],0)
        end

        #check the number of tiles in the zone
        #if it is one with tileid=0 mean none of the tiles in zone and hence skip
        if(listing["TILEID"][iz]==[0])
            println("POSTP Warning: No tile in zone $(iz), skipping")
            continue
        end

        #Number of objects in this zone
        nobj_zone=jld_order["nobj"][iz]

        #poiting to the group in the file for this zone
        tg=jld_order["zone_order"][iz]

        zone_imin=jld_order["beg_index"][iz]
        zone_imax=zone_imin+jld_order["nobj"][iz]-1
        zone_nobj=jld_order["nobj"][iz]

        #indices of objects in the zone
        tracer_targets["index_group"][1:nobj_zone] .= collect(zone_imin:zone_imax)

        #collect different columns directly available in JLD2 zones file
        for col in columns
            #get the columns based on the need
            #tval=load(jldfile,"$(tg)$(col)")
            tracer_targets[col][1:nobj_zone] .= load(jldfile,"$(tg)$(col)")

            #if(col in keys(tracer_targets))
            #    tracer_targets[col]=[tracer_targets[col];tval]
            #else
            #    tracer_targets[col][1:nobj_zone]=tval
            #end
        end

        #collect the columns from FBA procedure by reading the tile files
        fill_FBA_results!(config,tiles_pass,tracer,listing["TILEID"][iz],tracer_targets,nobj_zone,
            num_fba,nbitweight,nbit_pack)

        #columns to write
        col_write=[]
        array_write=[]
        tmp_dic=Dict{String,Array}()
        for col in keys(tracer_targets)
            #commented for debugging
            if(col in ["index_group","Count_R1"]) #columns to skip
                continue
            end

            #append!(col_write,[col])
            #append!(array_write,view(tracer_targets[col],1:nobj_zone))
            #println(col,nobj_zone,size(tracer_targets[col]))
            tmp_dic[col]=view(tracer_targets[col],1:nobj_zone)
        end

        #fits file to write this zone
        outfits=ZoneFITS_FileName(config,iz,tracer,group)
        FITS(outfits,"w") do fout
            #write(fout,tracer_targets)
            #write(fout,col_write,array_write)
            write(fout,tmp_dic) #Here we write to the fits file
        end
        println(now(UTC)," written: $(nobj_zone)",outfits)
        #break

    end #end of the iz loop
end

function fill_FBA_results!(config,tiles_pass,tracer,tile_list,tracer_targets,nobj_zone,
            num_fba,nbitweight,nbit_pack)
    #index to sort the objects in tile by index field
    zone_isort=sortperm(tracer_targets["index_group"][1:nobj_zone])

    #To account for how many times an object is assigned for realization 1
    fill!(tracer_targets["Count_R1"],0)

    for tile_id in tile_list
        if(tile_id==0) #Not an actuall tile used to fill values
            continue
        end
        #FBA jldfile for the tile
        fba_jldfile=TileFBA_FileName(config,tile_id)

        index_pre_tile=load(fba_jldfile,"$(tracer)_index")
        bool_ass=load(fba_jldfile,"$(tracer)_Bool_ass")
        location_ass=load(fba_jldfile,"$(tracer)_location_ass")

        #make sure the tile was assigned with consistent num_fba
        @assert num_fba==size(bool_ass,2) "inconsistent num_fba \nfor tile: $(col_tile) in $(fba_jldfile)"

        #find the sorting permutation
        isort_perm=sortperm(index_pre_tile)

        #rearrange with sorted index for this tile
        index_pre_tile=index_pre_tile[isort_perm]
        bool_ass=bool_ass[isort_perm,:]
        location_ass=location_ass[isort_perm,:]

        #find the pass for this tile
        #@show tile_id, tiles_pass[tile_id]
        tile_pass=tiles_pass[tile_id]

        #append the info from FBA #follow the logic of sorted indices
        ctile=1
        for ti in 1:size(index_pre_tile,1)
            while(index_pre_tile[ti]>tracer_targets["index_group"][zone_isort[ctile]])
                ctile +=1
                if(ctile>nobj_zone)
                    break
                end
            end

            if(ctile>nobj_zone)
                break
            end

            if(index_pre_tile[ti]==tracer_targets["index_group"][zone_isort[ctile]])
                ind_zone=zone_isort[ctile]

                #Now fill the needed information for first realization
                if(bool_ass[ti,1])
                    #track the counts for first realization
                    #println(typeof(ind_zone),typeof(tracer_targets["index_group"]))
                    tracer_targets["Count_R1"][ind_zone] +=1
                    tracer_targets["FBA_PASS_R1"][ind_zone]= tile_pass
                    tracer_targets["FBA_LOCATION_R1"][ind_zone]=location_ass[ti,1]
                end
                #Now fill the FBA_BITS
                PackBITS!(nbitweight,nbit_pack,bool_ass[ti,:],tracer_targets,ind_zone)
            end
        end

    end
    #Sanity checks
    if(maximum(tracer_targets["Count_R1"]) > config["target"][tracer]["Num_obs"])
        nmore=size(findall(tracer_targets["Count_R1"] .> config["target"][tracer]["Num_obs"]),1)
        msg_warn="$(tracer): $(nmore) objects have more than $(config["target"][tracer]["Num_obs"]) assignment for realization 1"
        println("PostP Warning: $(msg_warn)")
    end
    #println("counts (min,max):",minimum(tracer_targets["Count_R1"]),' ',maximum(tracer_targets["Count_R1"]))
end

"""Converts and array of boolean to UInt64 representation
example: value_out =UInt64(0)
bool_in=ones(Bool,32)
bool_in[1]=false #[false,false,false]#tar[1,1:4]
#@show value_out
value_out=Bool2UInt(bool_in,value_out)
@show bool_in
@show value_out,typeof(value_out)
@show bitstring(value_out)
"""
function Bool2UInt(bool_in,value_out::UInt64)
    #@show value_out
    value_out=UInt64(0)
    for ii in 1:size(bool_in,1)
        if(bool_in[ii])
            value_out += UInt64(2)^(ii-1)
            #println(ii,' ',value_out,' ',typeof(value_out))
        end
    end
    return value_out
end

"""Converts and array of boolean to UInt32 representation
example:
value_out =UInt32(0)
bool_in=ones(Bool,32)
bool_in[1]=false #[false,false,false]#tar[1,1:4]
#@show value_out
value_out=Bool2UInt(bool_in,value_out)
@show bool_in
@show value_out,typeof(value_out)
@show bitstring(value_out)
"""
function Bool2UInt(bool_in,value_out::UInt32)
    value_out=UInt32(0)
    #@show value_out
    for ii in 1:size(bool_in,1)
        if(bool_in[ii])
            value_out += UInt32(2)^(ii-1)
            #println(ii,' ',value_out,' ',typeof(value_out))
        end
    end
    return value_out
end

"""
Transfers the bollean array to array of UInt32 or UInt64
nbitweight: Number of bitweights columns
nbit_pack: 32 or 64 the unitsize of packing
bool_in: Input boolean array of atleast nbitweight*nbit_pack size
pack_out: output array of size nbitweight and type UInt32 or UInt64 based on nbit_pack
"""
function PackBITS!(nbitweight,nbit_pack,bool_in,tracer_targets,ind_zone)
    for ii in 1:nbitweight
        imin=((ii-1)*nbit_pack)+1
        imax=imin+nbit_pack-1
        tracer_targets["FBA_BITS$(ii)"][ind_zone]=Bool2UInt(bool_in[imin:imax],
            tracer_targets["FBA_BITS$(ii)"][ind_zone])
        #println(ii,' ',bool_in[imin:imax])
        #println(bitstring(pack_out[ii]))
    end
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
using DataStructures

@dist labeled_categorical(labels, probs) = labels[Gen.categorical(probs)]

export labeled_categorical

Rot2(a::Real) = Rotations.RotMatrix{2,Float64}(a)

function floatingPosesOf(sg)
    fp = S.floatingPosesOf(sg)
    collect(values(fp))
end

POSE_DISTANCE_FUNC_RADIUS = 1.0

function euclidean_distance(pt1::SVector{3,<:Real}, pt2::SVector{3,<:Real})
    return sqrt(sum((x2 - x1)^2 for (x1, x2) in zip(pt1, pt2)))
end
function angle_distance(rot1::R.Rotation{3}, rot2::R.Rotation{3})
    q1 = PC.componentsWXYZ(Rotations.UnitQuaternion(rot1))
    q2 = PC.componentsWXYZ(Rotations.UnitQuaternion(rot2))
    return acos(min(1.0, 2 * sum(q1 .* q2)^2 - 1))
end
function pose_distance(p1::Pose, p2::Pose; r::Float64 = POSE_DISTANCE_FUNC_RADIUS)
    euclidean_distance(p1.pos, p2.pos) + r * angle_distance(p1.orientation, p2.orientation)
end

export Rot2, floatingPosesOf

function normalize_log_weights(log_weights)
    log_total_weight = Gen.logsumexp(log_weights)
    normalized_weights =  exp.(log_weights .- log_total_weight)
    normalized_weights
end

normalized_logprob(log_p) = exp.(log_p .- maximum(log_p)) ./ sum(exp.(log_p .- maximum(log_p)))

sample_n_from_categorical(weights, n) = [Gen.random(categorical,weights) for _ in 1:n]

export normalized_logprob, normalize_log_weights, sample_n_from_categorical

using Statistics

function discretize(cloud, resolution)
    round.(cloud ./ resolution) * resolution
end

function voxelize(cloud, resolution)
    cloud = round.(cloud ./ resolution) * resolution
    idxs = unique(i -> cloud[:,i], 1:size(cloud)[2])
    cloud[:, idxs]
end

function min_max(cloud::Matrix{<:Real})
    minimum(cloud;dims=2), maximum(cloud;dims=2)
end

function min_max(cloud::Vector{<:Real})
    minimum(cloud), maximum(cloud)
end

function center_cloud(cloud)
    cloud .- mean(cloud, dims=2)
end

function opposite_face(face::Symbol)::Symbol
    pairs = [[:top, :bottom],[:left, :right],[:front, :back]]
    for p in pairs
        if face ∈ p
            idx = findfirst(p .== face)
            return p[((idx%2) +1)]
        end
    end
    return nothing
end

function get_unexplained_points(full_cloud, sub_cloud; radius=1.0)
    obs_tree = NearestNeighbors.KDTree(full_cloud)

    valid = fill(true, size(full_cloud)[2])
    idxs = NearestNeighbors.inrange(obs_tree, sub_cloud, radius)
    valid[unique(vcat(idxs...))] .= false
    valid
end

function rotation_between_two_vectors(u,v)
    u = u ./ LinearAlgebra.norm(u)
    v = v ./ LinearAlgebra.norm(v)
    half = (u+v) ./ LinearAlgebra.norm(u+v)
    q = R.QuatRotation(LinearAlgebra.dot(u,half),LinearAlgebra.cross(u,half)...)
end

export voxelize, min_max, discretize, center_cloud, pose_distance


text:
import ..Datasets: dataset, datafile, reader, Mode, Test

wnli_init() = register(DataDep(
    "GLUE-WNLI",
    """
    Winograd NLI (WNLI) task (GLUE version)
    """,
    "https://dl.fbaipublicfiles.com/glue/data/WNLI.zip",
    "ae0e8e4d16f4d46d4a0a566ec7ecceccfd3fbfaa4a7a4b4e02848c0f2561ac46";
    post_fetch_method = fn -> begin
      mv(fn, "WNLI.zip")
      DataDeps.unpack("WNLI.zip")
      innerdir = "WNLI"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct WNLI <: Dataset end

function dataset(::Type{M}, d::WNLI) where M <: Mode
    ds = reader(datafile(M, d))
    header = split(take!(ds), '\t')
    field_num = length(header)
    needed_field = M == Test ? (2,3) : (2,3,4)
    rds = get_channels(String, length(needed_field))

    task = @async begin
        for s ∈ ds
            s = split(s, '\t')
            length(s) != field_num && continue
            for (i, j) ∈ enumerate(needed_field)
                put!(rds[i], s[j])
            end
        end
    end
    for rd ∈ rds
        bind(rd, task)
    end

    rds
end

function trainfile(::WNLI)
    datadep"GLUE-WNLI/train.tsv"
end

function devfile(::WNLI)
    datadep"GLUE-WNLI/dev.tsv"
end

function testfile(::WNLI)
    datadep"GLUE-WNLI/test.tsv"
end

get_labels(::WNLI) = ("0", "1")


text:
using ADSeismic
using ADCME
using PyPlot
using DelimitedFiles
using JLD2
using Statistics
mpi_init()

M = N = Int(round(sqrt(mpi_size())))
n = 100
@info M
param = MPIAcousticPropagatorParams(NX = n*M, NY = n*N, n = n,
     Rcoef=0.2, DELTAX=10, DELTAY=10, DELTAT=0.05, NSTEP = 2000)
compute_PML_Params!(param)

srci = [n*param.M÷5]
srcj = [n*param.N÷2]
srcv = hcat([Ricker(param, 100.0, 500.0)]...)
src = MPIAcousticSource(param, srci,srcj,srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
receiver =  MPIAcousticReceiver(param, rcvi, rcvj)

c = 1000*ones(n, n)
center = n*M÷2
width = n*M ÷ 8
for i = 1:n
    for j = 1:n
        ii = i + (param.II-1) * n
        jj = j + (param.JJ-1) * n
        x = i*param.DELTAX + (param.II-1) * n * param.DELTAX
        y = j*param.DELTAY + (param.JJ-1) * n * param.DELTAY
        if center -width <= ii <= center + width && center -width <=jj<=center +width
            c[i, j] = 2000.
        end
    end
end

propagator = MPIAcousticPropagatorSolver(param, src, c)
MPISimulatedObservation!(propagator, receiver)

sess = Session()

rcvv = run(sess, receiver.rcvv)
@save "data/dat$(mpi_rank())-$(mpi_size()).jld2" rcvv

mpi_finalize()


text:
"""
    new_linear_constraint(com::CS.CoM, func::SAF{T}, set) where {T<:Real}

Create a new linear constraint and return a `LinearConstraint` with already a correct index
such that it can be simply added with [`add_constraint!`](@ref)
"""
function new_linear_constraint(com::CS.CoM, func::SAF{T}, set) where {T<:Real}
    lc = get_linear_constraint(func, set)
    lc.idx = length(com.constraints) + 1
    return lc
end

function get_linear_constraint(func::SAF{T}, set) where {T<:Real}
    func = remove_zero_coeff(func)

    indices = [v.variable_index.value for v in func.terms]

    return LinearConstraint(0, func, set, indices)
end

function remove_zero_coeff(func::MOI.ScalarAffineFunction)
    terms = [term for term in func.terms if term.coefficient != 0]
    return MOI.ScalarAffineFunction(terms, func.constant)
end

"""
    get_indices(func::VAF{T}) where {T}

Get indices from the VectorAffineFunction
"""
function get_indices(func::VAF{T}) where {T}
    return [v.scalar_term.variable_index.value for v in func.terms]
end

"""
    get_inner_constraint

Create the inner constraint of a reified or indicator constraint
"""
function get_inner_constraint(com, func::VAF{T}, set::Union{ReifiedSet, IndicatorSet}, inner_set::MOI.AbstractVectorSet) where {T<:Real}
    f = MOIU.eachscalar(func)
    inner_internals = ConstraintInternals(
        0,
        f[2:end],
        set.set,
        get_indices(f[2:end]),
    )
    return init_constraint_struct(com, set.set, inner_internals)
end

"""
    get_inner_constraint(com, func::VAF{T}, set::Union{ReifiedSet, IndicatorSet}, inner_set::Union{ReifiedSet{A}, IndicatorSet{A}, MOI.IndicatorSet{A}}) where {A,T<:Real}

Create the inner constraint when the inner constraint is a reified or indicator constraint as well.
"""
function get_inner_constraint(com, func::VAF{T}, set::Union{ReifiedSet, IndicatorSet}, inner_set::Union{ReifiedSet{A}, IndicatorSet{A}, MOI.IndicatorSet{A}}) where {A,T<:Real}
    f = MOIU.eachscalar(func)
    inner_internals = ConstraintInternals(
        0,
        f[2:end],
        set.set,
        get_indices(f[2:end]),
    )

    inner_constraint = get_inner_constraint(com, f[2:end], inner_set, inner_set.set)
    complement_inner = get_complement_constraint(com, inner_constraint)
    indices = inner_internals.indices
    activator_internals = get_activator_internals(A, indices)
    if inner_set isa ReifiedSet
        constraint =
            ReifiedConstraint(inner_internals, activator_internals, inner_constraint, complement_inner)
    else
        constraint =
            IndicatorConstraint(inner_internals, activator_internals, inner_constraint)
    end
    return constraint
end

function get_inner_constraint(com, func::VAF{T}, set::Union{ReifiedSet, IndicatorSet, MOI.IndicatorSet}, inner_set::MOI.AbstractScalarSet) where {T<:Real}
    inner_terms = [v.scalar_term for v in func.terms if v.output_index == 2]
    inner_constant = func.constants[2]
    inner_set = set.set

    inner_func = MOI.ScalarAffineFunction{T}(inner_terms, inner_constant)

    return get_linear_constraint(inner_func, inner_set)
end

function get_inner_constraint(com, vars::MOI.VectorOfVariables, set::Union{ReifiedSet, IndicatorSet}, inner_set)
    inner_internals = ConstraintInternals(
        0,
        MOI.VectorOfVariables(vars.variables[2:end]),
        set.set,
        Int[v.value for v in vars.variables[2:end]],
    )
    return init_constraint_struct(com, set.set, inner_internals)
end

function get_activator_internals(A, indices)
    ActivatorConstraintInternals(A, indices[1] in indices[2:end], false, 0)
end

"""
    saf_is_svf(saf::MOI.ScalarAffineFunction)

Checks if a `ScalarAffineFunction` can be represented as a `SingleVariable`.
This can be used for example when having a `VectorAffineFunction` in `AllDifferentSet` constraint when the decision
has to be made whether a new constraint + variable has to be created.
"""
function is_svf(saf::MOI.ScalarAffineFunction)
    !iszero(saf.constant) && return false
    length(saf.terms) != 1 && return false
    !isone(saf.terms[1].coefficient) && return false
    return true
end

function get_extrema(model::Optimizer, saf::MOI.ScalarAffineFunction{T}) where T
    min_val = saf.constant
    max_val = saf.constant
    for term in saf.terms
        if term.coefficient < 0
            min_val += term.coefficient*model.variable_info[term.variable_index.value].upper_bound
            max_val += term.coefficient*model.variable_info[term.variable_index.value].lower_bound
        else
            min_val += term.coefficient*model.variable_info[term.variable_index.value].lower_bound
            max_val += term.coefficient*model.variable_info[term.variable_index.value].upper_bound
        end
    end
    return min_val, max_val
end

function is_discrete_saf(saf::MOI.ScalarAffineFunction{T}) where T
    !isapprox(saf.constant, round(saf.constant)) && return false, saf.constant
    for term in saf.terms
        !isapprox(term.coefficient, round(term.coefficient)) && return false, term.coefficient
    end
    return true, 0
end


text:
# This implements some of the interface of the Graphs.jl package:
# http://graphsjl-docs.readthedocs.org/en/latest/interface.html

using Traits
using Graphs

@traitdef Graphs_Basic{G} begin
    vertex_type(G) -> DataType
    edge_type(G) -> DataType
    # newly defined below:
    vertex_type(Type{G}) -> DataType
    edge_type(Type{G}) -> DataType
    is_directed(G) -> Bool
end

@traitdef Vertex_List{G} <: Graphs_Basic{G} begin
    # associated type
    I = Base.return_types(vertices, (G,))[1]

    num_vertices(G) -> Integer
    vertices(G) -> I

    @constraints begin
        # make sure that I is indeed iterable:
        istrait(IsIterable{I})
    end
end

@traitdef Edge_List{G} <: Graphs_Basic{G} begin
    # associated type
    V = vertex_type(G)
    E = edge_type(G)
    I = Base.return_types(edges, (G,))[1]

    # returns the number of edges contained in g.
    num_edges(G) -> Integer
    # returns an iterable view/container of all edges.
    edges(G) -> I
    # returns the source vertex of an edge e in graph g.
    source(E, G) -> V
    # returns the target vertex of an edge e in graph g.
    target(E, G) -> V

    @constraints begin
        # make sure that I is indeed iterable:
        istrait(IsIterable{I})
    end
end

@traitdef Edge_Map{G} <: Graphs_Basic{G} begin
    E = edge_type(G)

    # returns the index of a vertex e in graph g.
    edge_index(E, G) -> Integer
end

# # and more:
# vertex_map
# edge_map
# adjacency_list
# incidence_list
# bidirectional_adjacency_list
# bidirectional_incidence_list

#update vertex_type and edge_type to work on Type{}
Graphs.vertex_type{T<:AbstractGraph}(g::Type{T}) = g.parameters[1]
Graphs.edge_type{T<:AbstractGraph}(g::Type{T}) = g.parameters[2]

# Definitions from Graphs.jl/test/edgelist.jl
pairs = [(1,2), (1,3), (2,3), (2,4), (3,5), (4,5), (2,5)]
eds = Edge{Int}[Edge(i,p[1],p[2]) for (i,p) in enumerate(pairs)]

gd = simple_edgelist(5, eds)
gu = simple_edgelist(5, eds; is_directed=false)

# check traits:
GD = typeof(gd)
GU = typeof(gu)
E = Edge{Int}

@assert istrait(Graphs_Basic{GD})
@assert istrait(Graphs_Basic{GU})
@assert istrait(Vertex_List{GD})
@assert istrait(Vertex_List{GU})
@assert istrait(Edge_List{GD})
@assert istrait(Edge_List{GU})

# use for dispatch:
@traitfn function f{G, E; Edge_List{G}, Vertex_List{G}}(g::G, e::E)
    num_edges(g) + num_vertices(g) + edge_index(e)
end

try
    f(2, 3, 5) # doesn't work as Int are not graphs
catch e
    println(e) # MethodError(f,(2,3,5)
end

f(gd, eds[5]) # -> 17
# note that f will work with any types which implement Edge_List{G, E}, Vertex_List{G}


text:
"""
Contains functions for reading and writing data from files in LIBSVM-format.
"""
module LIBSVMFileIO

using SparseArrays
using Formatting

export libsvmsize, libsvmread, libsvmwrite

"""
    libsvmsize(file)

Read the number of examples/data points and the biggest non-zero feature index
from <file>.

# Example
```julia-repl
julia> nEx, nFeat = libsvmsize("path/to/libsvm/file/a1a")
(1605, 119)
```

**Note:** The number of features is not necessarily the same the number of
features given in the data set documentation. Trailing features that are zero
for all data points are not counted (and can not be counted). For instance, the
'a1a' data set have 123 features but libsvmsize return 119.

See also: [`libsvmread`](@ref), [`libsvmwrite`](@ref), [`libsvmsize`](@ref)
"""
function libsvmsize(file)
	nEx = 0
	nFeat = 0
	for ex in eachline(file)
		nEx += 1

		offset = 1
		feat_idx_rx = r"[0-9]+ *:"

		while true
			m_idx = match(feat_idx_rx, ex, offset)
			m_idx == nothing && break

			offset = m_idx.offset + length(m_idx.match)
			idx = parse(Int, m_idx.match[1:(end-1)])
			nFeat = max(nFeat, idx)
		end
	end
	return nEx, nFeat
end

"""
    libsvmread(file; <keyword arguments>)
Read the data and labels from <file>. Return data as a `Vector` of
`SparseVector` and labels a `Vector`.

# Example
```julia-repl
julia> data, labels = libsvmread("path/to/libsvm/file/a1a");

julia> typeof(data)
Array{SparseArrays.SparseVector{Float64,Int64},1}

julia> typeof(labels)
Array{Int64,1}

julia> length(data)
1605

julia> data, labels = libsvmread("path/to/libsvm/file/a1a", selection=11:20);

julia length(data)
10
```

# Keyword Arguments
- `labeltype::Number`: Number type to parse labels as (default: `Int`).
- `valuetype::Real`: Number type to parse data values as (default: `Float64`).
- `size::Tuple{Integer,Integer}`: Size of data set given as `(<nbr of data
  points>, <nbr of examples>)`. If given, no parsing of the data set size is
  made. Useful for speeding up repeated partial reads.
- `selection`: Collection of indices. If given, only load the data points with
  indices in `selection`.
- `multilabel::Bool`: The data points of the data set have multiple labels
  (default: false). For multi-label data the labels are returned as a `Vector`
  of `NTuples` of <labeltype>.
- `dense::Bool`: Return data as ordinary dense vectors instead of
  `SparseVector` (default: false)

See also: [`libsvmread`](@ref), [`libsvmwrite`](@ref), [`libsvmsize`](@ref)
"""
function libsvmread(
		file;
		labeltype=Int, valuetype=Float64,
		size=nothing, selection=nothing, multilabel=false, dense=false)

	size == nothing ? ( (nEx,nFeat) = libsvmsize(file)) : ( (nEx,nFeat) = size)
	selection == nothing && (selection = 1:nEx)

	function parseexample(ex)
		entries = split(ex, ' ')

		# Labels
		if multilabel
			m = match(r" *[-+,\.0-9]+ *", ex)
			l = Tuple(parse.(labeltype, split(m.match, ',')))
		else
			m = match(r" *[-+\.0-9]+ *", ex)
			l = parse(labeltype, m.match)
		end

		# Features
		d = (dense ? zeros(valuetype,nFeat) : spzeros(valuetype,nFeat))

		offset = 1
		idx_rx = r"[0-9]+ *:"
		val_rx = r": *[+-\.0-9]+"

		while true
			m_idx = match(idx_rx, ex, offset)
			m_val = match(val_rx, ex, offset)
			(m_idx == nothing || m_val == nothing) && break
			offset = m_val.offset + length(m_val.match)

			idx = parse(Int, m_idx.match[1:(end-1)])
			val = parse(valuetype, m_val.match[2:end])
			d[idx] = val
		end
		return l, d
	end

	data_container = (dense ? Vector{valuetype} : SparseVector{valuetype,Int})
	label_container = (multilabel ? (NTuple{N,labeltype} where N) : labeltype)

	data = Vector{data_container}(undef,length(selection))
	labels = Vector{label_container}(undef, length(selection))

	open(file, "r") do f
		i_store = 1
		for i_ex = 1:nEx
			ex = readline(f)
			if i_ex in selection
				labels[i_store], data[i_store] = parseexample(ex)
				i_store += 1
			end
		end
	end

	return data, labels
end

"""
    libsvmwrite(data, labels, file; precision=6)

Write <data> and <labels> to <file> in LIBSVM-format.

# Example
```julia-repl
julia> data = [randn(5) for _ = 1:10];

julia> labels = [(rand(1:2), randn()) for _ = 1:10];

julia> libsvmwrite(data, labels, "example")

julia> for l in eachline("example") println(l) end
1,1.374311 1:0.447671 2:0.151464 3:-0.251898
1,-0.455168 1:-0.782467 2:0.160167 3:-0.850168
1,1.618833 1:-0.946307 2:-1.073485 3:0.504319
2,-1.141834 1:0.545098 2:-0.15826 3:-1.271802
1,0.70919 1:0.414993 2:0.232563 3:1.5618
```

# Arguments
- `data` and `labels` are iterable collections containing iterable collections of
  numbers.
- `precision` is the maximal number of decimal places to write to `file`.
  Trailing zeros are never written.

See also: [`libsvmread`](@ref), [`libsvmwrite`](@ref), [`libsvmsize`](@ref)
"""
function libsvmwrite(data, labels, file; precision=6)
	length(data) != length(labels) &&
		error("Number of examples and labels not the same.")

	formatnumber(x::Integer) = string(x)
	function formatnumber(x)
		fmt_spec = FormatSpec(".$(precision)f")
		str = fmt(fmt_spec, x)
		# Prune trailing zeros and decimal periods
		str = str[1:findlast(c->c!='0',str)]
		str[end] == '.' && (str = str[1:(end-1)])
		return str
	end

	function writelabels(labels, f)
		s = ""
		for l in labels
			s *= formatnumber(l)*","
		end
		write(f, s[1:(end-1)])
	end

	function writefeatures(data, f)
		s = ""
		for (i,d) in enumerate(data)
			val = formatnumber(d)
			parse(Float64, val) == 0.0 && continue # Ignore zeros
			s *= formatnumber(i)*":"*val*" "
		end
		write(f, s[1:(end-1)])
	end

	open(file, "w") do f
		for i = 1:length(data)
			writelabels(labels[i], f)
			write(f, " ")
			writefeatures(data[i], f)
			write(f, "\n")
		end
	end
end

end # Module


text:
module MiraculousNN

export fit

using CUDA
using ProgressMeter

include("metric.jl")
include("forprop.jl")
include("backprop.jl")
include("loss.jl")
include("utils.jl")

using .metric
using .forprop
using .backprop
using .loss
using .utils

CUDA.allowscalar(false)

fit(x, y, epochs::Int; params...) = fit(x, y; epochs = epochs, params...)

function fit(
    x,
    y;
    epochs = 10,
    batch_size = 16,
    learning_rate = 0.1,
    wd = 0.01,
    verbose = true,
    use_gpu = true,
)
    input_size = size(x, 1)
    output_size = size(y, 1)

    x, y = if use_gpu && CUDA.functional()
        CuArray(x), CuArray(y)
    else
        Array(x), Array(y)
    end

    nn_params = _init_layers(typeof(x), [input_size, 200, 200, output_size])

    nlayers = get_nlayers(nn_params)
    num_batches = div(size(x, 2), batch_size) - 1

    train_accuracy = 0
    nn_values = Dict()
    grads = Dict()

    for i = 1:epochs
        p = Progress(
            num_batches;
            desc = "Epoch $i",
            barglyphs = BarGlyphs("|-> |"),
            enabled = verbose,
        )

        for b_idx = 1:num_batches
            b = b_idx * batch_size
            x_batch = view(x, :, b:(b + batch_size - 1))
            y_batch = view(y, :, b:(b + batch_size - 1))

            nn_values = _forward(x_batch, nn_params)

            grads = _backward(nn_params, nn_values, x_batch, y_batch)
            _update_weights!(nn_params, grads, learning_rate, wd)

            nn_values = _forward(x, nn_params)
            pred_proba = nn_values["A$nlayers"]

            cost = _compute_cost(y, pred_proba, wd, nn_params)
            train_accuracy = metric.accuracy(y, pred_proba)

            ProgressMeter.next!(
                p;
                showvalues = [(:cost, cost), (:train_accuracy, train_accuracy)],
            )
        end
    end
end

function _init_layers(T::Type, layer_sizes)::Dict{String,T}
    params = Dict()

    for i = 2:length(layer_sizes)
        W = rand(Float64, layer_sizes[i], layer_sizes[i - 1]) * 0.01
        B = zeros(Float64, layer_sizes[i], 1)
        params["W$(i - 1)"] = T(W)
        params["B$(i - 1)"] = T(B)
    end
    return params
end

_init_layers(layer_sizes) = _init_layers(Matrix{Float64}, layer_sizes)

function _forward(x, params::Dict{String,T})::Dict where {T}
    nlayers = get_nlayers(params)
    nn_values = Dict{String,T}()

    last_output = x

    for i = 1:nlayers
        nn_values["Z$i"] = params["W$i"] * last_output .+ params["B$i"]

        nn_values["A$i"] = if i == nlayers
            forprop.softmax(nn_values["Z$i"])
        else
            forprop.relu.(nn_values["Z$i"])
        end

        last_output = nn_values["A$i"]
    end
    return nn_values
end

function _backward(
    params::Dict{String,T},
    nn_values::Dict{String,T},
    x,
    y,
)::Dict{String,T} where {T}
    nlayers = get_nlayers(params)
    num_samples = size(y, 2)
    grads = Dict()

    dA = nn_values["A$nlayers"]
    dZ = backprop.softmax.(dA, y)

    for i = nlayers:-1:1
        if i != nlayers
            dA = transpose(params["W$(i + 1)"]) * dZ
            dZ = backprop.relu.(dA, nn_values["A$i"])
        end

        next_output = i == 1 ? x : nn_values["A$(i - 1)"]

        grads["W$i"] = dZ * transpose(next_output) / num_samples
        grads["B$i"] = sum(dZ, dims = 2) / num_samples
    end

    return grads
end

get_nlayers(params::Dict)::Int = (keys(params) |> length) / 2

function _update_weights!(
    params::Dict{String,T},
    grads::Dict{String,T},
    learning_rate,
    wd,
)::Nothing where {T}
    nlayers = get_nlayers(params)

    for i = 1:nlayers
        params["W$i"] .-= learning_rate * (grads["W$i"] + 2 * wd .* params["W$i"])
        params["B$i"] .-= learning_rate * grads["B$i"]
    end
end

function _compute_cost(y, pred_proba, wd, params)
    weight_penalty = sum(sum(params[k] .^ 2) for k in keys(params) if k[1] == 'W')
    num_samples = size(y, 2)

    return (sum(loss.logistic.(y, pred_proba)) + wd * weight_penalty) / num_samples
end

end

max_stars_repo_name: prenc/MiraculousNN

max_issues_repo_name: prenc/MiraculousNN

max_forks_repo_name: prenc/MiraculousNN


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools

pdb_filepath = "/mnt/g/Python Programme/BALL.jl/1,2-Benzodiazepine.pdb"

comp = System(pdb_filepath)
comp = DataFrameSystem(comp)

print(@benchmark SSSR(comp))


text:
"""
    hamiltonian(skt::SlaterKosterTable)

Return a function equivalent to
`(args...) -> hamiltonian(skt, args...)`.

"""
hamiltonian(skt::SlaterKosterTable) = (args...) -> hamiltonian(skt, args...)

"""
    hamiltonian(skt::SlaterKosterTable, r::NTuple{3,<:Real}, ϕ1, ϕ2)

Return the hamiltonian matrix element corresponding to hopping from
`ϕ2` to `ϕ1`, where `ϕ1` is located at `r` relative to `ϕ2`, given in
Angstrom.

`ϕ1`, `ϕ2` are pairs of `element_symbol => orbital_symbol`, where
`element_symbol` is a valid `element_symbol` from `elementsymbols(skt)`,
and `orbital_symbol` are tesseral spherical harmonic orbitals.

Accepted symbols for the orbitals are:

:s, :px, :py, :pz

!!! note
    The following orbitals are not yet supported: :dx2y2, :dxy, :dxz,
    :dyz, :dz2, :f...

"""
function hamiltonian(skt::SlaterKosterTable, r, ϕ1, ϕ2)
    orbital_symbol1, orbital_symbol2 = last(ϕ1), last(ϕ2)
    if _orbitalrank(orbital_symbol2) < _orbitalrank(orbital_symbol1)
        ϕ1, ϕ2 = ϕ2, ϕ1
        r = -r
    end
    table = skt.data[first(ϕ1), first(ϕ2)]
    return hamiltonian(table, r, last(ϕ1), last(ϕ2))
end

"""
    hamiltonian(skt::SlaterKosterTable, ϕ)

Return the diagonal hamiltonian matrix element corresponding to `ϕ`.
`ϕ` is a pair of `element_symbol => orbital_symbol`, where
`element_symbol` is a valid `element_symbol` from
`elementsymbols(skt)`, and `orbital_symbol` is a tesseral spherical
harmonic orbital.

Accepted symbols for the orbitals are:

:s, :px, :py, :pz

!!! note
    The following orbitals are not yet supported: :dx2y2, :dxy, :dxz,
    :dyz, :dz2, :f...

"""
hamiltonian(skt::SlaterKosterTable, ϕ) =
    hamiltonian(skt.data[first(ϕ), first(ϕ)], last(ϕ))

function hamiltonian(skt::HomoNuclearTable, orbital_symbol::Symbol)
    orbital = Symbol(first(string(orbital_symbol)))
    df = skt.energy_table
    E = orbital === :s ? df.Es :
        orbital === :p ? df.Ep :
        orbital === :d ? df.Ed :
        orbital === :f ? df.Ef :
        error("invalid orbital symbol: $orbital_symbol")
    return electronvolt*E[]
end

function hamiltonian(skt::PrimitiveTable, r, orbital_symbol1::Symbol, orbital_symbol2::Symbol)
    # Convert distance to atomic units for indexing
    r = r/a0
    @assert length(r) == 3 "`length(r) = $(length(r)) != 3`: `r` is not a valid 3-dimensional vector: $r"
    min_dist, max_dist = extrema(skt.integral_table[!, :dist])
    norm_r = norm(r)
    @assert min_dist <= norm_r <= max_dist "`norm(r) = $norm_r` is outside table bounds: [$min_dist, $max_dist]"
    grid_dist = skt.grid_dist
    lo = floor(Int, norm_r/grid_dist)
    hi = norm_r == max_dist ? lo : lo + 1
    # To support automatic differentiation, we cannot use the
    # remainder operator (`%`) here, but subtract instead
    rem = norm_r/grid_dist - lo
    # Interpolate the table for the given distance
    slice = max(1, lo-1):min(skt.n_grid_points, hi+1)
    iprev, ilo = lo > 1 ? (1, 2) : (1, 1)
    ihi, inext = hi < skt.n_grid_points ? (ilo + 1, ilo + 2) : (ilo + 1, ilo + 1)
    table = DataFrame(Dict(
        name => spline_catmullrom(rem, col[iprev], col[ilo], col[ihi], col[inext]) for (name,col) in collect(pairs(eachcol(skt.integral_table[slice, :])))
    ))
    # Normalize the separation vector
    d = r ./ norm_r
    return electronvolt*hamiltonian(table, d, Val(orbital_symbol1), Val(orbital_symbol2))[]
end

function hamiltonian(df::DataFrame, d, ::Val{:s}, ::Val{:s})
    return df.Hss0
end

for (i, x) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
    @eval @inline function hamiltonian(df::DataFrame, d, ::Val{:s}, ::Val{Symbol("p", $x)})
        return d[$i]*df.Hsp0
    end

    for (j, y) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
        if i == j
            @eval @inline function hamiltonian(df::DataFrame, d, ::Val{Symbol("p", $x)}, ::Val{Symbol("p", $x)})
                return d[$i]^2*df.Hpp0 + (1 - d[$i]^2)*df.Hpp1
            end
        else
            @eval @inline function hamiltonian(df::DataFrame, d, ::Val{Symbol("p", $x)}, ::Val{Symbol("p", $y)})
                return d[$i]*d[$j]*(df.Hpp0 - df.Hpp1)
            end
        end
    end
end


text:
### InclusionConstraint ###
struct InclusionConstraint{SubSetType, SupSetType, KWT} <: SetConstraint
    subset::SubSetType
    supset::SupSetType
    kws::KWT
end
function need_variablify(c::InclusionConstraint)
    return need_variablify(c.subset) || need_variablify(c.supset)
end
function variablify(c::InclusionConstraint)
    return JuMP.build_constraint(error, variablify(c.subset),
                                 PowerSet(variablify(c.supset)); c.kws...)
end
function clear_spaces(c::InclusionConstraint)
    clear_spaces(c.subset)
    clear_spaces(c.supset)
end
function Sets.perspective_variable(c::InclusionConstraint)
    return synchronize_perspective(Sets.perspective_variable(c.subset),
                                   Sets.perspective_variable(c.supset))
end
function create_spaces(c::InclusionConstraint, spaces::Spaces)
    sub = create_spaces(c.subset, spaces)
    sup = create_spaces(c.supset, spaces)
    return merge_spaces(spaces, sub, sup)
end

JuMP.function_string(print_mode, c::InclusionConstraint) = string(c.subset)
function JuMP.in_set_string(print_mode, c::InclusionConstraint)
    string(print_mode == JuMP.IJuliaMode ? "\\subseteq" : "⊆", " ", c.supset)
end
struct PowerSet{S}
    set::S
end

# Fallback, might be because `subset` or `sup_powerset` is a `SetVariableRef` or
# a `Polyhedron` (which is handled by `JuMP.add_constraint`).
function JuMP.build_constraint(_error::Function, subset, sup_powerset::PowerSet;
                               kws...)
    InclusionConstraint(subset, sup_powerset.set, kws)
end

### InclusionConstraint for sets ###

## Set in Set ##
# See [LTJ18] B. Legat, P. Tabuada and R. M. Jungers.
# *Computing controlled invariant sets for hybrid systems with applications to model-predictive control*.
# 6th IFAC Conference on Analysis and Design of Hybrid Systems ADHS 2018, **2018**.
function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage,
                                                       <:Sets.LinearImage})
    return set_space(space, DualSpace)
end
# AS ⊆ S <=> S ⊆ A^{-1}S so PrimalSpace work
# AS ⊆ S <=> A^{-T}S∘ ⊆ S∘ so DualSpace work
function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage,
                                                       <:SetVariableRef})
    return space
end
# We can always transform  an ellipsoid to primal or dual space so we can handle
# any space
function set_space(space::Space,
                   ::InclusionConstraint{<:SetVariableRef,
                                         <:Sets.AbstractEllipsoid{T}}) where T<:Number
    return space
end

# S-procedure: Q ⊆ P <=> xQx ≤ 1 => xPx ≤ 1 <=> xPx ≤ xQx <=> Q - P is PSD
function JuMP.build_constraint(_error::Function,
                               subset::Sets.Ellipsoid,
                               sup_powerset::PowerSet{<:Sets.Ellipsoid};
                               S_procedure_scaling = nothing)
    @assert S_procedure_scaling === nothing || isone(S_procedure_scaling)
    Q = subset.Q
    P = sup_powerset.set.Q
    return psd_constraint(Symmetric(Q - P))
end

# S-procedure: Q ⊆ P <=> q(x) ≤ 1 => p(x) ≤ 1 <=> p(x) ≤ q(x) <= q - p is SOS
function JuMP.build_constraint(_error::Function,
                               subset::Union{Sets.PolySet,
                                             Sets.ConvexPolySet},
                               sup_powerset::PowerSet{<:Union{Sets.PolySet,
                                                              Sets.ConvexPolySet}};
                               S_procedure_scaling = nothing)
    @assert S_procedure_scaling === nothing || isone(S_procedure_scaling)
    q = subset.p
    p = sup_powerset.set.p
    JuMP.build_constraint(_error, q - p, SOSCone())
end

# S-procedure: Q ⊆ P <=> q - p is SOS
function s_procedure(model, subset, supset; S_procedure_scaling = nothing)
    q = subset.p
    p = supset.p
    if S_procedure_scaling === nothing
        S_procedure_scaling = @variable(model)
        # We want to avoid creating a non-convex problem. If one of `p` and `q`
        # is a polynomial with constant coefficients, we multiply the variable
        # by this one
        if MP.coefficienttype(q) <: Number
            s = S_procedure_scaling * q - p
        else
            s = q - S_procedure_scaling * p
        end
    else
        s = q - S_procedure_scaling * p
    end
    return JuMP.build_constraint(error, s, SOSCone())
end

# S-procedure: Q ⊆ P <=> q - p is SOS
function JuMP.add_constraint(model::JuMP.Model,
                             constraint::InclusionConstraint{<:Sets.Householder,
                                                             <:Sets.Householder},
                             name::String = "")
    return JuMP.add_constraint(model,
                               s_procedure(model, constraint.subset,
                                           constraint.supset;
                                           constraint.kws...), name)
end

_add_constraint_or_not(model, ::Nothing) = nothing
_add_constraint_or_not(model, con) = JuMP.add_constraint(model, con)

function _preprocess_domain(domain)
    # TODO `detecthlinearity!` fails to ignore `[1e-17, 0]` coming from
    #      `zero_eliminate` of `[1e-17, 0, 1]` so we remove duplicates to drop it
    domain = polyhedron(
        removeduplicates(hrep(domain), Polyhedra.default_solver(domain)),
        library(domain)
    )
    detecthlinearity!(domain)
    return domain
end
function _linear_part(model, domain)
    Λ = [zero(JuMP.AffExpr) for i in 1:fulldim(domain)]
    for (i, h) in enumerate(halfspaces(domain))
        iszero(h.β) || error("only cones are supported")
        λ = @variable(model, lower_bound = 0.0, base_name = "λ[$i]")
        Λ = MA.broadcast!(MA.add_mul, Λ, λ, h.a)
    end
    return Λ
end
function lin_in_domain(model, h::Vector, domain::Polyhedra.HRep)
    domain = _preprocess_domain(domain)
    dim(domain) <= 0 && return
    a = h + _linear_part(model, domain)
    if hashyperplanes(domain)
        # If we are in lower dimension,
        # we can reduce the size of `a` so that the linear constraint has smaller size.
        V = _linspace(domain)
        a = V' * a
    end
    return JuMP.build_constraint(error, a, MOI.Zeros(length(a)))
end
function add_constraint_inclusion_domain(
    model::JuMP.Model,
    subset::Sets.PolarPoint,
    supset::Sets.PolarPoint,
    domain::Polyhedra.Polyhedron)
    return _add_constraint_or_not(model, lin_in_domain(model, (subset.a - supset.a), domain))
end

function _quadratic_part(model, domain)
    Λ = [zero(JuMP.AffExpr) for i in 1:fulldim(domain), j in 1:fulldim(domain)]
    for (i, hi) in enumerate(halfspaces(domain))
        for (j, hj) in enumerate(halfspaces(domain))
            i <= j && break
            (iszero(hi.β) && iszero(hj.β)) || error("only cones are supported")
            A = hi.a * hj.a' + hj.a * hi.a'
            λ = @variable(model, lower_bound = 0.0, base_name = "λ[$i,$j]")
            Λ = MA.broadcast!(MA.add_mul, Λ, λ, A)
        end
    end
    return Λ
end
function _linspace(domain)
    L = Matrix{Polyhedra.coefficient_type(domain)}(undef, nhyperplanes(domain), fulldim(domain))
    for (i, h) in enumerate(hyperplanes(domain))
        iszero(h.β) || error("only cones are supported")
        L[i, :] = h.a
    end
    return LinearAlgebra.nullspace(L)
end
function psd_in_domain(model, Q::Symmetric, domain::Polyhedra.HRep)
    domain = _preprocess_domain(domain)
    dim(domain) <= 0 && return
    A = Q - _quadratic_part(model, domain)
    if hashyperplanes(domain)
        # If we are in lower dimension,
        # we can reduce the size of `A` so that the PSD constraint has smaller size.
        V = _linspace(domain)
        A = V' * A * V
    end
    return psd_constraint(Symmetric(A))
end
function lifted_psd_in_domain(model, Q::Symmetric, domain)
    # TODO `detecthlinearity!` fails to ignore `[1e-17, 0]` coming from
    #      `zero_eliminate` of `[1e-17, 0, 1]` so we remove duplicates to drop it
    domain = polyhedron(
        removeduplicates(hrep(domain), Polyhedra.default_solver(domain)),
        library(domain)
    )
    detecthlinearity!(domain)
    dim(domain) <= 0 && return
    Λ = [zero(JuMP.AffExpr) for i in 1:fulldim(domain)]
    for (i, hi) in enumerate(halfspaces(domain))
        λ = @variable(model, lower_bound = 0.0, base_name = "λ[$i]")
        Λ = MA.broadcast!(MA.add_mul, Λ, λ, hi.a)
    end
    off = Q[2:end, 1] - Λ
    A = [Q[1, 1] off'
         off     Q[2:end, 2:end] - _quadratic_part(model, domain)]
    if hashyperplanes(domain)
        V = _linspace(domain)
        V = [1 zeros(1, size(V, 2))
             zeros(size(V, 1), 1) V]
        A = V' * A * V
    end
    return psd_constraint(Symmetric(A))
end
function add_constraint_inclusion_domain(
    model::JuMP.Model,
    subset::Sets.Ellipsoid,
    supset::Sets.Ellipsoid,
    domain::Polyhedra.Polyhedron)
    return _add_constraint_or_not(model, psd_in_domain(model, Symmetric(subset.Q - supset.Q), domain))
end
function JuMP.add_constraint(model::JuMP.Model,
                             constraint::InclusionConstraint{<:Sets.Piecewise,
                                                             <:Sets.Piecewise},
                             name::String = "")
    subset = constraint.subset
    supset = constraint.supset
    for (i, si) in enumerate(subset.sets)
        for (j, sj) in enumerate(supset.sets)
            add_constraint_inclusion_domain(model, si, sj, subset.pieces[i] ∩ supset.pieces[j])
        end
    end
end

# S ⊆ T <=> T* ⊇ S*
function JuMP.build_constraint(_error::Function,
                               subset::Sets.HouseDualOf,
                               sup_powerset::PowerSet{<:Sets.HouseDualOf};
                               kws...)
    S = subset
    T = sup_powerset.set
    JuMP.build_constraint(_error, Sets.perspective_dual(T), PowerSet(Sets.perspective_dual(S)); kws...)
end

# S ⊆ T <=> polar(T) ⊆ polar(S)
function JuMP.build_constraint(_error::Function,
                               subset::Sets.Polar,
                               sup_powerset::PowerSet{<:Sets.Polar}; kws...)
    S = subset
    T = sup_powerset.set
    JuMP.build_constraint(_error, Polyhedra.polar(T), PowerSet(Polyhedra.polar(S));
                          kws...)
end

# See [LTJ18]
function JuMP.build_constraint(_error::Function,
                               subset::Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}},
                               sup_powerset::PowerSet{<:Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}}};
                               kws...)
    dim = Sets.dimension(subset)
    DynamicPolynomials.@polyvar x[1:dim]
    JuMP.build_constraint(_error, apply_map(subset, x),
                          PowerSet(apply_map(sup_powerset.set, x)); kws...)
end
function JuMP.build_constraint(_error::Function,
                               subset::Sets.AbstractSet,
                               sup_powerset::PowerSet{<:Sets.LinearPreImage{<:Sets.AbstractSet}};
                               kws...)
    x = Sets.space_variables(subset)
    JuMP.build_constraint(_error, subset,
                          PowerSet(apply_map(sup_powerset.set, x)); kws...)
end
function JuMP.build_constraint(_error::Function,
                               subset::Sets.LinearImage{<:Sets.AbstractSet},
                               sup_powerset::PowerSet{<:Sets.AbstractSet};
                               kws...)
    JuMP.build_constraint(
        _error, subset.set,
        PowerSet(Sets.LinearPreImage(sup_powerset.set, subset.A)); kws...)
end
function JuMP.build_constraint(_error::Function,
                               subset::Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}},
                               sup_powerset::PowerSet{<:Sets.AbstractSet};
                               kws...)
    x = Sets.space_variables(sup_powerset.set)
    JuMP.build_constraint(_error, apply_map(subset, x), sup_powerset; kws...)
end

## Set in Polyhedron ##
function set_space(space::Space, ::InclusionConstraint{<:SetVariableRef,
                                                       <:Polyhedra.Rep})
    return set_space(space, DualSpace)
end

# Ellipsoid #
function JuMP.add_constraint(model::JuMP.Model,
                             constraint::InclusionConstraint{<:Sets.AbstractSet,
                                                             <:Polyhedra.Rep},
                             name::String = "")
    ◯ = constraint.subset
    □ = constraint.supset
    for hp in hyperplanes(□)
        @constraint(model, ◯ ⊆ hp)
    end
    for hs in halfspaces(□)
        @constraint(model, ◯ ⊆ hs)
    end
end

#  S∘ ⊆ [⟨a,   x⟩ ≤ β]
#  S∘ ⊆ [⟨a/β, x⟩ ≤ 1]
# a/β ∈ S
function JuMP.build_constraint(_error::Function, subset::Sets.Polar,
                               sup_powerset::PowerSet{<:Polyhedra.HyperPlane})
    @assert iszero(sup_powerset.set.β) # Otherwise it is not symmetric around the origin
    JuMP.build_constraint(_error, Line(sup_powerset.set.a),
                          Polyhedra.polar(subset))
end
function JuMP.build_constraint(_error::Function, subset::Sets.Polar,
                               sup_powerset::PowerSet{<:Polyhedra.HalfSpace})
    JuMP.build_constraint(_error, ScaledPoint(sup_powerset.set.a, sup_powerset.set.β), Polyhedra.polar(subset))
end

# τ^{-1}(τ(S)*) ⊆ [⟨a, x⟩ ≤ β]
#         τ(S)* ⊆ [⟨(β, -a), (z, x)⟩ ≥ 0]
#       (β, -a) ∈ τ(S)
function JuMP.build_constraint(_error::Function, subset::Sets.PerspectiveDual,
                               sup_powerset::PowerSet{<:Polyhedra.HyperPlane})
    JuMP.build_constraint(_error, SymScaledPoint(-sup_powerset.set.a, sup_powerset.set.β), Polyhedra.polar(subset))
end
function JuMP.build_constraint(_error::Function, subset::Sets.PerspectiveDual,
                               sup_powerset::PowerSet{<:Polyhedra.HalfSpace})
    JuMP.build_constraint(_error, ScaledPoint(-sup_powerset.set.a, sup_powerset.set.β), Sets.perspective_dual(subset))
end

## Polyhedron in Set ##
function set_space(space::Space, ::InclusionConstraint{<:Polyhedra.Rep,
                                                       <:SetVariableRef})
    return set_space(space, PrimalSpace)
end

function JuMP.add_constraint(
    model::JuMP.Model,
    constraint::InclusionConstraint{
        <:Polyhedra.VRep,
        <:Union{Sets.AbstractSet{<:JuMP.AbstractJuMPScalar},
                Polyhedra.HRepElement{<:JuMP.AbstractJuMPScalar}}},
    name::String = ""
)
    □ = constraint.subset
    ◯ = constraint.supset
    for line in lines(□)
        @constraint(model, line in ◯)
    end
    for ray in rays(□)
        @constraint(model, ray in ◯)
    end
    for point in points(□)
        @constraint(model, point in ◯)
    end
end

# Approach 1:
#     [x'Qx ≤ 1] ⊆ [⟨a, x⟩ ≤ β]
# <=> [x'Qx ≤ 1] ⊆ [⟨a/β, x⟩ ≤ 1]
# <=> [x'Qx ≤ 1] ⊆ [(⟨a/β, x⟩)^2 ≤ 1]
# <=> [x'Qx ≤ 1] ⊆ [x'(aa' / β^2)x ≤ 1]
# <=> Q ⪰ aa' / β^2
# Schur's Lemma
# [β^2 a']
# [a   Q ]
# Approach 2:
#     [x'Qx ≤ 1] ⊆ [⟨a, x⟩ ≤ β]
# <=> [x'Qx ≤ 1] ⊆ [⟨-a, x⟩ ≤ β]
# <=> [x'Qx - z^2 ≤ 0] ⊆ [⟨-a, x⟩z - βz^2 ≤ 0]
# <=> λx'Qx - λz^2 ≥ ⟨-a, x⟩z - βz^2
# <=> λx'Qx - λz^2 ≥ ⟨-a, x⟩z - βz^2
# <=> λx'Qx + ⟨a, x⟩z + (β - λ) z^2 ≥ 0
# [β-λ a'/2]
# [a/2 λQ  ]
# TODO how to prove that they are equivalent ?
# Approach 1 is better if β is constant
# Approach 2 is better if Q is constant
_unstable_constantify(x) = x
function _unstable_constantify(x::JuMP.GenericAffExpr)
    return isempty(JuMP.linear_terms(x)) ? x.constant : x
end
function _psd_matrix(subset::Sets.Ellipsoid, hs::HalfSpace)
    return Symmetric([
        _unstable_constantify(hs.β)^2 hs.a'
        hs.a subset.Q
    ])
end
function JuMP.build_constraint(
    _error::Function,
    subset::Sets.Ellipsoid,
    sup_powerset::PowerSet{<:HalfSpace}
)
    return psd_constraint(Symmetric(_psd_matrix(subset, sup_powerset.set)))
end

function add_constraint_inclusion_domain(
    model::JuMP.Model,
    subset::Sets.Ellipsoid,
    supset::HalfSpace,
    domain
)
    return _add_constraint_or_not(
        model,
        lifted_psd_in_domain(model, _psd_matrix(subset, supset), domain)
    )
end

function JuMP.add_constraint(
    model::JuMP.Model,
    constraint::InclusionConstraint{
        <:Sets.Piecewise,
        <:HalfSpace},
    name::String = ""
)
    subset = constraint.subset
    supset = constraint.supset
    for (piece, set) in zip(subset.pieces, subset.sets)
        add_constraint_inclusion_domain(model, set, constraint.supset, piece)
    end
end

#     [p(x) ≤ 1] ⊆ [⟨a, x⟩ ≤ β]
# <= λ(1 - p(x)) ≤ β - ⟨a, x⟩ (necessary if p is SOS-convex)
#              0 ≤ λ(p(x) - 1) - ⟨a, x⟩ + β
# Set `x = 0`: 0 ≤ λ(p(0) - 1) + β
# hence if p(0) ≤ 1 (i.e. 0 ∈ S), then λ ≤ β / (1 - p(0))
# Homogeneous case: λ ≤ β
# Use build_constraint when SumOfSquares#66 if λ = β (e.g. homogeneous)
function JuMP.add_constraint(model::JuMP.Model,
                             constraint::InclusionConstraint{<:Union{Sets.ConvexPolySet,
                                                                     Sets.ConvexPolynomialSet},
                                                             <:HalfSpace},
                            name::String = "")
    p = Sets.gauge1(constraint.subset)
    h = constraint.supset
    x = Sets.space_variables(constraint.subset)
    hs = dot(h.a, x) - h.β
    if MP.coefficienttype(p) <: Number
        λ = @variable(model, lower_bound=0.0, base_name = "λ")
        cref = @constraint(model, λ * (p - 1) - hs in SOSCone())
    elseif Polyhedra.coefficient_type(h) <: Number
        λ = @variable(model, lower_bound=0.0, base_name = "λ")
        cref = @constraint(model, (p - 1) - λ * hs in SOSCone())
    else
        β = _unstable_constantify(h.β)
        if β isa Number
            # TODO what is a good value for β ???
            λ = β/2
        else
            λ = @variable(model, lower_bound=0.0, base_name = "λ")
        end
        cref = @constraint(model, λ * (p - 1) - hs in SOSCone())
    end
    return cref
end
function JuMP.build_constraint(_error::Function,
                               subset::Sets.Householder,
                               sup_powerset::PowerSet{<:HalfSpace})
    # 0 ≤ βz + ⟨-a, x⟩
    x = [sup_powerset.set.β; -sup_powerset.set.a]
    H = Sets._householder(subset.h)
    # The householder transformation is symmetric and orthogonal so no need to
    # worry about whether we should invert or transpose H
    y = H * x
    JuMP.build_constraint(_error, subset.set,
                          PowerSet(HalfSpace(y[2:end], -y[1])))
end


text:
type Cluster <: FeatureSpace
  vectors::Dict{Any, FeatureVector}
  vector_sum::FeatureVector
  mdata::Any
  Cluster() = new(Dict{Any,FeatureVector}(),FeatureVector())
end

# returns FeatureVector indexed by [key]
function getindex(c::Cluster, key)
  return c.vectors[key]
end

# maps a [key] to a FeatureVector [fv]
function setindex!(c::Cluster, fv::FeatureVector, key)
  if haskey(c, key)
    subtract!(c.vector_sum,c[key])
  end
  add!(c.vector_sum,fv)
  c.vectors[key] = fv
end

function haskey(c::Cluster, key)
  return Base.haskey(c.vectors, key)
end

function length(c::Cluster)
  return Base.length(c.vectors)
end

# returns all keys in the cluster
function keys(c::Cluster)
  return Base.keys(c.vectors)
end

# returns all FeaturVectors in the cluster
function values(c::Cluster)
  return Base.values(c.vectors)
end

# check to see if the Cluster is empty.
function isempty(c::Cluster)
  return Base.isempty(c.vectors)
end

# returns cluster centroid
function centroid(c::Cluster)
  return c.vector_sum/length(c)
end

# returns the distance between the centroids of the provided clusters
function distance(c1::Cluster, c2::Cluster, dist::Function = dist_cos)
  return dist(centroid(c1),centroid(c2))
end

# average distance to centroid
function dist_centroid(cluster::Cluster, dist_func::Function=dist_euclidean, norm::Bool=true)
  features = values(cluster.vectors)
  cluster_avg_dist = 0
  cent = centroid(cluster)
  for fv in features
    if dist_func == dist_euclidean
      cluster_avg_dist += dist_func(cent,fv,norm)^2
    else
      cluster_avg_dist += dist_func(cent,fv,norm)
    end
  end
  return (cluster_avg_dist/length(features))
end

# creates distance matrix for a cluster
function dist_matrix(cluster::Cluster, dist_func::Function=dist_euclidean, norm::Bool=true, write_csv::Bool=false)
  len = length(cluster)
  dmatrix = Array(Number,(len,len))
  fvs = collect(keys(cluster))
  for i in [1:len]
    for j in [i:len]
      if i == j
        dmatrix[i,j] = 0
      else
        dist = dist_func(cluster[fvs[i]],cluster[fvs[j]],norm)
        dmatrix[i,j] = dist
        dmatrix[j,i] = dist
      end
    end
  end

  if write_csv
    f = open(join(["matrix_",string(dist_cos),".cvs"],""), "w")
    write(f, string(dist_func))
    for fv in fvs[1:len]
      write(f,join(["\t",string(fv)],""))
    end
    write(f,"\n")
    for i in [1:len]
      write(f,string(fvs[i]))
      for j in [1:len]
        write(f,join(["\t",string(dmatrix[i,j])],""))
      end
      write(f,"\n")
    end
    close(f)
  end

  return dmatrix
end

# prints out the cluster cleanly.
function display(cl::Cluster)
  print_with_color(:white,string(typeof(cl)))
  k = length(cl)

  print_with_color(:white," with $k FeatureVectors")
  if k > 0
    print_with_color(:white,":\n")
    if k > 5
      k = 5
    end

    for key in collect(keys(cl))[1:k]
      l = length(cl[key])
      if length(key) > 40
        key = key[1:40]
        key = key*"…"
      end
      print_with_color(:white,"    $key with $l features\n")
    end
    if length(cl) > 5
      print_with_color(:white,"\t\t\t⋮\n")
    end
  end
end

function show(io::IO, cl::Cluster)
  print(io,string(typeof(cl)))
  k = length(cl.vectors)
  print(io," with $k FeatureVectors")
end


text:
using Metrics
using Metrics: precision

@testset "Classification" begin

    yb_true = [1, 1, 0, 0]
    yb_pred = [0.9, 0.1, 0.1, 0.1]
    y_true = Metrics.onehot_encode([1, 2, 1, 3, 0], 0:3)
    y_sparse = [1, 2, 1, 3, 0]
    y_pred = [.1 .8 .05 .05; .05 .5 .35 .1; .1 .6 .2 .1; .0 .1 .6 .3; .85 .1 .0 .05]'
    weights = [0.3, 0.1, 0.4, 0.2]

    @test binary_accuracy(yb_pred, yb_true) == 0.75
    @test categorical_accuracy(y_pred, y_true) == 0.6
    @test sparse_categorical(y_pred, y_sparse) == 0.6

    @test top_k_categorical(y_pred, y_true; k=2) == 1.0
    @test top_k_sparse_categorical(y_pred, y_sparse; k=2) == 1.0

    @test precision(y_pred, y_true) == 0.4999999999999999
    @test precision(y_pred, y_true; avg_type="micro") == 0.4
    @test precision(y_pred, y_true; avg_type="weighted") == 0.09999999999999998
    @test precision(y_pred, y_true; avg_type="weighted", sample_weights=weights) == 0.14999999999999997

    @test recall(y_pred, y_true) == 0.33333333333333326
    @test recall(y_pred, y_true; avg_type="micro") == 0.4
    @test recall(y_pred, y_true; avg_type="weighted") == 0.06666666666666665
    @test recall(y_pred, y_true; avg_type="weighted", sample_weights=weights) == 0.08333333333333331

    @test f_beta_score(y_pred, y_true) == 0.3999999999999998
    @test f_beta_score(y_pred, y_true; β=2) == 0.4545454545454544

    @test specificity(y_pred, y_true) == 0.825
    @test specificity(y_pred, y_true; avg_type="micro") == 0.8
    @test specificity(y_pred, y_true; avg_type="weighted") == 0.19
    @test specificity(y_pred, y_true; avg_type="weighted", sample_weights=weights) == 0.22249999999999998

    @test cohen_kappa(y_pred, y_true) == 0.23051948051948054
end


text:
# Solves the SVM classification problem using the barrier method.

# The main function.
# This implementation is deterministic for illustration purpose.
# In practice we may want to use random starting point.
function svmbarrier(X::AbstractMatrix,
                    Y::AbstractVector,
                    C::Real,
                    ϵ::Real;
                    t0::Real=1.,
                    μ::Real=100.,
                    α::Real=0.01, β::Real=0.05)
  w = zeros(size(X, 2))
  z = ones(Y)*2
  x0 = [w; z]
  A = [-Y.*X (-eye(size(Y, 1))); zeros(X) (-eye(size(Y, 1)))]
  b = [-1.*ones(Y); zeros(Y)]
  wsize = size(w, 1)
  f(x) = sum(x[1: wsize].^2)/2 + C*sum(x[wsize+1: end])
  ∇f(x) = [x[1: size(w, 1)]; C*ones(Y)]
  Hf(x) = Diagonal([ones(wsize); zeros(Y)])
  x, αdual, numstepsarray = barriermethod(
    f, ∇f, Hf, A, b, x0; ϵ=ϵ, t0=t0, μ=μ, α=α, β=β)
  x[1: wsize], αdual[1: size(Y, 1)], numstepsarray
end

# Implements the barrier method given the objective function and its
# first and second derivatives.
# The problem is supposed to have a linear constraint Ax ≤ b.
function barriermethod(f, ∇f, Hf,
                       A::AbstractMatrix,
                       b::AbstractVector,
                       x0::AbstractVector;
                       ϵ::Real=1e-3,
                       t0::Real=1.,
                       μ::Real=100.,
                       α::Real=0.01, β::Real=0.05)
  (all(A*x0.<b)
   || throw(DomainError("x0 must be strictly feasible.")))
  m = float(size(Y, 1))
  x, t = x0, t0
  numstepsarray::Vector{Int} = []
  while true
    x, numsteps = centeringstep(f, ∇f, Hf, A, b, x, t; α=α, β=β)
    push!(numstepsarray, numsteps)
    m/t < ϵ && break
    t *= μ
  end
  αdual = 1./(t*(b-A*x))
  x, αdual, numstepsarray
end

# Carries out a single centering step of the barrier method.
# The problem is supposed to have a linear constraint Ax ≤ b.
function centeringstep(f, ∇f, Hf,
                       A::AbstractMatrix,
                       b::AbstractVector,
                       x0::AbstractVector,
                       t::Real;
                       ϵ::Real=1e-8,
                       α::Real=0.01, β::Real=0.05)
  x = copy(x0)
  numsteps = 0
  obj(x) = t*f(x)-sum(log.(b-A*x))
  while true
    d = 1./(b-A*x)
    ∇ = t*∇f(x)+A'*d
    H = t*Hf(x)+A'*Diagonal(d)^2*A
    Δxnt = -H\∇
    λ² = -∇⋅Δxnt
    λ²/2 < ϵ && break
    tbacktrace = 1
    xtest = x + tbacktrace*Δxnt
    while(any(A*xtest .≥ b) || obj(xtest) > obj(x)-tbacktrace*α*λ²)
      tbacktrace *= β
      xtest = x + tbacktrace*Δxnt
    end
    # numeric limit
    obj(xtest) == obj(x) && break
    x = xtest
    numsteps += 1
  end
  x, numsteps
end


text:
export test_cell

test_cell(p=u0(), l=8, d=6, h=3, t=0.03, w2wr=0.9;
        location=p, length=l, depth=d, height=h,
        wall_thickness=t, slab_thickness=t, roof_thickness=t,
        window_to_wall_ratio=w2wr) =
  let path = rectangular_path(location, length, depth),
      w2wr1 = (1-window_to_wall_ratio)/2
    with_slab_family(thickness=slab_thickness) do
      slab(path)
    end
    with_wall_family(thickness=wall_thickness) do
      with_window_family(width=window_to_wall_ratio*length, height=window_to_wall_ratio*height) do
        w = wall(path, top_level=height)
        add_window(w, add_xy(location, length*w2wr1, height*w2wr1))
      end
    end
    with_roof_family(thickness=roof_thickness) do
      roof(path, level=h)
    end
  end


text:
using IntcodeVM, Sockets, LightGraphs

const port = 60015
IntcodeVM.run_async("input.txt", port)
socket = connect(port)

const directions = [[0,1],[0,-1],[-1,0],[1,0]]
const unwind = [2, 1, 4, 3]
distances = Dict([0,0] => 0.0)
oxygen_location = [Inf, Inf]
xmin = xmax = ymin = ymax = 0

function dfs(point)
    global xmin = min(xmin, point[1])
    global xmax = max(xmax, point[1])
    global ymin = min(ymin, point[2])
    global ymax = max(ymax, point[2])
    # We are at some location (x,y).
    # We need to recursively explore in every direction.
    # Don't go back the same way you came.
    # You can tell where you can from by distance (parent is closer to start).
    for i = 1:4
        # Possible direction north/south/west/east.
        dst = directions[i] + point

        # An infinite distance to this point means we know there is a wall
        # at this location.
        get(distances, dst, NaN) == Inf && continue

        # A distance <= our current position means we have explored this
        # point before, and our current position in the DFS tree can do
        # no better.
        get(distances, dst, Inf) <= distances[point] && continue

        # We don't know what this point is. Let's try it.
        println(socket, i)
        response = parse(Int, readline(socket))

        # Hit a wall. Mark this point as infinite distance and move on.
        if response == 0
            distances[dst] = Inf
        # Repair droid moved in the requested direction.
        # Update the distance to this point and recursively continue DFS.
        else
            response == 2 && println("Found oxygen at $(dst)")
            response == 2 && global oxygen_location = dst
            distances[dst] = distances[point] + 1
            dfs(dst)

            # In this recursive algorithm, the position of the robot does not
            # automatically unwind. We have to explicitly bring him back to
            # the starting position.
            println(socket, unwind[i])

            # We expect that the robot goes back to where it came from. If it
            # doesn't then there must be a flaw in our algorithm.
            parse(Int, readline(socket)) != 0 || error("Robot should be able to return")
        end
    end
end

function print_map()
    # count backwards for y direction because your terminal prints top to bottom.
    for y=ymax:-1:ymin
        for x=xmin:xmax
            if [x,y] == oxygen_location
                print('O')
            elseif [x,y] == [0,0]
                print('@')
            elseif get(distances, [x,y], Inf) == Inf
                print('#')
            else
                print('.')
            end
        end
        println()
    end
end

dfs([0,0])
print_map()
println("Day 15 Part 1: $(distances[oxygen_location])")
close(socket)

# I need to do some more thinking about this one. Fun graph problem!
# I am tempted to try to solve this with Floyd-Warshall, but with about 800
# non-wall tiles I think this approach might not be feasible.

v = sort(collect(keys(filter(isfinite ∘ last, distances))))
G = Graph(length(v))
for i=1:length(v)
    for j=1:4
        # Interestingly, a binary search through the sorted array makes no real difference.
        # Reminds me of Bjarne Stroustrup's comments on the performance of linear searches
        # against vectors (https://www.youtube.com/watch?v=YQs6IC-vgmo).
        if length(searchsorted(v, v[i] + directions[j])) > 0
            add_edge!(G, i, searchsortedfirst(v, v[i] + directions[j]))
        end
        #if v[i] + directions[j] in v
        #    add_edge!(G, i, findfirst(==(v[i] + directions[j]), v))
        #end
    end
end

o2 = findfirst(==(oxygen_location), v)
print("Day 15 Part 2: ")
print(findmax(floyd_warshall_shortest_paths(G).dists[o2,:])[1])
println(" (using Floyd-Warshall!)")

# In retrospect, this was super overkill. When I started thinking about part 2
# I incorrectly assumed we needed the longest of all optimal paths in the network
# and somehow got mentally stuck on All-Pairs. Really you just need to run any
# SPF algorithm against the oxygen location. I think someone else figured out
# an even more clever way to do this but I have not looked at the Reddits.


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ d3fd503e-ee83-4784-9e1e-48ce486e96d7
using SparkSQL

# ╔═╡ 0cddf535-f9ba-4671-852e-7bfa6c027259
using DataFrames, Decimals, Dates

# ╔═╡ a2916264-5da9-4fee-b528-d39f6ac3dff1
SparkSQL.initJVM()

# ╔═╡ bb8df7fc-c778-4091-b809-dfad5a6be414
sprk = SparkSQL.SparkSession("spark://localhost:7077", "Julia Lang Spark App");

# ╔═╡ b44aef6c-5a0c-4cae-896b-2f2ce4bd403a
############################################################################
#
# Lesson 1:
#
# Convert a CSV file to Parquet format with Apache Spark from Julia.
#
############################################################################

# ╔═╡ d45d1b01-404a-40b3-894b-e37407b464e8
# Step 1:  Read the CSV file using SQL. Cast (convert) the values to the correct data types. The CSV file is also missing column headers so we assign column names.

# ╔═╡ 113e3688-d972-46dd-b324-5ed9a71853ca
CSVdata = sql(sprk, "SELECT _c0 AS TICKER, to_date(_c1, 'yyyymmdd') AS TRADING_DATE, CAST(_c5 AS double) AS CLOSE FROM CSV.`/tmp/sampleData.csv`;");

# ╔═╡ 6179c4a4-5922-493e-b44b-96497fe66c35
# Step 2: Create a temporary database view in order to be able refer to the data using SQL.  In this example, the view name will be "stocks".

# ╔═╡ 3a091694-5c2a-4546-883f-0de0d0700736
createOrReplaceTempView(CSVdata, "stocks");

# ╔═╡ 343a7cd4-0ed3-475b-bfae-9d0d36e7b477
# Step 3: Use SQL to create a new Parquet file in the user specified location.

# ╔═╡ 190db402-8012-4873-b2ce-e48750bb0a8f
createFile = sql(sprk, "CREATE TABLE IF NOT EXISTS stocksFile USING PARQUET LOCATION '/tmp/stocks.PARQUET' AS SELECT * FROM stocks;")

# ╔═╡ 45a048a5-811f-43a0-b0ff-451b060022eb
############################################################################
#
# Lesson 2:
#
# Using Apache Spark and Julia data together.
#
############################################################################

# ╔═╡ 2639ef0d-53aa-4775-972f-991341f6a528
# How to query Apache Spark from Julia

# ╔═╡ 0f10e1fe-4d94-4656-a833-b78ca1bb809e
sparkQuery = sql(sprk, "SELECT * FROM stocks;")

# ╔═╡ e47fa0c2-1d26-41a8-9852-02976aa30590
# How to move Apache Spark DataSets to a Julia DataFrames:

# Use the SparkSQL.jl "toJuliaDF" function to bring Spark data into a Julia DataFrame. SparkSQL.jl creates the Julia DataFrame with the correct data types and column names automatically.

# ╔═╡ 833d323d-3944-4368-ac17-8bdd25692130
fromSparkToJulia = toJuliaDF(sparkQuery)

# ╔═╡ 8f9ee109-cfc7-4191-8cf4-611434895bde
# Create a Julia DataFrame

# ╔═╡ efcb9465-9fdd-4a87-9375-435ff2df2512
JuliaDataFrame = DataFrame(tickers = ["CRM", "IBM"])

# ╔═╡ 00062b98-5a98-45b6-9985-ed8016a01505
# How to move Julia DataFrames to Apache Spark DataSets.

# Use the SparkSQL.jl "toSparkDS" function to move that Julia DataFrame to a new Apache Spark DataSet.

# ╔═╡ 525b9215-c6c0-4079-9239-71c19bf96b61
movetoSpark = toSparkDS(sprk, JuliaDataFrame)

# ╔═╡ fc0c9c17-9a7a-435e-87ac-53bed8015b7f
# Create a database view in Spark for Julia DataFrame Data.

# ╔═╡ bade3add-f080-4dd2-b3e4-726e30d060d5
createOrReplaceTempView(movetoSpark, "julia_data")

# ╔═╡ 40680e7f-24d0-417b-9f2a-78086d84404c
# How to query Apache Spark and Julia data using SQL.

# Example query shows how to filter the stocks view with a subselect on julia data.

# ╔═╡ 0c2abeda-d5dc-4f59-a1c0-337f33396a47
chainSparkAndJuliaData = sql(sprk, "SELECT * FROM stocks WHERE TICKER IN (SELECT split(value, ',' )[0] AS tickers FROM julia_data)")

# ╔═╡ 02758b94-b949-4fc0-8ce9-5205ace25096
# Bring the results back from Spark to a Julia DataFrame.

# ╔═╡ 0a2eb12b-8680-49e7-9772-3b5ba49a2a2b
chainedData = toJuliaDF(chainSparkAndJuliaData)

# ╔═╡ ea3b1579-1e75-4c95-ba2e-9374f7265d35
# The results show only CRM and IBM ticker data.  Signifying the subquery worked.

# ╔═╡ e78c7b62-3dbd-43a0-972c-bbde2e3051e1
describe(chainedData)

# ╔═╡ Cell order:
# ╠═d3fd503e-ee83-4784-9e1e-48ce486e96d7
# ╠═a2916264-5da9-4fee-b528-d39f6ac3dff1
# ╠═bb8df7fc-c778-4091-b809-dfad5a6be414
# ╠═b44aef6c-5a0c-4cae-896b-2f2ce4bd403a
# ╠═d45d1b01-404a-40b3-894b-e37407b464e8
# ╠═113e3688-d972-46dd-b324-5ed9a71853ca
# ╠═6179c4a4-5922-493e-b44b-96497fe66c35
# ╠═3a091694-5c2a-4546-883f-0de0d0700736
# ╠═343a7cd4-0ed3-475b-bfae-9d0d36e7b477
# ╠═190db402-8012-4873-b2ce-e48750bb0a8f
# ╠═45a048a5-811f-43a0-b0ff-451b060022eb
# ╠═0cddf535-f9ba-4671-852e-7bfa6c027259
# ╠═2639ef0d-53aa-4775-972f-991341f6a528
# ╠═0f10e1fe-4d94-4656-a833-b78ca1bb809e
# ╠═e47fa0c2-1d26-41a8-9852-02976aa30590
# ╠═833d323d-3944-4368-ac17-8bdd25692130
# ╠═8f9ee109-cfc7-4191-8cf4-611434895bde
# ╠═efcb9465-9fdd-4a87-9375-435ff2df2512
# ╠═00062b98-5a98-45b6-9985-ed8016a01505
# ╠═525b9215-c6c0-4079-9239-71c19bf96b61
# ╠═fc0c9c17-9a7a-435e-87ac-53bed8015b7f
# ╠═bade3add-f080-4dd2-b3e4-726e30d060d5
# ╠═40680e7f-24d0-417b-9f2a-78086d84404c
# ╠═0c2abeda-d5dc-4f59-a1c0-337f33396a47
# ╠═02758b94-b949-4fc0-8ce9-5205ace25096
# ╠═0a2eb12b-8680-49e7-9772-3b5ba49a2a2b
# ╠═ea3b1579-1e75-4c95-ba2e-9374f7265d35
# ╠═e78c7b62-3dbd-43a0-972c-bbde2e3051e1

max_stars_repo_name: propelledanalytics/Tutorials

max_issues_repo_name: propelledanalytics/Tutorials

max_forks_repo_name: propelledanalytics/Tutorials