text:
import Base: show, getindex, iterate, HasLength, HasEltype, length, sort

function sort(dag::DAG, bs::Vector{Vector{Symbol}})

  #topological_sort!(dag)

  dct = OrderedDict()
  for (i, s) in enumerate(dag.vars)
    dct[s] = i
  end

  # Arrange first sym < second sym

  for (i, s) in enumerate(bs)
    if isless(dct[s[2]], dct[s[1]])
      tmp = s[1]
      bs[i][1] = s[2]
      bs[i][2] = tmp
    end
  end
  fset = Symbol[]
  for (i, s) in enumerate(bs)
    push!(fset, s[1])
  end
  sort!(unique!(fset))
  bs2 = sort(bs; by = x -> x[1])

  # Now sort on second symbol within first symbol

  bs3 = Vector{Symbol}[]
  for s in fset
    indx = filter(x -> x[1] == s, bs2)
    sort!(indx; by = x -> dct[x[2]])
    for i in indx
      push!(bs3, i)
    end
  end
  bs3
end

function d_sep_combinations(dag::DAG, bs::Vector{Symbol})
  cs = Vector{Symbol}[Symbol[]]
  if length(bs) > 2
    for c in combinations(bs[3:end])
      push!(cs, c)
    end
  end
  bs_new = Vector{Symbol}[]
  for c in cs
    if length(c) == 0
      if d_separation(dag, bs[1], bs[2])
        push!(bs_new, [bs[1], bs[2]])
      end
    else
      if d_separation(dag, bs[1], bs[2]; c=c)
        push!(bs_new, [bs[1], bs[2], c...])
      end
    end
  end
  bs_new
end

struct BasisSet
  dag::DAG
  bs::Vector{Vector{Symbol}}

  # Inner constructor
  #BasisSet(dag, bs) = new(dag, sort(bs; by = x -> x[1]))

  function BasisSet(dag, bs)
    # Update,create and return the sorted BasisSet
    bs2 = sort(dag, bs)
    bs3 = Vector{Symbol}[]
    for b in bs2
      bs3 = vcat(bs3, d_sep_combinations(dag, b))
    end
    new(dag, bs3)
  end

end

iterate(b::BasisSet, state=1) =
  state > length(b.bs) ? nothing : (b.bs[state], state+1)

getindex(b::BasisSet, i::Int) = b.bs[i]

HasLength(b::BasisSet) = length(b.bs)

HasEltype(b::BasisSet) = eltype(b.bs)

length(b::BasisSet) = length(b.bs)

function bs_show(io::IO, bs::BasisSet)
  println("BasisSet[")
  for ci in bs.bs
      show(ConditionalIndependency(ci))
  end
  println("]")
end

show(io::IO, bs::BasisSet) = bs_show(io, bs)

export
  BasisSet


text:
mutable struct Rake
    stop_words::Vector{String}
    min_char_length::Int
    max_words_length::Int
    min_keyword_frequency::Int
    min_words_length_adj::Int
    max_words_length_adj::Int
    min_phrase_freq_adj::Int
end

struct Score
    rank::Float64
    frequency::Int
end

function Rake(stop_words; min_char_length=3, max_words_length=3, min_keyword_frequency=10,
              min_words_length_adj=1, max_words_length_adj=1, min_phrase_freq_adj=2)
    Rake(stop_words, min_char_length, max_words_length, min_keyword_frequency,
         min_words_length_adj, max_words_length_adj, min_phrase_freq_adj)
end

import Base.run

function run(self::Rake, text)

    # Convert all multiple whitespace to single whitespace
    cleaned_text = remove_redundant_whitespace(text)

    # Convert to lower case
    cleaned_text = lowercase(cleaned_text)

    # Split on punctuation
    split_on = ['.', '?', '!', '\n', ';', ',', ':', '\u2019', '\u2013']
    sentences = split(cleaned_text, split_on)

    # Split on stopwords
    phrases = split_at_stopwords(sentences, Rake.stop_words)

    # Find all possible keywords
    all_possible_keyphrases = find_all_possible_keyphrases(phrases)

    # Score all possible keywords
    keyphrase_with_score = calculate_keyphrase_scores(all_possible_keyphrases)

    # Function to compare our (keyword, score) tuples
    #tuple_compare(x,y) = x[2].rank > y[2].rank
    #sorted_keywords = sort(keyword_candidates, lt=tuple_compare)

    return keyphrase_with_score
end

ext: jl

max_stars_repo_name: EtymoIO/RAKE

max_issues_repo_name: EtymoIO/RAKE


text:
using MetabolicSim
using Test
using DiffEqBase, OrdinaryDiffEq, DiffEqCallbacks

#generate parameters
N,M = 2,2
Rm,ρ = [1.0, 1.0],[0.1, 0.1]
u = [1.0 0.2 ; 0.2 1.0]
l = [0.5 0.0 ; 0.0 0.2]
p = make_parameters(N,M,u,Rm,l,ρ)

#generate integrator
prob = ODEProblem(MetabolicSim.dx!, ones(N+M), (0,1.0), p)
integrator = DiffEqBase.init(prob, Tsit5())

#test default functions
@testset begin
    @test typeof(MetabolicSim.biomass_func_default(integrator)) == Float64

    new_p = MetabolicSim.new_parameters_default(integrator)
    @test size(new_p[1]) == (N+1, M)
    @test typeof(new_p[2]) == Float64
end

f_biomass = MetabolicSim.biomass_func_default
f_new_p = MetabolicSim.new_parameters_default
#test callback function assertion errors
@testset begin
    f_wrong_args(a) = Float64(1.0)
    f_wrong_result(integrator) = [1.0]

    @test_throws AssertionError MetabolicSim.assert_assembly_callbacks(f_wrong_args, f_new_p)
    @test_throws AssertionError MetabolicSim.assert_assembly_callbacks(f_biomass, f_wrong_args)
    @test_throws AssertionError MetabolicSim.assert_assembly_callbacks(f_wrong_result, f_new_p)
    @test_throws AssertionError MetabolicSim.assert_assembly_callbacks(f_biomass, f_wrong_result)
end

#test consumer_equilibirum assembly function
@testset begin
    #test function
    @test_nowarn add_at_equilibrium(f_biomass,f_new_p, tmax = 100.0, abstol = 1e-8, reltol = 1e-6)

    #test abs_warning
    abstol = 1.0
    abs_warn = "abstol is set greater than `1e-5` (`$abstol`). Are you sure you dont want tighter tolerances?"
    @test_logs (:warn, abs_warn) add_at_equilibrium(f_biomass,f_new_p, tmax = 100.0, abstol = abstol, reltol = 1e-6)

    #test rel_warning
    reltol = 1.0
    rel_warn = "reltol is set greater than `1e-5` (`$reltol`). Are you sure you dont want tighter tolerances?"
    @test_logs (:warn, rel_warn) add_at_equilibrium(f_biomass,f_new_p, tmax = 100.0, abstol = 1e-8, reltol = reltol)

    #test simulalation
    tmax = 500.0
    cb = add_at_equilibrium(f_biomass,f_new_p, tmax = tmax, abstol = 1e-6, reltol = 1e-6)

    @test cb[1] == :equilibrium
    @test isa(cb[2],SciMLBase.DECallback)

    sol = simulate(ones(N+M), p, cb; t_start=0.0, t_end= 1000.0)

    @test_nowarn sol = simulate(ones(N+M), p, cb; t_start=0.0, t_end= 1000.0)
    @test length(sol.u[end]) == 5
    #test assembly after tmax
    tmax_ind = findfirst(sol.t .> tmax)
    @test all(length.(sol.u[tmax_ind:end]) .== length(sol.u[tmax_ind]))
end

#test add_at_t assembly function
@testset begin
    t_stops = [1.0, 20.0, 30.0]
    @test_nowarn add_at_t(f_biomass,f_new_p,t_stops)

    #test t_tstops are the same
    cb = add_at_t(f_biomass,f_new_p,t_stops)
    cb[1] == :at_time_stop
    isa(cb[2], SciMLBase.DECallback)
    cb[3] == t_stops

    sol = simulate(ones(N+M), p, cb; t_start=0.0, t_end= 100.0)

    #test species are added at correct times
    @test length.(sol.(t_stops)) == (M + N .+ collect(0:2)) #at t_stop
    @test length.(sol.(t_stops .+ 0.1)) == (M + N .+ collect(1:3)) #after t_stops

    #test sol parameters are updated
    p_ = sol.prob.p
    @test p_.N == N+3
    @test p_.M == M
    @test size(p_.u) == (N+3, M)
    @test length(p_.Rm) == N+3

end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
__precompile__()
module vOptGeneric
import JuMP, MathOptInterface, Combinatorics
const MOI = MathOptInterface

export vModel,
    getvOptData,
    printX_E,
    getY_N,
    vSolve,
    @addobjective,
    parseMOP,
    writeMOP

include("MOP.jl")
include("algorithms.jl")

mutable struct vOptData
    objs::Vector{JuMP.GenericAffExpr}          #Objectives
	objSenses::Vector{MOI.OptimizationSense}   #Objective Senses
    Y_N::Vector{Vector{Float64}}               #Objective values for each point
    X_E::Vector{Vector{Float64}}               #Variable values for each point
end
vOptData() = vOptData([], [], [], [])

function JuMP.copy_extension_data(data::vOptData, new_model::JuMP.AbstractModel, model::JuMP.AbstractModel)
    new_objs = [JuMP.copy(obj, new_model) for obj in data.objs]
    return vOptData(new_objs, copy(data.objSenses), [], [])
end

function getvOptData(m::JuMP.Model)
    !haskey(m.ext, :vOpt) && error("This model wasn't created with vOptGeneric")
    return m.ext[:vOpt]::vOptData
end

function vModel(optimizer_factory; args...)
    m = JuMP.Model(optimizer_factory ; args...)
    m.optimize_hook = vSolve
    # m.printhook = printhook
    m.ext[:vOpt] = vOptData()
    return m
end

function vModel(; args...)
    m = JuMP.Model(; args...)
    m.optimize_hook = vSolve
    # m.printhook = printhook
    m.ext[:vOpt] = vOptData()
    return m
end

function vSolve(m::JuMP.Model ; relax=false, method=nothing, step = 1., round_results = false, verbose = true, args...)

    vd = getvOptData(m)

    if relax != false
        @warn "linear relaxation not yet implemented"
    end

    if method == :epsilon
        solve_eps(m, step, round_results, verbose ; relaxation=relax, args...)
    elseif method == :dicho || method == :dichotomy
        solve_dicho(m, round_results, verbose ; relaxation=relax, args...)
    elseif method == :Chalmet || method == :chalmet
        solve_Chalmet(m, step, verbose ; relaxation=relax, args...)
    elseif method == :lex || method == :lexico
        solve_lexico(m, verbose ; relaxation=relax, args...)
    else
        @warn("use solve(m, method = :(epsilon | dichotomy | chalmet | lexico) )")
    end

end

# function printhook(io::IO, m::Model)
#     vd = getvOptData(m)
#     for i = 1:length(vd.objs)
#         println(vd.objSenses[i] == :Min ? "Min " : "Max ", vd.objs[i])
#     end
#     str = JuMP.model_str(JuMP.REPLMode, m)
#     index = findfirst("Subject to", str)
#     index !== nothing && print(str[first(index):end])
# end

macro addobjective(m, args...)
    if length(args) != 2
        error("in @addobjective: needs three arguments: model,
                objective sense (Max or Min) and linear expression.")
    end
    m = esc(m)
    args[1] != :Min && args[1] != :Max && error("in @addobjective: expected Max or Min for objective sense, got $(args[1])")
    sense = args[1] == :Min ? MOI.MIN_SENSE : MOI.MAX_SENSE
    expr = esc(args[2])
    return quote
        f = JuMP.AffExpr() + JuMP.@expression($m, $expr)
        !isa(f, JuMP.GenericAffExpr) && error("in @addobjective : vOptGeneric only supports linear objectives")
        vd = $m.ext[:vOpt]
        push!(vd.objSenses, $(esc(sense)))
        push!(vd.objs, f)
        f
    end
end

function printX_E(m::JuMP.Model)
    vd = getvOptData(m)
    for i = 1:length(vd.Y_N)
        print(vd.Y_N[i]," : ")
        for var in JuMP.all_variables(m)
            val = JuMP.value(var, i)
            if val != 0
                if JuMP.is_binary(var) || JuMP.is_integer(var)
                    print(JuMP.name(var), "=", round(Int, val), " ")
                else
                    print(JuMP.name(var)," =", val," ")
                end
            end
        end
        println()
   end
end

function JuMP.value(v::JuMP.VariableRef, i::Int)
    vd = getvOptData(v.model)
    return vd.X_E[i][v.index.value]
end

function getY_N(m::JuMP.Model)
    return getvOptData(m).Y_N
end

end

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
function _test_saturation(data; bitchip::Int=14, bitpix::Int=32, bzero::Int=0, numkin::Int=1)

    bitchip_sat::Int = (2^bitchip - 1) * numkin
    bitpix_sat::Int = 2^bitpix - bzero

    data = convert(Array{Int,3},data)

    o1::Array{Int,1} = []
    o2::Array{Int,1} = []

    [append!(o1, CamiXon.find_all(vec(data[:,:,i]), bitchip_sat; count=true)[1] for i ∈ axes(data,3))]
    [append!(o2, CamiXon.find_all(vec(data[:,:,i]), bitpix_sat; count=true)[1] for i ∈ axes(data,3))]

    sum(o1) > 0 ? println("Warning: physical saturation detected ($(sum(o1)) times)") : false
    sum(o2) > 0 ? println("Warning: $bitpix bit saturation detected ($(sum(o2)) times)") : false

end


text:
using Base.Test
using XSV

xsv_string, xsv_data = genxsv(num_rows=100, num_cols=10)

@test length(xsv_data)==100
@test length(xsv_data[1])==10

parsed_data = readxsv(xsv_string)

for i=1:100
	for j=1:10
		@test parsed_data[i][j]==xsv_data[i][j]
	end
end


text:
abstract type Face{Df,Dp} end

struct Segment{D,T}<:Face{1,D}
  vertices::Tuple{Point{D,T},Point{D,T}}
end

struct Triangle{D,T}<:Face{2,D}
  vertices::Tuple{Point{D,T},Point{D,T},Point{D,T}}
end

struct Tetrahedron{D,T}<:Face{3,D}
  vertices::Tuple{Point{D,T},Point{D,T},Point{D,T},Point{D,T}}
end

struct Plane{D,T}
  origin::Point{D,T}
  normal::VectorValue{D,T}
end

struct CartesianPlane{D,T}
  d::Int8
  value::T
  positive::Bool
end

# Constructors

Segment(a::Point...) = Segment(a)

Triangle(a::Point...) = Triangle(a)

Tetrahedron(a::Point...) = Tetrahedron(a)

function Plane(a::Face{Df,Dp}) where {Df,Dp}
  o = center(a)
  n = normal(a)
  Plane(o,n)
end

function CartesianPlane(p::Point{D},d::Integer,orientation::Integer) where D
  @notimplementedif abs( orientation ) ≠ 1
  T = eltype(p)
  CartesianPlane{D,T}(Int8(d),p[d],orientation>0)
end

# Getters

origin(a::Plane) = a.origin

function normal(a::CartesianPlane{D,T}) where {D,T}
  n = zero( VectorValue{D,T} )
  data = Tuple(n)
  v = a.positive ? 1 : -1
  data = Base.setindex(data,v,a.d)
  VectorValue( data )
end

function Base.zero(::Type{<:Plane{D,T}}) where {D,T}
  Plane(zero(Point{D,T}),zero(Point{D,T}))
end

function Base.getindex(::Face)
  @abstractmethod
end

Base.getindex(a::Segment,i::Integer) = a.vertices[i]

Base.getindex(a::Triangle,i::Integer) = a.vertices[i]

Base.getindex(a::Tetrahedron,i::Integer) = a.vertices[i]

Base.lastindex(a::Face) = num_vertices(a)

function get_polytope(::Face)
  @abstractmethod
end

get_polytope(::Segment) = SEGMENT

function get_polytope(::Triangle)
  TRI.face_orientations[3] = -1
  TRI
end

get_polytope(::Tetrahedron) = TET

num_dims(::Face{Df}) where Df = Df

num_point_dims(::Face{Df,Dp}) where {Df,Dp} = Dp

face_dim(::Face{D}) where D = D

function num_vertices(a::Face)
  p = get_polytope(a)
  num_vertices(p)
end

function num_edges(a::Face)
  p = get_polytope(a)
  num_edges(p)
end

function num_facets(a::Face)
  p = get_polytope(a)
  num_facets(p)
end

function num_faces(a::Face,d::Integer)
  p = get_polytope(a)
  num_faces(p,d)
end

function num_faces(a::Face)
  p = get_polytope(a)
  num_faces(p)
end

get_vertex(a::Face,vertex::Integer) = get_dface(a,vertex,Val{0}())

get_edge(a::Face,edge::Integer) = get_dface(a,edge,Val{1}())

get_facet(a::Face{D},facet::Integer) where D = get_dface(a,facet,Val{D-1}())

get_cell(a::Face) = a

@inline function get_dface(a::Face,dface::Integer,::Val{d}) where d
  @notimplementedif !is_simplex(a)
  p = get_polytope(a)
  dface_to_nodes = get_face_vertices(p,d)[dface]
  dface_vertices = ntuple( i -> a[dface_to_nodes[i]], Val{d+1}() )
  simplex_face( dface_vertices )
end

function get_dface(a::Face,dface::Integer,::Val{0})
  a[dface]
end

is_simplex(a::Face) = is_simplex(get_polytope(a))

is_n_cube(a::Face) = is_n_cube(get_polytope(a))

simplex_face(v::Point...) = simplex_face(v)

@inline function simplex_face(v::NTuple{N,<:Point}) where N
  if N == 1
    v[1]
  elseif N == 2
    Segment(v)
  elseif N == 3
    Triangle(v)
  elseif N == 4
    Tetrahedron(v)
  else
    @notimplemented
  end
end

function is_cartesian(f::Face)
  @notimplementedif !is_simplex(f)
  false
end

function center(a::Face)
  c = zero( typeof(a[1]) )
  for i in 1:num_vertices(a)
    c += a[i] / num_vertices(a)
  end
  c
end

center(a::Point) = a

center(a::Plane) = origin(a)

function normal(::Face)
  @notimplemented
end

normal(a::Plane) = a.normal

function normal(a::Face{1,2})
  v = a[2] - a[1]
  n = orthogonal(v)
  n / norm(n)
end

function normal(a::Face{2,3})
  @notimplementedif !is_simplex(a)
  v1 = a[2]-a[1]
  v2 = a[3]-a[1]
  n = orthogonal(v1,v2)
  _norm = norm(n)
  @assert _norm > 0
  n / _norm
end

function normal(::Face,::Integer)
  @abstractmethod
end

function normal(f::Face{D,D},facet::Integer) where D
  p = get_polytope(f)
  n = normal( get_facet(f,facet) )
  n * get_facet_orientations(p)[facet]
end

function normal(f::Face{1,D},facet::Integer) where D
  if facet == 1
    n = f[1]-f[2]
  elseif facet == 2
    n = f[2]-f[1]
  else
    @unreachable
  end
  n / norm(n)
end

function normal(f::Face{2,3},facet::Integer)
  p = get_polytope(f)
  n = normal(f)
  _f = get_facet(f,facet)
  v = _f[2]-_f[1]
  v /= norm(v)
  v *= get_facet_orientations(p)[facet]
  n × v
end

function measure(f::Face)
  @notimplementedif !is_simplex(f)
  simplex_measure(f)
end

function measure(f::Face{1})
  norm( f[2] -f[1] )
end

function measure(pmin::Point,pmax::Point)
  v = pmax - pmin
  prod(Tuple(v))
end

function simplex_measure(f::Face{Df,Dp}) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  @notimplementedif !is_simplex(f)
  factor = 1/factorial(Df)
  v = ntuple( i -> f[i+1]-f[1], Val{Df}() )
  norm( orthogonal(v...) ) * factor
end

function simplex_measure(f::Face{2,2})
  @notimplementedif !is_simplex(f)
  factor = 1 / 2
  v1 = f[2] - f[1]
  v2 = f[3] - f[1]
  abs( v1 × v2 ) * factor
end

function simplex_measure(f::Face{3,3})
  @notimplementedif !is_simplex(f)
  factor = 1 / 6
  v1 = f[2] - f[1]
  v2 = f[3] - f[1]
  v3 = f[4] - f[1]
  abs( (v1×v2) ⋅ v3 ) * factor
end

simplex_measure(f::Face{1}) = measure(f)

Base.length(f::Face{1}) = measure(f)

function surface(f::Face{Df,Dp}) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  measure(f)
end

volume(f::Face{D,D}) where D = measure(f)

function min_height(e::Face{1})
  length(e)
end

function min_height(f::Face{Df}) where Df
  @notimplementedif !is_simplex(f)
  min_dist = Inf
  for i in 1:num_vertices(f)
    max_dist = 0.0
    vertex = f[i]
    for j in 1:num_facets(f)
      facet = get_facet(f,j)
      dist = distance(facet,vertex)
      max_dist = max(dist,max_dist)
    end
    min_dist = min(max_dist,min_dist)
  end
  min_dist
end

get_bounding_box(p::Point) = p,p

function get_bounding_box(f::Face)
  if is_cartesian(f)
    f[1],f[end]
  else
    pmin = pmax = f[1]
    for i in 1:num_vertices(f)
      pmin = Point(min.(Tuple(pmin),Tuple(f[i])))
      pmax = Point(max.(Tuple(pmax),Tuple(f[i])))
    end
    pmin,pmax
  end
end

# Distances

function signed_distance(p::Point,Π::Plane)
  o = origin( Π )
  n = normal( Π )
  (p-o) ⋅ n
end

function signed_distance(point::Point{D},Π::CartesianPlane{D}) where D
  Π.positive ? point[Π.d] - Π.value : Π.value - point[Π.d]
end

distance(a::Point,b::Point) =  norm( a - b )

distance(s::Face{1},p::Point) = distance(p,s)

function distance(p::Point{D},s::Face{1,D}) where D
  s1_s2 = s[2] - s[1]
  l = norm(s1_s2)
  v = s1_s2 / l
  s1_p = p - s[1]
  s1_projection = ( s1_p ⋅ v ) * v
  s2_projection = s1_projection - s1_s2
  p_projection = s1_projection - s1_p
  if norm( s1_projection ) ≤ l && norm( s2_projection) ≤ l
    norm(p_projection)
  else
    min( distance(s[1],p), distance(s[2],p) )
  end
end

function distance(p::Point{Dp},f::Face{Df,Dp}) where {Df,Dp}
  if contains_projection(f,p)
    distance_to_infinite_face(f,p)
  else
    distance_to_boundary(f,p)
  end
end

distance(f::Face,p::Point) = distance(p,f)

function distance_to_infinite_face(a::Point{D},b::Point{D}) where D
  distance(a,b)
end

function distance_to_infinite_face(s::Face{1,3},p::Point{3})
  p′ = projection(p,s)
  distance(p,p′)
end

function distance_to_infinite_face(f::Face{Df,Dp},p::Point{Dp}) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  n = normal(f)
  c = center(f)
  abs( (p-c) ⋅ n )
end

distance_to_infinite_face(a::Face{D,D},b::Point{D}) where D = 0.0

function distance_to_boundary(a::Face,b::Point)
  min_dist = Inf
  for i in 1:num_facets(a)
    facet = get_facet(a,i)
    if isa(facet,Point) || measure(facet) ≠ 0
      dist = distance(facet,b)
      if dist < min_dist
        min_dist = dist
      end
    end
  end
  min_dist
end

# Projections

function projection(p::Point{D},q::Point{D}) where D
  q
end

function projection(p::Point{D},s::Face{1,D}) where D
  c = center(s)
  v = s[2] - s[1]
  v = v / norm(v)
  c + ( ( p - c ) ⋅ v ) * v
end

function projection(p::Point{Dp},f::Face{Df,Dp}) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  n = normal(f)
  c = center(f)
  p + ( ( c - p ) ⋅ n ) * n
end

function projection(p::Point,Π::Plane)
  n = normal(Π)
  c = origin(Π)
  p + ( ( c - p ) ⋅ n ) * n
end

function projection(p::Point,Π::CartesianPlane)
  Point( Base.setindex(Tuple(p),Π.value,Π.d) )
end

function contains_projection(f::Face,p::Point)
  for i in 1:num_facets(f)
    facet = get_facet(f,i)
    if isa(facet,Point) || measure(facet) ≠ 0
      c = center(facet)
      n = normal(f,i)
      if (p-c) ⋅ n > 0
        return false
      end
    end
  end
  true
end

# Voxel predicates

function voxel_intersection(f::Face{1,2},pmin::Point,pmax::Point,p::Polytope)
  D = num_point_dims(f)
  for i in 1:num_vertices(f)
    v = f[i]
    voxel_intersection(v,pmin,pmax) && return true
  end
  voxel_intersection(f,pmin,pmax)
end

function voxel_intersection(f::Face{2,3},pmin::Point,pmax::Point,p::Polytope)
  D = num_point_dims(f)
  for i in 1:num_vertices(f)
    v = f[i]
    voxel_intersection(v,pmin,pmax) && return true
  end
  for i in 1:num_edges(f)
    e = get_edge(f,i)
    voxel_intersection(e,pmin,pmax) && return true
  end
  n = normal(f)
  c = center(f)
  plane = Plane(c,n)
  v_to_dists = _compute_distances(plane,pmin,pmax,p)
  all( d -> d > 0, v_to_dists ) && return false
  all( d -> d < 0, v_to_dists ) && return false
  for e in 1:num_edges(p)
    if _is_edge_cut(p,v_to_dists,e)
      i = _intersection_point(p,pmin,pmax,v_to_dists,e)
      if contains_projection(f,i)
        return true
      end
    end
  end
  false
end

function voxel_intersection(p::Point,pmin::Point,pmax::Point)
  all( Tuple(p) .> Tuple(pmin) ) || return false
  all( Tuple(p) .< Tuple(pmax) ) || return false
  true
end

function voxel_intersection(e::Face{1,D},pmin::Point{D},pmax::Point{D}) where D
  t_min,t_max = 0.0,1.0
  for d in 1:D
    p_d = e[1][d]
    v_d = e[2][d] - e[1][d]
    if v_d < 0
      v_d = -v_d
      p_d = - p_d + pmin[d] + pmax[d]
    end
    if v_d != 0
      t = ( pmin[d] - p_d ) / v_d
      if t > t_min
        t_min = t
      end
      t = ( pmax[d] - p_d ) / v_d
      if t < t_max
        t_max = t
      end
    else
      if p_d < pmin[d]
        return false
      end
      if p_d > pmax[d]
        return false
      end
    end
  end
  t_min < t_max
end

function _compute_distances(
  Π::Plane,
  pmin::Point,
  pmax::Point,
  p::Polytope{D}) where D

  @assert is_n_cube(p)
  ntuple(i->signed_distance(_get_vertex(p,pmin,pmax,i),Π),Val{2^D}())
end

function _get_vertex(p::Polytope{D},pmin::Point,pmax::Point,i::Integer) where D
  Point(ntuple(d-> (i-1) & (1<<(d-1)) == 0 ? pmin[d] : pmax[d], Val{D}()))
end

function _is_edge_cut(p::Polytope,v_to_dists::Tuple,e::Integer)
  e_to_v = get_face_vertices(p,1)[e]
  sign(v_to_dists[e_to_v[1]]) ≠ sign(v_to_dists[e_to_v[2]])
end

function _intersection_point(
  p::Polytope,
  pmin::Point,
  pmax::Point,
  v_to_dists::Tuple,
  e::Integer)

  e_to_v = get_face_vertices(p,1)[e]
  d1 = v_to_dists[e_to_v[1]]
  d2 = v_to_dists[e_to_v[2]]
  v1 = _get_vertex(p,pmin,pmax,e_to_v[1])
  v2 = _get_vertex(p,pmin,pmax,e_to_v[2])
  (v1*d2-v2*d1)/(d2-d1)
end

# Compute planes

function bisector_plane(edge::Face{1,3},Π1::Plane,Π2::Plane)
  n1 = normal(Π1)
  n2 = normal(Π2)
  n1 ⋅ n2 ≉ -1 || error("Edge too sharp")
  n = n1-n2
  if norm(n) < 1
    v = edge[2]-edge[1]
    v /= norm(v)
    _n = n1+n2
    @assert norm(_n) > 1
    _n /= norm(_n)
    n = _n × v
    @assert norm(n) ≈ 1
  end
  n /= norm(n)
  @assert norm(n) ≈ 1
  c = center(edge)
  Plane(c,n)
end

function get_cell_planes(p::Polytope,pmin::Point,pmax::Point)
  @notimplementedif !is_n_cube(p)
  D = num_dims(p)
  N = 2*D
  @assert num_facets(p) == N
  Π_cell = lazy_map(
    i -> CartesianPlane(isodd(i)*pmin+iseven(i)*pmax,D-((i-1)÷2),1),
    1:N )
  Π_ids = - ( 1:N )
  Π_inout = lazy_map( iseven, 1:N )
  Π_cell, Π_ids, Π_inout
end

## Orthogonal

function orthogonal(a::VectorValue{2})
  VectorValue( -a[2], a[1] )
end

@generated function orthogonal(a::NTuple{N,VectorValue{D}}) where {N,D}
  entries = ""
  for i in 1:D
    data = ""
    for j in 1:D-1
      for k in 1:D
        if k != i
          data *= "a[$j][$k],"
        end
      end
    end
    if iseven(D+i)
      entries *= " + "
    else
      entries *= " - "
    end
    entries *= "det(TensorValue($data)),\n"
  end
  str = "VectorValue(\n$entries)"
  Meta.parse(str)
end

function orthogonal(a::VectorValue{D}...) where D
  length(a) == D-1 || error(
      "orthogonal(::VectorValue{D}...) only defined for D-1 VectorValues{D}'s")
  orthogonal(a,)
end


text:
function skew(v::Array{Float64,1})
    S = zeros(3,3)
    S[1,:] = [0., -v[3], v[2]]
    S[2,:] = [v[3],0.,-v[1]]
    S[3,:] = [-v[2],v[1],0.]
    return S
end

const FloatInt = Union{Float64,Int}
const VectorFloatInt = AbstractVector{<:Real} # Union{Vector{Float64},Vector{Int}}

mutable struct Quaternion
    s::Float64
    v::Array{Float64,1}
    Quaternion() = new()
    Quaternion(s::FloatInt) = new(1.0,zeros(3))
    Quaternion(s::FloatInt,v::VectorFloatInt) = new(s,v)
end

Quaternion(v::VectorFloatInt) = Quaternion(v[1],v[2:4])
Quaternion(v::NTuple{4,Float64}) = Quaternion(v[1],[v[2:4]...])

mutable struct AngleAxis
    theta::Float64
    ax::Array{Float64,1}
    AngleAxis() = new()
    AngleAxis(s::FloatInt) = new(0,[1;0;0])
    AngleAxis(s::FloatInt,v::VectorFloatInt) = new(s,v)
end

const AxisAngle = AngleAxis

mutable struct SO3
    R::Array{Float64,2}
    SO3() = new()
    SO3(dummy::FloatInt) = new(Matrix{Float64}(LinearAlgebra.I, 3,3))
    SO3(r::Array{Float64,2}) = new(r)
end

# Only I UniformScaling is applicable, eg 3*I : [3 0 0; 0 3 0; 0 0 3] does not exist in SO3
SO3(us::UniformScaling) = SO3(Matrix{Float64}(I, 3,3))

mutable struct so3
    S::Array{Float64,2}
    so3() = new()
    so3(s::FloatInt) = new(zeros(3,3))
    so3(v::VectorFloatInt) = new(skew(v))
    so3(S::Array{Float64,2}) = new(S)
end

function setso3!(s::so3, v::Union{Array, SubArray})
  s.S[1,2] = -v[3]
  s.S[1,3] = v[2]
  s.S[2,3] = -v[1]
  s.S[2,1] = v[3]
  s.S[3,1] = -v[2]
  s.S[3,2] = v[1]
  nothing
end

mutable struct Euler
    R::Float64
    P::Float64
    Y::Float64
    # these will be inconsistent with R, P, Y -- not for public consumption...
    fastconvert::Quaternion
    Euler() = new()
    Euler(s::FloatInt) = new(0.0,0.0,0.0, Quaternion(0))
    Euler(r::FloatInt,p::FloatInt,y::FloatInt) = new(r,p,y, Quaternion(0))
    Euler(r::FloatInt,p::FloatInt,y::FloatInt,q::Quaternion) = new(r,p,y, q)
    Euler(v::VectorFloatInt) = new(v[1],v[2],v[3], Quaternion(0))
    Euler(v::VectorFloatInt, q::Quaternion) = new(v[1],v[2],v[3], q)
end

#explicit convertion constructors
SO3(ss::so3) = convert(SO3, ss)
SO3(A::AngleAxis) = convert(SO3, A)
SO3(q::Quaternion) = convert(SO3, q)
SO3(E::Euler) = convert(SO3, E)

Quaternion(ss::so3) = convert(Quaternion, ss)
Quaternion(R::SO3) = convert(Quaternion, R)
Quaternion(A::AngleAxis) = convert(Quaternion, A)
Quaternion(E::Euler) = convert(Quaternion, E)

AngleAxis(R::SO3) = convert(AngleAxis, R)
AngleAxis(q::Quaternion) = convert(AngleAxis, q)

Euler(SO3::SO3) = convert(Euler, SO3)
Euler(q::Quaternion) = convert(Euler, q)

mutable struct SE3
  R::SO3
  t::Vector{Float64}
  SE3() = new()
  SE3(dummy::FloatInt) = new(SO3(0.0), zeros(3))
  SE3(t::VectorFloatInt, r::SO3) = new(r,t)
  SE3(v::VectorFloatInt, E::Euler) = new(convert(SO3,E), v[1:3] )
  SE3(v::VectorFloatInt, aa::AngleAxis) = new(convert(SO3,aa), v[1:3] ) # maybe replace with promote via dispatch
  SE3(v::VectorFloatInt, q::Quaternion) = new(convert(SO3,q), v[1:3] ) # maybe replace with promote via dispatch
  SE3(M::Array{Float64,2}) = new(SO3(M[1:3,1:3]), vec(M[1:3,4]) )
end

SE3(us::UniformScaling) = SE3(zeros(Float64,3), SO3(I))
SE3(t::VectorFloatInt, us::UniformScaling) = SE3(t, SO3(I))

function fast_norm_TU(u)
  # dest[1] = ...
  n = length(u)
  T = eltype(u)
  s = zero(T)
  @fastmath @inbounds @simd for i in 1:n
      s += u[i]^2
  end
  @fastmath @inbounds return sqrt(s)
end

function normalize!(q::Quaternion, tol=1e-6)
  # Quaternion(s, [x,y,z])
  s = q.s^2
  @fastmath @inbounds @simd for i in 1:3
      s += q.v[i]^2
  end
  @fastmath mag1 = 1.0/sqrt(s)
  # @fastmath @inbounds for i in 1:3
  #   q.v[i] *= mag1
  # end
  @fastmath @inbounds q.v .*= mag1
  @fastmath q.s *= mag1
  nothing
end

function normalize(q::Quaternion, tol=0.00001)
  qq = deepcopy(q)
  normalize!(qq,tol)
  return qq
end

function normalize(v::Array{Float64,1})
  return v / fast_norm_TU(v)
end

function matrix!(M::Array{Float64,2}, a::SE3)
  M[1:3,1:3] = a.R.R
  M[1:3,4] = a.t
  nothing
end
# Return 4x4 matrix version of SE3 transform, also see matrix! for in place version
function matrix(a::SE3)
  T = Matrix{Float64}(LinearAlgebra.I,4,4)
  matrix!(T, a)
  return T
end

function q_conj!(q::Quaternion)
    normalize!(q)
    q.v = -q.v
    nothing
end
function q_conj(q::Quaternion)
    qq = deepcopy(q)
    q_conj!(qq)
    return qq
end

transpose(a::SO3) = SO3(transpose(a.R)[:,:]) # TODO use adjoint instead of collect
adjoint(a::SO3) = SO3(a.R'[:,:])
inverse(a::SO3) = transpose(a)

# put result back in a
# function A_invB!(a::SE3, b::SE3)
#  = SE3( ( matrix(b)' \ (matrix(a)') )' )
# end

function *(a::SO3, b::SO3)
  return SO3(a.R*b.R)
end

function *(a::SE3, b::SE3)
  return SE3(vec(a.R.R*b.t + a.t), a.R*b.R)
end

function *(T::SE3, p::Vector{Float64})
  @fastmath return T.R.R*p + T.t
end

function *(q1::Quaternion, q2::Quaternion)
  ee  = [q1.s; q1.v] * [q2.s; q2.v]'
  w = ee[1,1] - ee[2,2] - ee[3,3] - ee[4,4]
  x = ee[1,2] + ee[2,1] + ee[3,4] - ee[4,3]
  y = ee[1,3] - ee[2,4] + ee[3,1] + ee[4,2]
  z = ee[1,4] + ee[2,3] - ee[3,2] + ee[4,1]
  if (w < 0.0)
    return Quaternion(-w, -[x; y; z])
  end
  return Quaternion(w, [x; y; z])
end

#mangled type products, return first type or nearest Group type (doesn't return an Algebra)
*(a::SO3, bq::Quaternion) = a*convert(SO3,bq)
*(a::SO3, b::so3) = a*convert(SO3,b)
*(a::SO3, baa::AngleAxis) = a*convert(SO3,baa)

*(aq::Quaternion, b::SO3) = aq*convert(Quaternion,b)
*(aq::Quaternion, b::so3) = aq*convert(Quaternion, convert(SO3,b) )
*(aq::Quaternion, baa::AngleAxis) = aq*convert(Quaternion, baa)

*(a::AngleAxis, b::AngleAxis) = convert(AngleAxis,convert(Quaternion,a)*convert(Quaternion,b))
*(a::AngleAxis, bq::Quaternion) = convert(AngleAxis, convert(Quaternion,a)*bq)
*(a::AngleAxis, b::SO3) = convert(AngleAxis, convert(Quaternion,a)*b )
*(aa::AngleAxis, b::so3) = aa*convert(SO3,b)

*(a::so3,b::so3) = convert(SO3,a)*convert(SO3,b)
*(a::so3, b::SO3) = convert(SO3,a)*b
*(a::so3, bq::Quaternion) = convert(Quaternion, convert(SO3,a) )*bq
*(a::so3, b::AngleAxis) = convert(AngleAxis, convert(SO3,a))*b

inverse(a::SE3) = SE3( matrix(a) \ I )
# TODO -- optimize this, and abstraction is wrong here
# Xj = Xi ⊕ ΔX
# Xj ⊖ ΔX = Xi
# Xi \ Xj = ΔX   # ⊖ Xi ⊕ Xj = ΔX
A_invB(a::SE3, b::SE3) = SE3( collect( ( matrix(b)' \ (matrix(a)') )' )  )
ominus(xi::SE3, xj::SE3) = A_invB(xi,xj)
oplus(xi::SE3, xj::SE3) = xi*xj
⊖(xi::SE3, xj::SE3) = A_invB(xi,xj)
⊕(xi::SE3, Δx::SE3) = xi*Δx
# SE3(0) ⊖ xi ⊕ xj = Δx
\(xi::SE3, xj::SE3) = SE3( matrix(xi) \ matrix(xj) )

\(qi::Quaternion,qj::Quaternion) = q_conj(qi) * qj

function doinversetransform!(dst::Vector{Float64}, aTb::SE3, src::Vector{Float64})::Nothing
  mul!(dst, transpose(aTb.R.R), src)  # At_mul_B!(dst, aTb.R.R, src)
  @fastmath @inbounds for i in 1:3, j in 1:3
    dst[i] -= aTb.R.R[j,i]*aTb.t[i]
  end
  nothing
end
# bR

# comparison functions

compare(a::SO3, b::SO3; tol::Float64=1e-14) = norm(a.R*transpose(b.R) - I) < tol

# function compare(a::SE3, b::SE3; tol::Float64=1e-14)
#   norm(a.t-b.t) < tol ? nothing : return false
#   return compare(a.R,b.R, tol=tol)
# end
function compare(a::SE3,b::SE3; tol::Float64=1e-10)
  norm(a.t-b.t)<tol && TransformUtils.compare(a.R, b.R, tol=tol)
end
function compare(a::Quaternion, b::Quaternion; tol::Float64=1e-14)
  qiq = a*q_conj(b)
  return tol <= qiq.s <= 1.0+tol && norm(qiq.v) < tol
end
function compare(a::AngleAxis,b::AngleAxis; tol::Float64=1e-14)
  aTb = q_conj(convert(Quaternion,a))*b
  return compare(Quaternion(0),aTb, tol=tol)
end

# convert functions

function convert(::Type{Quaternion}, v::VectorFloatInt)
    return Quaternion(v[1],v[2:4])
end

function convert(::Type{Quaternion}, aa::AngleAxis)
    v = normalize(aa.ax)
    theta = aa.theta/2
    w = cos(theta)
    x = v[1] * sin(theta)
    y = v[2] * sin(theta)
    z = v[3] * sin(theta)
    return Quaternion(w, [x, y, z])
end

function convert(::Type{AngleAxis}, q::Quaternion)
    theta = acos(q.s) * 2.0
    if norm(q.v)==0
      return AngleAxis(theta, [1.,0,0])
    end
    return AngleAxis(theta, normalize(q.v))
end

function convert!(R::SO3, q::Quaternion)
  # q.s = q.s;
  # x = q.v[1];
  # y = q.v[2];
  # z = q.v[3];
  @inbounds begin
    if q.s < 0.0
      q.s = -q.s
      q.v[1:3] = -q.v[1:3]
    end
    # nrm = sqrt(q.s^2+sum(q.v[1:3].^2))
    # if (nrm < 0.999)
    #   println("q2C -- not a unit quaternion nrm = $(nrm)")
    #   R = Matrix{Float64}(LinearAlgebra.I, 3,3)
    # else
      normalize!(q)
      # nrm = 1.0/nrm
      # q.s *= nrm
      # q.v .*= nrm
      # # x = x*nrm
      # # y = y*nrm
      # # z = z*nrm
      w2, x2, y2, z2 = q.s*q.s, q.v[1]*q.v[1], q.v[2]*q.v[2], q.v[3]*q.v[3]
      # y2 = y*y
      # z2 = z*z
      # w2 =
      xy = 2.0*q.v[1]*q.v[2]
      xz = 2.0*q.v[1]*q.v[3]
      yz = 2.0*q.v[2]*q.v[3]
      wx = 2.0*q.s*q.v[1]
      wy = 2.0*q.s*q.v[2]
      wz = 2.0*q.s*q.v[3]
      #build matrix
      # R = zeros(3,3)
      R.R[1,1] = w2+x2-y2-z2
      R.R[1,2] = xy-wz
      R.R[1,3] = xz+wy
      R.R[2,1] = xy+wz
      R.R[2,2] = w2-x2+y2-z2
      R.R[2,3] = yz-wx
      R.R[3,1] = xz-wy
      R.R[3,2] = yz+wx
      R.R[3,3] = w2-x2-y2+z2
    # end
  end
  # return SO3(R)
  nothing
end
function convert(::Type{SO3}, q::Quaternion)
  R = SO3(0)
  convert!(R,q)
  return R
end

function convert(::Type{Quaternion}, S::SO3)
  @inbounds begin
    rot = S.R
    tr = rot[1,1] + rot[2,2] + rot[3,3];

    if (tr > 0.0)
      S = sqrt(tr+1.0) * 2.0; # S=4*qw
      qw = 0.25 * S;
      qx = (rot[3,2] - rot[2,3]) / S;
      qy = (rot[1,3] - rot[3,1]) / S;
      qz = (rot[2,1] - rot[1,2]) / S;
    else
      if ((rot[1,1] > rot[2,2]) && (rot[1,1] > rot[3,3]))
        S = sqrt(1.0 + rot[1,1] - rot[2,2] - rot[3,3]) * 2.0; # S=4*qx
        qw = (rot[3,2] - rot[2,3]) / S;
        qx = 0.25 * S;
        qy = (rot[1,2] + rot[2,1]) / S;
        qz = (rot[1,3] + rot[3,1]) / S;
      else
        if (rot[2,2] > rot[3,3])
          S = sqrt(1.0 + rot[2,2] - rot[1,1] - rot[3,3]) * 2.0 # S=4*qy
          qw = (rot[1,3] - rot[3,1]) / S;
          qx = (rot[1,2] + rot[2,1]) / S;
          qy = 0.25 * S;
          qz = (rot[2,3] + rot[3,2]) / S;
        else
          S = sqrt(1.0 + rot[3,3] - rot[1,1] - rot[2,2]) * 2.0 # S=4*qz
          qw = (rot[2,1] - rot[1,2]) / S
          qx = (rot[1,3] + rot[3,1]) / S
          qy = (rot[2,3] + rot[3,2]) / S
          qz = 0.25 * S;
        end
      end
    end
  end

  if (qw >= 0)
    q = Quaternion(qw,[qx;qy;qz])
  else
    q = Quaternion(-qw,-[qx;qy;qz])
  end
  return q
end

function convert(::Type{SO3}, aa::AngleAxis)
  convert(SO3,convert(Quaternion, aa))
end

function convert(::Type{AngleAxis}, s::SO3)
  convert(AngleAxis, convert(Quaternion, s))
end

function convert!(Eu::Euler, q::Quaternion)
  #Function to convert quaternion to Euler angles, from Titterton and Weston
  #04, following quaternion convention of Titterton

  # a = q.s
  # b = q.v[1]
  # c = q.v[2]
  # d = q.v[3]

  Eu.R = atan(2.0*(q.v[2]*q.v[3] + q.v[1]*q.s), q.s*q.s - q.v[1]^2 - q.v[2]^2 + q.v[3]^2)
  # Eu.R = -atan(2.0*(q.v[2]*q.v[3] - q.v[1]*q.s), q.s*q.s - q.v[1]^2 - q.v[2]^2 + q.v[3]^2)
  #-atan(2.0*(q.v[2]*q.v[3] - q.v[1]*q.s),1.0-2.0*(q.v[1]^2+q.v[2]^2)); # -atan(2.0*(q.v[2]*q.v[3] - q.v[1]*q.s), q.s*q.s - q.v[1]^2 - q.v[2]^2 + q.v[3]^2) #numerically more stable
  Eu.P = asin(2.0*(q.v[2]*q.s - q.v[1]*q.v[3]));
  Eu.Y =  atan(2.0*(q.v[1]*q.v[2] + q.v[3]*q.s), q.s^2 + q.v[1]^2 - q.v[2]^2 - q.v[3]^2)
  # Eu.Y = -atan(2.0*(q.v[1]*q.v[2] - q.v[3]*q.s),1.0-2.0*(q.v[2]^2+q.v[3]^2)); # -atan(2.0*(q.v[1]*q.v[2] - q.v[3]*q.s), q.s^2 + q.v[1]^2 - q.v[2]^2 - q.v[3]^2)

  nothing
end
function convert(::Type{Euler}, q::Quaternion)
  Eu = Euler(0)
  convert!(Eu, q)
  return Eu
end

function convert(::Type{so3}, G::SO3)
  so3(logmap(G))
end

function convert(::Type{so3}, q::Quaternion)
  so3(logmap(convert(SO3,q)))
end

function convert(::Type{Quaternion}, alg::so3)
  convert(Quaternion, convert(SO3, alg))
end

function convert(::Type{Euler}, R::SO3)
  convert(Euler, convert(Quaternion, R))
end

function convert!(q::Quaternion, E::Euler)::Nothing

  @fastmath halfroll = 0.5*E.R;
  @fastmath halfpitch = 0.5*E.P;
  @fastmath halfyaw = 0.5*E.Y;
  @fastmath sin_r2 = sin(halfroll);
  @fastmath sin_p2 = sin(halfpitch);
  @fastmath sin_y2 = sin(halfyaw);
  @fastmath cos_r2 = cos(halfroll);
  @fastmath cos_p2 = cos(halfpitch);
  @fastmath cos_y2 = cos(halfyaw);

  @fastmath q.s = cos_r2 * cos_p2 * cos_y2 + sin_r2 * sin_p2 * sin_y2;
  @fastmath @inbounds q.v[1] = sin_r2 * cos_p2 * cos_y2 - cos_r2 * sin_p2 * sin_y2;
  @fastmath @inbounds q.v[2] = cos_r2 * sin_p2 * cos_y2 + sin_r2 * cos_p2 * sin_y2;
  @fastmath @inbounds q.v[3] = cos_r2 * cos_p2 * sin_y2 - sin_r2 * sin_p2 * cos_y2;

  # Enforce positive scalar quaternions following conversion from Euler angles
  @inbounds if (q.s<0.0)
    q.s = -q.s
    q.v[1:3] = -q.v[1:3];
  end
  # q = q./norm(q);
  normalize!(q)
  # return Quaternion(q[1],q[2:4])
  nothing
end
  # mutable struct HalfAngles
  #   halfang::Vector{Float64}
  #   halfsin::Vector{Float64}
  #   halfcos::Vector{Float64}
  #   HalfAngles() = new()
  #   HalfAngles(::Int) = new(zeros(3), zeros(3), zeros(3))
  # end
  # unsafe_reuseE2Q::Vector{HalfAngles}
  # thr_reuse = E.unsafe_reuseE2Q[Threads.threadid()]
  # thr_reuse.halfang[1] = 0.5*E.R;
  # thr_reuse.halfang[2] = 0.5*E.P;
  # thr_reuse.halfang[3] = 0.5*E.Y;
  # thr_reuse.halfang[1] = 0.5*E.R;
  # thr_reuse.halfang[2] = 0.5*E.P;
  # thr_reuse.halfang[3] = 0.5*E.Y;
  # thr_reuse.halfsin[1] = sin(thr_reuse.halfang[1]);
  # thr_reuse.halfsin[2] = sin(thr_reuse.halfang[2]);
  # thr_reuse.halfsin[3] = sin(thr_reuse.halfang[3]);
  # thr_reuse.halfcos[1] = cos(thr_reuse.halfang[1]);
  # thr_reuse.halfcos[2] = cos(thr_reuse.halfang[2]);
  # thr_reuse.halfcos[3] = cos(thr_reuse.halfang[3]);
  # @fastmath q.s = thr_reuse.halfcos[1] * thr_reuse.halfcos[2] * thr_reuse.halfcos[3] + thr_reuse.halfsin[1] * thr_reuse.halfsin[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[1] = thr_reuse.halfsin[1] * thr_reuse.halfcos[2] * thr_reuse.halfcos[3] - thr_reuse.halfcos[1] * thr_reuse.halfsin[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[2] = thr_reuse.halfcos[1] * thr_reuse.halfsin[2] * thr_reuse.halfcos[3] + thr_reuse.halfsin[1] * thr_reuse.halfcos[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[3] = thr_reuse.halfcos[1] * thr_reuse.halfcos[2] * thr_reuse.halfsin[3] - thr_reuse.halfsin[1] * thr_reuse.halfsin[2] * thr_reuse.halfcos[3];
  # @fastmath q.s = thr_reuse.halfcos[1] * thr_reuse.halfcos[2] * thr_reuse.halfcos[3]
  # @fastmath @inbounds q.v[1] = thr_reuse.halfsin[1] * thr_reuse.halfcos[2] * thr_reuse.halfcos[3]
  # @fastmath @inbounds q.v[2] = thr_reuse.halfcos[1] * thr_reuse.halfsin[2] * thr_reuse.halfcos[3]
  # @fastmath @inbounds q.v[3] = thr_reuse.halfcos[1] * thr_reuse.halfcos[2] * thr_reuse.halfsin[3]
  # @fastmath q.s += thr_reuse.halfsin[1] * thr_reuse.halfsin[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[1] -= thr_reuse.halfcos[1] * thr_reuse.halfsin[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[2] += thr_reuse.halfsin[1] * thr_reuse.halfcos[2] * thr_reuse.halfsin[3];
  # @fastmath @inbounds q.v[3] -= thr_reuse.halfsin[1] * thr_reuse.halfsin[2] * thr_reuse.halfcos[3];

function convert(::Type{Quaternion}, E::Euler)
  q = Quaternion(0)
  convert!(q, E)
  return q
end

function convert(::Type{SO3}, E::Euler)
  return convert(SO3,convert(Quaternion, E))
end

function convert!(R::SO3, E::Euler)
  # TODO -- refactor to inplace operations
  convert!(E.fastconvert, E)
  convert!(R, E.fastconvert)
  nothing
end

# type fastconvert!
#
# end

function rotate!(q1::Quaternion, v1::Array{Float64,1})
    #v = (q1*Quaternion(0.0,v1)*q_conj(q1)).v
    R = convert(SO3, q1);
    v = R.R * v1
    for i = 1:3 v1[i] = v[i] end
    nothing
end
function rotate(q1::Quaternion, v1::Array{Float64,1})
    vv1 = deepcopy(v1)
    rotate!(q1, vv1)
    return vv1
end

function expmOwn(alg::so3)
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  # I = Matrix{Float64}(LinearAlgebra.I, 3,3)
  if (v_norm>1e-6)
    #1E-6 is chosen because double numerical LSB is around 1E-18 for nominal values [-pi..pi] and (1E-7)^2 is at 1E-14, but 1E-7 rad/s is 0.02 deg/h
    R = I + sin(v_norm)/v_norm*alg.S + (1.0-cos(v_norm))/(v_norm^2)*(alg.S^2)
  else
    R = I + alg.S + 0.5*(alg.S^2)
  end
  return R
end

function logmap(grp::SO3)
  R = grp.R
  tra = LinearAlgebra.tr(R)
  if abs(tra+1.0) < 1e-10
    if abs(R[3,3]+1.0) > 1e-10
      ret = (pi / sqrt(2.0+2.0*R[3,3] )) * [R[1,3]; R[2,3]; (1.0+R[3,3])]
    else
      if abs(grp.R[1,1]+1.0) > 1e-10
        ret = (pi / sqrt(2.0+2.0*R[2,2])) * [R[1,2]; (1.0+R[2,2]); R[3,2]]
      else
        ret = (pi / sqrt(2.0+2.0*R[1,1])) * [(1.0+R[1,1]); R[2,1]; R[3,1]]
      end
    end
  else
    tra = round(tra, digits=12)
    tr_3 = tra-3.0
    # TODO: this sign might have been the wrong way round, changed on 10/10/2018 for Taylor series when tr_3 small
    ## also introduced abs(tr_3) rather than tr_3 < 1e-7, which looks to be a bug
    if (abs(tr_3) < 1e-7)
      magnitude = 0.5 - tr_3*tr_3/12.0
    else
      theta = acos((tra-1.0)/2.0)
      magnitude = theta/(2.0*sin(theta)) # 0.5/sinc(theta)
    end
    ret = magnitude*[ (R[3,2]-R[2,3]); (R[1,3]-R[3,1]); (R[2,1]-R[1,2])]
  end
  return ret
end

function logmap(q::Quaternion)
  logmap(convert(SO3, q))
end

# Chikjain's book p??
function rightJacExmap(alg::so3)
  Gam = alg.S
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  # I = Matrix{Float64}(LinearAlgebra.I, 3,3)
  if (v_norm>1e-7)
    Jr = I + (v_norm-sin(v_norm))/(v_norm^3)*(Gam^2) - (1.0-cos(v_norm))/(v_norm^2+0.0)*Gam
  else
    Jr = I
  end
  return Jr
end

function rightJacExmapinv(alg::so3)
  Gam = alg.S
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  # I = Matrix{Float64}(LinearAlgebra.I, 3,3)
  if (v_norm>1e-7)
    brace = (  1.0/(v_norm^2) - (1.0+cos(v_norm))/(2.0*v_norm*sin(v_norm))  )*Gam
    Jr_inv = I + 0.5*Gam + brace*Gam
  else
    Jr_inv = I
  end
  return Jr_inv
end

# function computeRatesFromRPY(Q, prevQ, dt)
#   if (dt<1):
#     # Assume oneside rectangular derivative, we could use a trapezoidal
#     # approach to obtain a more accurate estimate of the rotation rate
#     # TODO: Improve derivative scheme
#     S = (    np.dot(q2R(Q),np.transpose(q2R(prevQ)))  - np.identity(3)    ) / dt;
#     bw = np.matrix([[-S[2,1].item()],[-S[0,2].item()],[-S[1,0].item()]])
#   else
#     bw = zeros(3,1);
#   end
#   return bw
# end

# See Julia's implementation the matrix exponential function -- expm!()
# https://github.com/acroy/julia/blob/0bce8951f19e8f47d361e73b5b5bd6283926c01b/base/linalg/dense.jl
expmOwnT(M::Array{Float64,2}) = (Matrix{Float64}(LinearAlgebra.I,size(M,1),size(M,1)) + 0.5*M)*((Matrix{Float64}(LinearAlgebra.I,size(M,1),size(M,1)) - 0.5*M) \ Matrix{Float64}(LinearAlgebra.I,size(M,1),size(M,1)));
expmOwn1(M::Array{Float64,2}) = I + M;
expmOwn2(M::Array{Float64,2}) = I + M + 0.5*(M^2);
function expmOwn3(M::Array{Float64,2})
  M052 = 0.5*M^2;
  return I + M + M052 + M052*M/3.0;
end
function expmOwn4(M::Array{Float64,2})
  M052 = 0.5*M^2;
  return I + M + M052 + M052*M/3.0 + M052*M052/6.0
end

function convert(::Type{SO3}, alg::so3)
  v = vee(alg)
  nv = norm(v)
  if nv < 1e-3
    return SO3(exp(alg.S))
  else
    invnv = 1.0/nv
    return SO3(I + invnv*(sin(nv)*alg.S + invnv*(1.0-cos(nv))*(alg.S^2) ) )
  end
end

# manifold distance, add, and subtract
function logmap_SO2(Rl::Matrix{<:Real})
  ct = abs(Rl[1,1]) > 1.0 ? 1.0 : Rl[1,1]  # reinserting the sign below
  -sign(Rl[2,1])*acos(ct)
end

function wrapRad(th::Real)
  #rem is not the best function here so using rem2pi(x, RoundNearest) returns in interval [-π, π]
  #TODO is [-π, π] a problem or should it strictly be [-pi, pi)
  return rem2pi(th, RoundNearest)
end
# function wrapRad(th::Float64)
  # th = rem(th, 2pi) # returns in range (-2pi,2pi)
  # now move to range [-pi, pi)
  # if th >= pi
  #   th -= 2.0*pi
  # elseif th < -pi
  #   th += 2.0*pi
  # end
  # return th
# end
# rem(pi, 2pi)
# rem(4pi, 2pi)
# rem(4.5pi, 2pi)
# rem(-4.5pi, 2pi)
# rem(7pi, 2pi)
# rem(-7pi, 2pi)
# rem(-pi, 2pi)

R(th::Real) = [[cos(th);-sin(th)]';[sin(th);cos(th)]'];
R(;x::Real=0.0,y::Real=0.0,z::Real=0.0) = convert(SO3, so3([x,y,z]))

function vee!(rv::Vector{<:Real}, alg::so3)
  rv[1] = -alg.S[2,3]
  rv[2] = alg.S[1,3]
  rv[3] = -alg.S[1,2]
  nothing
end

function vee(alg::so3)
   rv = zeros(3)
   @inbounds vee!(rv, alg)#[-alg.S[2,3], alg.S[1,3], -alg.S[1,2]]
   return rv
end

function vee!(rv::Vector{<:Real}, q::Quaternion)
  rv[1] = q.s
  rv[2:4] = q.v[1:3]
  nothing
end
function vee(q::Quaternion)
   rv = zeros(4)
   @inbounds vee!(rv, q)
   return rv
end

# vectorize SE3 group to translation and AngleAxis numbers
# dim7 out: [xyz, theta, axyz]
function veeAngleAxis(G::SE3)
  aa = convert(AngleAxis, G.R)
  v = zeros(7)
  v[1:3] = G.t
  v[4] = aa.theta
  v[5:7] = aa.ax[:]
  return v
end

function vec(q::Quaternion)
  v = zeros(4)
  v[1] = q.s
  v[2:4] = q.v[:]
  return v
end

# vectorize SE3 group to translation and Quaternion numbers
# dim7 out: [xyz, cos(th/2), sin(th/2)xyz]
function veeQuaternion(G::SE3)
  q = convert(Quaternion, G.R)
  v = zeros(7)
  v[1:3] = G.t
  v[4:7] = vee(q)
  # v[4] = q.s
  # v[5:7] = q.v[:]
  return v
end

# vectorize SE3 group to translation and Quaternion numbers
# dim7 out: [xyz, cos(th/2), sin(th/2)xyz]
function veeEuler(G::SE3)
  E = convert(Euler, G.R)
  v = zeros(6)
  v[1:3] = G.t
  v[4] = E.R
  v[5] = E.P
  v[6] = E.Y

  return v
end

function deltaso3vee(aQb::Quaternion,aQc::Quaternion)
  dq = aQb\aQc
  s = logmap(dq)
end

#


text:
# BitStreams are streams of spikes. The spikes are generated several ways:
# 1) From virtual neurons
# 2) Interneurons (IN)s
# 3) external source. Ex images
# 4) Poisson generators

# BitStream as Stimulus is in motion:

# A stream has an optional input and an output.
# Represents a stream of spikes.
abstract type AbstractBitStream end

# -------------------------------------------------------------------------
# Axons
# -------------------------------------------------------------------------

# Axons have only one input and one output, and can have an optional
# delay.

# A stream's output is routed to a Axon's input.
# On each pass the input is "shifted" to the Axon's output.
# Some Axons have zero delay so the input/outputs are the same on
# each pass (aka DirectAxon)
#
# A simulation pass has a "pre" and post step.
# The "pre" step takes the output of a stream and places on the input.
# The "post" step shifts the input through the delay towards the output.
#
# A Axon's delay is a simple representation of distance.
abstract type AbstractAxon end

abstract type AbstractDelay end

abstract type AbstractCell end
abstract type AbstractSoma end
abstract type AbstractDendrite end
abstract type AbstractCompartment end
abstract type AbstractSynapse end


text:
# This file was generated, do not modify it.

using Pkg # hideall
Pkg.activate("_literate/EX-breastcancer/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

using UrlDownload
using DataFrames
using PrettyPrinting
using PyPlot
using MLJ
ioff() # hide
MLJ.color_off(); # hide

RANDOM_SEED = 42;

url = "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data";
feature_names = ["ID", "Class", "mean radius", "mean texture", "mean perimeter", "mean area", "mean smoothness", "mean compactness", "mean concavity", "mean concave points", "mean symmetry", "mean fractal dimension", "radius error", "texture error", "perimeter error", "area error", "smoothness error", "compactness error", "concavity error", "concave points error", "symmetry error", "fractal dimension error", "worst radius", "worst texture", "worst perimeter", "worst area", "worst smoothness", "worst compactness", "worst concavity", "worst concave points", "worst symmetry", "worst fractal dimension"]
data = urldownload(url, true, format = :CSV, header = feature_names);

figure(figsize=(8, 6))
hist(data.Class)
xlabel("Classes")
ylabel("Number of samples")
plt.savefig(joinpath(@OUTPUT, "Target_class.svg")); # hide

df = DataFrame(data)[:, 2:end];

pprint(first(df,10))

pprint(describe(df))

pprint(schema(df))

coerce!(df, :Class => OrderedFactor{2});

y, X = unpack(df, ==(:Class),name->true, rng = RANDOM_SEED);

transformer_instance = Standardizer()
transformer_model = machine(transformer_instance, X)
fit!(transformer_model)
X = MLJ.transform(transformer_model, X);

train, test = partition(eachindex(y), 0.8, shuffle=true, rng=RANDOM_SEED);

for m in models(matching(X, y))
    println("Model name = ",m.name,", ","Prediction type = ",m.prediction_type,", ","Package name = ",m.package_name);
end

model_names=Vector{String}();
loss_acc=[];
loss_ce=[];
loss_f1=[];

figure(figsize=(8, 6))
for m in models(matching(X, y))
    if m.prediction_type==Symbol("probabilistic") && m.package_name=="ScikitLearn" && m.name!="LogisticCVClassifier"
        #Excluding LogisticCVClassfiier as we can infer similar baseline results from the LogisticClassifier

        #Capturing the model and loading it using the @load utility
        model_name=m.name
        package_name=m.package_name
        eval(:(clf = @load $model_name pkg=$package_name verbosity=1))

        #Fitting the captured model onto the training set
        clf_machine = machine(clf(), X, y)
        fit!(clf_machine, rows=train)

        #Getting the predictions onto the test set
        y_pred = MLJ.predict(clf_machine, rows=test);

        #Plotting the ROC-AUC curve for each model being iterated
        fprs, tprs, thresholds = roc(y_pred, y[test])
        plot(fprs, tprs,label=model_name);

        #Obtaining different evaluation metrics
        ce_loss = mean(cross_entropy(y_pred,y[test]))
        acc = accuracy(mode.(y_pred), y[test])
        f1_score = f1score(mode.(y_pred), y[test])

        #Adding the different obtained values of the evaluation metrics to the respective vectors
        push!(model_names, m.name)
        append!(loss_acc, acc)
        append!(loss_ce, ce_loss)
        append!(loss_f1, f1_score)
    end
end

#Adding labels and legend to the ROC-AUC curve
xlabel("False Positive Rate")
ylabel("True Positive Rate")
legend(loc="best", fontsize="xx-small")
title("ROC curve")
plt.savefig(joinpath(@OUTPUT, "breastcancer_auc_curve.svg")) # hide

model_info=DataFrame(ModelName=model_names,Accuracy=loss_acc,CrossEntropyLoss=loss_ce,F1Score=loss_f1);

pprint(sort!(model_info,[:CrossEntropyLoss]));

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Random

const chars = ['a':'z'; 'A':'Z'; '0':'9'; ['/','+','-','=']]

"Return a Base64 encoded random passwort of length `len`."
genpasswd(len) = randstring(chars, len)

"Return an available port number ≥ `start`."
function getPort(start::Integer)
    s = UDPSocket()
    port = start
    while true
        bind(s, Sockets.localhost, port) && break
        port += 1
        port > 65535 && throw(SystemError("no port available"))
    end
    close(s)
    port
end

"Return the `Sockets.InetAddr` of a socket `sock`."
function getHostPort(sock::UDPSocket)
    s = UDPSocket()
    p = getPort(1000)
    @assert bind(s, Sockets.localhost, p) "No port available"
    send(sock, Sockets.localhost, p, serialize(0))
    hp, _ = recvfrom(s)
    hp
end

max_stars_repo_name: pbayer/erjulix

max_issues_repo_name: pbayer/erjulix

max_forks_repo_name: pbayer/erjulix


text:
using Documenter, PoincareInvariants

DocMeta.setdocmeta!(
    PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation.PaduaTransforms,
    :DocTestSetup,
    quote
        using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation.PaduaTransforms
    end
)

makedocs(
    sitename = "PoincareInvariants.jl",
    modules=[PoincareInvariants],
    pages = [
        "Home" => "index.md",
        "First Poincaré Invariants" => "first_poincare_invariants.md",
        "Second Poincaré Invariants" => [
            "Guide" => "second_poincare_invariants.md",
            "Chebyshev Implementation" => "chebyshev_implementation.md",
            "Padua Transforms" => "padua_transforms.md"
        ],
        "Canonical Symplectic Structures" => "canonical_symplectic_structures.md",
        "Reference" => "reference.md"
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true")
)

# deploydocs(
#     repo   = "github.com/JuliaGNI/PoincareInvariants.jl"
# )


text:
module TapirBenchmarks

export avgfilter1d_seq!,
    avgfilter1d_setup,
    avgfilter1d_constprop_seq!,
    avgfilter1d_constprop_setup,
    avgfilter1d_constprop_tapir_dac!,
    avgfilter1d_constprop_tapir_seq!,
    avgfilter1d_constprop_threads!,
    avgfilter1d_tapir_dac!,
    avgfilter1d_tapir_seq!,
    avgfilter1d_threads!,
    avgfilter2d_seq!,
    avgfilter2d_setup,
    avgfilter2d_tapir_dac!,
    avgfilter2d_tapir_seq!,
    avgfilter2d_threads!,
    divide_at_mean_seq,
    divide_at_mean_tapir,
    divide_at_mean_threads,
    meanvar_seq,
    meanvar_tapir,
    meanvar_threads

using Base: Tapir
using Base.Experimental: Const
using StaticArrays: SVector
using Statistics: mean, var

macro grainsize(n::Integer)
    Expr(:loopinfo, (Symbol("tapir.loop.grainsize"), Int(n)))
end

let configpath = joinpath(@__DIR__, "config.jl"),
    defaultpath = joinpath(@__DIR__, "default-config.jl")

    if !isfile(configpath)
        cp(defaultpath, configpath)
    end
end
include("config.jl")

if USE_ALIASSCOPE
    const var"@maybe_aliasscope" = Base.Experimental.var"@aliasscope"
else
    macro maybe_aliasscope(ex)
        esc(ex)
    end
end

include("folds.jl")

include("avgfilter1d.jl")
include("avgfilter2d.jl")
include("avgfilter1d_constprop.jl")
include("meanvar.jl")
include("divide_at_mean.jl")
include("divide_at_mean_with_map.jl")

end


text:
cachetypeof(::Type{<:ShardedRecycler}) = ShardedCache

# Workaround for avoid storing types as `DataType``
ShardedRecycler(::Type{Factory}, cache::ShardedCache{T}) where {Factory,T} =
    ShardedRecycler{T,Type{Factory}}(Factory, cache)
ShardedRecycler{T}(::Type{Factory}, cache::ShardedCache{T}) where {Factory,T} =
    ShardedRecycler{T,Type{Factory}}(Factory, cache)

ShardedRecycler{T}(factory::Factory, cache::ShardedCache{T}) where {Factory,T} =
    ShardedRecycler{T,Factory}(factory, cache)

aslimit(limit::Integer) = convert(UInt, limit)
function aslimit(limit::Float64)
    isequal(limit, Inf) && return typemax(UInt)
    throw(ArgumentError("expected `Inf` or a positive integer"))
end

function ShardedCache{T}(; limit::Union{Integer,Float64} = 8) where {T}
    limit = aslimit(limit)
    pools = [T[] for _ in 1:Threads.nthreads()]
    return ShardedCache(limit, Val(T), pools)
end

function ShardedCache{T}(objects; kwargs...) where {T}
    cache = ShardedCache{T}(kwargs...)
    pools = cache.pools
    i = 1
    for x in objects
        push!(pools[i], x)
        length(pools[i]) == cache.limit && break
        i = i == length(pools) ? 1 : i + 1
    end
    return cache
end

function Recyclers.unsafe_init!(cache::ShardedCache)
    pools = cache.pools
    sizehint!(empty!(pools), Threads.nthreads())
    for _ in 1:Threads.nthreads()
        push!(pools, eltype(cache)[])
    end
    return cache
end

function Recyclers.maybeget!(cache::ShardedCache{T}) where {T}
    pool = cache.pools[Threads.threadid()]
    isempty(pool) && return nothing
    return Some{T}(pop!(pool))
end

function Recyclers.recycle!(cache::ShardedCache{T}, object::T) where {T}
    pool = cache.pools[Threads.threadid()]
    if length(pool) < cache.limit
        push!(pool, object)
        return true
    else
        return false
    end
end

function Recyclers.unsafe_empty!(cache::ShardedCache)
    foreach(empty!, cache.pools)
    return cache
end

function Recyclers.unsafe_destroy!(cache::ShardedCache)
    empty!(cache.pools)
    return cache
end


text:
using CSV
using DataFrames
using ScottishTaxBenefitModel
using .Definitions

household_name = "model_households"
people_name = "model_people"

hh_dataset = CSV.File("$(MODEL_DATA_DIR)/$(household_name).tab", delim='\t' ) |> DataFrame
people_dataset = CSV.File("$(MODEL_DATA_DIR)/$(people_name).tab", delim='\t') |> DataFrame

dropmissing!(people_dataset,:data_year) # kill!!! non hbai kids

scottish_hhlds = hh_dataset[(hh_dataset.region .== 299999999),:]
scottish_people = semijoin(people_dataset, scottish_hhlds,on=[:hid,:data_year])

shs_matches = CSV.File( "$(MATCHING_DIR)/mapped_from_shs_vars.tab" )|>DataFrame
#
# mh = innerjoin(scottish_hhlds,shs_matches,on=[:data_year,:hid],makeunique=true)
# @assert size(mh)[1] == size( scottish_hhlds )[1]
# match in previously calculated councils - no need for actual join
#
scottish_hhlds.council = shs_matches.council
scottish_hhlds.nhs_board = shs_matches.nhs_board
scottish_hhlds.bedrooms = shs_matches.bedrooms

CSV.write("$(MODEL_DATA_DIR)model_households_scotland.tab", scottish_hhlds, delim = "\t")
CSV.write("$(MODEL_DATA_DIR)model_people_scotland.tab", scottish_people, delim = "\t")


text:
#"""
#Return a vector of vectors indicating the indices of identical values in a
#`Column` or identical rows of a `Table` using a sorting algorithm (accepting
#keyword arguments of Julia's `sort` interface).
#"""
"""
Return a vector of vectors indicating the indices of identical values in a
`Column` or identical rows of a `Table` using a hashing algorithm.
"""
function groupinds(C::Union{Column,Table}; kwargs...) # TODO this method name might be added to Base in 0.5 (https://github.com/JuliaLang/julia/pull/15503/)
    out = Dict{eltype(C),Vector{Int}}()
    i = 1
    for i = 1:length(C)
        if !in(C.data[i], keys(out))
            out[C.data[i]] = [i]
        else
            push!(out[C.data[i]],i)
        end
    end
    collect(values(out))
end

"""
Return the indices of the first unique elements in a `Column` or a `Table`.
"""
function uniqueind(x::Column)
    out = Vector{Int}()
    seen = Set{eltype(x)}()
    i = 1
    for y in x
        if !in(y, seen)
            push!(seen, y)
            push!(out, i)
        end
        i += 1
    end
    out
end

function uniqueind(x::Table)
    out = Vector{Int}()
    seen = Set{eltype(x)}()
    i = 1
    for y in x
        if !in(y.data, seen)
            push!(seen, y.data)
            push!(out, i)
        end
        i += 1
    end
    out
end

function Base.unique(x::Union{Column,Table})
    x[uniqueind(x)]
end

"""
Mutating form of `unique`.
"""
function unique!(x::Column)
    idx = uniqueind(x)
    x.data[1:length(idx)] == x.data[idx]
    resize!(x.data, length(idx))
end

function unique!(x::Table)
    idx = uniqueind(x)
    for i = 1:ncol(x)
        x.data[i][1:length(idx)] == x.data[i][idx]
        resize!(x.data[i], length(idx))
    end
end

Base.union(cell::Cell) = Column(cell)
Base.union(col::Column) = unique(col)
Base.union{Name}(col1::Union{Cell{Name},Column{Name}}, col2::Union{Cell{Name},Column{Name}}) =
    Column{Name}(union(col1.data, col2.data))
Base.union{Name}(col1::Union{Cell{Name},Column{Name}}, col2::Union{Cell{Name},Column{Name}}, cols::Union{Cell,Column}...) =
    union(union(col1, col2), cols...)

Base.intersect{Name}(col1::Union{Cell{Name},Column{Name}}, col2::Union{Cell{Name},Column{Name}}) =
    Column{Name}(intersect(col1.data, col2.data))

Base.setdiff{Name}(col1::Union{Cell{Name},Column{Name}}, col2::Union{Cell{Name},Column{Name}}) =
    Column{Name}(setdiff(col1.data, col2.data))

Base.union(row::Row) = vcat(row)
Base.union(table::Table) = unique(table)
function Base.union{Names}(table1::Union{Row{Names},Table{Names}}, table2::Union{Row{Names},Table{Names}})
    seen = Set{eltypes(typeof(table1))}()
    idx1 = Vector{Int}()
    for i = 1:length(table1)
        if !in(table1[i].data, seen)
            push!(seen, table1[i].data)
            push!(idx1, i)
        end
    end
    idx2 = Vector{Int}()
    for i = 1:length(table2)
        if !in(table2[i].data, seen)
            push!(seen, table2[i].data)
            push!(idx2, i)
        end
    end

    vcat(table1[idx1], table2[idx2])
end
function Base.union(table1::Union{Row,Table}, table2::Union{Row,Table}, tables::Union{Row,Table}...)
    union(union(table1, table2), tables...) # TODO a bit inefficient to build the hash multiple times...
end

function Base.union!{Names}(table1::Table{Names}, table2::Union{Row{Names},Table{Names}})
    seen = Set{eltypes(typeof(table1))}()
    repeats = Vector{Int}()
    for i = 1:length(table1)
        if !in(table1[i].data, seen)
            push!(seen, table1[i].data)
        else
            push!(repeats, i)
        end
    end
    deleteat!(table1, repeats)
    for i = 1:length(table2)
        if !in(table2[i].data, seen)
            push!(seen, table2[i].data)
            push!(table1, table2[i])
        end
    end
end
function Base.union!(table1::Union{Row,Table}, table2::Union{Row,Table}, tables::Union{Row,Table}...)
    union!(union!(table1, table2), tables...) # TODO a bit inefficient to build the hash multiple times...
end

function Base.intersect{Names}(table1::Union{Row{Names},Table{Names}}, table2::Union{Row{Names},Table{Names}})
    seen1 = Set{eltypes(typeof(table1))}()
    for i = 1:length(table1)
        if !in(table1[i].data, seen1)
            push!(seen1, table1[i].data)
        end
    end
    seen2 = Set{eltypes(typeof(table2))}()
    idx2 = Vector{Int}()
    for i = 1:length(table2)
        if in(table2[i].data, seen1)
            if !in(table2[i].data, seen2)
                push!(seen2, table2[i].data)
                push!(idx2, i)
            end
        end
    end

    table2[idx2]
end

function Base.setdiff{Names}(table1::Union{Row{Names},Table{Names}}, table2::Union{Row{Names},Table{Names}})
    idx1 = Dict{eltypes(typeof(table1)),Int}()
    for i = 1:length(table1)
        if !in(table1[i].data, keys(idx1))
            idx1[table1[i].data] = i
        end
    end
    for i = 1:length(table2)
        if in(table2[i].data, keys(idx1))
            delete!(idx1,table2[i].data)
        end
    end

    table1[collect(values(idx1))]
end


text:
#
# Warn about the order of the ITensor after contractions
#

const default_warn_order = 14

const warn_order =
  Ref{Union{Int, Nothing}}(default_warn_order)

"""
    ITensors.get_warn_order()

Return the threshold for the order of an ITensor above which
ITensors will emit a warning.

You can set the threshold with the function `set_warn_order!(N::Int)`.
"""
get_warn_order() = warn_order[]

"""
    ITensors.set_warn_order(N::Int)

After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the order `N`.

This function returns the initial warning threshold (what it was
set to before this function was called).

You can get the current threshold with the function `ITensors.get_warn_order(N::Int)`. You can reset to the default value with
`ITensors.reset_warn_order()`.
"""
function set_warn_order(N::Union{Int, Nothing})
  N_init = get_warn_order()
  warn_order[] = N
  return N_init
end

"""
    ITensors.reset_warn_order()

After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the default order
$default_warn_order.

This function returns the initial warning threshold (what it was
set to before this function was called).
"""
reset_warn_order() = set_warn_order(default_warn_order)

"""
    ITensors.disable_warn_order()

After this is called, ITensor will not warn about ITensor
contractions that result in large ITensor orders.

This function returns the initial warning threshold (what it was
set to before this function was called).
"""
disable_warn_order() = set_warn_order(nothing)

"""
    @disable_warn_order

Disable warning about the ITensor order in a block of code.

# Examples
```julia
A = ITensor(IndexSet(_ -> Index(1), Order(8)))
B = ITensor(IndexSet(_ -> Index(1), Order(8)))
A * B
@disable_warn_order A * B
@reset_warn_order A * B
@set_warn_order 17 A * B
@set_warn_order 12 A * B
```
"""
macro disable_warn_order(block)
  quote
    local old_order = disable_warn_order()
    r = $(esc(block))
    set_warn_order(old_order)
    r
  end
end

"""
    @set_warn_order

Temporarily set the order threshold for warning about the ITensor
order in a block of code.

# Examples
```julia
@set_warn_order 12 A * B

@set_warn_order 15 begin
  C = A * B
  E = C * D
end
```
"""
macro set_warn_order(new_order, block)
  quote
    local old_order = set_warn_order($(esc(new_order)))
    r = $(esc(block))
    set_warn_order(old_order)
    r
  end
end

"""
    @reset_warn_order

Temporarily sets the order threshold for warning about the ITensor
order in a block of code to the default value $default_warn_order.

# Examples
```julia
@reset_warn_order A * B
```
"""
macro reset_warn_order(block)
  quote
    local old_order = reset_warn_order()
    r = $(esc(block))
    set_warn_order(old_order)
    r
  end
end

#
# A global timer used with TimerOutputs.jl
#

using NDTensors: timer

#
# Get the current number of BLAS threads
# For VERSION >= v"1.6" this will become
# using LinearAlgebra; BLAS.get_num_threads()
#

using NDTensors: blas_get_num_threads

#
# Block sparse multithreading
#

"""
$(NDTensors.enable_threaded_blocksparse_docstring(@__MODULE__))
"""
using_threaded_blocksparse() = NDTensors._using_threaded_blocksparse[]

"""
$(NDTensors.enable_threaded_blocksparse_docstring(@__MODULE__))
"""
enable_threaded_blocksparse() = NDTensors._enable_threaded_blocksparse()

"""
$(NDTensors.enable_threaded_blocksparse_docstring(@__MODULE__))
"""
disable_threaded_blocksparse() = NDTensors._disable_threaded_blocksparse()

#
# Turn enable or disable combining QN ITensors before contracting
#

const _using_combine_contract = Ref(false)

using_combine_contract() = _using_combine_contract[]

function enable_combine_contract()
  _using_combine_contract[] = true
  return nothing
end

function disable_combine_contract()
  _using_combine_contract[] = false
  return nothing
end

#
# Turn debug checks on and off
#

using_debug_checks() = false

macro debug_check(ex)
  quote
    if using_debug_checks()
      $(esc(ex))
    end
  end
end

function enable_debug_checks()
  if !getfield(@__MODULE__, :using_debug_checks)()
    Core.eval(@__MODULE__, :(using_debug_checks() = true))
  end
end

function disable_debug_checks()
  if getfield(@__MODULE__, :using_debug_checks)()
    Core.eval(@__MODULE__, :(using_debug_checks() = false))
  end
end

#
# Turn contraction sequence optimizations on and off
#

const _using_contraction_sequence_optimization = Ref(false)

using_contraction_sequence_optimization() = _using_contraction_sequence_optimization[]

function enable_contraction_sequence_optimization()
  _using_contraction_sequence_optimization[] = true
  return nothing
end

function disable_contraction_sequence_optimization()
  _using_contraction_sequence_optimization[] = false
  return nothing
end


text:
function getHeaders(this::GDAXUser,
                    method::String,
                    request_path::String,
                    body::Dict)::Dict{String, String}
    timestamp = Int64(round(Dates.datetime2unix(now(Dates.UTC))))
    body_str  = isempty(body) ? "" : JSON.json(body)
    prehash_signature = string(timestamp, uppercase(method), request_path, body_str)
    h = HMACState("sha256", base64decode(this.secret))
    Nettle.update!(h, prehash_signature)
    signature = base64encode(Nettle.digest!(h))

    return Dict("CB-ACCESS-KEY" => this.key,
                    "CB-ACCESS-SIGN" => signature,
                    "CB-ACCESS-TIMESTAMP" => string(timestamp),
                    "CB-ACCESS-PASSPHRASE" => this.passphrase,
                    "Content-Type" => "application/json")
end

function send_request(this::GDAXUser,
                        method::String,
                        request_path::String;
                        query = Dict(),
                        json = Dict(),
                        timeout = 1.0)

    query_str = Requests.format_query_str(query)
    query_str = "?" * query_str
    while !isempty(query_str) && last(query_str) == '?'
        query_str = chop(query_str)
    end

    handle = eval(parse("Requests." * lowercase(method)))

    return handle(this.rest_url * request_path,
                 query = query,
                 json = isempty(json) ? nothing : json,
                 timeout = timeout,
                 headers = getHeaders(this, method, request_path * query_str, json))
end

getProducts(this::GDAXUser) = send_request(this, "get", "/products")
getServerTime(this::GDAXUser) = send_request(this, "get", "/time")
getProductOrderBook(this::GDAXUser, product_id::String; level::Int64 = 2) = send_request(this, "get", "/products/" * product_id * "/book", query = Dict("level" => level))
getProductTicker(this::GDAXUser, product_id::String) = send_request(this, "get", "/products/" * product_id * "/ticker")
getProductTrades(this::GDAXUser, product_id::String) = send_request(this, "get", "/products/" * product_id * "/trades")
get24HourStats(this::GDAXUser, product_id::String) = send_request(this, "get", "/products/" * product_id * "/stats")
getCurrencies(this::GDAXUser) = send_request(this, "get", "/currencies")
getPosition(this::GDAXUser) = send_request(this, "get", "/position")
getAccounts(this::GDAXUser) = send_request(this, "get", "/accounts")
getAccountHistory(this::GDAXUser, account_id::String) = send_request(this, "get", "/accounts/" * account_id * "/ledger")
getAccountHolds(this::GDAXUser, account_id::String) = send_request(this, "get", "/accounts/" * account_id * "/holds")

"""
Response Items

Each bucket is an array of the following information:

    * time bucket start time
    * low lowest price during the bucket interval
    * high highest price during the bucket interval
    * open opening price (first trade) in the bucket interval
    * close closing price (last trade) in the bucket interval
    * volume volume of trading activity during the bucket interval
"""
function getHistoricRates(this::GDAXUser,
                            product_id::String,
                            from::DateTime,
                            to::DateTime,
                            granularity::Int64)
    send_request(this, "get", "/products/" * product_id * "/candles", query = Dict("start" => from,
                                                                  "end" => to,
                                                                  "granularity" => granularity))
end

"""
acc_type    fill or account
start_date 	Starting date for the report (inclusive)
end_date 	Ending date for the report (inclusive)
product_id 	ID of the product to generate a fills report for. E.g. BTC-USD. Required if type is fills
account_id 	ID of the account to generate an account report for. Required if type is account
format 	pdf or csv (defualt is pdf)
email 	Email address to send the report to (optional)
"""
function getReport(this::GDAXUser,
                    report_type::String,
                    start_date::DateTime,
                    end_date::DateTime,
                    product_id::String,
                    account_id::String,
                    format::String,
                    email::String)

    request_path = "/reports"
    data = Dict("type" => report_type,
                "start_date" => Dates.format(start_date, "yyyy-mm-ddTHH:MM:SS\\Z"),
                "end_date" => Dates.format(end_date, "yyyy-mm-ddTHH:MM:SS\\Z"),
                "product_id" => product_id,
                "account_id" => account_id,
                "format" => format,
                "email" => email)

    return send_request(this, "post", request_path, json = data)
end

function placeOrder(this::GDAXUser,
                    side::String,
                    product_id::String,
                    lots::Float64,
                    price::Float64;
                    post_only::Bool = true,
                    stp::String = "cb",
                    time_in_force::String = "GTC")
    order = Dict{String, Any}("side" => side,
                              "product_id" => product_id,
                              "price" => string(signif(price, 6)),
                              "size" => string(signif(lots, 8)),
                              "post_only" => post_only,
                              "client_oid" => string(Base.Random.uuid4()),
                              "type" => "limit",
                              "stp" => stp,
                              "time_in_force" => time_in_force)

    send_request(this, "post", "/orders", json = order)
end

function cancelOrder(this::GDAXUser, order_id::String)
    send_request(this, "delete", "/orders/" * order_id)
end

function cancelAll(this::GDAXUser)
    send_request(this, "delete", "/orders")
end

function listOrders(this::GDAXUser; product_id::String = "", status::String = "all")
    query = Dict("status" => status)
    !isempty(product_id) && (query["product_id"] = product_id)
    return send_request(this, "get", "/orders", query = query)
end

max_stars_repo_name: DatName/gdax-julia

max_issues_repo_name: DatName/gdax-julia

max_forks_repo_name: DatName/gdax-julia


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Elefridge
using JLD

path = "/Users/milan/cams/unstructured"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)
# ndims = ["lon","lat","vert"]
ndims = ["vert","hori"]
nbits = 32

IC = fill(0.0,n,length(ndims),nbits)

for (i,file) in enumerate(filelist)
    varname = split(split(file,"_")[end],".")[1]
    varnames[i] = varname
    println("---")
    println("Reading $varname")

    # vert x lat x lon
    X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    # convert biased exponent to signed exponent
    signed_exponent!(X)

    # permute dimensions to have lon x lat x vert, lat x lon x vert, and original
    # IC[i,1,:] = bitinformation(permutedims(X,[3,2,1]))
    # IC[i,2,:] = bitinformation(permutedims(X,[2,3,1]))
    # IC[i,3,:] = bitinformation(X)

    IC[i,1,:] = bitinformation(X)
    IC[i,2,:] = bitinformation(copy(X'))

    @save "/Users/milan/cams/entropy/information_all_unstructured.jld" varnames IC
end


text:
@testset "Fibers" begin
    @polyvar x y z
    f = x^2+y^2+1

    AF = AmoebaFiber2D(f)
    PolynomialAmoebas.update_fiber!(AF, (0.0, 0.0))

    @test_throws AssertionError AmoebaFiber2D(x^2+1)
    @test_throws AssertionError AmoebaFiber2D(x^2+z+y)

    CF = CoamoebaFiber2D(f)
    PolynomialAmoebas.update_fiber!(CF, (0.0, 0.0))
    @test_throws AssertionError CoamoebaFiber2D(x^2+1)
    @test_throws AssertionError CoamoebaFiber2D(x^2+z+y)

    CF2 = ContourFiber2D(f)
    PolynomialAmoebas.update_fiber!(CF2, 0.0)
    PolynomialAmoebas.fix_axis!(CF2, :w2)
    PolynomialAmoebas.update_fiber!(CF2, 1.0)
    @test_throws AssertionError ContourFiber2D(x^2+1)
    @test_throws AssertionError ContourFiber2D(x^2+z+y)

    I = ImaginaryFiber2D(f)
    @test_throws AssertionError ImaginaryFiber2D(x^2+1)

    g = x^2+y^2+z^2+1
    AG = AmoebaFiber3D(g)
    PolynomialAmoebas.update_fiber!(AG, (0.0, 0.0, 0.0))
    @test_throws AssertionError AmoebaFiber3D(x^2+1)
    @test_throws AssertionError AmoebaFiber3D(x^2+z)

    CG = CoamoebaFiber3D(g)
    PolynomialAmoebas.update_fiber!(CG, (0.0, 0.0, 1.0))
    @test_throws AssertionError CoamoebaFiber3D(x^2+1)
    @test_throws AssertionError CoamoebaFiber3D(x^2+z)

    IG = ImaginaryFiber3D(g)
    @test_throws AssertionError ImaginaryFiber3D(x^2+1)
end