text:
# Title:  Variable Capture
# Artist: Lyndon White (@oxinabox)
# Format: Literate.jl script
# License: MIT
# ---

# # What are we doing?
#
# We would like to capture all the local variables from methods we call into the context's metadata.

using Cassette

Cassette.@context CaptureCtx

# helpers for getting  the name of a slot given its slot number
slotname(ir::Core.CodeInfo, slotnum::Integer) = ir.slotnames[slotnum]
slotname(ir::Core.CodeInfo, slotnum) = slotname(ir, slotnum.id)

# inserts insert `ctx.metadata[:x] = x`
function record_slot_value(ir, metadata_slot, slotnum)
    name = slotname(ir, slotnum)

    #call_ast(:(Base.setindex!(dest_slot, lhs, QuoteNode(name))))
    return Expr(
        :call,
        Expr(:nooverdub, GlobalRef(Base, :setindex!)),
        metadata_slot,
        slotnum,
        QuoteNode(name)
    )
end

# After every assigment: `x = foo`, insert `ctx.metadata[:x] = x`
function instrument_assignments!(ir, metadata_slot)
    is_assignment(stmt) = Base.Meta.isexpr(stmt, :(=))
    stmtcount(stmt, i) = is_assignment(stmt) ? 2 : nothing
    function newstmts(stmt, i)
        lhs = stmt.args[1]
        record = record_slot_value(ir, metadata_slot, lhs)
        return [stmt, record]
    end
    Cassette.insert_statements!(ir.code, ir.codelocs, stmtcount, newstmts)
end

# At the start of the method we need to capture all  the argument values
function instrument_arguments!(ir, method, metadata_slot)
    # start from 2 to skip #self
    arg_names = Base.method_argnames(method)[2:end]
    arg_slots = 1 .+ (1:length(arg_names))
    @assert(
        ir.slotnames[arg_slots] == arg_names,
        "$(ir.slotnames[arg_slots]) != $(arg_names)"
    )

    Cassette.insert_statements!(
        ir.code, ir.codelocs,
        # Only add to line 1 as doing this at the begining
        # add 1 because we need to put the original line 1 back in
        (stmt, i) -> i == 1 ? length(arg_slots) + 1 : nothing,

        (stmt, i) -> [
            map(arg_slots) do slotnum
                slot = Core.SlotNumber(slotnum)
                record_slot_value(ir, metadata_slot, slot)
            end;
            stmt
        ]
    )
end

"""
     create_slot!(ir, namebase="")
Adds a slot to the IR with a name based on `namebase`.
It will be added at the end
returns the new slot number
"""
function create_slot!(ir, namebase="")
    slot_name = gensym(namebase)
    push!(ir.slotnames, slot_name)
    push!(ir.slotflags, 0x00)
    slot = Core.SlotNumber(length(ir.slotnames))
    return slot
end

"""
    put_metadata_in_its_slot!(ir, metadata_slot)

Attaches the cassette metadata object to the slot given.
This will be added as the very first statement to the ir.
"""
function put_metadata_in_its_slot!(ir, metadata_slot)
    getmetadata = Expr(
        :call,
        Expr(
            :nooverdub,
            GlobalRef(Core, :getfield)
        ),
        Expr(:contextslot),
        QuoteNode(:metadata)
    )

    Cassette.insert_statements!(
        ir.code, ir.codelocs,
        (stmt, i) -> i == 1 ? 2 : nothing,
        (stmt, i) -> [Expr(:(=), metadata_slot, getmetadata), stmt]
    )
end

# ## The main method of our actual pass
# This has to actually setup the capturing of the info we need and make sure it goes into metadata
function instrument_variables!(::Type{<:CaptureCtx}, reflection::Cassette.Reflection)
    ir::Core.CodeInfo = reflection.code_info
    # Create a slot that we will put the metadata in, make it the last one
    metadata_slot = create_slot!(ir, "metadata")

    # Now the real part where we determine about assigments
    instrument_assignments!(ir, metadata_slot)

    # record all the initial values so we get the parameters
    instrument_arguments!(ir, reflection.method, metadata_slot)

    # insert the initial `metadata slot` assignment into the IR.
    # Do this last so it doesn't get caught in our assignment catching
    put_metadata_in_its_slot!(ir, metadata_slot)

    return ir
end

const variable_capture_pass = Cassette.@pass instrument_variables!

# ## Define what to do when our function calls a function.
# We are going to add that dictionary as an  element of our dictionary
function Cassette.overdub(ctx::CaptureCtx, f, args...)
    if Cassette.canrecurse(ctx, f, args...)
        _ctx = Cassette.similarcontext(ctx, metadata = Dict())
        try
            return Cassette.recurse(_ctx, f, args...)
        finally
            # If the called function had any local variables save them
            # into a sub-dictionary
            if length(_ctx.metadata) > 0
                ctx.metadata[gensym(Symbol(f))] = _ctx.metadata
            end
        end
    else
        return Cassette.fallback(ctx, f, args...)
    end
end

# ## Hook it all up

function capture_variables(f, args...)
    ctx = CaptureCtx(;metadata=Dict(), pass=variable_capture_pass)
    Cassette.recurse(ctx, f, args...)
    return ctx.metadata
end

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

# # Demo
using Test
@testset "Basic local variable capture" begin
    function foo(a)
        b = a+10
        c = b+10
        return (a,b,c)
    end

    vars = capture_variables(foo, 10)

    @testset "Assignments" begin
        @test vars[:b] == 20
        @test vars[:c] == 30
    end

    @testset "arguments" begin
        @test vars[:a] == 10
    end

    display(vars)
end

# ### Should output something like:
# ```julia
# Dict{Any,Any} with 5 entries:
#   :a                => 10
#   :b                => 20
#   :c                => 30
#   Symbol("##+#379") => Dict{Any,Any}(:y=>10,:x=>10)
#   Symbol("##+#380") => Dict{Any,Any}(:y=>10,:x=>20)\
# ```


text:
abstract AbstractPoint{T <: Real}

immutable Point{T} <: AbstractPoint
	x::T
	y::T
end

x(p::Point) = p.x
y(p::Point) = p.y
Base.eltype{T}(::Type{Point{T}}) = T

Base.start(p::Point) = 1
Base.next(p::Point, state) = (getfield(p,fieldname(Point,state)), state + 1)
Base.done(p::Point, state) = state > nfields(Point)
Base.length(p::Point) = nfields(Point)
Base.size(p::Point) = (nfields(Point),)

Base.getindex(p::Point,i) = getfield(p,fieldname(Point,i))
Base.endof(p::Point)  = length(p)

immutable Line{T <: AbstractPoint}
	p1::T
	p2::T
end

function perpendicular(l1::Line, l2::Line)
	m1 = (y(l1.p2) - y(l1.p1)) / (x(l1.p2) - x(l1.p1))
	m2 = (y(l2.p2) - y(l2.p1)) / (x(l2.p2) - x(l2.p1))

	return isapprox(m1,-(1/m2))
end

⟂(l1,l2) = perpendicular(l1,l2)

isleft(p::AbstractPoint, l::Line) = ((x(l.p2) - x(l.p1)) * (y(p) - y(l.p1) ) - (x(p) - x(l.p1)) * (y(l.p2) - y(l.p1)))

Base.norm(l::Line) = norm([x(l.p2) - x(l.p1), y(l.p2) - y(l.p1)])

abstract AbstractPolygon{T <: AbstractPoint}

immutable Polygon{T} <: AbstractPolygon{T}
	vertexs::Vector{T}
	function Polygon(v::Vector)
		if v[end] != v[1]
			push!(v,v[1])
		end
		return new(v)
	end
end

Polygon{T}(v::Vector{T}) = Polygon{T}(v)

Base.start(p::Polygon) = start(p.vertexs)
Base.next(p::Polygon, state) = next(p.vertexs, state)
Base.done(p::Polygon, state) = done(p.vertexs, state)
Base.length(p::Polygon) = length(p.vertexs)
Base.size(p::Polygon) = size(p.vertexs)

Base.getindex(p::Polygon,i) = p.vertexs[i]
Base.endof(p::Polygon)  = endof(p.vertexs)
Base.eltype(p::Polygon) = eltype(p.vertexs)

function boundingbox(poly::AbstractPolygon)
	minx = mapreduce(x,min,poly.vertexs)
	miny = mapreduce(y,min,poly.vertexs)
	maxx = mapreduce(x,max,poly.vertexs)
	maxy = mapreduce(y,max,poly.vertexs)

	(Point(minx,miny), Point(maxx,maxy))
end

function Base.in(p::AbstractPoint, poly::AbstractPolygon)
	wn = 0
	for idx in 1:(length(poly) - 1)
		V1 = poly[idx]
		V2 = poly[idx + 1]

		if y(V1) <= y(p)
			if y(V2) > y(p)
				if isleft(p, Line(V1,V2)) > 0
					wn+=1
				end
			end
		else
			if y(V2) <= y(p)
				if isleft(p, Line(V1,V2)) < 0
					wn-=1
				end
			end
		end
	end
	return wn != 0
end

function Base.rand{T}(rng, S::AbstractPolygon{T}, n::Int64)
	pT = eltype(T)
	ret = Point{pT}[]

	(minp,maxp) = boundingbox(S)

	bw = x(maxp) - x(minp)
	bh = y(maxp) - y(minp)

	while length(ret) < n
		rnd = [x(minp),y(minp)] .+ (rand(rng,pT,2) .* [bw,bh])

		newPoint = Point{pT}(rnd...)

		if newPoint in S
			push!(ret,newPoint)
		end
	end
	ret
end

Base.rand(rng,S::AbstractPolygon) = rand(rng,S,1)

Base.rand(S::AbstractPolygon, n) = rand(Base.Random.GLOBAL_RNG,S,n)

Base.rand(S::AbstractPolygon) = rand(Base.Random.GLOBAL_RNG,S,1)[1]

ext: jl


text:
module TestVectors

using BangBang: push!!
using MicroCollections
using Test

const EMPTY_ITR = (UNKNOWN for _ in 1:0)
@assert Base.IteratorEltype(EMPTY_ITR) isa Base.EltypeUnknown

const UNSTABLE_ONE = Ref{Any}(1)
const SINGLETON_ITR = (UNSTABLE_ONE[] for _ in 1:1)
@assert Base.IteratorEltype(SINGLETON_ITR) isa Base.EltypeUnknown

@testset "empty" begin
    @test EmptyVector() === emptyshim(Vector)
    @test eltype(EmptyVector()) === Union{}
    @test EmptyVector(Int[]) === emptyshim(Vector, Int)::AbstractVector{Int}
    @test EmptyVector(()) === emptyshim(Vector)
    @test EmptyVector{Symbol}(Int[]) === emptyshim(Vector, Symbol)
    @test EmptyVector{Symbol}(EMPTY_ITR) === emptyshim(Vector, Symbol)
    @test push!!(EmptyVector(), 0) == [0]
end

@testset "singleton" begin
    @test SingletonVector([0]) === singletonshim(Vector, 0)
    @test SingletonVector([0])::AbstractVector{Int} == [0]
    @test SingletonVector([nothing])::AbstractVector{Nothing} == [nothing]
    @test vcat(SingletonVector([0]), [0.5])::Vector{Float64} == [0.0, 0.5]
    @test SingletonVector(SINGLETON_ITR) === singletonshim(Vector, 1)
    @test SingletonVector{Float64}(SINGLETON_ITR) === singletonshim(Vector, 1.0)
    @test SingletonVector(Real[0])::AbstractVector{Real} == [0]
    @test SingletonVector{Int}(Real[0])::AbstractVector{Int} == [0]
    @test SingletonVector((0,))::AbstractVector{Int} == [0]
    @test push!!(SingletonVector([0]), 0.5) == [0.0, 0.5]
end

@testset "singleton BitVector" begin
    ys = singletonshim(BitVector, true)
    @test ys[1] === true
    @test vcat(ys, [false]) == [true, false]
    @test vcat(ys, [2]) == [1, 2]
    @test_throws InexactError push!!(singletonshim(BitVector, 2), false)
end

end  # module


text:
function tiddly_greet()
    return("tiddlywinks is working")
end


text:
function sortedintersect(A::Array, B::Array)
    ABintersect = typeof(A)()
    temp_index_A = 1
    temp_index_B = 1

    while (temp_index_A <= length(A)) & (temp_index_B <= length(B))
        if A[temp_index_A] == B[temp_index_B]
            push!(ABintersect,A[temp_index_A])
            temp_index_A += 1
            temp_index_B += 1
        elseif A[temp_index_A] < B[temp_index_B]
            temp_index_A += 1
        else
            temp_index_B += 1
        end
    end

    return ABintersect
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using ies
using Test

@testset "ies.jl" begin
    # Write your own tests here.
end


text:
using Gtk

const win = GtkWindow("Get Mouse Position", 600, 800)
const butn = GtkButton("Click Me Somewhere")
push!(win, butn)

callback(b, evt) = set_gtk_property!(win, :title, "Mouse Position: X is $(evt.x), Y is $(evt.y)")
signal_connect(callback, butn, "button-press-event")

showall(win)

c = Condition()
endit(w) = notify(c)
signal_connect(endit, win, :destroy)
wait(c)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"A location in 3D space"
Base.@kwdef struct Coordinate
    x::Float64
    y::Float64
    z::Float64
end

"Parametrizes a scene"
Base.@kwdef struct Video_Params
    p_objects::Float64 = 0.9
    n_possible_objects = 8
    possible_objects::Vector{Int64} = collect(1:n_possible_objects)
    probs_possible_objects = collect(ones(n_possible_objects)./n_possible_objects)
    x_min::Float64 = -5
    y_min::Float64 = 0
    z_min::Float64 = -5
    x_max::Float64 = 5
    y_max::Float64 = 3
    z_max::Float64 = 5
    delta = 1 #sd for gaussian #midpoint of logistic curve for prior over distance between objects
    #k = 10. #slope of logistic curve for prior over distance between objects
    num_receptive_fields::Int64 = 1
    #camera parameters that don't change
    image_dim_x::Int64 = 256
    image_dim_y::Int64 = 256
    horizontal_FoV::Float64 = 55
    vertical_FoV::Float64 = 55
end

"Parametrizes the state of the camera, which changes within a scene"
Base.@kwdef struct Camera_Params
    camera_location::Coordinate
    camera_focus::Coordinate
end

Base.@kwdef struct Receptive_Field
    "upper left"
    p1::Tuple{Int64, Int64}
    "lower right"
    p2::Tuple{Int64, Int64}
end

"""
Parametrizes a line segment in 3d space.

``L = \\mathrm{startpoint} + (a, b, c)t``, where ``t\\in [0, 1]``.
``\\mathrm{endpoint} = \\mathrm{startpoint} + (a, b, c)``.
"""
Base.@kwdef struct Line_Segment
    start::Coordinate
    endpoint::Coordinate
    a::Float64
    b::Float64
    c::Float64
end

"""
Stores a probability distribution over the object categories.

Used for defining the proposal function.
"""
Base.@kwdef struct Perturb_Params
    probs_possible_objects::Vector{Float64}
end

export Coordinate
export Video_Params
export Camera_Params
export Permanent_Camera_Params
export Receptive_Field
export Line_Segement
export Vector
export Perturb_Params

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
include("test.jl")
include("test_showvalues.jl")

println("")
println("All tests complete")


text:
export iid

iid(n::Integer...) = dist -> iid(dist, n...)

iid(dist, n...) = dist ^ n

iid(dist::Dists.Distribution, n...) = MeasureBase.powermeasure(dist, n)


text:
export WeightedMeasure

"""
    struct WeightedMeasure{R,M} <: AbstractMeasure
        logweight :: R
        base :: M
    end

"""
struct WeightedMeasure{R,M} <: AbstractMeasure
    logweight :: R
    base :: M
end

logweight(μ::WeightedMeasure) = μ.logweight

function Base.show(io::IO, μ::WeightedMeasure)
    io = IOContext(io, :compact => true)
    print(io, exp(μ.logweight), " * ", μ.base)
end

function Base.show_unquoted(io::IO, μ::WeightedMeasure, indent::Int, prec::Int)
    io = IOContext(io, :compact => true)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

function logdensity(sm::WeightedMeasure{R,M}, x::X) where {X, R, M <: AbstractMeasure}
    logdensity(sm.base, x) + sm.logweight
end

function Base.:*(k::T, m::AbstractMeasure) where {T <: Number}
    logk = log(k)
    return WeightedMeasure{typeof(logk), typeof(m)}(logk,m)
end

Base.:*(m::AbstractMeasure, k::Real) = k * m

≪(::M, ::WeightedMeasure{R,M}) where {R,M} = true
≪(::WeightedMeasure{R,M}, ::M) where {R,M} = true

basemeasure(μ::WeightedMeasure) = μ.base

TV.as(μ::WeightedMeasure) = TV.as(μ.base)

sampletype(μ:: WeightedMeasure) = sampletype(μ.base)


text:
using PartedArrays
using StaticArrays
using ForwardDiff
using BenchmarkTools
using Rotations
using MeshCatMechanisms
using MeshCat
using CoordinateTransformations
using LinearAlgebra
include("doggo_model.jl")

vis = Visualizer()
open(vis)
setobject!(vis,Triad(1))

# Position of leg frame
d1 = [2,-1,0]
d2 = [2,1,0]
d3 = [-2,-1,0]
d4 = [-2,1,0]
d_legs = (d1,d2,d3,d4)

# Orientation of leg frames
q1 = Quaternion(RotX(pi/2))
q2 = Quaternion(RotZ(pi)*RotX(pi/2))
q3 = q1
q4 = q2
q_legs = (q1,q2,q3,q4)
q1
T_legs = [compose(Translation(d_legs[i]), LinearMap(Quat(q_legs[i]))) for i = 1:4]
settransform!(vis,T_legs[1])
T_legs

# Doggo
include("doggo_default.jl")

# Inputs
r0 = @SVector [0.1, 0.2, 0.3]
q0 = Quaternion(RotY(-pi/2))

q = [deg2rad(45), deg2rad(-45),
     deg2rad(45), deg2rad(-45),
     deg2rad(45), deg2rad(-45),
     deg2rad(45), deg2rad(-45)]

# Forward Kinematics
function leg2Dto3D(r)
    z = zeros(7)
    z[1:2] = r[1:2]
    q = Quat(RotZ(r[3]))
    z[4] = q.w
    z[5] = q.x
    z[6] = q.y
    z[7] = q.z
    return z
end

function transform_leg(q)
    r = leg2Dto3D.(fk(leg, q[1:2]))
    pos = [q_legs[1]*r[i][1:3] + d_legs[1] for i = 1:4]
    rot = [q_legs[1]*Quat(r[i][4:7]...) for i = 1:4]
    return pos, rot
end

function leg2Dto3D_2(r)
    T = Translation(r[1],r[2],0) ∘ LinearMap(RotZ(r[3]))
end
function transform_leg2(q)
    r = leg2Dto3D_2.(fk(leg, q[1:2]))
    T = [T_legs[1] ∘ r[k] for k = 1:4]
    rot = [trans.linear for trans in T]
    pos = [trans.translation for trans in T]
    T
end
T1 = transform_leg(q)
T2 = transform_leg2(q)
T1 .≈ T2
# @btime transform_leg($q)
# @btime transform_leg2($q)

# First Leg, First Link
function fk11(q)
    leg_num = 1
    link_num = 3

    # Leg kinematics
    q_2d = fk(leg, q[8:9])[link_num]

    # Convert to 3D coordinate frame
    t11 = leg2Dto3D_2(q_2d)

    # Move to body frame
    tb1 = T_legs[leg_num] ∘ t11

    # Create world transformation
    r0 = q[1:3]
    q0 = q[4:7]
    twb = Translation(r0) ∘ LinearMap(Quat(q0...,true))

    # Transform to World
    tw1 = twb ∘ tb1
    [tw1.translation; SVector(Quaternion(tw1.linear))]
end

function fk11_(q)
    leg_num = 1
    link_num = 3

    # Leg kinematics
    q_2d = fk(leg, q[8:9])[link_num]

    # Convert to 3D coordinate frame
    t11 = leg2Dto3D_2(q_2d)

    # Move to body frame
    tb1 = T_legs[leg_num] ∘ t11

    # Create world transformation
    r0 = q[1:3]
    q0 = Quaternion(q[4:7])

    # Transform to World
    r = q0*tb1.translation + r0
    q = q0*Quaternion(tb1.linear)
    [r; SVector(q)]
end

function link_jacobian(q)
    leg_num = 1
    link_num = 3

    # Get leg kinematics
    r = fk(leg, q[8:9])[link_num]
    th = r[3]

    r0 = q[1:3]
    q0 = Quaternion(q[4:7])
    qw1 = q0*q_legs[leg_num]

    # Get translation from link to body frame
    t11 = leg2Dto3D_2(r)
    tb1 = T_legs[leg_num] ∘ t11
    rb1 = tb1.translation

    # Initialize Jacobian
    jac = zeros(7,9)

    # Jacobian of position wrt body
    jac[1:3,1:3] = Diagonal(I,3)
    jac[1:3,4:7] = deriv_conj(q0, rb1)

    # Get leg Jacobian and convert to 3D
    J1 = jacobian(leg, q[8:9])[link_num]
    dqdth = repeat([-sin(th/2)/2, 0, 0, cos(th/2)/2],1,2)
    J3d = [J1[1:2,:]; zeros(1,2); dqdth.*J1[3:3,:]]
    jac[1:3,8:9] = conj(qw1)*J3d[1:3,:]

    # Jacobian of orientation
    jac[4:7,4:7] = Rmult(Quaternion(tb1.linear))
    jac[4:7,8:9] = Lmult(qw1)*J3d[4:7,:]
    return jac
end

T_legs

r0 = @SVector [0, 0, 0.0]
q0 = Quaternion(RotY(-pi/2))
q0 = Quaternion(ones(Quat))
q = [r0; SVector(q0); deg2rad(40); deg2rad(-45)]
t11 = fk11(q)
T = Translation(t11[1:3]) ∘ LinearMap(Quat(t11[4:7]...))
settransform!(vis, T)
fk11_(q) ≈ t11

jac = link_jacobian(q)
jac_fd = ForwardDiff.jacobian(fk11,q)
ForwardDiff.jacobian(fk11_,q) ≈ jac

@code_warntype link_jacobian(q)

@btime ForwardDiff.jacobian($fk11_,$q)
@btime link_jacobian($q)

settransform!(vis,t11)

# World Transformation
Tw = [Twb ∘ T for T in T2]

settransform!(vis,T_legs[1])

# Test Quaternion stuff
q1 = Quat(q1)
q2 = Quat(q2)
q1_ = Quaternion(q1)
q2_ = Quaternion(q2)
q1_*q2_
Lmult(q1_)*q2_ == Rmult(q2_)*q1_

a = rand(3)
q1_*Quaternion(0.,a)*inv(q1_) == q1_*a
(Lmult(q1_)*Rmult(inv(q1_))*Quaternion(0.,a))[2:4] ≈ q1_*a
(Rmult(inv(q1_))*Lmult(q1_)*Quaternion(0.,a))[2:4] ≈ q1_*a
conj(q1_)*a ≈ q1_*a

r = a
s,v = scalar(q1_), vec(q1_)
j2 = 2(s*I + skew(v))*r
j3 = v*r' + I*v'r - 2s*skew(r) - skew(v)*skew(r) - skew(skew(v)*r)
jac = [zeros(1,4); j2 j3]

q = rand(Quat)
q = [q.w, q.x, q.y, q.z]

myfun(q) = Quaternion(q)*a
ForwardDiff.jacobian(myfun,q)
deriv_conj(Quaternion(q),a)

Quat(q...)*a ≈ Quaternion(q)*a
myfun2(q) = Quat(q[1],q[2],q[3],q[4])*a
ForwardDiff.jacobian(myfun2,q)

a = @SVector rand(3)
a = rand(3)
@btime $q1_*$a
@btime $q1*$a
q1 = Quat(q1)
Quat(q1)*a
Quat(q1)*Quat(q2)


text:
# -*- coding: utf-8 -*-
adparam = SNN.ADEXParameter(;a = 6.0,
    b = 7.3,
    cm = 803.1,
    v_rest= -63.2,
    tau_m=19.7,
    tau_w=351.0,
    v_thresh=-39.2,
    delta_T=6.4,
    v_spike = -36.6,
    v_reset = -59.2,
    spike_delta = 16.3)

E = SNN.AD(;N = 1, param=adparam)
E.I = [795.57128906]
SNN.monitor(E, [:v,:I])
SNN.sim!([E],[], dt=0.1ms, duration=2000ms)
# # Julia SNN Implementation of AdExp Neuron.
# [Adaptive_exponential_integrate and fire neuron](http://www.scholarpedia.org/article/Adaptive_exponential_integrate-and-fire_model)
# Dr. Wulfram Gerstner
# Romain Brette, Ecole Normale Supérieure, Paris, France

max_stars_repo_name: russelljjarvis/backupSNN

max_issues_repo_name: russelljjarvis/backupSNN

max_forks_repo_name: russelljjarvis/backupSNN


text:
# Aiyagari Model Test
# Michele Fornino
#---------------------------------------------------#
#					Aiyagari MODEL 					#
#---------------------------------------------------#

"""
	Aiyagari_params{F<:Real, I<:Integer}

	Composite type containing all parameters of the standard Aiyagari Model.

"""
@with_kw struct Aiyagari_params{F<:Real, I<:Integer} @deftype F

	# Parameters
	σ = 1.0
	w = 1.0
	r = log(1 + 0.01)
	ρ = r
	σz = 0.2
	θz = 0.1
	zd = 1.0

	# Derived Parameters
	γ = 1.0 / σ

	# Grids Parameters
	Na::I = 50
	Nz::I = 20
	zGrid_lower_quantile = 1e-2
	zGrid_upper_quantile = 1 - 1e-2
	aGrid_lower = 0.
	aGrid_upper = 20.

	# Assertions
	@assert σ > 0
	@assert r > 0
	@assert w > 0
	@assert ρ > 0
	@assert θz > 0
	@assert σz > 0
	@assert zd > 0

	# Check values for grids make sense
	@assert Na > 1
	@assert Nz > 0
	@assert zGrid_upper_quantile > zGrid_lower_quantile
	@assert zGrid_lower_quantile > 0.0
	@assert zGrid_upper_quantile < 1.0
	@assert aGrid_upper > aGrid_lower

	# Enforce definitions (prevents accidental creation of object with wrong derived parameters)
	@assert γ == 1 / σ

end

"""
	function HamiltonJacobiBellmanEquations.HJBEquation(params::Aiyagari_params)

Overloaded constructor for input of type Aiyagari_params.

"""
function HamiltonJacobiBellmanEquations.HJBEquation(params::Aiyagari_params)

	# Use UnPack.jl to localize all members of params
	@unpack_Aiyagari_params params

	# Utility Function and derivatives.
	function u(c::T)::T where T<:Real
		if σ == 1.0
			return log(c)
		else
			return (c ^ (1 - γ) - 1) / (1 - γ)
		end
	end

	function u′(c::T)::T where T<:Real
		return c ^ -γ
	end

	function u′⁻¹(x::T)::T where T<:Real
		return x ^ (- 1 / γ)
	end

	# function u(c::T)::T where T<:Real
	# 	return - 1 / γ * exp(-γ * c)
	# end

	# function u′(c::T)::T where T<:Real
	# 	return exp( - γ * c)
	# end

	# function u′⁻¹(x::T)::T where T<:Real
	# 	return - 1 / γ * log(x)
	# end

	# Returns Function
	function flow_returns(state::AbstractVector{T}, c::Union{T, AbstractVector{T}})::T where T<:Real
		return u(c[1])
	end

	# Policy Function in term of marginal value function
	function optimal_controls(state::AbstractVector{T}, ∂V∂a::Union{T, AbstractVector{T}})::T where T<:Real
		return u′⁻¹(∂V∂a[1])
	end

	# State equations
	function state_equations(state::AbstractVector{T}, c::Union{T, AbstractVector{T}})::T where T<:Real
		return w * state[2] + r * state[1] - c[1]
	end

	# State constraint
	function state_constraints(state::AbstractVector{T})::T where T<:Real
		return u′(w * state[2] + r * state[1])
	end

	# Tuple of endogenous state grids
	aGrid = Vector(range(aGrid_lower, aGrid_upper, length = Na))
	state_grids_endog = (aGrid, )

	# Exogenous shock
	exog_shock = ContinuousTimeMarkovChain(params)

	# Allocate hjb object
	return HJBEquation(flow_returns,
					   state_equations,
					   optimal_controls,
					   state_constraints,
					   state_grids_endog,
					   exog_shock,
					   ρ)

end

"""
	function MarkovChains.ContinuousTimeMarkovChain(params::Aiyagari_params)

Overloaded constructor for input of type Aiyagari_params.

"""
function MarkovChains.ContinuousTimeMarkovChain(params::Aiyagari_params)

	@unpack_Aiyagari_params params

	# Construct exog_shock ContinuousTimeMarkovChain object (if problem is stochastic).
	if Nz > 1

		# Obtain discretize diffusion process object as a MarkovChains.ContinuousTimeMarkovChain.
		μ = z -> -θz .* (z .- zd)
		σ = z -> sqrt(2 * θz * σz^2 / zd) .* sqrt.(z)
		dp = ItoDiffusionProcess(μ, σ)

		# Theoretical benchmark for CIR stationary distribution, used for zGrid.
		shape = zd^2 / σz^2
		rate = zd / shape
		Fz_theory = Gamma(shape, rate)
		zGrid_lower = quantile(Fz_theory, zGrid_lower_quantile)
		zGrid_upper = quantile(Fz_theory, zGrid_upper_quantile)

		# Allocate Grid
		zGrid = Vector(range(zGrid_lower, zGrid_upper, length = Nz))

		# Exogenous Shock: ContinuousTimeMarkovChain object
		exog_shock = ContinuousTimeMarkovChain(dp, zGrid)

	elseif Nz == 1

		exog_shock = ContinuousTimeMarkovChain(reshape([0], 1, 1), [zd])

	end

	return exog_shock

end


text:
using UltraDark
using Test
using MPI

resol = 16

output_dir = "output"
output_times = 0.0:0.002:0.01

output_config = OutputConfig(output_dir, output_times; box=false)
options = Config.SimulationConfig()

for grid_type in [Grids, PencilGrids]
    grids = grid_type(1.0, 16)
    grids.ψx .= (grids.x.^2 .+ grids.y.^2 .+ grids.z.^2).^0.5 ./ 1e9 # Set ψx to something non-zero
    @test simulate(grids, options, output_config) == nothing
end


text:
# https://github.com/rdeits/MeshCat.jl/pull/99#issuecomment-509067108

@testset "wait" begin
    @testset "with setobject" begin
        vis = Visualizer()
        setobject!(vis, Triad(1.0))
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
    end

    @testset "without setobject" begin
        vis = Visualizer()
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
    end
end


text:
using Test
using SumOfSquares
using DynamicPolynomials
using SemialgebraicSets

# Adapted from:
# SOSDEMO4 --- Matrix Copositivity
# Section 3.4 of SOSTOOLS User's Manual
#
# See also (3.41) of [BPT12].
#
# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

function horn_test(optimizer,
                   config::MOIT.TestConfig,
                   cone::SumOfSquares.PolyJuMP.PolynomialSet)
    atol = config.atol
    rtol = config.rtol

    # Horn matrix
    H = [1 -1  1  1 -1;
        -1  1 -1  1  1;
         1 -1  1 -1  1;
         1  1 -1  1 -1;
        -1  1  1 -1  1]

    @polyvar x[1:5]

    model = _model(optimizer)
    cref = @constraint(model, x' * H * x in CopositiveInner(cone))
    optimize!(model)

    @test termination_status(model) == MOI.INFEASIBLE

    model = _model(optimizer)
    orthant = @set x[1] ≥ 0 && x[2] ≥ 0 && x[3] ≥ 0 && x[4] ≥ 0 && x[5] ≥ 0
    cref = @constraint(model, sum(x) * x' * H * x in cone, domain = orthant)
    optimize!(model)

    @test termination_status(model) == MOI.INFEASIBLE

    model = _model(optimizer)
    orthant = orthant ∩ basicsemialgebraicset(
        algebraicset(polynomialtype(x[1])[]), polynomial.(monomials(x, 3)))
    cref = @constraint(model, sum(x) * x' * H * x in cone, domain = orthant)
    optimize!(model)

    if cone isa SOSCone
        @test termination_status(model) == MOI.OPTIMAL

        @test primal_status(model) == MOI.FEASIBLE_POINT
        @test isempty(certificate_monomials(cref))
        @test isempty(getmat(gram_matrix(cref)))
    else
        @test termination_status(model) == MOI.INFEASIBLE
    end
end
sos_horn_test(optimizer, config)   = horn_test(optimizer, config, SOSCone())
sdsos_horn_test(optimizer, config) = horn_test(optimizer, config, SDSOSCone())
dsos_horn_test(optimizer, config)  = horn_test(optimizer, config, DSOSCone())


text:
module IBMQClient

export AccountInfo, JobInfo, RemoteJob

using HTTP
using JSON
using UUIDs
using URIs
using Dates
using ConfParser
using Configurations
using REPL.TerminalMenus
using Crayons.Box

include("schema.jl")
include("adapters.jl")

include("ibmq.jl")
include("menu.jl")

end


text:
# Draw load adaption plot
# Author: Thomas Minier
using Gadfly
using RDatasets
using CalliPlots

custom_theme = Theme(
    key_position = :none,
    bar_spacing = 5px,
    major_label_font_size=15px,
    major_label_font="Computer Modern",
    major_label_color=colorant"#000000",
    minor_label_color=colorant"#000000",
    key_label_color=colorant"#000000",
    bar_highlight=colorant"#000000"
)

Gadfly.push_theme(custom_theme)

function colors()
 return Scale.color_discrete_manual(colorant"#990000", colorant"#FF5733", colorant"#FFC300", colorant"#DAF7A6")
end

serversNames = Dict(
    "http://curiosiphi:4000/watDiv_100" => "S1",
    "http://curiosiphi:4001/watDiv_100" => "S2",
    "http://curiosiphi:4002/watDiv_100" => "S3",
    "http://curiosiphi:4003/watDiv_100" => "S4",
    "http://curiosiphi:4004/watDiv_100" => "S5",
    "http://curiosiphi:3000/watDiv_100" => "S1",
    "http://curiosiphi:3001/watDiv_100" => "S2",
    "http://curiosiphi:3002/watDiv_100" => "S3",
    "http://curiosiphi:3003/watDiv_100" => "S4",
    "http://curiosiphi:3004/watDiv_100" => "S5",
)

function process(filename)
    file = readtable(filename)
    #  rename servers and convert timestamp in second
    for row in eachrow(file)
        row[:target] = serversNames[row[:target]]
        row[:timestamp] = floor(row[:timestamp]/1000)
    end
    # gather number of HTTP calls per server per timestamp
    df = by(file, [:target, :timestamp]) do x
        return nrow(x)
    end
    # gather timestamp by period (0-4, 5-19,20+)
    df = by(df, [:target, :timestamp]) do x
        timeType = "0s to 19s"
        if x[:timestamp][1] >= 20 #&& x[:timestamp][1] <= 30
            timeType = "20s to end"
        # elseif x[:timestamp][1] >= 31
        #     timeType = "31s+"
        end
        return DataFrame(x1 = x[:x1][1], x2 = timeType)
    end
    # sum http calls per server per period
    df = by(df, [:target, :x2]) do x
        return Float64(sum(x[:x1]))
    end
    return df
end

files = [
    "scripts/adaptivity/load_increase_3servers_run1.csv",
    "scripts/adaptivity/load_increase_3servers_run2.csv",
    "scripts/adaptivity/load_increase_3servers_run3.csv"
]
load_3servers = meanDF(:x1, map(process, files)...)

p = plot(load_3servers, xgroup=:x2, x=:target, y=:x1, color=:target,
Geom.subplot_grid(Geom.bar, Scale.x_discrete, free_x_axis=true),
Guide.xlabel("Elapsed time (seconds)"), Guide.ylabel("# HTTP requests"), colors())
draw(PDF("scripts/adaptivity/load_adaptivity.pdf", 4inch, 3inch), p)
draw(PNG("scripts/adaptivity/load_adaptivity.png", 7inch, 3inch), p)

max_stars_repo_name: Callidon/ulysses-tpf

max_issues_repo_name: Callidon/ulysses-tpf

max_forks_repo_name: Callidon/ulysses-tpf


text:
# -------- Policies ---------------------------------------------------------- #

uniform_policy(length) = ones(Float32, length) / length

function random_policy(length)
  policy = zeros(Float32, length)
  policy[rand(1:length)] = 1.
  policy
end

# -------- Glueing Single Outputs -------------------------------------------- #

function cat_outputs(outputs)

  n = length(outputs)
  v, p, f = first(outputs)

  vs = zeros(Float32, n)
  ps = zeros(Float32, length(p), n)
  fs = zeros(Float32, length(f), n)

  for (i, o) in enumerate(outputs)
    vs[i]   = o[1]
    ps[:,i] = o[2]
    fs[:,i] = o[3]
  end

  vs, ps, fs

end

# -------- Dummy Model ------------------------------------------------------- #

"""
This model returns value 0 and a uniform policy vector for each game state.
"""
struct DummyModel <: Model{Game, false} end

(m :: DummyModel)(g :: Game, args...) = 0f0, uniform_policy(policy_length(g)), Float32[]

(m :: DummyModel)(g :: Vector{<:Game}, args...) = cat_outputs(m.(g, args...))

Base.copy(m :: DummyModel) = m

# -------- Random Model ------------------------------------------------------ #

"""
This model returns value 0 and a randomly drawn policy vector distribution for
each game state.
"""
struct RandomModel <: Model{Game, false} end

(m :: RandomModel)(g :: Game, args...) = 0f0, random_policy(policy_length(g)), Float32[]

(m :: RandomModel)(g :: Vector{<: Game}, args...) = cat_outputs(m.(g, args...))

Base.copy(m :: RandomModel) = m

# -------- Rollout Model ----------------------------------------------------- #

"""
This model executes random moves until the game is over. The result is returned
as the value of the game state. It always proposes a uniform policy vector.
By that, the classical rollout step of MCTS is implemented when this model
is used for the tree search.
"""
struct RolloutModel <: Model{Game, false} end

function (m :: RolloutModel)(g :: Game, args...)

  result = random_playout(g)
  value = status(result) * current_player(g)

  Float32(value), uniform_policy(policy_length(g)), Float32[]

end

(m :: RolloutModel)(g :: Vector{<: Game}, args...) = cat_outputs(m.(g, args...))

Base.copy(m :: RolloutModel) = m

max_stars_repo_name: roSievers/jtac


text:
using Test
using RealNeuralNetworks.Utils.FakeSegmentations
using RealNeuralNetworks.Neurons
using RealNeuralNetworks.NodeNets
using RealNeuralNetworks.SWCs
using RealNeuralNetworks.Utils.SynapseTables

using CSV
using SparseArrays
using LinearAlgebra
using DataFrames

const NEURON_ID = 77625
const ASSET_DIR = joinpath(@__DIR__, "../asset")
const SWC_BIN_PATH = joinpath(ASSET_DIR, "$(NEURON_ID).swc.bin")
const ARBOR_DENSITY_MAP_VOXEL_SIZE = (2000,2000,2000)

@testset "test data structure transformation..." begin
    swc = SWCs.load(joinpath(ASSET_DIR, "example.swc"))
    nodeNet = NodeNet(swc)
    neuron = Neuron(nodeNet)
    swc = SWC(neuron)
    nodeNet = NodeNet(neuron)
    @test length(swc) == length(nodeNet) == Neurons.get_num_nodes(neuron)
end

@testset "test Neurons" begin
    println("\nload swc of a real neuron...")
    swc = SWCs.load_swc_bin( SWC_BIN_PATH )
    neuron = Neuron( swc )

    println("\nreset root node...")
    @time Neurons.reset_root(neuron, (258150.0f0, 135510.0f0, 778175.0f0))

    #neuron = Neurons.resample(neuron, Float32(40))
    println("\nget node list ...")
    @time nodeList = Neurons.get_node_list(neuron)
    @test !isempty(nodeList)
    println("\nget edge list ...")
    @time edgeList = Neurons.get_edge_list(neuron)
    @test !isempty(edgeList)
    println("\nget segment order list...")
    @time segmentOrderList = Neurons.get_segment_order_list( neuron )
    @test !isempty(segmentOrderList)

    println("\nclean up the neuron ...")
    println("original number of segments: $(Neurons.get_num_segments(neuron))")
    println("remove subtree in soma ...")
    @time neuron = Neurons.remove_subtree_in_soma(neuron)
    println("after remove subtree in soma: $(Neurons.get_num_segments(neuron))")

    println("\nremove hair segments ...")
    @time neuron = Neurons.remove_hair(neuron)
    println("after remove hair: $(Neurons.get_num_segments(neuron))")

    println("\nremove terminal blobs ...")
    @time neuron = Neurons.remove_terminal_blobs(neuron)
    println("after remove terminal blobs: $(Neurons.get_num_segments(neuron))")

    println("\nget principle direction based on the furthest terminal node pair...")
    @time vec, maxDistance = Neurons.get_furthest_terminal_node_pair_direction(neuron)
    @test maxDistance > 0

    println("\nremove redundent nodes...")
    @time neuron = Neurons.remove_redundent_nodes(neuron)
    println("after remove redundent nodes: $(Neurons.get_num_segments(neuron))")

    println("\nresampling ...")
    nodeNum1 = Neurons.get_node_num(neuron)
    @time neuron = Neurons.resample(neuron, Float32(400))
    nodeNum2 = Neurons.get_node_num(neuron)
    println("resampling of 400 nm reduced node number form ", nodeNum1, " to ", nodeNum2)
    @test nodeNum2 < nodeNum1

    println("\nsmooth ...")
    nodeNum1 = Neurons.get_node_num(neuron)
    @time neuron = Neurons.smooth(neuron)
    nodeNum2 = Neurons.get_node_num(neuron)
    println("smoothing the neuron should not change node number: ", nodeNum1, " => ", nodeNum2)
    @test nodeNum1 == nodeNum2

    @test Neurons.get_num_segments(neuron) > 0
    #println("get fractal dimension ...")
    #@time fractalDimension, _,_ = Neurons.get_fractal_dimension( neuron )
    #@show fractalDimension

    println("\nget surface area which is frustum based")
    @test Neurons.get_surface_area(neuron) > 0

    println("\nget frustum based volume")
    @test Neurons.get_volume(neuron) > 0

    println("\nget typical radius ...")
    @test Neurons.get_typical_radius( neuron ) > 0

    println("\nget asymmetry ...")
    @test Neurons.get_asymmetry( neuron ) > 0

    println("\nget mass center ...")
    @show Neurons.get_mass_center( neuron )

    println("\nget branching angle ...")
    @test Neurons.get_branching_angle( neuron, 5 ) > 0

    println("\nget path to root length ...")
    @test Neurons.get_path_to_soma_length( neuron, 5; nodeId=4 ) > 0

    println("\nsholl analysis ...")
    @time shollNumList = Neurons.get_sholl_number_list(neuron, 10000 )
    @test !isempty(shollNumList)

    println("\nget segment path length list ...")
    @time segmentPathLengthList = Neurons.get_segment_path_length_list( neuron )
    @show length( segmentPathLengthList )
    @test length( segmentPathLengthList ) == Neurons.get_num_segments(neuron)

    println("\nget terminal segment index list...")
    @time terminalSegmentIdList = Neurons.get_terminal_segment_id_list( neuron )
    @test !isempty( terminalSegmentIdList )

    println("\nget terminal node list ...")
    @time terminalNodeList = Neurons.get_terminal_node_list( neuron )
    @test !isempty( terminalNodeList )

    println("\ntest split a tree...")
    @time tree1, tree2 = split(neuron, 4; nodeIdInSegment = 5)
    @test !isempty(tree1)
    @test !isempty(tree2)
end

@testset "test synapse attaching functions" begin
    println("\nload swc of a real neuron...")
    @time swc = SWCs.load_swc_bin( SWC_BIN_PATH )
    @time neuron = Neuron( swc )

    println("\ndownsample node number...")
    @time Neurons.downsample_nodes(neuron)

    println("\nattaching presynapses in DataFrame...")
    # the default reading of CSV is immutable and the dropmissing will fail!
    # https://github.com/JuliaData/DataFrames.jl/issues/1393#issuecomment-500019773
    preSynapses = CSV.read( joinpath(ASSET_DIR, "$(NEURON_ID).pre.synapses.csv");  copycols=true)
    @time preSynapses = SynapseTables.preprocess(preSynapses, (5,5,45))
    println("get ", DataFrames.nrow(preSynapses), " synapses.")
    @time Neurons.attach_pre_synapses!(neuron, preSynapses)
    numPreSynapses = Neurons.get_num_pre_synapses(neuron)
    println("have ", numPreSynapses, " synapses after attachment.")
    @test numPreSynapses <= DataFrames.nrow(preSynapses)
    @test numPreSynapses > 2

    println("\nattaching postsynapses in DataFrame...")
    postSynapses = CSV.read( joinpath(ASSET_DIR, "$(NEURON_ID).post.synapses.csv"); copycols=true)
    @time postSynapses = SynapseTables.preprocess(postSynapses, (5,5,45))
    println("get ", DataFrames.nrow(postSynapses), " synapses.")
    @time Neurons.attach_post_synapses!(neuron, postSynapses)
    numPostSynapses = Neurons.get_num_post_synapses(neuron)
    println("have ", numPostSynapses, " synapses after attachment.")
    @test numPostSynapses <= DataFrames.nrow(postSynapses)
    @test numPostSynapses > 2

    println("\nattaching presynapse list...")
    preSynapseList = Neurons.get_all_pre_synapse_list(neuron)
    println("get ", length(preSynapseList), " synapses.")
    neuronx = Neuron(swc)
    @time Neurons.attach_pre_synapses!(neuronx, preSynapseList)
    numPreSynapses = Neurons.get_num_pre_synapses(neuronx)
    println("have ", numPreSynapses, " synapses after attachment.")
    @test numPreSynapses > 2

    println("\nattaching postsynapse list...")
    postSynapseList = Neurons.get_all_post_synapse_list(neuron)
    println("get ", length(postSynapseList), " synapses.")
    @time Neurons.attach_post_synapses!(neuronx, postSynapseList)
    numPostSynapses = Neurons.get_num_post_synapses(neuronx)
    println("have ", numPostSynapses, " synapses after attachment.")
    @test numPostSynapses > 2

    println("\nresample neuron and double check the synapses.")
    neuron2 = Neurons.resample(neuron, Float32(1000))
    @time numPostSynapses2 = Neurons.get_num_post_synapses(neuron2)
    @time numPreSynapses2 = Neurons.get_num_pre_synapses(neuron2)
    @test numPostSynapses >= numPostSynapses2
    @test numPreSynapses >= numPreSynapses2

    println("\npostprocessing neuron and make sure that the synapses are still there.")
    neuron2 = Neurons.postprocessing(neuron)
    @time numPostSynapses2 = Neurons.get_num_post_synapses(neuron2)
    @time numPreSynapses2 = Neurons.get_num_pre_synapses(neuron2)
    @test numPostSynapses >= numPostSynapses2
    @test numPreSynapses >= numPreSynapses2

    println("\nadjust segment class...")
    @time Neurons.adjust_segment_class!(neuron)

    preSynapseList = Neurons.get_all_pre_synapse_list(neuron)
    postSynapseList = Neurons.get_all_post_synapse_list(neuron)
    @test !isempty(preSynapseList)
    @test !isempty(postSynapseList)

    println("\nget synapse to soma path length list...")
    @time preSynapseToSomaPathLengthList  = Neurons.get_pre_synapse_to_soma_path_length_list( neuron )
    @time postSynapseToSomaPathLengthList = Neurons.get_post_synapse_to_soma_path_length_list( neuron )
    @test !isempty(preSynapseToSomaPathLengthList)
    @test !isempty(postSynapseToSomaPathLengthList)
end

@testset "test Neuron IO and resampling " begin
    println("\nload swc of a real neuron...")
    @time swc = SWCs.load_swc_bin( SWC_BIN_PATH )
    neuron = Neuron( swc )

    println("\ndownsample node number...")
    Neurons.downsample_nodes(neuron)

    println("\ncompute arbor density map...")
    @time arborDensityMap = Neurons.get_arbor_density_map(neuron,
                                                ARBOR_DENSITY_MAP_VOXEL_SIZE, 8.0)
    #@test norm(arborDensityMap[:]) ≈ Neurons.get_total_path_length(neuron)
    #@test norm(arborDensityMap[:]) ≈ 1.0
    neuron1 = neuron
    densityMap1 = Neurons.translate_soma_to_coordinate_origin(neuron1, arborDensityMap,
                                                             ARBOR_DENSITY_MAP_VOXEL_SIZE)
    @show axes(densityMap1)
    println("compute arbor density map distance...")
    @time d = Neurons.get_arbor_density_map_distance(densityMap1, densityMap1)
    @test d == 0.0

    fileName = joinpath(dirname(SWC_BIN_PATH), "$(NEURON_ID).swc.bin")
    neuron2 = Neuron( SWCs.load_swc_bin(fileName) )
    densityMap2 = Neurons.get_arbor_density_map(neuron, ARBOR_DENSITY_MAP_VOXEL_SIZE, 8.0)
    #@test norm(densityMap2[:]) ≈ 1.0
    densityMap2 = Neurons.translate_soma_to_coordinate_origin(neuron2, densityMap2,
                                                              ARBOR_DENSITY_MAP_VOXEL_SIZE)
    @show axes(densityMap2)
    println("compute arbor density map distance...")
    @time d = Neurons.get_arbor_density_map_distance(densityMap1, densityMap2)
    println("arbor density map distance: $d")
    #@test d > 0.0 && d < 2.0

    Neurons.save(neuron, "/tmp/neuron.swc")
    #Neurons.save_swc(neuron2, "/tmp/neuron2.swc")
    rm("/tmp/neuron.swc")
    #rm("/tmp/neuron2.swc")
end

@testset "test fake segmentation skeletonization" begin
    println("create fake cylinder segmentation...")
    @time seg = FakeSegmentations.broken_cylinder()
    println("skeletonization to build a Neuron ...")
    @time neuron = Neuron(seg)
    @test !isempty(Neurons.get_node_list(neuron))

    println("create fake ring segmentation ...")
    seg = FakeSegmentations.broken_ring()
    neuron = Neuron(seg)
    @test !isempty(Neurons.get_node_list(neuron))

    println("transform to SWC structure ...")
    @time swc = SWCs.SWC( neuron )
    tempFile = tempname() * ".swc"
    SWCs.save(swc, tempFile)
    rm(tempFile)
end


text:
using Gnuplot
@gp 1:10
#save("plt2_ex1.gp") # hide
#save(term="cairolatex standalone", output="plt2_ex1.tex") # hide
#save(term="pdf", output="plt2_ex1.pdf") # hide
save(term="png", output="./code/plt2_ex1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# A simple Unicode raster 8x8 font
# The raw font data was ported from a Rust crate:
# https://crates.io/crates/font8x8/0.2.3

export
    show_character, show_string, get_glyph, get_glyph_width,
    load_font!

# The pixel data for this font is in assets/raster_font.data
mutable struct RasterFont
    pixels::Array{Array{UInt8},1}
    glyphs::Dict{Char,Int64}

    function RasterFont()
        new([], Dict{Char,Int64}())
    end
end

function load_font!(font::RasterFont, file::String)
    # The format of the file is formatted as follows:
    #
    #      This is the pixel data for the character
    #   |-------------------------------------|
    # ! 0x18 0x3C 0x3C 0x18 0x18 0x00 0x18 0x00
    # ^
    # |
    #  \- This is the character = "!"

    try
        idx = 1
        lines = readlines(file)
        for line in lines
            ele = split(line, " ")
            # Add character to glyph dictionary
            font.glyphs[ele[1][1]] = idx
            idx += 1
            # Add data to raw data array
            e2 = parse(UInt8, ele[2])
            e3 = parse(UInt8, ele[3])
            e4 = parse(UInt8, ele[4])
            e5 = parse(UInt8, ele[5])
            e6 = parse(UInt8, ele[6])
            e7 = parse(UInt8, ele[7])
            e8 = parse(UInt8, ele[8])
            e9 = parse(UInt8, ele[9])
            push!(font.pixels, [e2, e3, e4, e5, e6, e7, e8, e9])
        end
    catch ex
        println("******************************************************")
        println(ex)
        println("Exception loading: ", file)
        println("******************************************************")
        return false
    end

    true
end

function get_glyph(font::RasterFont, char::Char)
    glyph = font.glyphs[char]
    font.pixels[glyph]
end

function get_glyph_width()
    8
end

# Shows character at the console
# example:
# "a" =
#░░░░░░░░
#░░░░░░░░
#░████░░░
#░░░░██░░
#░█████░░
#██░░██░░
#░███░██░
#░░░░░░░░
# Helpful debugging functions
function show_character(char::Char)
    row = get_glyph(char)
    for c in row
        # println(bitstring(c))
        for shift in [0,1,2,3,4,5,6,7]
            bit = (c >> shift) & 1
            if bit == 0
                print("░")
            else
                print("█")
            end
        end
        println("")
    end
end

function show_string(text::String)
    for c in text
        show_character(c)
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using .JSort
import Base: isless

struct GammaCascade
    gate::AbstractGate
    e::Vector{Float32}
    regions::Vector{Tuple{Float32, Float32}}
    keywords::Dict{Symbol, Any}
    #function GammaCascade(gate::Gate, e::AbstractVector, regions::AbstractVector)
        #new(gate, sort(e) |> reverse)
    #end
end
GammaCascade(gate::AbstractGate, e::AbstractVector, regions::AbstractVector) = GammaCascade(gate, e, regions, Dict(:align => Dict()))

function extend!(base::GammaCascade, others::AbstractArray{GammaCascade})
    lower = [c for c in sort(others) if c < base]
    for other in lower
        extend!(base, other)
    end
end

function extend!(base::GammaCascade, other::GammaCascade; width=3)
    from = center(base.gate)
    to = center(other.gate)
    for e in base.e
        if abs(from - e - to) < width
            append!(base.e, other.e)
            break
        end
    end
end

isless(c1::GammaCascade, c2::GammaCascade) = isless(c1.gate, c2.gate)
center(c::GammaCascade) = center(c.gate)

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# these types are mandatory indepedent of used solver##############################
#
"""
```
NMRinfo :: Type
```
This type contains the information for `NMRtype`. In this point of the project it contains two fields: dist, atom1 and atom2. Essentially, dist is the distance between atom1 and atom2 which are strings related to the name of atoms.
"""
struct NMRInfo
	dist::Float64
	atom1::String
	atom2::String
end

"""
```
NMRType :: Type
```
This type is used as input data. Essentially, this type carry itself fields for internal solvers. Although it is possible to input data directly into the type, it is highly recommended that this type be constructed from the function `preprocessing`, like in the following example

## Example

```julia-repl
nmr_data_example = preprocessing("1A57.nmr")
```
As return an element of `NMRtype` is created. We are assuming that a `1A57.nmr` file is provided.
"""
struct NMRType
	virtual_path :: Vector{Int64}
	additional_distance::Vector{Vector{Int64}}
	info :: SparseMatrixCSC{NMRInfo,Int64}
	dim :: Int64
end

"""
```
preprocessing :: Function
```
This function is able to read a `.nmr` file and generate a `NMRType`, which is an internal type with fields that can be used by internals solvers. Given a list in `.nmr` format, let us say the `1A57.nmr` file, the following example show how simple is to use this function.

## Example
```julia-repl
preprocessing("1A57.nmr")
```
"""
function preprocessing(file::String,opt="read")
	if opt == "read"
		nmrfile = readdlm("$(file)")
		I = Int64.(nmrfile[:,1])
		J = Int64.(nmrfile[:,2])
		n = last(I)
		vpath = generate_virtual_path(I,J)
		lenI = length(I)
		vadd = Vector{Vector{Int64}}(undef,I[lenI])
		for i=1:I[lenI]
			vadd[i]=[]
		end
		V = Vector{NMRInfo}(undef,lenI)
		i = 1
		while i<=lenI
			k = 3
			while i+k<=lenI && I[i]==I[i+k]
				push!(vadd[I[i]],J[i+k])
				k = k+1
			end
			i = i+k
		end

		for i=1:lenI
			V[i]=NMRInfo(nmrfile[i,5],nmrfile[i,7],nmrfile[i,8])
		end
		for i=1:lenI
			if I[i]!=J[i]
				push!(I,J[i])
				push!(J,I[i])
				push!(V,NMRInfo(V[i].dist,V[i].atom2,V[i].atom1))
			end
		end
		nmrt = NMRType(vpath,vadd,sparse(I,J,V),n)
	else
		error("Unidentified option or file")
	end
	return nmrt
end

"""
```
generate_virtual_path :: Function
```
This function is an auxiliary function used to define a useful vector called in our context as virtual path. This vector allows to handle with re-order approach.
"""
function generate_virtual_path(NMRdatavertex1::Vector{Int64},NMRdatavertex2::Vector{Int64})
	#	D = [NMRdata.vertex1 NMRdata.vertex2]
	virtual_path = [1,2,3,4]
	k = 5
	li = 4
	while k <= last(NMRdatavertex1)
		ind = findall(x->x==k,NMRdatavertex1)
		la = 0
		if NMRdatavertex2[ind[1:3]] ==  virtual_path[li-2:li]
			la+=1
		else
			if NMRdatavertex2[ind[1:2]] == virtual_path[li-2:li-1]
				la+=2
				push!(virtual_path,NMRdatavertex2[ind[3]])
			elseif NMRdatavertex2[ind[1:2]] == virtual_path[li-1:li]
				la+=2
				push!(virtual_path,NMRdatavertex2[ind[3]])
			elseif NMRdatavertex2[ind[2:3]] == virtual_path[li-2:li-1]
				la+=2
				push!(virtual_path,NMRdatavertex2[ind[1]])
			elseif NMRdatavertex2[ind[2:3]] == virtual_path[li-2:li]
				la+=2
				push!(virtual_path,NMRdatavertex2[ind[1]])
			else
				la+=4
				append!(virtual_path,NMRdatavertex2[ind[1:3]])
			end

		end
		li = li+la
		push!(virtual_path,k)
		k += 1
	end
	return virtual_path
end

# these types are mandatory indepedent of used solver##############################
"""
```
ConformationSetup :: Type
```
This type is used to define the main options of the function used to conformation. Essentially, this type has three fields: precision, solver, allsolutions.

precision ::  is associated to current internal solver which depends on some level do a pruning test.

solver :: is the name of function which define the solver.

allsolution :: is a bool variable used to define the behavior of the method. Is set as true, the internal solver try to find all possible solutions.

## Example

```julia-repl
options = ConformationSetup(0.00001,classicBP,true)
```
As return an element options was created an it will used in conformation function. It is a mandatory definition before run conformation function.

"""
struct ConformationSetup
	precision :: Float64
	virtual_precision :: Float64
	solver :: Function
	evalLDE :: Bool
	allsolutions :: Bool
	function ConformationSetup(a,c,f,d=true,b=1.0e-6)
		return new(a,b,c,d,f)
	end
end

"""
```
AtomType :: Type
```
It is the most primitive type and it is used to store the spatial positions of atoms. It is a mutable type. This type has three field: .x, .y and .z
"""
mutable struct AtomType
	x :: Float64
	y :: Float64
	z :: Float64
	element :: String
	function AtomType(a,b,c,d="None")
		return new(a,b,c,d)
	end
end

"""
```
MoleculeType :: Type
```
It is used to store a (just one) solution find by conformation process. It is a mutable struct and has two fields: .atoms (used to store a vector of AtomType) and .lde (used to store the lde of an molecule)
"""
mutable struct MoleculeType
	atoms :: Vector{AtomType}
	lde :: Float64
end

function writefile(A::MoleculeType,file="SpatialPosition";format = "xyz")
	if isa(A,MoleculeType)
		B = [length(A.atoms)," "]
		for i=1:length(A.atoms)
			push!(B,[A.atoms[i].element,A.atoms[i].x,A.atoms[i].y,A.atoms[i].z])
		end
		open("$(file).$(format)",io) do io
			writedlm(io,B)
		end
	else
		error("Not recognize format")
	end
end

"""
```
Counter :: Type
```
It is a mutable type for counting operations in solvers. This can handle with counters for nodes, virtual positions, ddf, branchs and prunes.
"""
mutable struct Counter
	node :: Vector{Int64}
	virtual_path :: Vector{Int64}
	ddf :: Vector{Int64}
	branch :: Int64
	prune :: Int64
end

"""
```
ConformationOutput :: Type
```
It is a mutable type used to store the output provided by conformation function. With this type is possible to handle with some importants elements given by the following fields:.number,.molecules,.elapsedtime,.bytes and .gctime
"""
mutable struct ConformationOutput
	solver :: Function
	number :: Int64
	molecules :: Dict{Int64,MoleculeType}
	nop :: Counter
end

###################################################################################
# these functions are used by ClassicBP  ##########################################
"""
```
bondangle :: Function
```
This is an auxiliary function used by ClassicBP solver in order to compute the bond angle by cosine rule. As output the cosine and sine are given.
"""
function bondangle(d23,d24,d34)
	c = (-d24^2 + d34^2 + d23^2)/(2.0*d34*d23)
	#	c = (-D[i-2,i]^2+ D[i-1,i]^2+ D[i-2,i-1]^2)/(2.0*D[i-1,i]*D[i-2,i-1])
	if c<-1.0
		c=-1.0
	end
	if c>1.0
		c=1.0
	end
	s = sqrt(1.0-c^2)
	@debug "value of cθ and sθ" c,s
	return c,s
end

"""
```
badtorsionangle :: Function
```
This is an auxiliary function used by classicBP solver in order to compute the torsion angle. As output the cosine and sine of the torsion angle are given.
"""
function badtorsionangle(d12,d13,d14,d23,d24,d34)#i=4,...,n
	#d12=D[i-3,i-2]
	#d13=D[i-3,i-1]
	#d14=D[i-3,i]
	#d23=D[i-2,i-1]
	#d24=D[i-2,i]
	#d34=D[i-1,i]
	a = d12*d12 + d24*d24 - d14*d14
	a = a/(2.0*d12*d24)
	b = d24*d24 + d23*d23 - d34*d34
	b = b/(2.0*d24*d23)
	c = d12*d12 + d23*d23 - d13*d13
	c = c/(2.0*d12*d23)
	e = 1.0 - b^2;
	f = 1.0 - c^2;
	if (e < 0.0 || f < 0.0)
		@debug "some problem in torsion angle"
		return -2
	end
	ef = sqrt(e*f)
	valc = (a - b*c)/(ef)
	if (valc < -1.0)
		valc = -1.0
	end
	if (valc >  1.0)
		valc =  1.0
	end
	vals=sqrt(1.0-valc^2)
	@debug "value of cω and sω" valc,vals
	return valc,vals
	#number of operations [10,20,3,2]
end

"""
```
torsionmatrix :: Function
```
This is an auxiliary function  used by classicBp solver to compute the torsion array.
"""
function torsionmatrix(cosθ,sinθ,cosω,sinω,d34,B,sign::Bool)
	if sign == true

		B=zeros(4,4)
		B[1,1] = -cosθ
		B[1,2] = -sinθ
		B[1,4] = -d34*cosθ
		B[2,1] = sinθ*cosω
		B[2,2] = -cosθ*cosω
		B[2,3] = -sinω
		B[2,4] = d34*B[2,1]
		B[3,1] = sinθ*sinω
		B[3,2] = -cosθ*sinω
		B[3,3] = cosω
		B[3,4] = d34*B[3,1]
		B[4,4] = 1
	else

		B[2,3] = -B[2,3]
		B[3,1] = -B[3,1]
		B[3,2] = -B[3,2]
		B[3,4] = -B[3,4]
	end
	return B
end

"""
```
pruningtest :: Function
```
This functions is an auxiliary function used to test if some molecule is
feasible or not.
"""
function pruningtest(v::MoleculeType,i::Int64,D::NMRType,ε::Float64,count::Vector{Int64})
	if isempty(D.additional_distance[i])

		@debug "result of prunning 1"
		return 1,count
	else
		for j in D.additional_distance[i]
			dij =  (v.atoms[i].x-v.atoms[j].x)^2+(v.atoms[i].y-v.atoms[j].y)^2 +(v.atoms[i].z-v.atoms[j].z)^2
			count += [5,3,0,0]
			if (D.info[i,j].dist^2 -dij)^2 >ε
				count += [1,1,0,0]
				@debug "result of prunning 0"
				return 0,count
			end
		end
	end

	@debug "result of prunning 1"
	return 1,count
end

####################################################################################
# these functions can be used by any solver#########################################
"""
```
LDE :: Function
```
This function is an auxiliary function used to compute the LDE of some molecule.
This function change the .lde field of the molecule.
"""
function LDE(v::MoleculeType,D::NMRType)
	dij=0.0
	nne = findnz(D.info)
	num_nne = length(nne[1])
	maxε = 0.0
	for k=1:num_nne
		i=nne[1][k]
		j=nne[2][k]
		aux = (D.info[i,j].dist^2-((v.atoms[i].x-v.atoms[j].x)^2+(v.atoms[i].y-v.atoms[j].y)^2 +(v.atoms[i].z-v.atoms[j].z)^2))^2
		maxε = max(maxε,aux)
		dij = dij+(aux)/D.info[i,j].dist
	end
	#	println("max ε = $(maxε)")
	#	if nad>0.0
	v.lde = dij/(2.0*num_nne)
	#	else
	#		return dij
	#	end
end

"""
```
outputfilter
```
This function is especially useful when we need to filter an information about a solution in ConformationOutput. For example, let us suppose that we want to compute the worst LDE measure. So, we need to type:

julia> a = conformation(data,option)

julia> outputfilter(a,"lde")
"""
function outputfilter(a::ConformationOutput, option = "lde")
	mol = a.molecules
	num = a.number
	if option == "lde"
		vlde = zeros(num)

		for i=1:num
			vlde[i]=mol[i].lde
		end
		return minimum(vlde)
	end
	if option == "xyz"
		nl = length(a.molecules[1].atoms)
		A=Array{Vector{Float64},2}(undef,nl,num)
		for j=1:num
			for i=1:nl
				A[i,j]=[a.molecules[j].atoms[i].x, a.molecules[j].atoms[i].y,
					a.molecules[j].atoms[i].z]
			end
		end
		return A
	end
end

"""
```
build_distance_matrix
```
This function can be useful to visualize a the distance array associated to an
MoleculeType solution

## Example

```julia-repl
v = conformation(data,opt)
build_distance_array(v.molecules[1].atoms)

```
"""
function build_distance_matrix(v::Array{AtomType,1})
	n=length(v)
	D=zeros(n,n)
	for i=1:n
		for j=1:n
			D[i,j] = sqrt((v[i].x-v[j].x)^2+(v[i].y-v[j].y)^2+(v[i].z-v[j].z)^2)
		end
	end
	return D
end

function Base.copy(A::AtomType)
	return AtomType(A.x,A.y,A.z,A.element)
end

function Base.copy(M::MoleculeType)
	nl = length(M.atoms)
	vm = Vector{AtomType}(undef,nl)
	for i=1:nl
		vm[i]=copy(M.atoms[i])
	end
	return MoleculeType(vm,M.lde)
end

function Base.:≈(A::MoleculeType,B::MoleculeType)
	na = length(A.atoms)
	nb = length(B.atoms)
	v = zeros(3)
	if na == nb
		for i=1:na
			v = [A.atoms[i].x-B.atoms[i].x,A.atoms[i].y-B.atoms[i].y,A.atoms[i].z-B.atoms[i].z]
			if norm(v)>1.0e-4
				return false
			end
		end
	end
	return true
end
####################################################################################
function symsparse(I,J,v)
	for i=1:length(I)
		if I[i]!=J[i]
			push!(I,J[i])

			push!(J,I[i])
			push!(v,v[i])
		end
	end
	return sparse(I,J,v)
end

# these functions are used by quaternion_bp#########################################
function qbondangle(d23,d24,d34)
	c = (-d24^2+ d34^2+ d23^2)/(4.0*d34*d23)
	if c<-0.5
		c=-0.5
	end
	if c>0.5
		c=0.5
	end
	#cm = c/2.0
	return sqrt(0.5 + c),sqrt(0.5 - c)
	# number of operations =  [4,5,1,2]
end

function qtorsionangle(d12,d13,d14,d23,d24,d34)
	############################################

	d12m = d12*d12
	d23m = d23*d23
	d24m = d24*d24
	d34m = d34*d34
	dtil3 = d12m + d23m - d13*d13
	dtil2 = d24m + d23m - d34m
	valc = d23m*(d12m + d24m - d14^2) - 0.5*dtil3*dtil2
	valc = valc/sqrt((4.0*d12m*d23m - dtil3^2)*(4.0*d23m*d24m - dtil2^2))
	#	if (valc < -1.0)
	#		valc = -1.0
	#	end
	#	if (valc >  1.0)
	#		valc =  1.0
	#	end
	#	vals=sqrt(1.0-valc^2)

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

	#	a = d12*d12 + d24*d24 - d14*d14
	#	a = a/(2.0*d12*d24)
	#	b = d24*d24 + d23*d23 - d34*d34
	#	b = b / (2.0*d24*d23)
	#	c = d12*d12 + d23*d23 - d13*d13
	#	c = c / (2.0*d12*d23)
	#	e = 1.0 - b^2;
	#	f = 1.0 - c^2;
	#	if (e < 0.0 || f < 0.0)
	#		return -2
	#	end
	#	ef = 2.0*sqrt(e*f)
	#	valc = (a - b*c)/(ef)
	if (valc < -0.5)
		valc = -0.5
	end
	if (valc >  0.5)
		valc =  0.5
	end
	return sqrt(0.5+valc),sqrt(0.5-valc)
	# number of operations = [11,16,1,3]
end

# Quaternion small library
mutable struct Quaternion
	s :: Float64
	v1 :: Float64
	v2 :: Float64
	v3 :: Float64
end

#function qsign(q::Quaternion)
#	return Quaternion(-q.s, -q.v1, -q.v2, -q.v3)
#end

#function qsum(q::Quaternion,w::Quaternion)
#	return Quaternion(q.s + w.s, q.v1 + w.v1, q.v2 + w.v2, q.v3 + w.v3)
#end

#function qminus(q::Quaternion,w::Quaternion)
#	return Quaternion(q.s - w.s, q.v1 - w.v1, q.v2 - w.v2, q.v3 - w.v3)
#end

function qprod(q::Quaternion,w::Quaternion)
	return  Quaternion(q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3,
			   q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2,
			   q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1,
			   q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s)
	# number of operations = [12,16,0,0]
end

#function conj(q::Quaternion)
#	return  Quaternion(q.s, -q.v1, -q.v2, -q.v3)
#end

function rotopt(Q::Quaternion,t::Float64)
	sl = 2.0*t
	p1 = sl*(Q.s^2+Q.v1^2-0.5)
	p2 = sl*(Q.v2*Q.v1 + Q.v3*Q.s)
	p3 = sl*(Q.v3*Q.v1 - Q.v2*Q.s)
	return Quaternion(0.0,p1,p2,p3)
	# number of operations = [4,10,0,0]
end
# to be fair with memory acess in comparations
function prodmatrix(A::Array{Float64,2},B::Array{Float64,2})
	C=zeros(4,4)
	C[4,4] = 1.0
	for i=1:3
		C[i,1] = A[i,1]*B[1,1] + A[i,2]*B[2,1] + A[i,3]*B[3,1]
		C[i,2] = A[i,1]*B[1,2] + A[i,2]*B[2,2] + A[i,3]*B[3,2]
		C[i,3] = A[i,2]*B[2,3] + A[i,3]*B[3,3]
		C[i,4] = A[i,1]*B[1,4] + A[i,2]*B[2,4] + A[i,3]*B[3,4] + A[i,4]
	end
	return C
end
#[24,33,0,0]
#[48,64,0,0]
"""
	convert_to_dataframe

A function used to convert an array of AtomType elements in a dataframe

# Examples
```
julia-repl
julia> sol = conformation(data,options) # assuming all solutions founded

julia> convert_to_dataframe(sol.molecules[1].atoms)

returns a dataframe object.
```
"""
function convert_to_dataframe(A::Array{AtomType,1})
	n = length(A)
	df = DataFrame()
	df.atom = []
	df.x = []
	df.y = []
	df.z = []
	for k = 1:n
		push!(df.atom,A[k].element)
		push!(df.x,A[k].x)
		push!(df.y,A[k].y)
		push!(df.z,A[k].z)
	end
	return df
end


text:
using PkgBenchmark
curr = benchmarkpkg("..")
prev = benchmarkpkg("..", "HEAD~")
judge(curr, prev).benchmarkgroup


text:
using Distributed
@everywhere using Dates
@everywhere using SparseArrays
@everywhere using LinearAlgebra

include("utils.jl")
include("Corpora.jl")

include_everywhere("GloVe.jl")
include_everywhere("Bias.jl")

PRINT_EVERY = 1_000

function pre_compute(M, X, word_indices)
    num_words = length(word_indices)
    H = Dict{Int64,Array{Float64,2}}()
    G = Dict{Int64,Array{Float64,1}}()
    for i in word_indices
        H[i] = inv(GloVe.∇²Li(M.U, X, i))
        G[i] = GloVe.∇Li(M.W, M.b_w, M.U, M.b_u, X, i)
    end
    return H, G
end

function make_job_list(corpus; first=1, last=2^62)
    # Eventually should allow prcessing subsets
    return first:min(corpus.num_documents, last)
end

function queue_jobs(corpus, job_channel, job_list)
    println("$(now()) - Started queuing jobs.")
    flush(stdout)
    open(corpus.corpus_path) do f
        for doc_num in job_list
            # println("made $doc_num")
            doc = Corpora.get_text(f, corpus, doc_num)
            put!(job_channel, (doc_num, doc))
        end
    end
    println("All jobs queued.")
end

function handle_results(out_file, results, num_jobs)
    println("$(now()) - Started handling results.")
    println("Outfile: $(abspath(out_file))")
    i::Int64 = 1
    num_biases = length(Bias.WEAT_WORD_SETS)
    headers = vcat(["pid", "doc_num"], ["ΔBIF_$x" for x in 1:num_biases])
    open(out_file, "w") do f
        println(f, join(headers, ", "))
        while (i <= num_jobs)
            # println("handled $i")
            pid, doc_num, ΔBIF = take!(results)
            println(f, "$pid, $doc_num, " * join(ΔBIF, ", "))
            if (i == 1 || i % PRINT_EVERY == 0)
                println("$(now()) - Results from job $(i)")
                flush(stdout)
                flush(f)
            end
            i += 1
        end
    end
    println("Handled all results.")
    flush(stdout)
end

@everywhere function percent_diff_bias(document, M, X, inv_hessians, gradients, weat_idx_sets,
                           effect_sizes, doc_num)
    # Make the IF approximation
    target_indices = unique([i for set in weat_idx_sets for inds in set for i in inds])
    deltas = GloVe.my_compute_IF_deltas(document, M, X, target_indices, inv_hessians, gradients, "results/diff_bias-C0-V20-W8-D25-R0.05-E15-S1.csv", doc_num)

    # println("size of deltas is:")
    # println(size(deltas))
    # M.W = M.W + deltas # maybe this should be minus? it is the change in the embeddings.
    # deltas = GloVe.compute_IF_deltas(document, M, X, target_indices, inv_hessians, gradients)

    # Compute the bias change
    B̃ = [Bias.effect_size(M.W, set, deltas) for set in weat_idx_sets]
    # this is the percent difference bias. I don't think this makes sense for us
    # this is actually telling us
    # let's just do it. I think these should be zeros now
    return [100 * (b - b̃) / b for (b, b̃) in zip(effect_sizes, B̃)]
end

@everywhere function run_worker(jobs, results, M, X, inv_hessians, gradients, weat_idx_sets,
                            effect_sizes)
    pid = myid()
    println(stdout, "Starting worker $(pid)")
    flush(stdout)
    while true
        if !isready(jobs)
            println(stdout, "Worker $(pid) waiting for jobs...")
            flush(stdout)
        end
        (doc_num, doc) = try take!(jobs)
        catch
            break # Channel has been closed
        end

        ΔBIF = percent_diff_bias(doc, M, X, inv_hessians, gradients, weat_idx_sets,
                                   effect_sizes, doc_num)

        put!(results, (pid, doc_num, ΔBIF))
    end
    println(stdout, "Ending worker $(pid)")
    flush(stdout)
end

function main()
    embedding_path = get(ARGS, 1, "embeddings/vectors-C0-V20-W8-D25-R0.05-E15-S1.bin")
    corpus_path = get(ARGS, 2, "corpora/simplewikiselect.txt")
    out_file = get(ARGS, 3, "results/diff_bias-C0-V20-W8-D25-R0.05-E15-S1-dummy-B2.csv")

    println("Computing Differential Bias")
    print("$(now()) - Loading model and corpus... ")
    M = GloVe.load_model(embedding_path)
    corpus = Corpora.Corpus(corpus_path)
    println("Done.")
    println("Vocabulary: $(M.vocab_path) ($(M.V) words)")
    println("Embedding: $(M.embedding_path) ($(M.D) dimensions)")
    println("Corpus: $(corpus.corpus_path) ($(corpus.num_words) tokens, $(corpus.num_documents) docs)")

    place_holder=M.W
    # WEAT BIAS
    weat_idx_sets = [Bias.get_weat_idx_set(set, M.vocab) for set in Bias.WEAT_WORD_SETS]
    all_weat_indices = unique([i for set in weat_idx_sets for inds in set for i in inds])
    effect_sizes = [Bias.effect_size(M.W, set) for set in weat_idx_sets]
    p_values = [Bias.p_value(M.W, set) for set in weat_idx_sets]
    println("WEAT effect sizes: $effect_sizes")
    println("WEAT p-values: $p_values")

    print("$(now()) - Loading coocurence matrix... ")
    cooc_path = "$(dirname(embedding_path))/cooc$(match(r"-C[0-9]+-V[0-9]+-W[0-9]+", embedding_path).match).bin"
    X = GloVe.load_cooc(cooc_path, M.V, all_weat_indices)
    println("Done.")
    println("Cooc: $(abspath(cooc_path)) ($(nnz(X)) nnz)")

    print("$(now()) - Precomputing Hessians and Gradients... ")
    inv_hessians, gradients = pre_compute(M, X, all_weat_indices)
    println("Done.")

    # println(all(M.W .== M.W))

    # Setup for parallel computations
    jobs = RemoteChannel(()->Channel{Tuple{Int64, String}}(100))
    results = RemoteChannel(()->Channel{Tuple}(100))
    job_list = make_job_list(corpus)

    @sync begin
        @async queue_jobs(corpus, jobs, job_list)

        for p in workers()
            @async remote_do(run_worker, p, jobs, results, M, X, inv_hessians, gradients, weat_idx_sets, effect_sizes)
        end

        # I don't need to write a file
        @async handle_results(out_file, results, length(job_list))
    end

    # here's where we would add the new model

    close(jobs)
    close(results)

    # Recalculating the WEAT bias after updating the offending word vectors
    # WEAT BIAS
    weat_idx_sets = [Bias.get_weat_idx_set(set, M.vocab) for set in Bias.WEAT_WORD_SETS]
    all_weat_indices = unique([i for set in weat_idx_sets for inds in set for i in inds])
    effect_sizes = [Bias.effect_size(M.W, set) for set in weat_idx_sets]
    p_values = [Bias.p_value(M.W, set) for set in weat_idx_sets]
    println("New WEAT effect sizes: $effect_sizes")
    println("New WEAT p-values: $p_values")

    println("$(now()) - Done.")
end

main()

max_stars_repo_name: hmcgovern/understanding-bias

max_issues_repo_name: hmcgovern/understanding-bias

max_forks_repo_name: hmcgovern/understanding-bias


text:
### TRIANGLE
area(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T} = norm(vector_area(v1, v2, v3))
centroid(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T} = (v1 + v2 + v3) * Float64(1/3)  # 4 times faster than dividing by 3
vector_area(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T} = cross(v2 - v1, v3 - v2) * 0.5
function triangle_area(v123::NTuple{3,SVector{3,T}}, n̂::SVector{3,T}) where {T}
	vec_area = vector_area(v123[1], v123[2], v123[3])
	return dot(n̂, vec_area)
end
triangleNormal(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T} = normalize(vector_area(v1, v2, v3))

for funName in (:area, :centroid, :triangleNormal, :triangleCross)
    @eval begin
        function $funName(sv::SVector{3,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3])
        end
    end
end

### TETRAHEDRON
function centroid(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}, v4::SVector{3,T}) where {T}
    return (v1 + v2 + v3 + v4) * 0.25
end

function volume(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}, v4::SVector{3,T}) where {T}
    # NOTE: This is an algebraic refactorization of a symbolic answer
    a1, a2, a3 = v1[1], v1[2], v1[3]
    b1, b2, b3 = v2[1], v2[2], v2[3]
    c1, c2, c3 = v3[1], v3[2], v3[3]
    d1, d2, d3 = v4[1], v4[2], v4[3]
    V =       (b1 - a1) * (c2 * d3 - c3 * d2)
    V = muladd(b2 - a2,    c3 * d1 - c1 * d3, V)
    V = muladd(b3 - a3,    c1 * d2 - c2 * d1, V)
    V = muladd(c1 - d1,    a3 * b2 - a2 * b3, V)
    V = muladd(c2 - d2,    a1 * b3 - a3 * b1, V)
    V = muladd(c3 - d3,    a2 * b1 - a1 * b2, V)
    return V * Float64(1/6)
end

for funName in (:centroid, :volume)
    @eval begin
        function $funName(sv::SVector{4,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3], sv[4])
        end
    end
end


text:
using PyPlot, Dierckx

cubehermitev1(x) = 2*x^3-3*x^2+1
cubehermitem1(x) = x^3-2*x^2+x
cubehermitev2(x) = -2*x^3+3*x^2
cubehermitem2(x) = (x-1)*x^2

@vectorize_1arg Number cubehermitev1
@vectorize_1arg Number cubehermitem1
@vectorize_1arg Number cubehermitev2
@vectorize_1arg Number cubehermitem2

function _plothermite()
    clf()
    x = linspace(0,1,300)
    plot(x, cubehermitev1(x), label="v1", color="r")
    plot(x, cubehermitem1(x), label="m1", color="g")
    plot(x, cubehermitev2(x), label="v2", color="orange")
    plot(x, cubehermitem2(x), label="m2", color="b")
    legend(loc="center right")
end

"""A callable object that acts as a monotone cubic spline, with a
choice of behaviors at the boundaries."""
type MonotoneSpline
    x      ::Vector{Real}
    y      ::Vector{Real}
    m      ::Vector{Real}
    dx     ::Vector{Real}
    bc     ::String
    function MonotoneSpline(x::Vector, y::Vector, m::Vector,
        dx::Vector; bc="error")
        new(x, y, m, dx, bc)
    end
end

function MonotoneSpline(x::Vector, y::Vector; bc="error")

    n=length(x)
    @assert n==length(y)

    secantslopes = zeros(Float64, n-1)
    for i=1:n-1
        secantslopes[i] = (y[i+1]-y[i])/(x[i+1]-x[i])
    end
    dx = x[2:end] .- x[1:end-1]
    m = zeros(Float64, n)
    m[2:n-1] = 0.5*(secantslopes[2:end] + secantslopes[1:end-1])
    m[1] = secantslopes[1]
    m[n] = secantslopes[n-1]

    for k=1:n-1
        if secantslopes[k] == 0
            m[k] = m[k+1] = 0
        else
            alpha = m[k] / secantslopes[k]
            beta = m[k+1] / secantslopes[k]
            if alpha < 0 || beta < 0
                throw(ArgumentError)
            end
            rsq = alpha^2 + beta^2
            if rsq>9
                tau = sqrt(9.0/rsq)
            else
                tau = 1.0
            end
            m[k] = tau*alpha*secantslopes[k]
            m[k+1] = tau*beta*secantslopes[k]
        end
    end
    return MonotoneSpline(x, y, m, dx; bc=bc)
end

function Base.call(spl::MonotoneSpline, z::Number)
    if z < spl.x[1]
        if spl.bc=="error"
            println("Oh no! 1 spl.bc=$(spl.bc)")
            throw(ArgumentError)
        elseif spl.bc=="constant"
            return spl.y[1]
        else
            return spl.y[1]+spl.m[1]*(z-spl.x[1])
        end
    elseif z > spl.x[end]
        if spl.bc=="error"
            println("Oh no! 2 spl.bc=$(spl.bc)")
            throw(ArgumentError)
        elseif spl.bc=="constant"
            return spl.y[end]
        else
            return spl.y[end]+spl.m[end]*(z-spl.x[end])
        end
    end

    bin = searchsortedfirst(spl.x, z)-1
    if bin<1; bin=1; end
    t = (z-spl.x[bin])/spl.dx[bin]
    return (cubehermitev1(t)*spl.y[bin] + cubehermitem1(t)*spl.m[bin]*spl.dx[bin] +
    cubehermitev2(t)*spl.y[bin+1] + cubehermitem2(t)*spl.m[bin+1]*spl.dx[bin])
end
Base.call(spl::MonotoneSpline, z::AbstractArray) = map(spl, z)

"""A callable object that acts as a monotone cubic spline in log-log space,
with a choice of behaviors at the boundaries."""
type MonotoneSplineLogLog
    x      ::Vector{Real}
    y      ::Vector{Real}
    ms     ::MonotoneSpline
    function MonotoneSplineLogLog(x::Vector, y::Vector; bc="error")
        ms = MonotoneSpline(log(x), log(y); bc=bc)
        new(x, y, ms)
    end
end

Base.call(spl::MonotoneSplineLogLog, z::Number) = exp(spl.ms(log(z)))
Base.call(spl::MonotoneSplineLogLog, z::AbstractArray) = map(spl, z)

"""Return the MonotoneSplineLogLog that has the same knots as c1 but passes
through the control points that result from applying c2 on top of c1."""
function compose_splinelog(c1::MonotoneSplineLogLog, c2::MonotoneSplineLogLog)
    x = c1.x
    y = c2(c1(x))
    MonotoneSplineLogLog(x, y; bc="extrapolate")
end

"""Return the MonotoneSpline that has the same knots as c1 but passes
through the control points that result from applying c2 on top of c1."""
function compose_spline(c1::MonotoneSpline, c2::MonotoneSpline)
    x = c1.x
    y = c2(c1(x))
    MonotoneSpline(x, y; bc="extrapolate")
end

function test1()
    x = collect(0:9)
    y = [0,1,3,10,10,11,12,14,14,16]
    f = MonotoneSpline(x, y; bc="extrapolate")
    f2 = Dierckx.Spline1D(float(x), float(y))
    a = linspace(-2,10,500)
    clf()
    plot(x,y,"ok",label="Input data")
    plot(a, f(a), "r", label="Monotone spline")
    plot(a, evaluate(f2,a), "b", label="Dierckx.Spline1D")
    legend(loc="lower right")
    return f
end

function test2()
    model(x) = 1000*((x./1000).^0.6 ) + x
    x = collect(linspace(0,8000,12))
    y = model(x)
    f = MonotoneSpline(x, y; bc="extrapolate")
    g = MonotoneSplineLogLog(x, y; bc="extrapolate")
    s(z) = evaluate(Dierckx.Spline1D(x, y), z)
    clf()
    a = linspace(0,8000,600)
    subplot(211)
    plot(a, model(a), "k", lw=2)
    plot(a, f(a), "r")
    plot(a, g(a), "b")
    plot(a, s(a), "g")
    plot(x,y,"ok")

    subplot(212)
    plot(a, model(a)-f(a), "r", label="Monotone spline")
    plot(a, model(a)-g(a), "b", label="Log-space mono spline")
    plot(a, model(a)-s(a), "g", label="Dierckx.Spline1D")
    legend()
end


text:
module JBEM

include("BEMMesh.jl")

end # module


text:
for i = true : 5
	println("yolo")
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# 24x3 array
data = [ 0.837698   0.49452   2.54352
        -0.294096  -0.39636   0.728619
        -1.62089   -0.44919   1.20592
        -1.06458   -0.68214  -1.12841
         2.14341    0.7309    0.644968
        -0.284139  -1.133     1.98615
         1.19879    2.55633  -0.526461
        -0.032277   0.11701  -0.249265
        -1.02516   -0.44665   2.50556
        -0.515272  -0.578     0.515139
         0.259474  -1.24193   0.105051
         0.178546  -0.80547  -0.016838
        -0.607696  -0.21319  -1.40657
         0.372248   0.93341  -0.667086
        -0.099814   0.52698  -0.253867
         0.743166  -0.79375   2.11131
         0.109262  -1.28021  -0.415184
         0.499346  -0.95897  -2.24336
        -0.191825  -0.59756  -0.63292
        -1.98255   -1.5936   -0.935766
        -0.317612   1.33143  -0.46866
         0.666652  -0.81507   0.370959
        -0.761136   0.10966  -0.997161
        -1.09972    0.28247  -0.846566 ]

wt = [ 0.4 1.0 0.6 ]

uniform_weights = normalize(ones(size(data,1)))

rwt = [ 91
        65
        57
        34
        83
        70
        74
        26
        37
        21
        62
        22
        91
        83
        89
        39
        20
        83
        36
        75
        43
        78
        83
        17 ]
w = normalize(convert(Array{Float64}, rwt))

# 50x25 array
bigdata = [
    0.231  0.059  0.957  0.172  0.779  0.093  0.029  0.643  0.104  0.468  0.315  0.209  0.582  0.044  0.309  0.902  0.528  0.974  0.117  0.468  0.319  0.706  0.648  0.368  0.447
    0.828  0.791  0.316  0.88   0.293  0.359  0.984  0.109  0.439  0.495  0.088  0.13   0.563  0.448  0.728  0.683  0.655  0.26   0.915  0.43   0.756  0.157  0.616  0.758  0.613
    0.949  0.904  0.083  0.591  0.185  0.942  0.215  0.504  0.745  0.219  0.899  0.61   0.335  0.85   0.521  0.266  0.933  0.124  0.401  0.558  0.737  0.011  0.384  0.785  0.031
    0.203  0.038  0.285  0.13   0.174  0.344  0.646  0.065  0.384  0.84   0.844  0.235  0.959  0.648  0.537  0.514  0.423  0.19   0.855  0.72   0.195  0.009  0.775  0.565  0.177
    0.91   0.336  0.797  0.165  0.498  0.488  0.649  0.571  0.31   0.189  0.189  0.519  0.949  0.993  0.372  0.202  0.231  0.392  0.631  0.53   0.75   0.561  0.36   0.698  0.857
    0.525  0.846  0.231  0.301  0.38   0.469  0.934  0.422  0.732  0.166  0.528  0.176  0.402  0.654  0.513  0.164  0.27   0.123  0.187  0.937  0.188  0.93   0.828  0.521  0.86
    0.552  0.597  0.343  0.952  0.993  0.749  0.086  0.549  0.295  0.344  0.723  0.431  0.33   0.14   0.091  0.336  0.419  0.262  0.071  0.368  0.172  0.305  0.989  0.425  0.649
    0.4    0.265  0.412  0.844  0.923  0.83   0.99   0.815  0.942  0.028  0.655  0.593  0.981  0.36   0.753  0.863  0.782  0.558  0.352  0.818  0.438  0.989  0.609  0.926  0.696
    0.961  0.658  0.202  0.362  0.153  0.443  0.48   0.297  0.075  0.381  0.923  0.496  0.068  0.752  0.12   0.802  0.279  0.979  0.61   0.408  0.076  0.709  0.816  0.97   0.872
    0.915  0.179  0.402  0.696  0.003  0.332  0.33   0.91   0.318  0.573  0.635  0.291  0.8    0.515  0.305  0.531  0.607  0.133  0.041  0.725  0.911  0.413  0.256  0.463  0.506
    0.643  0.276  0.788  0.951  0.208  0.804  0.263  0.426  0.575  0.851  0.347  0.902  0.969  0.775  0.628  0.054  0.382  0.478  0.25   0.532  0.432  0.578  0.11   0.203  0.479
    0.479  0.387  0.059  0.287  0.065  0.295  0.527  0.003  0.605  0.187  0.794  0.18   0.396  0.228  0.363  0.099  0.145  0.703  0.894  0.603  0.141  0.099  0.172  0.89   0.546
    0.18   0.921  0.625  0.656  0.865  0.296  0.124  0.349  0.551  0.193  0.494  0.325  0.013  0.14   0.514  0.251  0.115  0.09   0.307  0.015  0.315  0.373  0.71   0.193  0.85
    0.39   0.978  0.906  0.319  0.873  0.368  0.919  0.004  0.045  0.74   0.806  0.453  0.76   0.219  0.869  0.164  0.809  0.472  0.975  0.201  0.454  0.954  0.429  0.569  0.913
    0.094  0.622  0.213  0.197  0.817  0.941  0.056  0.605  0.122  0.435  0.636  0.823  0.038  0.388  0.592  0.27   0.418  0.862  0.964  0.803  0.466  0.49   0.592  0.831  0.905
    0.886  0.169  0.444  0.011  0.663  0.293  0.613  0.511  0.694  0.286  0.401  0.584  0.194  0.882  0.541  0.228  0.053  0.724  0.039  0.82   0.537  0.955  0.467  0.31   0.867
    0.735  0.282  0.53   0.397  0.961  0.297  0.065  0.201  0.083  0.932  0.09   0.363  0.53   0.33   0.967  0.6    0.794  0.36   0.065  0.297  0.619  0.954  0.553  0.313  0.543
    0.695  0.247  0.7    0.268  0.211  0.198  0.155  0.567  0.876  0.997  0.103  0.228  0.304  0.698  0.898  0.771  0.83   0.663  0.541  0.94   0.252  0.88   0.049  0.247  0.3
    0.917  0.214  0.162  0.285  0.929  0.178  0.501  0.321  0.606  0.621  0.238  0.745  0.494  0.774  0.928  0.66   0.193  0.917  0.213  0.497  0.061  0.204  0.326  0.371  0.751
    0.733  0.791  0.396  0.067  0.905  0.481  0.246  0.542  0.633  0.352  0.133  0.227  0.1    0.214  0.115  0.474  0.307  0.53   0.132  0.684  0.168  0.813  0.476  0.865  0.988
    0.771  0.507  0.799  0.744  0.971  0.706  0.95   0.714  0.663  0.708  0.025  0.568  0.683  0.689  0.755  0.746  0.759  0.404  0.564  0.201  0.59   0.59   0.289  0.593  0.2
    0.796  0.09   0.129  0.104  0.283  0.998  0.079  0.86   0.734  0.525  0.934  0.178  0.436  0.996  0.149  0.808  0.822  0.961  0.287  0.428  0.597  0.661  0.263  0.429  0.887
    0.237  0.238  0.366  0.796  0.427  0.849  0.899  0.814  0.571  0.992  0.521  0.669  0.77   0.975  0.465  0.313  0.848  0.527  0.859  0.307  0.005  0.756  0.79   0.383  0.141
    0.109  0.418  0.831  0.863  0.221  0.705  0.91   0.731  0.132  0.332  0.689  0.233  0.695  0.174  0.216  0.781  0.189  0.301  0.842  0.423  0.322  0.16   0.754  0.573  0.243
    0.412  0.59   0.222  0.151  0.653  0.928  0.858  0.292  0.592  0.235  0.574  0.279  0.579  0.595  0.258  0.37   0.948  0.313  0.046  0.341  0.72   0.715  0.298  0.086  0.598
    0.278  0.276  0.966  0.497  0.41   0.979  0.923  0.249  0.326  0.544  0.38   0.505  0.369  0.329  0.761  0.957  0.757  0.632  0.884  0.537  0.435  0.054  0.676  0.145  0.459
    0.369  0.706  0.034  0.595  0.788  0.652  0.15   0.219  0.576  0.196  0.079  0.212  0.465  0.636  0.758  0.244  0.403  0.022  0.748  0.634  0.484  0.466  0.324  0.083  0.291
    0.105  0.179  0.385  0.055  0.829  0.311  0.593  0.911  0.761  0.23   0.597  0.851  0.586  0.133  0.034  0.831  0.378  0.342  0.057  0.428  0.778  0.625  0.628  0.927  0.959
    0.853  0.543  0.402  0.721  0.109  0.947  0.053  0.577  0.623  0.133  0.598  0.72   0.7    0.759  0.763  0.556  0.825  0.112  0.835  0.14   0.427  0.173  0.025  0.09   0.126
    0.055  0.973  0.188  0.014  0.523  0.656  0.11   0.332  0.819  0.086  0.287  0.13   0.917  0.3    0.519  0.171  0.646  0.014  0.376  0.462  0.865  0.771  0.741  0.08   0.105
    0.118  0.822  0.953  0.399  0.244  0.026  0.438  0.3    0.765  0.179  0.621  0.259  0.781  0.454  0.347  0.121  0.203  0.487  0.677  0.982  0.453  0.035  0.164  0.378  0.828
    0.809  0.225  0.838  0.32   0.225  0.672  0.209  0.71   0.06   0.474  0.746  0.126  0.63   0.928  0.003  0.156  0.073  0.194  0.833  0.235  0.035  0.94   0.858  0.785  0.351
    0.025  0.172  0.691  0.962  0.701  0.943  0.601  0.332  0.191  0.201  0.764  0.517  0.127  0.706  0.473  0.579  0.136  0.579  0.184  0.961  0.601  0.697  0.12   0.26   0.506
    0.395  0.715  0.395  0.307  0.16   0.497  0.876  0.507  0.055  0.787  0.353  0.367  0.691  0.003  0.887  0.942  0.642  0.815  0.004  0.102  0.942  0.5    0.609  0.694  0.698
    0.367  0.631  0.486  0.475  0.931  0.44   0.816  0.493  0.231  0.12   0.637  0.291  0.321  0.839  0.689  0.822  0.444  0.433  0.11   0.327  0.858  0.758  0.986  0.519  0.794
    0.717  0.631  0.143  0.525  0.314  0.283  0.324  0.595  0.398  0.617  0.996  0.915  0.604  0.721  0.178  0.345  0.149  0.383  0.215  0.052  0.727  0.651  0.894  0.549  0.51
    0.16   0.045  0.436  0.674  0.467  0.038  0.242  0.899  0.526  0.097  0.871  0.057  0.845  0.486  0.378  0.48   0.532  0.735  0.967  0.4    0.853  0.578  0.937  0.432  0.793
    0.342  0.509  0.557  0.241  0.151  0.405  0.668  0.255  0.3    0.937  0.192  0.204  0.516  0.244  0.379  0.473  0.34   0.213  0.809  0.335  0.536  0.777  0.167  0.724  0.253
    0.516  0.925  0.96   0.712  0.57   0.793  0.528  0.912  0.739  0.477  0.365  0.306  0.969  0.368  0.18   0.703  0.425  0.525  0.889  0.553  0.325  0.612  0.264  0.607  0.185
    0.691  0.251  0.734  0.521  0.615  0.529  0.091  0.928  0.123  0.465  0.38   0.598  0.34   0.67   0.192  0.786  0.585  0.743  0.487  0.445  0.38   0.193  0.036  0.992  0.99
    0.158  0.038  0.621  0.153  0.682  0.747  0.205  0.717  0.209  0.265  0.106  0.626  0.288  0.762  0.877  0.023  0.046  0.312  0.455  0.861  0.118  0.476  0.647  0.03   0.682
    0.546  0.542  0.957  0.541  0.629  0.477  0.341  0.254  0.772  0.263  0.797  0.235  0.875  0.524  0.107  0.669  0.63   0.083  0.399  0.064  0.268  0.917  0.129  0.12   0.113
    0.858  0.371  0.171  0.868  0.271  0.934  0.72   0.675  0.017  0.299  0.666  0.47   0.541  0.627  0.364  0.579  0.876  0.223  0.441  0.503  0.979  0.832  0.567  0.075  0.037
    0.587  0.633  0.244  0.791  0.257  0.782  0.66   0.319  0.227  0.371  0.817  0.782  0.938  0.831  0.831  0.135  0.504  0.437  0.785  0.128  0.594  0.799  0.985  0.756  0.471
    0.827  0.408  0.469  0.521  0.712  0.337  0.199  0.127  0.034  0.884  0.003  0.045  0.604  0.433  0.69   0.234  0.429  0.749  0.045  0.663  0.066  0.714  0.152  0.92   0.354
    0.157  0.239  0.192  0.041  0.342  0.585  0.033  0.937  0.682  0.354  0.978  0.282  0.952  0.363  0.572  0.652  0.06   0.748  0.916  0.734  0.073  0.325  0.386  0.576  0.664
    0.184  0.877  0.863  0.514  0.115  0.169  0.713  0.359  0.171  0.755  0.054  0.069  0.925  0.904  0.541  0.359  0.714  0.011  0.15   0.278  0.726  0.72   0.591  0.132  0.195
    0.236  0.436  0.73   0.179  0.953  0.533  0.69   0.28   0.713  0.353  0.596  0.49   0.774  0.343  0.712  0.128  0.026  0.231  0.018  0.053  0.566  0.107  0.821  0.966  0.516
    0.74   0.989  0.963  0.24   0.602  0.354  0.901  0.454  0.671  0.556  0.565  0.513  0.255  0.036  0.055  0.918  0.142  0.973  0.382  0.155  0.346  0.701  0.201  0.26   0.149
    0.468  0.181  0.096  0.156  0.14   0.302  0.87   0.352  0.162  0.01   0.87   0.005  0.207  0.693  0.74   0.992  0.017  0.349  0.962  0.575  0.539  0.812  0.868  0.886  0.35
]

w1 = ones(size(data,1))
w1[7] = 3

repl = [ 0.837698   0.49452   2.54352
        -0.294096  -0.39636   0.728619
        -1.62089   -0.44919   1.20592
        -1.06458   -0.68214  -1.12841
         2.14341    0.7309    0.644968
        -0.284139  -1.133     1.98615
         1.19879    2.55633  -0.526461
         1.19879    2.55633  -0.526461
         1.19879    2.55633  -0.526461
        -0.032277   0.11701  -0.249265
        -1.02516   -0.44665   2.50556
        -0.515272  -0.578     0.515139
         0.259474  -1.24193   0.105051
         0.178546  -0.80547  -0.016838
        -0.607696  -0.21319  -1.40657
         0.372248   0.93341  -0.667086
        -0.099814   0.52698  -0.253867
         0.743166  -0.79375   2.11131
         0.109262  -1.28021  -0.415184
         0.499346  -0.95897  -2.24336
        -0.191825  -0.59756  -0.63292
        -1.98255   -1.5936   -0.935766
        -0.317612   1.33143  -0.46866
         0.666652  -0.81507   0.370959
        -0.761136   0.10966  -0.997161
        -1.09972    0.28247  -0.846566 ]

rw = ones(size(repl,1))

nan_test_data = [ 0.0  -1.0   1.0   0.0  -1.0  -1.0   1.0  -1.0  -1.0  -1.0
                  1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0   0.0   1.0   1.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  1.0  -1.0   0.0   1.0   0.0   1.0   0.0   0.0   0.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  0.0   1.0  -1.0   0.0  -1.0   1.0   1.0   1.0   1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  0.0   1.0   0.0   0.0   0.0   0.0   1.0   0.0  -1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                 -1.0  -1.0  -1.0   0.0   0.0   1.0   1.0   0.0   1.0  -1.0
                  0.0   0.0   1.0   1.0   1.0  -1.0   0.0  -1.0   0.0  -1.0
                 -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  0.0   1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0  -1.0  -1.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  0.0   0.0   0.0   0.0   0.0   0.0   1.0   1.0   1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                 -1.0   1.0   1.0   0.0  -1.0   1.0   0.0   1.0   0.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  1.0   0.0   0.0   1.0   0.0   1.0   0.0  -1.0   0.0  -1.0
                 -1.0   1.0   1.0  -1.0  -1.0  -1.0   1.0   0.0  -1.0   0.0
                  1.0  -1.0   0.0   1.0   0.0   1.0   0.0   0.0   0.0   0.0 ]

nan_test_wt = [0.02831793460992489,0.004744050454818514,0.02822019260230257,0.005236932428348762,0.0492706652801098,0.006454522058030116,0.12219757929079839,0.044007167153604775,0.1322092637943642,0.018847091843058467,0.03117301149825311,0.027270955184623146,0.04017412455447496,0.0005958465507276494,0.019037064090875612,0.050632029210891875,0.07630930567408355,0.009961180834579806,0.02822019260230257,0.041029418374897736,0.009204977191773241,0.08762854701121198,0.023362497542365402,0.07330510717709664,0.04259034298648226]


text:
function bfunc(x)
    return x^2
end


text:
### Exports ###

export SynapticConnection, SynapticInput

### Types ###

mutable struct SynapticConnection{D}
  uuid::UUID
  op::Symbol
  data::D
end

# FIXME: Support multiple dimensions
@with_kw mutable struct SynapticInput{Frontend,LearnAlg,Mod,NDims}
  inputSize::Int
  outputSize::Int

  condRate::Float32 = 0.1
  traceRate::Float32 = 0.5

  # FIXME: Don't hard-code 16
  frontend::Frontend = GenericFrontend(inputSize, 16)
  learn::LearnAlg = HebbianDecayLearn()
  modulator::Mod = nothing

  C::Array{Float32,NDims} = zeros(Float32, outputSize, inputSize)
  W::Array{Float32,NDims} = rand(Float32, outputSize, inputSize)
  M::Array{Bool,NDims} = ones(Bool, outputSize, inputSize)
  T::Array{Float32,NDims} = zeros(Float32, outputSize, inputSize)
  O::Vector{Float32} = zeros(Float32, outputSize)
end
# FIXME: Apply conf
function SynapticInput(dstcont::CPUContainer, outputSize, conf)
  dstnode = root(dstcont)
  if haskey(conf, :modulator)
    modulator = conf[:modulator]
  else
    modulator = nothing
  end
  SynapticInput(;
    inputSize=first(size(dstnode)),
    outputSize=outputSize,
    modulator=modulator
  )
end

### Methods ###

function reinit!(conn::SynapticInput)
  @unpack frontend, C, W, T = conn
  reinit!(frontend)
  fill!(C, 0f0)
  rand!(W, 0f0:0.01f0:1f0)
  fill!(T, 0f0)
end

### Includes ###

include("learn.jl")
include("modulate.jl")
include("generic.jl")
include("conv.jl")


text:
module utils

function to_categorical(T::Type, a::Vector, num_classes)::Matrix{T}
    num_instances = length(a)
    cat_a = zeros(T, num_classes, num_instances)

    for (class, col) in zip(a, eachcol(cat_a))
        col[class + 1] = 1  # stupid as hell, needs to be reworked
    end

    return cat_a
end

to_categorical(a, num_classes) = to_categorical(Float32, a, num_classes)

end

max_stars_repo_name: prenc/MiraculousNN

max_issues_repo_name: prenc/MiraculousNN

max_forks_repo_name: prenc/MiraculousNN


text:
using Documenter, ColorBlendModes

module CompositingExamples
    using ColorBlendModes, PNGFiles

    export load, save # just for examples
    export generate

    images = Dict{String, AbstractMatrix}()

    function load(filename)
        img = get(images, filename, nothing)
        if img === nothing
            img = PNGFiles.load(joinpath("assets", filename))
            images[filename] = img
        end
        img
    end

    function save(filename, image)
        PNGFiles.save(joinpath("assets", filename), image)
        nothing
    end

    function generate(bm::BlendMode)
        grad_256  = load("grad_256.png")
        juliadots = load("juliadots.png")
        wave_256  = load("wave_256.png")
        out = bm.(bm.(grad_256, juliadots), wave_256, opacity=0.4)
        save(keyword(bm) * ".png", out)
    end

    function generate(op::CompositeOperation, bm::BlendMode)
        blue  = load("blue.png")
        green = load("green.png")
        out = blend.(blue, green, mode=bm, op=op)
        save(keyword(op) * "_" * keyword(bm) * ".png", out)
    end
end

DocMeta.setdocmeta!(ColorBlendModes, :DocTestSetup,
                    :(using ColorBlendModes, ColorTypes, FixedPointNumbers);
                    recursive=true)

makedocs(
    clean = false,
    modules = [ColorBlendModes],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                             assets = ["assets/mix.css", "assets/favicon.ico"]),
    checkdocs = :exports,
    sitename = "ColorBlendModes",
    pages    = Any[
        "Introduction" => "index.md",
        "Blending and Compositing" => "blending-and-compositing.md",
        "Blend Modes" => "blend-modes.md",
        "Composite Operations" => "composite-operations.md",
        "Color Space Dependence" => "color-space-dependence.md",
        "Utility Functions" => "utility-functions.md",
        "Index" => "function-index.md",
        ]
    )

deploydocs(
    repo = "github.com/kimikage/ColorBlendModes.jl.git",
    target = "build")


text:
export SizedStream

mutable struct SizedStream <: IO
    source::IO
    size::UInt64
    nread::UInt64
    ondone::Function
    onstart::Function
end

function SizedStream(source, size; ondone::Function=NOOP, onstart::Function=NOOP)
    SizedStream(source, size, 0, ondone, onstart)
end

function Base.read(s::SizedStream, ::Type{UInt8})
    if eof(s)
        throw(EOFError())
    end
    if s.nread == 0
        s.onstart()
    end
    s.nread += 1
    b = read(s.source, UInt8)
    if s.nread == s.size
        s.ondone()
    end
    b
end

function Base.eof(s::SizedStream)
    return eof(s.source) || s.nread >= s.size
end

function Base.isopen(s::SizedStream)
    return isopen(s.source)
end

function Base.close(s::SizedStream)
    nothing
end