text:
export ⊕,depiece,pieces,PiecewiseSpace

## Interlace interator
# this interlaces the entries treating everything equally,
# but will stop interlacing finite dimensional when it runs
# out.
# Dimensions are either Int or ∞

immutable InterlaceIterator{DMS<:Tuple}
    dimensions::DMS
end

InterlaceIterator(V::Vector) = InterlaceIterator(tuple(V...))

Base.start(it::InterlaceIterator) = (findfirst(d-> d≠0,it.dimensions),1)
function Base.next(it::InterlaceIterator,st)
    if st[1] == length(it.dimensions)
        k = 1
        j = st[2]+1
    else
        k = st[1]+1
        j = st[2]
    end
    nst=(k,j)
    if done(it,nst)
        (st,nst)
    else
        # call next if we are not done
        (st,j > it.dimensions[k]?next(it,nst)[2]:nst)
    end
end

# are all Ints, so finite dimensional
function Base.done{N}(it::InterlaceIterator{NTuple{N,Int}},st)
    for k=1:length(it.dimensions)
        if st[2] ≤ it.dimensions[k] && k ≤ st[1]
            return false
        end
    end
    return true
end

# Are not all Ints, so ∞ dimensional
Base.done(it::InterlaceIterator,st) = false
Base.eltype(it::InterlaceIterator) = Tuple{Int,Int}

## SumSpace encodes a space that can be decoupled as f(x) = a(x) + b(x) where a is in S and b is in V

abstract DirectSumSpace{SV,T,DD,d} <: Space{T,DD,d}

dimension(sp::DirectSumSpace) = mapreduce(dimension,+,sp.spaces)

spaces(s::Space) = (s,)
spaces(sp::DirectSumSpace) = sp.spaces

InterlaceIterator(sp::DirectSumSpace) = InterlaceIterator(map(dimension,sp.spaces))
interlacer(sp::DirectSumSpace) = InterlaceIterator(sp)
interlacer(sp::Space) = InterlaceIterator(tuple(dimension(sp)))
cache(Q::InterlaceIterator) = CachedIterator(Q)

for TYP in (:SumSpace,:TupleSpace)
    @eval begin
        immutable $TYP{SV,T,DD,d} <: DirectSumSpace{SV,T,DD,d}
            spaces::SV
            $TYP(dom::Domain) = new(tuple(map(typ->typ(dom),SV.parameters)...))
            $TYP(sp::Tuple) = new(sp)
        end

        $TYP(sp::Tuple) = $TYP{typeof(sp),mapreduce(basistype,promote_type,sp),
                             typeof(domain(first(sp))),domaindimension(first(sp))}(sp)
    end
end

immutable PiecewiseSpace{SV,T,DD<:UnionDomain,d} <: DirectSumSpace{SV,T,DD,d}
    spaces::SV
    PiecewiseSpace(dom::AnyDomain)=new(tuple(map(typ->typ(dom),SV.parameters)...))
    PiecewiseSpace(dom::UnionDomain)=new(tuple(map((typ,dom)->typ(dom),SV.parameters,dom.domains)...))
    PiecewiseSpace(sp::Tuple)=new(sp)
end

function PiecewiseSpace(spin::Tuple)
    sp=tuple(union(spin)...)  # remove duplicates

    PiecewiseSpace{typeof(sp),mapreduce(basistype,promote_type,sp),
                               typeof(UnionDomain(map(domain,sp))),domaindimension(first(sp))}(sp)
end

for TYP in (:SumSpace,:TupleSpace,:PiecewiseSpace)
    @eval begin
        $TYP(A::$TYP,B::$TYP)=$TYP(tuple(A.spaces...,B.spaces...))

        $TYP(A::Space,B::$TYP)=$TYP(tuple(A,B.spaces...))
        $TYP(A::$TYP,B::Space)=$TYP(tuple(A.spaces...,B))
        $TYP(A::Space...)=$TYP(A)
        $TYP(sp::Array)=$TYP(tuple(sp...))

        canonicalspace(A::$TYP)=$TYP(sort([A.spaces...]))
    end
end

for TYP in (:SumSpace,:TupleSpace)
    @eval setdomain(A::$TYP,d::Domain)=$TYP(map(sp->setdomain(sp,d),A.spaces))
end

setdomain(A::PiecewiseSpace,d::UnionDomain)=PiecewiseSpace(map((sp,dd)->setdomain(sp,dd),A.spaces,d.domains))

function spacescompatible{S<:DirectSumSpace}(A::S,B::S)
    if length(A) != length(B)
        false
    else
        ret=true
        for k=1:length(A)
            if !spacescompatible(A[k],B[k])
                return false
            end
        end
        true
    end
end

function Base.promote_rule{SV,B,DD,d,V,T<:Number}(::Type{Fun{SumSpace{SV,B,DD,d},V}},::Type{T})
    for k=1:length(SV.parameters)
        pt=promote_type(Fun{SV.parameters[k],V},T)
        if pt != Fun
            return Fun{SumSpace{Tuple{SV.parameters[1:k-1]...,pt.parameters[1],SV.parameters[k+1:end]...},
                       B,DD,d},promote_type(V,T)}
        end
    end
    Fun
end

Base.promote_rule{SV,B,DD,d,T<:Number}(::Type{Fun{SumSpace{SV,B,DD,d}}},::Type{T})=promote_rule(Fun{SumSpace{SV,B,DD,d},Float64},T)

function Base.promote_rule{SV,B,DD,d,V,T<:Number}(::Type{Fun{PiecewiseSpace{SV,B,DD,d},V}},::Type{T})
    # if any doesn't support promoting, just leave unpromoted

    newfsp=map(s->promote_type(Fun{s,V},T),SV.parameters)
    if any(s->s==Fun,newfsp)
        Fun
    else
        newsp=map(s->s.parameters[1],newfsp)
        Fun{PiecewiseSpace{Tuple{newsp...},B,DD,d},promote_type(V,T)}
    end
end

Base.promote_rule{SV,B,DD,d,T<:Number}(::Type{Fun{PiecewiseSpace{SV,B,DD,d}}},::Type{T})=promote_rule(Fun{PiecewiseSpace{SV,B,DD,d},Float64},T)

for OP in (:(Base.getindex),:(Base.length),:(Base.next),:(Base.start),:(Base.done),:(Base.endof))
    @eval $OP(S::DirectSumSpace,k...)=$OP(S.spaces,k...)
end

#support tuple set
for OP in (:(Base.start),:(Base.done),:(Base.endof))
    @eval $OP{SS<:DirectSumSpace}(f::Fun{SS},k...)=$OP(space(f),k...)
end

Base.next{SS<:DirectSumSpace}(f::Fun{SS},k)=f[k],k+1

for TYP in (:SumSpace,:TupleSpace)
    @eval domain(A::$TYP)=domain(A.spaces[end])      # TODO: this assumes all spaces have the same domain
                                                     #        we use end to avoid ConstantSpace
end

Space(d::UnionDomain)=PiecewiseSpace(map(Space,d.domains))
domain(S::PiecewiseSpace)=UnionDomain(map(domain,S.spaces))

function spacescompatible(A::Tuple,B::Tuple)
    if length(A) != length(B)
        return false
    end
    #assumes domain doesn't impact sorting
    asort=sort([A...]);bsort=sort([B...])
    for k=1:length(asort)
        if !spacescompatible(asort[k],bsort[k])
            return false
        end
    end

    return true
end

function union_rule(A::SumSpace,B::Space)
    if !domainscompatible(A,B)
        NoSpace()
    else
        for sp in A.spaces
            if isconvertible(B,sp)
                return A
            end
        end
        SumSpace(A,B)
    end
end

## evaluate

for OP in (:(Base.last),:(Base.first))
    @eval begin
        $OP{SS<:SumSpace}(f::Fun{SS})=mapreduce($OP,+,vec(f))
        $OP{SS<:PiecewiseSpace}(f::Fun{SS})=$OP($OP(vec(f)))
        $OP{SS<:TupleSpace}(f::Fun{SS})=$OP(vec(f))
    end
end

evaluate(f::AbstractVector,S::SumSpace,x)=mapreduce(vf->evaluate(vf,x),+,vec(Fun(f,S)))

function evaluate(f::AbstractVector,S::PiecewiseSpace,x)
    d=domain(S)
    g=Fun(f,S)

#    ret=zero(promote_type(eltype(f),eltype(S)))
    for k=1:numpieces(d)
        if in(x,d[k])
            return g[k](x)
        end
    end
    return 0*first(g)
end

function evaluate(v::AbstractVector,S::PiecewiseSpace,x::AbstractVector)
    f=Fun(v,S)
    [f(xk) for xk in x]
end

function evaluate(v::AbstractVector,S::TupleSpace,x...)
    f=Fun(v,S)
    eltype(f)[f[k](x...) for k=1:length(f.space)]
end

## calculus
for TYP in (:SumSpace,:TupleSpace,:PiecewiseSpace)
    for OP in (:differentiate,:integrate)
        @eval function $OP{D<:$TYP,T}(f::Fun{D,T})
            fs=map($OP,vec(f))
            sp=$TYP(map(space,fs))
            Fun(interlace(fs,sp),sp)
        end
    end
    for OP in (:(Base.real),:(Base.imag),:(Base.conj))
        @eval begin
            $OP{SV,DD,d}(f::Fun{$TYP{SV,RealBasis,DD,d}}) = Fun($OP(f.coefficients),f.space)
            function $OP{SV,T,DD,d}(f::Fun{$TYP{SV,T,DD,d}})
                fs=map($OP,vec(f))
                sp=$TYP(map(space,fs))
                Fun(interlace(fs,sp),sp)
            end
        end
    end
end

for TYP in (:SumSpace,:TupleSpace)
    @eval function Base.cumsum{D<:$TYP,T}(f::Fun{D,T})
        fs=map(cumsum,vec(f))
        sp=$TYP(map(space,fs))
        Fun(interlace(fs,sp),sp)
    end
end

for TYP in (:SumSpace,:PiecewiseSpace), OP in (:(Base.sum),:linesum)
    @eval $OP{V<:$TYP}(f::Fun{V})=mapreduce($OP,+,vec(f))
end

function Base.cumsum{V<:PiecewiseSpace}(f::Fun{V})
    vf=pieces(f)
    r=zero(eltype(f))
    for k=1:length(vf)
        vf[k]=cumsum(vf[k]) + r
        r=last(vf[k])
    end
    depiece(vf)
end

Base.cumsum{V<:PiecewiseSpace}(f::Fun{V},d::Domain)=mapreduce(g->cumsum(g,d),+,pieces(f))

bilinearform{S<:PiecewiseSpace,V<:PiecewiseSpace}(f::Fun{S},g::Fun{V}) = sum(map(bilinearform,pieces(f),pieces(g)))
linebilinearform{S<:PiecewiseSpace,V<:PiecewiseSpace}(f::Fun{S},g::Fun{V}) = sum(map(linebilinearform,pieces(f),pieces(g)))

# assume first domain has 1 as a basis element

function Base.ones{T<:Number}(::Type{T},S::SumSpace)
    @assert length(S.spaces)==2
    if isconvertible(ConstantSpace(),S.spaces[1])
        ones(T,S[1])⊕zeros(T,S[2])
    else
        zeros(T,S[1])⊕ones(T,S[2])
    end
end

Base.ones(S::SumSpace)=ones(Float64,S)

Base.ones{T<:Number,SS,V}(::Type{T},S::PiecewiseSpace{SS,V}) = depiece(map(ones,S.spaces))
Base.ones(S::PiecewiseSpace)=ones(Float64,S)

identity_fun(S::PiecewiseSpace) = depiece(map(identity_fun,S.spaces))

# vec

function Base.getindex{DSS<:DirectSumSpace}(f::Fun{DSS},k::Integer)
    sp=space(f).spaces
    it=InterlaceIterator(map(dimension,sp))
    N=length(f.coefficients)
    d=dimension(sp[k])

    # preallocate: we know we have at most N coefficients
    ret=Array(eltype(f),N)
    j=1  # current coefficient
    p=0  # current length
    for (n,m) in it
        if j > N || m > d
            break
        end
        if n==k
            ret[m]=f.coefficients[j]
            p+=1
        end
        j+=1
    end
    resize!(ret,p)  # through out extra coefficients
    Fun(ret,sp[k])
end

# interlace coefficients according to iterator
function interlace{T,V<:AbstractVector}(::Type{T},v::AbstractVector{V},it::InterlaceIterator)
    ret=Vector{T}()
    N=mapreduce(length,max,v)
    for (n,m) in it
        if m > N
            break
        end

        if m ≤ length(v[n])
            push!(ret,v[n][m])
        else
            push!(ret,zero(T))
        end
    end
    ret
end

interlace{V<:AbstractVector}(v::AbstractVector{V},it::InterlaceIterator) =
    interlace(mapreduce(eltype,promote_type,v),v,it)

interlace{F<:Fun}(v::AbstractVector{F},sp::DirectSumSpace) =
    interlace(map(coefficients,v),InterlaceIterator(sp))

function interlace(v::Tuple,sp::DirectSumSpace)
    V=Array(Vector{mapreduce(eltype,promote_type,v)},length(v))
    for k=1:length(v)
        V[k]=coefficients(v[k])
    end
    interlace(V,InterlaceIterator(sp))
end

Base.vec(S::DirectSumSpace) = S.spaces
Base.vec{S<:DirectSumSpace}(f::Fun{S}) = Fun[f[j] for j=1:length(space(f).spaces)]

pieces{S<:PiecewiseSpace}(f::Fun{S}) = vec(f)

for (Dep,Sp) in ((:depiece,:PiecewiseSpace),(:detuple,:TupleSpace))
    @eval begin
        function $Dep{F<:Fun}(v::Vector{F})
            spaces=map(space,v)
            sp=$Sp(spaces)
            Fun(interlace(v,sp),sp)
        end
        function $Dep(v::Tuple)
            spaces=map(space,v)
            sp=$Sp(spaces)
            Fun(interlace(v,sp),sp)
        end

        $Dep(v::Vector{Any})=$Dep(tuple(v...))
    end
end

interlace{FF<:Fun}(f::AbstractVector{FF}) = detuple(f)

# convert a vector to a Fun with TupleSpace
Fun(v::Vector{Any},sp::TupleSpace) = detuple(map(Fun,v,sp.spaces))
Fun{F<:Fun}(v::Vector{F},sp::TupleSpace) = detuple(map(Fun,v,sp.spaces))

## transforms

function points(d::PiecewiseSpace,n)
   k=div(n,length(d))
    r=n-length(d)*k

    [vcat([points(d.spaces[j],k+1) for j=1:r]...);
        vcat([points(d.spaces[j],k) for j=r+1:length(d)]...)]
end

function transform(S::PiecewiseSpace,vals::Vector,plan...)
    n=length(vals)
    K=length(S)
    k=div(n,K)
    PT=coefficient_type(S,eltype(vals))
    if k==0
        ret=Array(PT,n)
        for j=1:n
            ret[j]=transform(S[j],[vals[j]])[1]
        end

        ret
    else
        r=n-K*k
        M=Array(PT,k+1,K)

        for j=1:r
            M[:,j]=transform(S[j],vals[(j-1)*(k+1)+1:j*(k+1)])
        end
        for j=r+1:length(S)
            M[1:k,j]=transform(S[j],vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k])
            M[k+1,j]=zero(PT)
        end

    vec(M.')
    end
end

itransform(S::PiecewiseSpace,cfs::Vector,plan...)=vcat([itransform(S.spaces[j],cfs[j:length(S):end]) for j=1:length(S)]...)

itransform(S::SumSpace,cfs,plan...)=Fun(cfs,S)(points(S,length(cfs)))
itransform!(S::SumSpace,cfs,plan...)=(cfs[:]=Fun(cfs,S)(points(S,length(cfs))))

## SumSpace{ConstantSpace}
# this space is special

union_rule(P::PiecewiseSpace,C::ConstantSpace{AnyDomain})=PiecewiseSpace(map(sp->union(sp,C),P.spaces))

## choosedomainspace

function choosedomainspace{T}(A::InterlaceOperator{T,1},sp::TupleSpace)
    # this ensures correct dispatch for unino
    sps = Vector{Space}(
        filter(x->!isambiguous(x),map(choosedomainspace,A.ops,sp.spaces)))
    if isempty(sps)
        UnsetSpace()
    else
        union(sps...)
    end
end


text:
using LibGit2
using Serialization

unicode_url = "https://github.com/unicode-table/unicode-table-data.git"
local_path = joinpath(@__DIR__, "unicode-data")
serialized_file = joinpath(@__DIR__, "unicode-dict.serial")

# Check if local path already exists. If so, delete is
ispath(local_path) && rm(local_path, recursive = true)
LibGit2.clone(unicode_url, local_path)

# Path to the actual data
symbols_path = joinpath(local_path, "loc", "en", "symbols")

function add!(dict, str::String)
    # Filter out empty lines
    isempty(str) && return

    # Get the caracters up until the first colon, this is the codepoint for this entry
    ind = findfirst(":", str)
    num = parse(UInt32, SubString(str, 1, prevind(str, first(ind))); base = 16)
    info = strip(SubString(str, nextind(str, last(ind)), lastindex(str)))

    # Add entry if it exists
    !isempty(info) && (dict[num] = info)
    return
end

function collect_codepoints(dir)
    dict = Dict{UInt32, String}()
    for (root, _, files) in walkdir(dir), file in files
        open(joinpath(root, file)) do f
            for ln in eachline(f)
                add!(dict, ln)
            end
        end
    end
    return dict
end

# Walk through and collect all of the code points.
codepoints = collect_codepoints(symbols_path)

# Save
open(io -> serialize(io, codepoints), serialized_file; write = true)

ext: jl


text:
module  TestManipulate

using PopGen
using DataFrames
using Test

cats = @nancycats;
sharks = @gulfsharks;

@testset "locations" begin
    x = rand(length(samples(cats))) ; y = rand(length(samples(cats)))
    locations!(cats, long = x, lat = y)
    @test cats.meta.longitude == Float32.(x)
    @test cats.meta.latitude == Float32.(y)
    @test locations(cats).longitude == Float32.(x)
    @test locations(cats).latitude == Float32.(y)
end

@testset "decimal-minutes locations" begin
    x = fill("11 22.33W", length(samples(cats))) ; y = fill("-41 31.52", length(samples(cats)))
    locations!(cats, long = x, lat = y)
    @test all(cats.meta.longitude .== Float32(-11.3722))
    @test all(cats.meta.latitude .== Float32(-41.5253))
end

@testset "loci and genotypes" begin
    @test length(loci(cats)) == 9
    @test get_genotype(cats, sample = "N115", locus = "fca8") == (135,135)
    N115 = get_genotypes(cats, "N115")
    @test length(N115) == 9
    @test typeof(N115) == Vector{Union{Missing, Tuple{Int16,Int16}}}
    @test typeof(get_genotypes(cats, name = "N115" , locus = "fca8")) <: SubDataFrame
    @test names(get_genotypes(cats, name = "N115" , locus = "fca8")) == ["name", "population", "locus", "genotype"]
    @test size(get_genotypes(cats, name = ["N115", "N7"] , locus = "fca8")) == (2,4)
    @test size(get_genotypes(cats, name = "N115" , locus = ["fca8", "fca37"])) == (2,4)
    @test size(get_genotypes(cats, name = ["N115", "N7"] , locus = ["fca8", "fca37"])) == (4,4)
    @test length(genotypes(sharks, "contig_475")) == 212
end

@testset "populations" begin
    @test length(populations(cats)) == 17
    @test typeof(populations(cats, counts = true).population) == Vector{String}
    @test typeof(populations(cats, counts = true).count) == Vector{Int}

    rn_dict = Dict("1" => "one", "2" => "two")
    populations!(cats, rn_dict)
    @test "one" ∈ cats.meta.population && "1" ∉ cats.meta.population
    @test "two" ∈ cats.meta.population && "2" ∉ cats.meta.population

    rn_vect = string.(1:17)
    populations!(cats, rn_vect)
    @test "one" ∉ cats.meta.population && "1" ∈ cats.meta.population
    @test "two" ∉ cats.meta.population && "2" ∈ cats.meta.population

    rn_vect_ii = ["N215", "N297"]
    rn_vect_iin = ["one", "seventeen"]
    populations!(cats, rn_vect_ii, rn_vect_iin)
    @test "one" ∈ cats.meta.population && "seventeen" ∈ cats.meta.population
    @test cats.meta[cats.meta.name .== "N215", :population] == ["one"]
    @test cats.meta[cats.meta.name .== "N297", :population] == ["seventeen"]
end

@testset "exclusion" begin
    tmp = exclude(@nancycats, name = "N100", population = ["1", "15"])
    @test length(samples(tmp)) == 215
    @test length(populations(tmp)) == 15

    tmp = exclude(@nancycats, name = "N102", locus= "fca8", population = "3")
    @test length(loci(tmp)) == 8
    @test length(populations(tmp)) == 16
    @test length(samples(tmp)) == 225
end

end # module


text:
@testset "estimate_img_properties" begin
    img = generate_surface(SynthSphere(50), 0.5, [0,0,1])
    ρ,I,σ,τ = estimate_img_properties(img)
    @test ρ ≈ 0.5369128178710199
    @test I ≈ [0.0, 0.0, 1.0]
    @test σ ≈ 0.0
    @test τ ≈ 1.5707963267948966
end


text:
# Benchmarks the version of transit_poly which uses a Type
# for holding all of the variables related to the transit
# computation.

using Statistics
using PyPlot

# Defines the function that computes the transit:
include("../../../src/transit_poly_struct.jl")

nu = [1,2,3,5,8,13,21,34,55,89,144]; #233,377,610];
nnu = length(nu)
nb = [100,316,1000,3160,10000,31600,100000];
nnb = length(nb)

timing_ratio = zeros(nnb,nnu)

# Define a function that loops over impact parameter
# and defines arrays to hold the results:

function profile_transit_poly(trans,nb)
fgrad = zeros(2+trans.n)
timing = zeros(length(nb))
tmean = zeros(5)
for j=1:length(nb)
  nb0 = nb[j]
  b = zeros(nb0)
  flux = zeros(nb0)
  for k=1:5
    elapsed = time_ns()
    for i=1:nb0
      b[i] = sqrt(((i-float(nb0)/2)*2/float(nb0)*(1.0+2.0*trans.r))^2)
      trans.b = b[i]
      flux[i] = transit_poly!(trans)
    end
    tmean[k] = (time_ns()-elapsed)*1e-9

    #plot(b,flux)
  end
  timing[j] = median(tmean)
end
return timing
end

# Instantiate a transit with r=0.1, b=0.5, u_n=[0.2,0.3]:
r = 0.1; b = 0.5
# Load in PyPlot (Julia matplotlib wrapper):
clf()
cmap = get_cmap("plasma")
for iu = 1:nnu
  u_n = ones(nu[iu])/float(nu[iu])
  trans = transit_init(r,b[1],u_n,true)
# Transform from u_n to c_n coefficients:
  timing = zeros(nnb)
# Call the function:
  @time timing =profile_transit_poly(trans,nb)
  if iu == 1
    timing_ratio[:,1] = timing
  else
    timing_ratio[:,iu] = timing./timing_ratio[:,1]
  end
  loglog(nb,timing,color=cmap(float(iu)/(1.25 * nnu)))
  loglog(nb,timing,"o",color=cmap(float(iu)/(1.25 * nnu)),label=string("n: ",nu[iu]))
end
xlabel("Number of points")
ylabel("Timing [sec]")
legend(loc="upper left",fontsize=8,ncol=3)
timing_ratio[:,1].=1.0
savefig("benchmark_transit_poly.pdf", bbox_inches="tight")

clf()
tmed = vec(median(timing_ratio;dims=1))
loglog(nu,tmed,"o", color="C0")
plot(nu,tmed,label="Measured",linewidth=2, color="C0")
alp = log(tmed[nnu])/log(nu[nnu])
plot(nu,nu.^0.2,linestyle="--",label=L"$n^{0.2}$",color="C1")
plot(nu,tmed[nnu]*(nu/nu[nnu]).^2,linestyle="-.",label=L"$n^2$",color="C1")
xlabel("Number of limb-darkening coefficients")
ylabel("Relative timing")
legend(loc="upper left")
axis([1,144,0.7,2e2])
savefig("benchmark_limbdark_timing.pdf", bbox_inches="tight")

#return
#end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
###
### ToricHomotopy
###
# This is part of the implementation of polyhedral homotopy, namely to reverse the toric
# degeneration from 0 to 1.
# This is only an internal homotopy.

const StructVectorComplexF64 = StructArrays.StructArray{
    Complex{Float64},
    1,
    NamedTuple{(:re, :im),Tuple{Vector{Float64},Vector{Float64}}},
    Int64,
}

Base.@kwdef struct ToricHomotopy{S} <: AbstractHomotopy
    system::ModelKit.CompiledSystem{S}
    system_coeffs::Vector{ComplexF64}
    weights::Vector{Float64}

    #t^w
    t_weights::Vector{Float64}
    #t^w for t complex or negative
    complex_t_weights::StructVectorComplexF64
    coeffs::Vector{ComplexF64}
    dt_coeffs::Vector{ComplexF64}
    x::Vector{ComplexF64}
    # cache value
    t_coeffs::Base.RefValue{ComplexF64}
    t_taylor_coeffs::Base.RefValue{ComplexF64}

    # taylor coefficients
    taylor_coeffs::TaylorVector{5,ComplexF64}
    tc3::TaylorVector{4,ComplexF64}
    tc2::TaylorVector{3,ComplexF64}
end

function ToricHomotopy(
    system::ModelKit.CompiledSystem,
    system_coeffs::Vector{Vector{ComplexF64}},
)
    m = ModelKit.nparameters(system)
    m1 = sum(length, system_coeffs)
    m == m1 ||
        throw(ArgumentError("System parameters and coefficients do not have the same size, got $m and $m1"))
    n = nvariables(system)
    taylor_coeffs = TaylorVector{5}(ComplexF64, m)
    ToricHomotopy(
        system = system,
        system_coeffs = reduce(vcat, system_coeffs),
        weights = zeros(m),
        t_weights = zeros(m),
        complex_t_weights = StructArrays.StructArray(zeros(ComplexF64, m)),
        coeffs = zeros(ComplexF64, m),
        dt_coeffs = zeros(ComplexF64, m),
        x = zeros(ComplexF64, n),
        t_coeffs = Ref(complex(0.0, 0.0)),
        t_taylor_coeffs = Ref(complex(0.0, 0.0)),
        taylor_coeffs = taylor_coeffs,
        tc3 = TaylorVector{4}(taylor_coeffs),
        tc2 = TaylorVector{3}(taylor_coeffs),
    )
end

Base.size(H::ToricHomotopy) = size(H.system)

function update_weights!(
    H::ToricHomotopy,
    support::AbstractVector{<:AbstractMatrix},
    lifting::AbstractVector{<:AbstractVector},
    cell::MixedSubdivisions.MixedCell;
    min_weight::Union{Nothing,Float64} = nothing,
    max_weight::Union{Nothing,Float64} = nothing,
)
    l = 1
    s_max, s_min = 0.0, Inf
    n = length(cell.normal)
    for (i, Aᵢ) in enumerate(support)
        wᵢ = lifting[i]
        βᵢ = cell.β[i]
        mᵢ = size(Aᵢ, 2)
        aᵢ, bᵢ = cell.indices[i]
        for j = 1:mᵢ
            if j == aᵢ || j == bᵢ
                H.weights[l] = 0.0
            else
                sᵢⱼ = wᵢ[j] - βᵢ
                for k = 1:n
                    @inbounds sᵢⱼ += Aᵢ[k, j] * cell.normal[k]
                end
                H.weights[l] = sᵢⱼ
                s_max = max(s_max, sᵢⱼ)
                s_min = min(s_min, sᵢⱼ)
            end
            l += 1
        end
    end

    if min_weight !== nothing
        λ = s_min / min_weight
        H.weights ./= λ
        s_min, s_max = min_weight, s_max / λ
    elseif max_weight !== nothing
        λ = s_max / max_weight
        H.weights ./= λ
        s_min, s_max = s_min / λ, max_weight
    end

    H.t_coeffs[] = NaN
    H.t_taylor_coeffs[] = NaN

    s_min, s_max
end

function evaluate_weights!(u::Vector{Float64}, weights::Vector{Float64}, t::Float64)
    s = log(t)
    n = length(weights)
    LoopVectorization.@avx for i = 1:n
        u[i] = exp(weights[i] * s)
    end
    u
end

function evaluate_weights!(
    u::Vector{Float64},
    v::Vector{Float64},
    weights::Vector{Float64},
    t::ComplexF64,
)
    s = log(t)
    n = length(weights)
    sr, si = reim(s)
    LoopVectorization.@avx for i = 1:n
        r = exp(weights[i] * sr)
        wsi = weights[i] * si
        # TODO: This should actually by s, c = sincos(wsi)
        # However, this it not yet supported by LoopVectorization
        s = sin(wsi)
        c = cos(wsi)
        u[i] = r * c
        v[i] = r * s
    end
end

function coeffs!(H::ToricHomotopy, t)
    # H.t_coeffs[] == t && return H.coeffs
    c = H.coeffs
    if !isreal(t) || real(t) < 0
        evaluate_weights!(
            H.complex_t_weights.re,
            H.complex_t_weights.im,
            H.weights,
            complex(t),
        )
        @inbounds for i = 1:length(c)
            uᵢ = H.system_coeffs[i]
            ctwᵢ = H.complex_t_weights[i]
            c[i] = uᵢ * ctwᵢ
        end
    elseif iszero(t)
        @inbounds for i = 1:length(c)
            wᵢ = H.weights[i]
            uᵢ = H.system_coeffs[i]
            c[i] = ifelse(wᵢ == 0, uᵢ, complex(0.0, 0.0))
        end
    else
        evaluate_weights!(H.t_weights, H.weights, real(t))
        @inbounds for i = 1:length(c)
            uᵢ = H.system_coeffs[i]
            twᵢ = H.t_weights[i]
            c[i] = uᵢ * twᵢ
        end
    end

    H.t_coeffs[] = t
    c
end

function dt_coeffs!(H::ToricHomotopy, t)
    if iszero(t)
        @inbounds for i = 1:length(H.dt_coeffs)
            wᵢ = H.weights[i]
            uᵢ = H.system_coeffs[i]
            H.dt_coeffs[i] = ifelse(wᵢ == 1, uᵢ, complex(0.0, 0.0))
        end
    else
        c = coeffs!(H, t)
        if isreal(t)
            s = real(t)
            @fastmath s⁻¹ = inv(s)
            for i = 1:length(c)
                @inbounds H.dt_coeffs[i] = H.weights[i] * c[i] * s⁻¹
            end
        else
            @fastmath t⁻¹ = inv(t)
            for i = 1:length(c)
                @inbounds H.dt_coeffs[i] = H.weights[i] * c[i] * t⁻¹
            end
        end
    end

    H.dt_coeffs
end

function evaluate!(u, H::ToricHomotopy, x::AbstractVector, t)
    ModelKit.evaluate!(u, H.system, x, coeffs!(H, t))
end

function evaluate_and_jacobian!(u, U, H::ToricHomotopy, x::AbstractVector, t)
    ModelKit.evaluate_and_jacobian!(u, U, H.system, x, coeffs!(H, t))
    nothing
end

function taylor!(u, v::Val{1}, H::ToricHomotopy, x::Vector, t)
    ModelKit.evaluate!(u, H.system, x, dt_coeffs!(H, t))
end

function taylor_coeffs!(H::ToricHomotopy, t::Real)
    H.t_taylor_coeffs[] != t || return H.taylor_coeffs

    tc = H.taylor_coeffs
    # c, c¹, c², c³, c⁴ = H.taylor_coeffs
    if iszero(t)
        for i = 1:length(tc)
            wᵢ = H.weights[i]
            uᵢ = H.system_coeffs[i]
            if wᵢ == 0
                tc[i, 1] = uᵢ
                tc[i, 2] = tc[i, 3] = tc[i, 4] = tc[i, 5] = 0.0
            elseif wᵢ == 1
                tc[i, 2] = uᵢ
                tc[i, 1] = tc[i, 3] = tc[i, 4] = tc[i, 5] = 0.0
            else
                tc[i, 1] = tc[i, 2] = tc[i, 3] = tc[i, 4] = tc[i, 5] = 0.0
            end
        end
    else
        if t !== H.t_coeffs[]
            evaluate_weights!(H.t_weights, H.weights, t)
        end
        t⁻¹ = inv(t)
        @inbounds for i = 1:length(tc)
            wᵢ = H.weights[i]
            uᵢ = H.system_coeffs[i]
            twᵢ = H.t_weights[i]
            tc[i, 1] = uᵢ * twᵢ
            tw1 = wᵢ * twᵢ * t⁻¹
            tc[i, 2] = uᵢ * tw1
            tw2 = 0.5 * (wᵢ - 1) * tw1 * t⁻¹
            tc[i, 3] = uᵢ * tw2
            tw3 = (wᵢ - 2) * tw2 * t⁻¹ / 3
            tc[i, 4] = uᵢ * tw3
            tw4 = 0.25 * (wᵢ - 3) * tw3 * t⁻¹
            tc[i, 5] = uᵢ * tw4
        end

    end

    H.t_taylor_coeffs[] = t

    H.taylor_coeffs
end

function taylor!(u, v::Val{2}, H::ToricHomotopy, tx::TaylorVector, t)
    taylor_coeffs!(H, real(t))
    ModelKit.taylor!(u, v, H.system, tx, H.tc2)
end
function taylor!(u, v::Val{3}, H::ToricHomotopy, tx::TaylorVector, t)
    taylor_coeffs!(H, real(t))
    ModelKit.taylor!(u, v, H.system, tx, H.tc3)
end
function taylor!(u, v::Val{4}, H::ToricHomotopy, tx::TaylorVector, t)
    taylor_coeffs!(H, real(t))
    ModelKit.taylor!(u, v, H.system, tx, H.taylor_coeffs)
end


text:
using MLStyle

export canonical
canonical(x) = x

# TODO: Make sure local variables are handled properly (e.g. local function args)
function canonical(expr :: Expr)
    # @show expr
    r = canonical
    @match expr begin
        :($x |> $f) => begin
            rf = r(f)
            rx = r(x)
            :($rf($rx)) |> r
        end

        Expr(:block, body...) => Expr(:block, canonical.(body)...)

        Expr(:do, :(For($x)), :($f)) => begin
            rf = r(f)
            rx = r(x)

            :(For($rf, $rx)) |> r
        end

        Expr(:do, :(For($(x...))), :($f)) => begin
            rf = r(f)
            rx = r.(x)

            rxtup = Expr(:tuple, rx...)
            :(For($rf, $rxtup)) |> r
        end

        :((iid($n))($dist)) => begin
            rn = r(n)
            rdist = r(dist)
            :(iid($rn,$rdist)) |> r
        end

        :($f($(args...))) => begin
            rf = r(f)
            rx = map(r,args)
            :($rf($(rx...)))
        end

        # TODO: This was intended to work around the closure issues by rewriting it as a local fucntion with all values passed explicitly as arguments. Doesn't seem to work, at least not yet
        # :($x -> begin $lnn; $fbody end) => begin
        #     vs = setdiff(variables(fbody), variables(x))
        #     @gensym f
        #     pars = astuple(variables(x))
        #     pars = :($(Expr(:parameters)))
        #     for v in vs
        #         pushfirst!(pars.args, :($v=$v))
        #     end
        #     quote
        #         $f($pars) = $fbody
        #         $f
        #     end
        # end |> r

        x => x
    end
end

function canonical(m :: Model)
    args = m.args :: Vector{Symbol}
    vals  = map(canonical, m.vals)
    dists = map(canonical, m.dists)
    retn = m.retn
    Model(getmodule(m), args, vals, dists, retn)
end

ex1 = :(map(1:10) do x x^2 end)

# canonical(ex1)
# canonical(linReg1D)


text:
module BisectPy

include("bisect.jl")
include("insort.jl")
include("find.jl")

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d47521c0-25d5-11eb-09e7-1934afa8689d
# I used these packages.

begin
	using PlutoUI
	using Images
end

# ╔═╡ 6ab4e9c0-114c-11eb-2d22-05776bd25df1
md"# Hello!
##### This is the second [Project Euler](https://projecteuler.net/archives) question where I find the sum of the even-valued terms of the Fibonacci Sequence under 4 million!"

# ╔═╡ 1524cf62-114d-11eb-3dc3-4debd151734b
md"# Use the slider to find Fibonacci Terms!"

# ╔═╡ 89b61e60-114d-11eb-0375-4bcab2fef870
md"##### If you slide the slider, you get the Fibonacci term!"

# ╔═╡ 22731370-25d9-11eb-049b-b546c56ec5ad
@bind fibonacciterm Slider(0:12 , show_value=true)

# ╔═╡ f9488b00-1157-11eb-22ed-bdc7f113fbe9
md"# How does it work?????"

# ╔═╡ 41c9f2c0-114d-11eb-0d12-11a7feced986
md"##### I am using Binet's Formula to calculate the terms in the Fibonacci Sequence!"

# ╔═╡ 168eb990-25db-11eb-28db-79e076dc7ea3
md"### Binet's Formula: "

# ╔═╡ f3cfe3e0-25dd-11eb-1fbc-9f9c70a63261
md"# $ \frac{\varphi ^ {n}}{\sqrt{5}} $
## $ \varphi $ is the golden ratio . "

# ╔═╡ 2bc5c920-25db-11eb-000b-bf9df6093842
md"##### I will expain the Golden Ratio in the next notebook. You just need to know that it is about 1.618 and instead of calling it Golden Ratio all the time, we use the Greek letter $ \varphi $."

# ╔═╡ e7abfa92-1147-11eb-2b44-159345577c2a
# The Golden Ratio!
φ = (1 + sqrt(5))/2

# ╔═╡ 3ee0f5de-1148-11eb-2039-93fa7b1be9af
# The Binet's Formula is :
# the golden ratio (~ 1.618) raised to the nth power and divide it by sqrt of 5 and rounding the answer to the nearest integer. The error approaches 0 as you approach the infinitely big terms!

begin
	a = (φ ^ fibonacciterm)/sqrt(5)
	b = round(Int, a)
end

# ╔═╡ 6ef159f2-1157-11eb-04fa-a9542f145bf1
# Grammar is IMPORTANT!!!

if fibonacciterm == 1
	md"##### The $fibonacciterm st term in the Fibonacci Sequence is $b."
elseif fibonacciterm == 2
	md"##### The $fibonacciterm nd term in the Fibonacci Sequence is $b."
elseif fibonacciterm == 3
	md"##### The $fibonacciterm rd term in the Fibonacci Sequence is $b."
else
	md"##### The $fibonacciterm th term in the Fibonacci Sequence is $b."
end

# ╔═╡ a3fef210-1b75-11eb-0311-5348902769b0
md"# I can calculate many Fibonacci Numbers at once!"

# ╔═╡ 96044e72-1b8f-11eb-0bd5-bf6b54c9464d
c = 99

# ╔═╡ a6b1de40-1b8f-11eb-0e31-f17fed3c148b
d = 0 : c

# ╔═╡ b05dc25e-1b8f-11eb-1a64-dd58a16b7dec
e = φ .^ d

# ╔═╡ b7dc3f30-1b8f-11eb-316d-af7911153b99
f = e/sqrt(5)

# ╔═╡ 83d045b0-1b8f-11eb-02e0-317499abf430
fibonacci_numbers = round.(f)

# ╔═╡ 40007ef0-2057-11eb-13c5-f7df94b8a45c
md"# I can just find the even-valued Fibonacci terms"

# ╔═╡ 5981e260-2057-11eb-14d6-eba9727f24a2
md"### Every third Fibonacci term is even!"

# ╔═╡ 651d1360-2057-11eb-3fbb-e7f27cb8358b
g = 100

# ╔═╡ cb4b9070-2058-11eb-3581-8577499eafc7
h = 1 : g

# ╔═╡ 6b0a8890-205a-11eb-2f7a-1fba9e2c34c0
i = φ .^ (3 * h)

# ╔═╡ 81033a70-205a-11eb-284a-3fa820bc460e
j = i / sqrt(5)

# ╔═╡ 9b287f00-205a-11eb-1a54-dba8016126e1
even_valued_fibonacci_numbers = round.(j)

# ╔═╡ c88ce7d0-1b85-11eb-292d-1da5aa66d749
md"# I control the size of the Fibonacci Numbers."

# ╔═╡ 09903500-1b88-11eb-0365-09cc57ba4c8c
md"##### I am just setting the domain for the Fibonacci Numbers because I do not want them to be too big (The even-valued Fibonacci terms cannot be greater than 4 million) ."

# ╔═╡ d1eb63b2-1b8f-11eb-085d-ebde562eb9db
# This will tell you which numbers are greater a certain number (4 million)in the terminal!

for k in 1 : g
	if even_valued_fibonacci_numbers[k] < 4000000
		println("The $(k) th term ($(even_valued_fibonacci_numbers[k])) is less than 4 million!")
		k + 1
	else
		println("The $(k) th term ($(even_valued_fibonacci_numbers[k])) is greater than 4 million!")
	end
end

# ╔═╡ cfaf4300-25da-11eb-04d8-0bfb97249bb2
md"# The results are in the terminal! "

# ╔═╡ eddc31c0-1b95-11eb-1282-892c325c95a6
# The terminal told me that the twelfth term is too big, but the eleventh one isn't.

upperbound = 11

# ╔═╡ 62313b10-1b96-11eb-08e9-d31d4ef0943f
l = 1 : upperbound

# ╔═╡ 6db45030-1b96-11eb-2c41-ed2e1c5441a6
m = φ .^ (3 * l)

# ╔═╡ 7f796a80-1b96-11eb-1582-7b4321ab93e7
n = m / sqrt(5)

# ╔═╡ 51c743f0-1b96-11eb-36c9-a9ce56cd8859
even_valued_fibonacci_terms_under_4_million = round.(n)

# ╔═╡ 1a2a1e20-205c-11eb-1462-c700d278e137
sum(even_valued_fibonacci_terms_under_4_million)

# ╔═╡ f295afd0-25d5-11eb-0e04-b386d7db338d
md"# The answer is 4,613,732! I got it right!"

# ╔═╡ 0d4225c0-25d6-11eb-1b80-cb6b323a7181
# I took a screenshot of Project Euler.

My_Image = load(download("https://i.imgur.com/V1pD1xR.png"))

# ╔═╡ db8154c0-25d5-11eb-0ddd-3b57dae3e708
md"# Extra Info"

# ╔═╡ af9236c0-25dc-11eb-0ffa-d724f99a5364
# The code is run in the REPL (Read, Eval, Print, Loop)

My_Image_2 = load(download("https://i.imgur.com/6Pb2xel.png"))

# ╔═╡ 08debbe0-25e7-11eb-1ec7-fb24a591077c
My_Image_2[532:1005, 545:1875]

# ╔═╡ 719b8920-25e1-11eb-2c6e-3f9ad3d236e4
md"# Epilogue"

# ╔═╡ 233dbcd0-25e1-11eb-26d3-f51fca7956d8
md"#### After doing this project, I developed a _completely autonomous algorithm_. I uploaded it on Google Classroom. Bye! "

# ╔═╡ Cell order:
# ╟─6ab4e9c0-114c-11eb-2d22-05776bd25df1
# ╟─1524cf62-114d-11eb-3dc3-4debd151734b
# ╟─89b61e60-114d-11eb-0375-4bcab2fef870
# ╠═22731370-25d9-11eb-049b-b546c56ec5ad
# ╠═6ef159f2-1157-11eb-04fa-a9542f145bf1
# ╟─f9488b00-1157-11eb-22ed-bdc7f113fbe9
# ╟─41c9f2c0-114d-11eb-0d12-11a7feced986
# ╟─168eb990-25db-11eb-28db-79e076dc7ea3
# ╠═f3cfe3e0-25dd-11eb-1fbc-9f9c70a63261
# ╟─2bc5c920-25db-11eb-000b-bf9df6093842
# ╠═e7abfa92-1147-11eb-2b44-159345577c2a
# ╠═3ee0f5de-1148-11eb-2039-93fa7b1be9af
# ╟─a3fef210-1b75-11eb-0311-5348902769b0
# ╠═96044e72-1b8f-11eb-0bd5-bf6b54c9464d
# ╠═a6b1de40-1b8f-11eb-0e31-f17fed3c148b
# ╠═b05dc25e-1b8f-11eb-1a64-dd58a16b7dec
# ╠═b7dc3f30-1b8f-11eb-316d-af7911153b99
# ╠═83d045b0-1b8f-11eb-02e0-317499abf430
# ╟─40007ef0-2057-11eb-13c5-f7df94b8a45c
# ╟─5981e260-2057-11eb-14d6-eba9727f24a2
# ╠═651d1360-2057-11eb-3fbb-e7f27cb8358b
# ╠═cb4b9070-2058-11eb-3581-8577499eafc7
# ╠═6b0a8890-205a-11eb-2f7a-1fba9e2c34c0
# ╠═81033a70-205a-11eb-284a-3fa820bc460e
# ╠═9b287f00-205a-11eb-1a54-dba8016126e1
# ╟─c88ce7d0-1b85-11eb-292d-1da5aa66d749
# ╟─09903500-1b88-11eb-0365-09cc57ba4c8c
# ╠═d1eb63b2-1b8f-11eb-085d-ebde562eb9db
# ╟─cfaf4300-25da-11eb-04d8-0bfb97249bb2
# ╠═08debbe0-25e7-11eb-1ec7-fb24a591077c
# ╠═eddc31c0-1b95-11eb-1282-892c325c95a6
# ╠═62313b10-1b96-11eb-08e9-d31d4ef0943f
# ╠═6db45030-1b96-11eb-2c41-ed2e1c5441a6
# ╠═7f796a80-1b96-11eb-1582-7b4321ab93e7
# ╠═51c743f0-1b96-11eb-36c9-a9ce56cd8859
# ╠═1a2a1e20-205c-11eb-1462-c700d278e137
# ╟─f295afd0-25d5-11eb-0e04-b386d7db338d
# ╠═0d4225c0-25d6-11eb-1b80-cb6b323a7181
# ╟─db8154c0-25d5-11eb-0ddd-3b57dae3e708
# ╠═d47521c0-25d5-11eb-09e7-1934afa8689d
# ╠═af9236c0-25dc-11eb-0ffa-d724f99a5364
# ╟─719b8920-25e1-11eb-2c6e-3f9ad3d236e4
# ╟─233dbcd0-25e1-11eb-26d3-f51fca7956d8


text:
using KitePodModels
using Test

KiteUtils.set_data_path("")
include("plot_step_response.jl")

@testset "KitePodModels.jl" begin
    kcu = KCU(se())
    @test get_depower(kcu) ≈ 0.236
    @test get_steering(kcu) ≈ 0.0
    set_depower_steering(kcu, 1.0, 1.0)
    KitePodModels.on_timer(kcu)
    @test get_depower(kcu) ≈ 0.23975
    @test get_steering(kcu) ≈ 0.01
    @test calc_alpha_depower(kcu, se().depower_offset/100.0) ≈ 0.0
    @test rad2deg(calc_alpha_depower(kcu, 0.37)) ≈ 20.499517119485063
end


text:
"""
abstract ParameterState{S<:ValueSupport, F<:VariateForm} <: VariableState{F}

Root of parameter state type hierarchy
"""
abstract type ParameterState{S<:ValueSupport, F<:VariateForm} <: VariableState{F} end

DiscreteParameterState{F<:VariateForm} = ParameterState{Discrete, F}
ContinuousParameterState{F<:VariateForm} = ParameterState{Continuous, F}

UnivariateParameterState{S<:ValueSupport} = ParameterState{S, Univariate}
MultivariateParameterState{S<:ValueSupport} = ParameterState{S, Multivariate}
MatrixvariateParameterState{S<:ValueSupport} = ParameterState{S, Matrixvariate}

ParameterStateVector{S<:ParameterState} = Vector{S}

value_support{S<:ValueSupport, F<:VariateForm}(::Type{ParameterState{S, F}}) = S
variate_form{S<:ValueSupport, F<:VariateForm}(::Type{ParameterState{S, F}}) = F

diagnostics(state::ParameterState) = Dict(zip(state.diagnostickeys, state.diagnosticvalues))

==(z::S, w::S) where {S<:ParameterState} = reduce(&, [getfield(z, n) == getfield(w, n) for n in fieldnames(S)])

isequal(z::S, w::S) where {S<:ParameterState} = reduce(&, [isequal(getfield(z, n), getfield(w, n)) for n in fieldnames(S)])


text:
using JSON
using AbstractTrees

mutable struct BinaryNode{T}
    data::T
    parent::BinaryNode{T}
    left::BinaryNode{T}
    right::BinaryNode{T}

    # Root constructor
    BinaryNode{T}(data) where T = new{T}(data)
    # Child node constructor
    BinaryNode{T}(data, parent::BinaryNode{T}) where T = new{T}(data, parent)
end

BinaryNode(data) = BinaryNode{typeof(data)}(data)

function leftchild(data, parent::BinaryNode)
    #!isdefined(parent, :left) || error("left child is already assigned")
    node = typeof(parent)(data, parent)
    parent.left = node
end
function rightchild(data, parent::BinaryNode)
    #!isdefined(parent, :right) || error("right child is already assigned")
    node = typeof(parent)(data, parent)
    parent.right = node
end

function AbstractTrees.children(node::BinaryNode)
    if isdefined(node, :left)
        if isdefined(node, :right)
            return (node.left, node.right)
        end
        return (node.left,)
    end
    isdefined(node, :right) && return (node.right,)
    return ()
end

AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)

Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()

function arrToTree(arr, parent = nothing)
    if length(arr) > 0
        if isnothing(parent)
            root = BinaryNode{Union{Int, Nothing}}(nothing)
        else
            root = parent
        end
        if typeof(arr[1]) == Int
            left = leftchild(arr[1], root)
        else
            left = leftchild(nothing, root)
            arrToTree(arr[1], left)
        end
        if typeof(arr[2]) == Int
            right = rightchild(arr[2], root)
        else
            right = rightchild(nothing, root)
            arrToTree(arr[2], right)
        end

        return root
    end
end

function readInput(filename)
    values = readlines(filename*".txt")
    numbers = []
    for line in values
        n = JSON.parse(line)
        numbers = push!(numbers, n)
    end
    return numbers
end

function explode(t)
    leaves = [l for l in AbstractTrees.Leaves(t)]
    #find deepest pair
    deepest, deepidx = deepestPair(t)
    #find left and right neighbor
    lidx = ridx = nothing
    if deepidx > 1
        lidx = deepidx - 1
    end
    if deepidx + 1 < length(leaves)
        ridx = deepidx + 2
    end
    lngb = rngb = nothing
    if isnothing(lidx) == false
        lngb = leaves[lidx].data
    end
    if isnothing(ridx) == false
        rngb = leaves[ridx].data
    end
    #compute new values for neighbors
    new_l = nothing
    if isnothing(lngb) == false
        new_l = deepest.left.data + lngb
    end
    new_r = nothing
    if isnothing(rngb) == false
        new_r = deepest.right.data + rngb
    end
    #change values of left and right neighbors

    if isnothing(new_l) == false
        lneigh = leaves[lidx]
        p = lneigh.parent
        if isleft(p, lneigh)
            lneigh = leftchild(new_l, p)
        else
            lneigh = rightchild(new_l,p)
        end
        if isdefined(p, :parent) == false
            t.left = p
        end
    end
    if isnothing(new_r) == false
        rneigh = leaves[ridx]
        p = rneigh.parent
        if isleft(p, rneigh)
            rneigh = leftchild(new_r, p)
        else
            rneigh = rightchild(new_r, p)
        end
        if isdefined(p, :parent) == false
            t.right = p
        end
    end
    #replace deepest pair root with 0
    if isleft(deepest.parent, deepest)
        deepest = leftchild(0, deepest.parent)
    else
        deepest = rightchild(0, deepest.parent)
    end

    return t

end

function countparents(t)
    count = 0
    while isdefined(t, :parent)
        count += 1
        t = t.parent
    end
    return count
end

function deepestPair(t)
    leaves = [l for l in AbstractTrees.Leaves(t)]
    max_depth = 0
    pair = t
    deepestindex = 0
    for (i,l) in enumerate(leaves)
        is_left = isleft(l.parent, l)
        if typeof(l.data) == Int
            if (is_left && typeof(l.parent.right.data) == Int) || (is_left == false && typeof(l.parent.left.data) == Int)
                d = countparents(l)
                if d > max_depth
                    pair = l.parent
                    deepestindex = i
                    max_depth = d
                end
            end
        end
    end

    return pair, deepestindex
end

function isleft(parent, child)
    if parent.left == child
        return true
    else
        return false
    end
end

function leftneighbor(t, pair)
    pairleft = pair.left.data
    leaves = [l.data for l in AbstractTrees.Leaves(t)]
    lidx = findfirst(x -> x == pairleft, leaves)
    if leaves[lidx + 1] != pair.right.data
        checked_leaves = copy(leaves)
        for i in 1:lidx
            checked_leaves[i] = -1
        end
        lidx = findfirst(x -> x == pairleft, checked_leaves)
    end
    if lidx == 1
        return nothing, 0
    else
        return leaves[lidx - 1], lidx - 1
    end
end

function rightneighbor(t, pair)
    pairright = pair.right.data
    leaves = [l.data for l in AbstractTrees.Leaves(t)]
    ridx = findfirst(x -> x == pairright, leaves)
    while ridx == 1 || leaves[ridx - 1] != pair.left.data
        checked_leaves = copy(leaves)
        for i in 1:ridx
            checked_leaves[i] = -1
        end
        ridx = findfirst(x -> x == pairright, checked_leaves)
    end
    if ridx == length(leaves)
        return nothing, 0
    else
        return leaves[ridx + 1], ridx + 1
    end
end

function getNeighbors(t, pair)
    pairleft = pair.left.data
    pairright = pair.right.data
    leaves = [l.data for l in AbstractTrees.Leaves(t)]
    lidx = ridx = nothing
    for i in 1:length(leaves)-1
        if leaves[i] == pairleft && leaves[i+1] == pairright
            if i > 1
                lidx = i-1
            end
            if i+1 < length(leaves)
                ridx = i+2
            end
            break
        end
    end
    return lidx, ridx
end

function treedepth(t)
    if typeof(t.data) == Int
        return 0
    end
    ld = treedepth(t.left)
    rd = treedepth(t.right)
    return 1 + max(ld, rd)

end

function needsEx(t)
    d = treedepth(t)
    if d > 4
        return true
    else
        return false
    end
end

function needsSplit(t)
    leaves = [l for l in AbstractTrees.Leaves(t)]
    for l in leaves
        if l.data >= 10
            return true
        end
    end
    return false
end

function split_sn(t)
    leaves = [l for l in AbstractTrees.Leaves(t)]
    for l in leaves
        if l.data >= 10
            p = l.parent
            new_l = floor(l.data/2)
            new_r = ceil(l.data/2)
            if isleft(p, l)
                sp = leftchild(nothing, p)
                p = p.left
                nl = leftchild(new_l, p)
                nr = rightchild(new_r,p)
            else
                sp = rightchild(nothing, p)
                p = p.right
                nl = leftchild(new_l, p)
                nr = rightchild(new_r,p)
            end
            break
        end
    end
    return t
end

function add(t1, t2)
    root = BinaryNode{Union{Int,Nothing}}(nothing)
    l = leftchild(nothing, root)
    r = rightchild(nothing, root)
    l.left = t1.left
    t1.left.parent = l
    l.right = t1.right
    t1.right.parent = l
    r.left = t2.left
    t2.left.parent = r
    r.right = t2.right
    t2.right.parent = r
    while needsEx(root) || needsSplit(root)
        while needsEx(root)
            root = explode(root)
        end
        if needsSplit(root)
            root = split_sn(root)
        end
    end
    return root
end

function getMagnitude(t)
    if typeof(t.left.data) == Int && typeof(t.right.data) == Int
        return 3*t.left.data + 2*t.right.data
    elseif typeof(t.left.data) == Int
        return 3*t.left.data + 2*getMagnitude(t.right)
    elseif typeof(t.right.data) == Int
        return 3* getMagnitude(t.left) + 2*t.right.data
    else
        return 3 * getMagnitude(t.left) + 2 * getMagnitude(t.right)
    end
end

function solve(filename)
    numbers = readInput(filename)
    trees = []
    for n in numbers
        tree = arrToTree(n)
        trees = push!(trees, tree)
    end
    t1 = deepcopy(trees[1])
    for (i, t) in enumerate(trees)
        if i != length(trees)
            t_next = deepcopy(trees[i+1])
            t1 = add(t1, t_next)
        end
    end
    println(getMagnitude(t1))

    maxMagn = 0
    max_i = 0
    max_j = 0
    for i in 1:length(trees)
        for j in 1:length(trees)
            t = deepcopy(trees)
            if i != j
                m = getMagnitude(add(t[i], t[j]))
                if m > maxMagn
                    maxMagn = m
                    max_i = i
                    max_j = j
                end
            end
        end
    end
    println("Max magnitude: "*string(maxMagn))
    println("Lines: "*string(max_i)*" and "*string(max_j))
end

#solve("example")
solve("input")


text:
#=
BFSParentsMultiSource:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Random

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSParentsMultiSource.jl")
dim = 10000
sourcdim = 3
matrix =  sprand(Bool, dim, dim, 0.1)
matrix = GBMatrix(matrix)
sources =  rand(0:0, sourcdim, dim)
sources[1, 1] = 1
sources[2, 2] = 2
sources[3, 3] = 3
sources = GBMatrix(sparse(sources))
@time ms_bfsp(matrix, sources, dim, sourcdim)


text:
#Mixture of t Factor Analyzers
function mtfa_msteps(tau, estep_estimates, Mquantities,lower_nu = 0.1,upper_nu = 200)
  #browser()
  #tau_ij  <- estep_estimates[,1]
    #xi_ij = estep_estimates[:,1]
    #zeta_ij = estep_estimates[:,2]
    #MAKE THIS USE GLOBAL VARIABLES
    Mquantities.mq.fnu = function f_nu_i(nu_i)
      #log(nu_i/2) + 1 - digamma(nu_i/2) + sum(tau_ij*(zeta_ij - xi_ij))/sum(tau_ij)
      return - (sum(tau)*(0.5*nu_i*log(0.5*nu_i) - log(gamma(0.5*nu_i))) + 0.5*nu_i*sum(tau .* (estep_estimates[:,2] .- estep_estimates[:,1])))
    end

    #function f_nu_i(nu_i)
      #log(nu_i/2) + 1 - digamma(nu_i/2) + sum(tau_ij*(zeta_ij - xi_ij))/sum(tau_ij)
    #  return log.(0.5*nu_i) .+ 1 .- digamma.(0.5*nu_i) .+ sum(tau .* (zeta_ij .- xi_ij))/sum(tau)
    #end

    #return [find_zero(f_nu_i, (lower_nu, upper_nu), Bisection())]

    Mquantities.mq.optim_min = optimize(Mquantities.mq.fnu, lower_nu, upper_nu)

    if !Optim.converged(Mquantities.mq.optim_min)
      print("Error estimating nu")
      return nothing
    end

    return [Optim.minimizer(Mquantities.mq.optim_min)]
end

#Mixture of Slash Factor Analyzers
function mslfa_msteps(tau, estep_estimates,Mquantities)
    return  [-sum(tau)/sum(tau.*estep_estimates[:,2])]
end

#Mixture of Contaminated Normal Factor Analyzers
function mcnfa_msteps(tau,estep_estimates,Mquantities)
    #tau_ij  = estep_estimates[:,1]
    #kappa_ij = estep_estimates[:,3]
    #delta_ij = estep_estimates[:,4]
    #p = estep_estimates[1,6] #This entry is just a column of p's
    #nu_i =
    #gamma_i =
  return [sum(tau .* estep_estimates[:,2])/sum(tau); estep_estimates[1,4]*sum(tau .* estep_estimates[:,2])/sum(tau .* estep_estimates[:,2] .* estep_estimates[:,3])]
end
#Mixture of Contaminated Normal Factor Analyzers

function mghfa_msteps(tau, estep_estimates, Mquantities, lower_lam = -100, upper_lam = 100, lower_omega = 1e-4, upper_omega = 500)
    #browser()
    #tau_ij  <- estep_estimates[,1]
    #tau_ij  <- estep_estimates[,1]
    #xi_ij <- estep_estimates[,1]
    #rho_ij <- estep_estimates[,2]
    ##zeta_ij <- estep_estimates[,3]
    #omega_k <- estep_estimates[1,4]
    #lambda_k <- estep_estimates[1,5]

  Mquantities.mq.zeta_bar = sum(estep_estimates[:,3].*tau)/sum(tau)
  Mquantities.mq.xi_bar = sum(estep_estimates[:,1].*tau)/sum(tau)
  Mquantities.mq.rho_bar = sum(estep_estimates[:,2].*tau)/sum(tau)

  Mquantities.mq.flam = function f_lam_i(lam_i)
      #log(nu_i/2) + 1 - digamma(nu_i/2) + sum(tau_ij*(zeta_ij - xi_ij))/sum(tau_ij)
      return - ( Mquantities.mq.zeta_bar*lam_i - log(besselk(lam_i,estep_estimates[1,4])))
  end

  Mquantities.mq.fomega = function f_omega_i(omega_i)
    #log(nu_i/2) + 1 - digamma(nu_i/2) + sum(tau_ij*(zeta_ij - xi_ij))/sum(tau_ij)
    return  log(besselk(estep_estimates[1,5],omega_i)) + 0.5*omega_i*( Mquantities.mq.xi_bar+ Mquantities.mq.rho_bar)
  end

  Mquantities.mq.min_lam = optimize(Mquantities.mq.flam, lower_lam, upper_lam)

  Mquantities.mq.min_omega = optimize(Mquantities.mq.fomega, lower_omega, upper_omega)

  if !Optim.converged(Mquantities.mq.min_lam) | !Optim.converged(Mquantities.mq.min_omega)
    print("Error estimating scaling density parameters")
    return nothing
  end

  return [Optim.minimizer(Mquantities.mq.min_omega);Optim.minimizer(Mquantities.mq.min_lam)]
end

function mbsfa_msteps(tau,estep_estimates,Mquantities)
#  #estep_order: tau, xi, rho
  #tau_ij  <- estep_estimates[,1]
  #xi_ij = estep_estimates[:,1]
  #rho_ij = estep_estimates[:,2]
  Mquantities.mq.xi_bar = sum(estep_estimates[:,1] .* tau)/sum(tau)
  Mquantities.mq.rho_bar = sum(estep_estimates[:,2] .* tau)/sum(tau)
  return [sqrt(Mquantities.mq.xi_bar + Mquantities.mq.rho_bar -2)]
end

function mlfa_msteps(tau,estep_estimates,Mquantities)
  #estep_order: tau, xi
    #au_ij  <- estep_estimates[,1]
    #xi_ij <- estep_estimates[:,1]
    Mquantities.mq.xi_bar = sum(estep_estimates[:,1] .* tau)/sum(tau)
    #if(((xi_bar-1)^2 + 8*xi_bar)<0){browser()}
    return [(-(Mquantities.mq.xi_bar-1) + sqrt((Mquantities.mq.xi_bar-1)^2 + 8*Mquantities.mq.xi_bar))/(2*Mquantities.mq.xi_bar)]
end

max_stars_repo_name: john-c-davey/FactorMixtures

max_issues_repo_name: john-c-davey/FactorMixtures

max_forks_repo_name: john-c-davey/FactorMixtures


text:
###Function to generate data for a spring-mass dynamic system
### t: timepoints to evaluate
### y0: initial conditions
### paras: values of the parameters (K,M)
function springmassdata(t::AbstractVector,y0::AbstractVector,paras::AbstractVector)
    a = sqrt(paras[1]/paras[2])
    y = [y0[1]*cos(a*t)+y0[2]/a*sin(a*t) -a*y0[1]*sin(a*t)+y0[2]*cos(a*t)]
end

###Function to generate noisy data for a spring-mass dynamic system
### t: timepoints to evaluate
### y0: initial conditions
### paras: values of the parameters (K,M)
### noisegenerator: an AbstractNoiseModel used to apply measurement noise to the dynamic system
springmassnoisydata(t::AbstractVector,y0::AbstractVector,paras::AbstractVector,noisegenerator::AbstractNoiseModel) = applynoise!(noisegenerator,springmassdata(t,y0,paras))

###ODE for the spring-mass dynamic system
###The second-order differentional equation a = -K/M * x
###is written as a system of 2 first-order differential equations
###Function arguments are as expected by the Sundials ODE solver package
### t: timepoints to evaluate the ODE
### y: the values of the variables (x,v=dx/dt)
### ydot: the derivate vaues (v=dx/dt,a=d2x/dt2)
### paras: the equation parameters (K,M)
function springmassode(t,y,ydot,paras)

    ydot[1] = y[2]
    ydot[2] = -paras[1]/paras[2]*y[1] #-K/M*X

    nothing
end

###Convenience function to create a spring-mass ODEModel
### t: timepoints to evaluate the ODE
### y0: initial conditions of the dynamic system
### measurementparas: the values for the parameters for which measurement data is generated (the "ground truth" of the experiment)
### variance: the variance of measurement noise
### paraminit: specification of the model parameters (K,M)
function springmassmodel(t::AbstractVector,y0::AbstractVector,measurementparas::AbstractVector,variance::AbstractVector,paraminit...)
    modelparameters = parameters([:K,:M],paraminit...)
    noisemodel = noise(:gaussian,variance)
    measurementdata = data(:function,t,springmassnoisydata,t,y0,measurementparas,noisemodel)
    model(:ode,modelparameters,measurementdata,noisemodel,springmassode,y0,2,[1,2];name = "Spring-Mass ODE")
end


text:
################################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
#   steelT3 model from
#   AMPL: A Modeling Language for Mathematical Programming, 2nd ed
#   by Robert Fourer, David Gay, and Brian W. Kernighan
#
#   Author : Louis Luangkesorn
#   Date: April 3, 2015
#
################################################################################

using JuMP

# Data

T = 4
prod = ["bands", "coils"]
area = Dict([("bands", ("east", "north")),
             ("coils", ("east", "west", "export"))])

avail = [40, 40, 32,40]

rate = Dict("bands" => 200,
             "coils" =>140)

inv0 = Dict("bands" => 10,
             "coils" => 0)

prodcost = Dict("bands" => 10,
                 "coils" => 11)
invcost = Dict("bands" => 2.5,
                "coils" => 3)

revenue = Dict("bands" => Dict("east" => [25.0 26.0 27.0 27.0],
                                 "north" => [26.5 27.5 28.0 28.5]),
                "coils" => Dict("east" =>[30    35    37    39],
                                "west" => [29    32    33    35],
                                "export" => [25    25    25    28])
              )

market = Dict("bands" => Dict("east" => [2000  2000  1500  2000],
                                "north" => [4000  4000  2500  4500]),
               "coils" => Dict("east" => [1000   800  1000  1100],
                                "west" => [2000  1200  2000  2300],
                                "export" => [1000   500   500   800])
              )

# Model

Prod = Model()
print(area["bands"])
print(market["coils"]["west"][1])

# Decision Variables

@variable(Prod, Make[p=prod, t=1:T] >= 0); # tons produced
@variable(Prod, Inv[p=prod, t=0:T] >= 0); # tons inventoried
@variable(Prod, market[p][a][t] >= Sell[p=prod, a=area[p],t=1:T] >= 0); # tons sold

@constraint(Prod, triconstr[p=prod, a=area[p], t=1:T],
               Sell[p, a, t] - market[p][a][t] <= 0)

@constraint(Prod, xyconst[t=1:T],
               sum{(1/rate[p]) * Make[p,t], p=prod} <= avail[t])

# Total of hours used by all products
# may not exceed hours available, in each week

@constraint(Prod, xyconst[p=prod],
               Inv[p,0] == inv0[p])
# Initial inventory must equal given value

@constraint(Prod, xyconst[p=prod, t=1:T],
               Make[p,t] + Inv[p, t-1] == sum{Sell[p,a,t], a=area[p]} + Inv[p,t])

# Tons produced and taken from inventory
# must equal tons sold and put into inventory

@objective(Prod, Max,
              sum{
                sum{revenue[p][a][t] * Sell[p, a, t] -
                      prodcost[p] * Make[p,t] -
                      invcost[p]*Inv[p,t],
                    a = area[p]},
                p=prod, t=1:T})
#maximize Total_Profit:
# Total revenue less costs for all products in all weeks

function PrintSolution(status, area, Make, Inventory, Sell, product, Time)
    println("RESULTS:")
    if status == :Optimal
      for p = product
        println("Make $(p)")
        for t = 1:T
          print("$(getvalue(Make[p,t]))\t")
        end
        println()
        println("Inventory $(p)")
        for t=1:T
          print("$(getvalue(Inventory[p,t]))\t")
        end
        println()
        for a = area[p]
          println("Sell $(p) $(a)")
          for t=1:T
            print("$(getvalue(Sell[p,a,t])) \t")
          end
        println()
        end
      end
    else
        println("  No solution")
    end
    println("")
end

status = solve(Prod)
PrintSolution(status, area, Make, Inv, Sell, prod, T)


text:
abstract type AbstractDatasetRegistry end

"""
    DatasetRegistry()

A store for information on datasets and dataset recipes for loading those datasets.
"""
Base.@kwdef struct DatasetRegistry <: AbstractDatasetRegistry
    datasets::Dict{String,Any} = Dict{String,Any}()
    recipes::Dict{String,Vector{DatasetRecipe}} = Dict{String,Vector{DatasetRecipe}}()
end

## Queries

"""
    listdatasources([registry])

List the dataset sources registered in `registry` (defaults to
`FastAI.defaultdataregistry()`).
"""
listdatasources(reg::DatasetRegistry) = collect(keys(reg.datasets))

"""
    datasetpath([registry], name)

(Down)load registered dataset source named `name` from the dataset registry.
Use [`listdatasources`](#) for a list of all dataset sources.
"""
datasetpath(reg::DatasetRegistry, name::String) = reg.datasets[name]()

"""
    finddatasets([registry]; name=nothing, blocks=Any)

Find preconfigured dataset recipes for datasets that match block
types `blocks` in all data sources (if `name == nothing`) or dataset source
with `name`. `blocks` can be given as a type or a nested tuple of block types.

Return a vector of `Pair`s `datasetname => recipe`

#### Examples

Loading a result

```julia
datasetname, recipe = finddatasets(blocks=(Image, Label))[1]
data, blocks = loadrecipe(recipe, datasetpath(datasetname))
```

Example searches

```
# Single-label image classification
finddatasets(blocks=(Image, Label))

# Single-label classification from any data
finddatasets(blocks=(Any, Label))

# Datasets with images as input data
finddatasets(blocks=(Image, Any))

# All ways to load `pascal2007`
finddatasets(name="pascal2007")
```
"""
function finddatasets(reg::DatasetRegistry; name=nothing, blocks=Any)::Vector{Pair{String,DatasetRecipe}}
    results = collect(Iterators.flatten(((k => v) for v in vs) for (k, vs) in reg.recipes))
    results = filter(((d, recipe),) -> typify(recipeblocks(recipe)) <: typify(blocks), results)
    if !isnothing(name)
        results = filter(((d, r),) -> d == name, results)
    end
    return results
end

"""
    loaddataset(name[, blocks = Any]) -> (data, blocks)

Load dataset `name` with a recipe that matches block types
`blocks`. The first matching recipe is selected and loaded.

## Examples

Load a data container suitable for single-label image classification:
```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
```

Load dataset with any recipe:
```julia
data, blocks = loaddataset(name)
```
"""
function loaddataset(reg::DatasetRegistry, name::String, blocks=Any)
    results = finddatasets(reg; name=name, blocks=blocks)
    isempty(results) && error("Could not find a recipe for dataset \"$name\" that matches block types `$blocks`")
    name_, recipe = results[1]
    @assert name == name_
    return loadrecipe(recipe, datasetpath(reg, name))
end

## Registration

"""
    registerdataset!([registry,] name, loadfn)

Register a dataset in `registry::DatasetRegistry` with a `name` and a
function `loadfn()` that downloads it and returns a path.
"""
function registerdataset!(reg::DatasetRegistry, dataset::String, loadfn)
    dataset ∈ keys(reg.datasets) && error(
        "Dataset name $(dataset) is already registered")
    reg.datasets[dataset] = loadfn
    return reg
end

"""
    registerrecipe!([registry,] name, recipe)

Register a recipe in `registry::DatasetRegistry` for dataset `name`.
"""
function registerrecipe!(reg::DatasetRegistry, dataset::String, recipe::DatasetRecipe)
    dataset ∈ keys(reg.datasets) || error("Dataset $dataset not found.")
    recipes = get!(reg.recipes, dataset, DatasetRecipe[])
    push!(recipes, recipe)
    return reg
end

function blocksmatch(recipe::DatasetRecipe, TBlocks)
    @show recipeblocks(recipe), TBlocks
	return recipeblocks(recipe) <: TBlocks
end


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

mutable struct OnErrorDeploymentExtended <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    deploymentName::Any # spec type: Union{ Nothing, String } # spec name: deploymentName

    function OnErrorDeploymentExtended(;provisioningState=nothing, type=nothing, deploymentName=nothing)
        o = new()
        validate_property(OnErrorDeploymentExtended, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(OnErrorDeploymentExtended, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(OnErrorDeploymentExtended, Symbol("deploymentName"), deploymentName)
        setfield!(o, Symbol("deploymentName"), deploymentName)
        o
    end
end # type OnErrorDeploymentExtended

const _property_map_OnErrorDeploymentExtended = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("type")=>Symbol("type"), Symbol("deploymentName")=>Symbol("deploymentName"))
const _property_types_OnErrorDeploymentExtended = Dict{Symbol,String}(Symbol("provisioningState")=>"String", Symbol("type")=>"String", Symbol("deploymentName")=>"String")
Base.propertynames(::Type{ OnErrorDeploymentExtended }) = collect(keys(_property_map_OnErrorDeploymentExtended))
Swagger.property_type(::Type{ OnErrorDeploymentExtended }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OnErrorDeploymentExtended[name]))}
Swagger.field_name(::Type{ OnErrorDeploymentExtended }, property_name::Symbol) =  _property_map_OnErrorDeploymentExtended[property_name]

const _allowed_OnErrorDeploymentExtended_type = ["LastSuccessful", "SpecificDeployment"]

function check_required(o::OnErrorDeploymentExtended)
    true
end

function validate_property(::Type{ OnErrorDeploymentExtended }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "OnErrorDeploymentExtended", :enum, val, _allowed_OnErrorDeploymentExtended_type)
    end
end


text:
#
# Example of a set of reference values that StateCheck will use to test current vlaues against.
# These are generated by the function scprintref() and parsed by the scdocheck()
#

#! format: off
varr = [
 [ "My gradients", "ν∇u[1]",  1.34348869323730468750e-04,  9.84732866287231445313e-01,  5.23545503616333007813e-01,  3.08209930764271777814e-01 ],
 [ "My gradients", "ν∇u[2]",  1.16317868232727050781e-01,  9.92088317871093750000e-01,  4.83800649642944335938e-01,  2.83350456014221541157e-01 ],
 [ "My gradients", "ν∇u[3]",  1.05845928192138671875e-03,  9.51775908470153808594e-01,  4.65474426746368408203e-01,  2.73615551085745090099e-01 ],
 [ "My gradients", "ν∇u[4]",  5.97668886184692382813e-02,  9.68048095703125000000e-01,  5.42618036270141601563e-01,  2.81570862027933854765e-01 ],
 [ "My gradients", "ν∇u[5]",  8.31030607223510742188e-02,  9.35931921005249023438e-01,  5.05405902862548828125e-01,  2.46073509972619536290e-01 ],
 [ "My gradients", "ν∇u[6]",  3.09681892395019531250e-02,  9.98341441154479980469e-01,  4.54375565052032470703e-01,  3.09461067853178561915e-01 ],
 [ "My gradients", "κ∇θ[1]",  8.47448110580444335938e-02,  9.94180679321289062500e-01,  5.27157366275787353516e-01,  2.92455951648181833313e-01 ],
 [ "My gradients", "κ∇θ[2]",  1.20514631271362304688e-02,  9.93527650833129882813e-01,  4.71063584089279174805e-01,  2.96449027197666359346e-01 ],
 [ "My gradients", "κ∇θ[3]",  8.14980268478393554688e-02,  9.55443382263183593750e-01,  5.05038917064666748047e-01,  2.77201022741208891187e-01 ],
 [    "My fields",   "u[1]",  3.53445491472876849315e-02,  9.73118774570230771204e-01,  4.53566376673364857197e-01,  3.28622448223506280485e-01 ],
 [    "My fields",   "u[2]",  4.23016659320296639635e-02,  9.67799553619200114696e-01,  5.40307230728526155517e-01,  2.94603153327737565803e-01 ],
 [    "My fields",   "θ[1]",  6.23675581701588210848e-02,  9.76550123041147521974e-01,  5.16046312418334207628e-01,  2.81983244164903890105e-01 ],
]
parr = [
 [ "My gradients", "ν∇u[1]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[2]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[3]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[4]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[5]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[6]",    16,     7,    16,     0 ],
 [ "My gradients", "κ∇θ[1]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[2]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[3]",    16,    16,    16,     0 ],
 [    "My fields",   "u[1]",    16,    16,    16,     0 ],
 [    "My fields",   "u[2]",    16,    16,    16,     0 ],
 [    "My fields",   "θ[1]",    16,    16,    16,     0 ],
]
#! format: on


text:
function constant_features_agent(a::Integer, mapf::MAPF)
    t0 = mapf.starting_times[a]
    s = mapf.sources[a]
    d = mapf.destinations[a]
    dist = mapf.distances_to_destinations[d][s]
    return Float64[a, t0, dist]
end

function path_features_agent(a::Integer, solution::Solution, mapf::MAPF)
    path = solution[a]
    g = mapf.graph
    duration = length(path)
    weight = path_weight(path, mapf)
    mean_outdegree = mean(outdegree(g, v) for (t, v) in path)
    mean_indegree = mean(indegree(g, v) for (t, v) in path)
    return Float64[duration, weight, mean_indegree, mean_outdegree]
end

function conflict_features_agent(a::Integer, solution::Solution, mapf::MAPF)
    A = nb_agents(mapf)
    conflicts0 = count_conflicts(a, solution, mapf; tol=0) / 1
    conflicts1 = count_conflicts(a, solution, mapf; tol=1) / 3
    conflicts2 = count_conflicts(a, solution, mapf; tol=2) / 5
    conflicts3 = count_conflicts(a, solution, mapf; tol=3) / 7
    crossings = count_conflicts(a, solution, mapf; tol=Inf)
    return (1 / A) .* Float64[conflicts0, conflicts1, conflicts2, conflicts3, crossings]
end

function all_features_agent(a::Integer, solution::Solution, mapf::MAPF)
    return vcat(
        constant_features_agent(a, mapf),
        path_features_agent(a, solution, mapf),
        conflict_features_agent(a, solution, mapf),
    )
end

function agents_embedding(mapf::MAPF)
    solution = independent_astar(mapf)
    x = reduce(hcat, all_features_agent(a, solution, mapf) for a in 1:nb_agents(mapf))
    s = std(x; dims=2)
    s[isapprox.(s, 0.0)] .= 1  # columns with zero variance
    x = (x .- mean(x; dims=2)) ./ s
    return x
end


text:
using Pkg
using ArrayAllocators
using ArrayAllocators.ByteCalculators
using OffsetArrays
using Test

# Load in subpackages
Pkg.develop(PackageSpec(path=joinpath(dirname(@__DIR__), "./NumaAllocators")))
Pkg.develop(PackageSpec(path=joinpath(dirname(@__DIR__), "./SafeByteCalculators")))

using NumaAllocators
using SafeByteCalculators

@testset "ArrayAllocators.jl" begin
    A = zeros(UInt8, 2048, 2048);
    B = Array{UInt8}(calloc, 2048, 2048);
    M = Array{UInt8}(malloc, 1024, 4096)
    Z = ArrayAllocators.zeros(UInt8, 2048, 2048)
    Z2 = ArrayAllocators.zeros(UInt8, (2048, 2048))
    @test A == B
    @test A == Z
    @test size(Z) == (2048, 2048)
    @test size(Z2) == (2048, 2048)
    @test size(M) == (1024, 4096)

    @test_throws OverflowError Array{UInt8}(calloc, 20480000, typemax(Int64))
    @test_throws OverflowError Array{UInt16}(calloc, 2, typemax(Int64)÷2)

    @test_throws OverflowError Array{UInt16}(MallocAllocator(), 2, typemax(Int64)÷2)
    @test_throws OverflowError Array{UInt16}(MallocAllocator{CheckedMulByteCalculator}(), 2, typemax(Int64)÷2)
    @test_throws OverflowError Array{UInt16}(MallocAllocator{WideningByteCalculator}(), 2, typemax(Int64)÷2)
    @test_throws OverflowError Array{UInt16}(MallocAllocator{SafeByteCalculator}(), 2, typemax(Int64)÷2)

    @static if Sys.iswindows()
        WV = Array{UInt8}(ArrayAllocators.Windows.virtual, 64, 1024);
        @test size(WV) == (64, 1024)
        @test WV == zeros(UInt8, 64, 1024)
    end
    @static if Sys.islinux() || Sys.iswindows()
        C = Array{UInt8}(NumaAllocator(0), 2048, 2048);
        @test A == C
        @test current_numa_node() isa Int
        @test highest_numa_node() isa Int
    end

    @static if Sys.isunix() || Sys.iswindows()
        D = Array{UInt8}(MemAlign(), 1024, 4096)
        @test size(D) == (1024, 4096)
        @test reinterpret(Int, pointer(D)) % ArrayAllocators.alignment(MemAlign()) == 0
        E = Array{UInt8}(MemAlign(2^16), 1024, 2048)
        @test size(E) == (1024, 2048)
        @test reinterpret(Int, pointer(E)) % 2^16 == 0
    end
end

@testset "Composition with OffsetArrays.jl" begin
    OA = OffsetArray{UInt8}(calloc, -1024:1023, -5:5)
    @test all(OA .== 0)
    @test size(OA) == (2048,11)
    @test OA[-1024,-3] == 0
    OA = OffsetArray{UInt8}(calloc, -5:3, Base.OneTo(9))
    @test all(OA .== 0)
    @test size(OA) == (9,9)
    @test OA[-1,5] == 0
    @static if Sys.islinux() || Sys.iswindows()
        numaOA = OffsetArray{UInt8}(numa(0), -5:5, 1:9)
        @test numaOA[-5, 1] == numaOA[1]
    end
end


text:
## ---

"""
```julia
vsort(A; dims, multithreaded=:auto)
```
Return a sorted copy of the array `A`, optionally along dimensions specified by `dims`.

If sorting over multiple `dims`, these dimensions must be contiguous (i.e.
`dims=(2,3)` but not `dims=(1,3)`). Note also that specifying `dims` other than
`:` creates `view`s, with some nonzero performance cost.

See also `vsort!` for a more efficient in-place version.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1, 3, 2, 4];

julia> vsort(A)
4-element Vector{Int64}:
 1
 2
 3
 4
"""
vsort(A; dims=:, multithreaded=:auto) = vsort!(copy(A), dims=dims, multithreaded=multithreaded)

"""
```julia
vsort!(A; dims, multithreaded=:auto)
```
Sort the array `A`, optionally along dimensions specified by `dims`.

If sorting over multiple `dims`, these dimensions must be contiguous (i.e.
`dims=(2,3)` but not `dims=(1,3)`). Note also that specifying `dims` other than
`:` creates `view`s, with some nonezero performance cost.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1, 3, 2, 4];

julia> vsort!(A)
4-element Vector{Int64}:
 1
 2
 3
 4
```
"""
function vsort!(A; dims=:, multithreaded=:auto)
    if (multithreaded===:auto && length(A) > 16383) || multithreaded===true
        _vtsort!(A, dims)
    else
        _vsort!(A, dims)
    end
end
export vsort!

# Sort linearly (reducing along all dimensions)
function _vsort!(A::AbstractArray, ::Colon)
    # IF there are NaNs, move them all to the end of the array
    A, iₗ, iᵤ = sortnans!(A)
    # Sort the non-NaN elements
    quicksort!(A, iₗ, iᵤ)
end

## --- as above, but multithreaded

# Sort linearly (reducing along all dimensions)
function _vtsort!(A::AbstractArray, ::Colon)
    # IF there are NaNs, move them all to the end of the array
    A, iₗ, iᵤ = sortnans!(A)
    # Sort the non-NaN elements
    quicksortt!(A, iₗ, iᵤ)
end

# Fall back to singlethreaded for multidimensional cases
_vtsort!(A::AbstractArray, dims) = _vsort!(A::AbstractArray, dims)

## -- multidimensional cases (singlethreaded only)

# Reducing / sorting over noncontiguous dims may be a problem
function iscontiguous(dims)
    for i = 2:length(dims)
        if !(dims[i] == dims[i-1] + 1 || dims[i] == dims[i-1] - 1)
            return false
        end
    end
    return true
end

# Sort one dim
_vsort!(A, dims::Int) = _vsort!(A, (dims,))

# Sort some dims
function _vsort!(A::AbstractArray{T,N}, dims::Tuple) where {T,N}
    iscontiguous(dims) || error("Only continuous `dims` are supported")
    __vsort!(A, dims)
end

# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_sort_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the sort over.
  Aᵥ = Expr(:call, :view, :A)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(Aᵥ[])
  inds = Vector{Symbol}(undef, N)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Aᵥ.args, :)
    else
      push!(nonreduct_inds, n)
      push!(Aᵥ.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices(A,$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices(A,$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(Aᵥ = $Aᵥ))
    push!(rblock.args, :(_vsort!(Aᵥ, :)))
    # Put it all together
    return quote
      @inbounds $loops
      return A
    end
  else
    return quote
      _vsort!(A, :)
      return A
    end
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_sort_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return __vsort!(A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:N, drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return __vsort!(A, $tc))))
      return q
    end
  end
  return staticdim_sort_quote(static_dims, N)
end

# Efficient @generated in-place sort
@generated function __vsort!(A::AbstractArray{T,N}, dims::D) where {T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_sort_quote(N, M, D)
end
@generated function __vsort!(A::AbstractArray{T,N}, dims::Tuple{}) where {T,N}
  :(return A)
end


text:
ENV["GKSwstype"]="100"

using DifferentiableNAS
using Flux
using Flux: logitcrossentropy, onecold, onehotbatch
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Plots
using PyPlot
gr()
include("CIFAR10.jl")
include("training_utils.jl")
include("visualize.jl")

trials = Dict(  #"test/models/osirim/darts_2021-01-13T14:17:49.647" => "DARTS_none",
                #"test/models/osirim/darts_2021-01-13T14:24:11.711"=> "DARTS_sans_none",
                #"test/models/osirim/admm_2021-01-13T17:43:21.59"=> "ADMM_sans_none",
                #"test/models/osirim/admm_2021-01-14T16:55:27.481"=> "ADMM_none",
                #"test/models/osirim/admm_2021-01-14T17:16:51.535"=> "ADMM_128_RTX",
                #"test/models/osirim/admm_2021-01-15T12:06:38.741"=> "ADMM_128_RTX_sans_none_disc",
                #"test/models/osirim/admm_6641599" => "ADMM_128_sliding_freq_sliding_disc_1",
                #"test/models/osirim/admm_6641601" => "ADMM_128_sliding_freq",
                #"test/models/osirim/admm_6641606" => "ADMM_sliding_freq_sliding_disc",
                #"test/models/osirim/darts_6641607" => "BNDARTS",
                #"test/models/osirim/bnadmm_6641778" => "BNADMM_sliding_freq_sliding_disc",
                #"test/models/osirim/bnadmm_6641847" => "BNADMM_none_fixed",
                #"test/models/osirim/bnadmm_6641882" => "BNADMM_progressive_1e-2",
                #"test/models/osirim/sigadmm_6642019" => "SIGADMM",
                #"test/models/osirim/bndarts_6642020" => "BNDARTS",
                #"test/models/osirim/bnadmm_6642126" => "BNADMM_prog_1e-1",
                #"test/models/osirim/bnadmm_6642153" => "BNADMM_prog_relutanh",
                #"test/models/osirim/bnadmm_6642226" => "BNADMM_rho_disc_1e-3",
                #"test/models/osirim/bnadmm_6642233" => "BNADMM_rho_disc_5e-3",
                #"test/models/osirim/bnadmm_6642280" => "BNADMM_eval",
                #"test/models/osirim/bnadmm_6642317" => "BNADMM_rho_disc_1e-3_2",
                #"test/models/osirim/bnadmm_6642319" => "BNADMM_rho_disc_alr9e-4",
                #"test/models/osirim/bnadmm_6642365" => "BNADMM_fix_data",
                #"test/models/osirim/bnadmm_6642368" => "BNADMM_eval",
                #"test/models/osirim/darts_6642375" => "DARTS",
                #"test/models/osirim/bnadmm_6642376" => "BNADMM_rho_1e-3",
                #"test/models/osirim/bnadmm_6642377" => "BNADMM_rho_5e-3",
                #"test/models/osirim/bnadmm_6642384" => "BNADMM_rho_1e-2",
                #"test/models/osirim/bnadmm_6642392" => "BNADMM_rho_1e-3",
                #"test/models/osirim/bnadmm_6642393" => "BNADMM_eval",
                #"test/models/osirim/darts_6642395" => "DARTS_none",
                #"test/models/osirim/bnadmm_6642397" => "BNADMM_alr_3e-3",
                #"test/models/osirim/bnadmm_6642406" => "BNADMM_eval",
                #"test/models/osirim/bnadmm_6642409" => "DARTS_eval",
                #"test/models/osirim/bnadmm_6642571" => "BNADMM_singleupdate",
                #"test/models/osirim/darts_6642654" => "DARTS_cosanneal",
                #"test/models/osirim/darts_6642774" => "DARTS_test30epochs_frac",
                #"test/models/osirim/darts_6642805" => "DARTS_test30epochs_full",
                #"test/models/osirim/darts_6642941" => "DARTS_test30check_small",
                #"test/models/osirim/darts_6642942" => "DARTS_test30check_big",
                #"test/models/osirim/bnadmm_6644989" => "BNADMM_speedup",
                #"test/models/osirim/admm_6645078" => "BNADMM_scheduled",
                #"test/models/osirim/darts_6645238" => "DARTS_speedup",
                #"test/models/osirim/darts_6645482" => "DARTS_32_dep",
                #"test/models/osirim/darts_6645483" => "DARTS_32_oldv",
                #"test/models/osirim/darts_6645488" => "DARTS_64_dep",
                #"test/models/osirim/darts_6645542" => "DARTS_64_dep_cuda11",
                "test/models/osirim/darts_6645746" => "DARTS64",
                "test/models/osirim/admm_6645848" => "BNADMM64scheduled",
                "test/models/osirim/bnadmm_6645849" => "BNADMM64",
                "test/models/osirim/evaldarts/eval_6645902" => "eval",
                #"test/models/olympe/bnadmm_528425" => "BNADMM_o1",
                #"test/models/olympe/bnadmm_528440" => "BNADMM_o2",
                #"test/models/olympe/bnadmm_528445" => "BNADMM_o3",
                #"test/models/olympe/bnadmm_528487" => "BNADMM_o4",
                #"test/models/olympe/bnadmm_529404" => "BNADMM_o5",
                #"test/models/olympe/bnadmm_529405" => "BNADMM_o6",
                #"test/models/olympe/bnadmm_529406" => "BNADMM_o7",
                #"test/models/olympe/bnadmm_530364" => "BNADMM_o8",
                #"test/models/olympe/bnadmm_530427" => "BNADMM_o9",
                #"test/models/olympe/bnadmm_532776" => "BNADMM_agpu",
                #"test/models/olympe/bnadmm_532778" => "BNADMM_acpu",
                "test/models/olympe/evaldarts/eval_540757" => "eval",
                "test/models/olympe/bnadmm_540761" => "BNADMM",
                "test/models/olympe/admm_540763" => "BNADMMscheduled",
                "test/models/olympe/darts_540764" => "DARTS",
                )

for folder_name in keys(trials)

    if ispath(joinpath(folder_name, "histbatch.bson"))
        @show file_name = string(folder_name, "/histbatch.bson")
        BSON.@load file_name histbatch

        if ispath(joinpath(folder_name, "args.bson"))
            BSON.@load joinpath(folder_name, "args.bson") args
            batchsize = args["batchsize"]
            trainval_fraction = args["trainval_fraction"]
        else
            BSON.@load joinpath(folder_name, "model.bson") argparams
            batchsize = argparams.batchsize
            trainval_fraction = 1.0
        end
        prim = [
            "none",
            "max_pool_3x3",
            "avg_pool_3x3",
            "skip_connect",
            "sep_conv_3x3",
            "sep_conv_5x5",
            #"sep_conv_7x7",
            "dil_conv_3x3",
            "dil_conv_5x5",
            #"conv_7x1_1x7"
        ]

        nodelabels = ["c_{k-2}", "c_{k-1}", "0", "1", "2", "3", "c_{k}"]

        connects = vcat([[(j,i) for j = 1:i-1] for i = 3:6]...)
        if typeof(histbatch)<:histories
            normal_ = histbatch.normal_αs
            reduce_ = histbatch.reduce_αs

            if length(normal_[1][1]) == 7
                prim = prim[2:8]
            end
            n_y_min = minimum([softmax(a[i])[j] for a in histbatch.normal_αs for i in 1:14 for j in 1:length(prim)])
            n_y_max = maximum([softmax(a[i])[j] for a in histbatch.normal_αs for i in 1:14 for j in 1:length(prim)])
            for i = 1:14
                p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(n_y_min,n_y_max), legend=false)
                for j = 1:length(prim)
                    plot!([softmax(a[i])[j] for a in histbatch.normal_αs], label=prim[j])
                end
            end
            plot(p..., layout = (2,7), size = (2200,600));
            savefig(string(folder_name, "/normal_plots.png"))

            r_y_min = minimum([softmax(a[i])[j] for a in histbatch.reduce_αs for i in 1:14 for j in 1:length(prim)])
            r_y_max = maximum([softmax(a[i])[j] for a in histbatch.reduce_αs for i in 1:14 for j in 1:length(prim)])
            p = Vector(undef, 14)
            for i = 1:14
                p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(r_y_min,r_y_max), legend=false)
                for j = 1:length(prim)
                    plot!([softmax(a[i])[j] for a in histbatch.reduce_αs], label=prim[j])
                end
            end
            plot(p..., layout = (2,7), size = (2200,600));
            savefig(string(folder_name, "/reduce_plots.png"))
        else
            normal_ = histbatch.normal_αs_sm
            reduce_ = histbatch.reduce_αs_sm
            batches = 1:length(normal_)
            epochs = batches ./ (1+(25000 ÷ batchsize)*trainval_fraction)
            if length(normal_[1][1]) == 7
                if occursin("DC",trials[folder_name])
                    prim = prim[1:7]
                else
                    prim = prim[2:8]
                end
            else
                fig, axs = plt.subplots(4, 5, sharex="col", sharey="row")
                for (i,con) in enumerate(connects)
                    if con[2] == 6
                        axs[con[2]-2,con[1]].set_xlabel(latexstring(nodelabels[con[1]]))
                    end
                    if con[1] == 1
                        axs[con[2]-2,con[1]].set_ylabel(string(nodelabels[con[2]]))
                    end
                    for j = length(prim):-1:2
                        axs[con[2]-2,con[1]].plot(epochs, [a[i][j] for a in histbatch.normal_αs_sm], label = prim[j])
                    end
                end
                for i in 1:3
                    for j in 2+i:5
                        axs[i,j].axis("off")
                    end
                end
                handles, labels = axs[1,1].get_legend_handles_labels()
                fig.legend(handles, labels, loc="upper right")
                PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-normal_plots_sans_none.png"))

                fig, axs = plt.subplots(4, 5, sharex="col", sharey="row")
                for (i,con) in enumerate(connects)
                    if con[2] == 6
                        axs[con[2]-2,con[1]].set_xlabel(string(nodelabels[con[1]]))
                    end
                    if con[1] == 1
                        axs[con[2]-2,con[1]].set_ylabel(string(nodelabels[con[2]]))
                    end
                    for j = length(prim):-1:2
                        axs[con[2]-2,con[1]].plot(epochs, [a[i][j] for a in histbatch.reduce_αs_sm], label = prim[j])
                    end
                    if con[1] in [1,2]
                        axs[con[2]-2,con[1]].plot(epochs, [a[i][findall(x->x=="skip_connect", prim)[1]] for a in histbatch.reduce_αs_sm], label = "fact_red", color = "tab:olive")
                    end
                end
                for i in 1:3
                    for j in 2+i:5
                        axs[i,j].axis("off")
                    end
                end
                handles, labels = axs[1,1].get_legend_handles_labels()
                fig.legend(handles, labels, loc="upper right")
                PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-reduce_plots_sans_none.png"))
            end

            fig, axs = plt.subplots(4, 5, sharex="col", sharey="row")
            for (i,con) in enumerate(connects)
                if con[2] == 6
                    axs[con[2]-2,con[1]].set_xlabel(latexstring(nodelabels[con[1]]))
                end
                if con[1] == 1
                    axs[con[2]-2,con[1]].set_ylabel(string(nodelabels[con[2]]))
                end
                for j = length(prim):-1:1
                    axs[con[2]-2,con[1]].plot(epochs, [a[i][j] for a in histbatch.normal_αs_sm], label = prim[j])
                end
            end
            for i in 1:3
                for j in 2+i:5
                    axs[i,j].axis("off")
                end
            end
            handles, labels = axs[1,1].get_legend_handles_labels()
            fig.legend(handles, labels, loc="upper right")
            PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-normal_plots.png"))

            fig, axs = plt.subplots(4, 5, sharex="col", sharey="row")
            for (i,con) in enumerate(connects)
                if con[2] == 6
                    axs[con[2]-2,con[1]].set_xlabel(latexstring(nodelabels[con[1]]))
                end
                if con[1] == 1
                    axs[con[2]-2,con[1]].set_ylabel(string(nodelabels[con[2]]))
                end
                for j = length(prim):-1:1
                    axs[con[2]-2,con[1]].plot(epochs, [a[i][j] for a in histbatch.reduce_αs_sm], label = prim[j])
                end
                if con[1] in [1,2]
                    axs[con[2]-2,con[1]].plot(epochs, [a[i][findall(x->x=="skip_connect", prim)[1]] for a in histbatch.reduce_αs_sm], label = "fact_red", color = "tab:olive")
                end
            end
            for i in 1:3
                for j in 2+i:5
                    axs[i,j].axis("off")
                end
            end

            handles, labels = axs[1,1].get_legend_handles_labels()
            fig.legend(handles, labels, loc="upper right")
            PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-reduce_plots.png"))
            if typeof(histbatch)<:historiessml
                fig, ax = plt.subplots(1)
                ax.set_title("Losses")
                ax.plot(epochs, histbatch.train_losses, label = "training loss")
                ax.plot(epochs, histbatch.val_losses, label = "validation loss")
                ax.legend()
                PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-loss.png"))
                fig, axs = plt.subplots(2)
                axs[1].plot(epochs, histbatch.train_losses, label = "training loss")
                axs[2].plot(epochs, histbatch.val_losses, label = "validation loss", color = "tab:orange")
                axs[1].legend()
                axs[2].legend()
                PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-loss_subplots.png"))
            end
        end

        visualize(normal_[length(normal_)], string(folder_name, "/", trials[folder_name], "-normal_graphviz"))
        visualize(reduce_[length(reduce_)], string(folder_name, "/", trials[folder_name], "-reduce_graphviz"))
        PyPlot.close("all")
    else
        @show file_name = string(folder_name, "/histeval.bson")
        BSON.@load file_name histeval
        fig, axs = plt.subplots(2)
        axs[1].plot(histeval.train_losses, label = "training loss")
        axs[2].plot(histeval.accuracies, label = "accuracy", color = "tab:green")
        axs[1].legend()
        axs[2].legend()
        PyPlot.savefig(string(folder_name, "/", trials[folder_name], "-plots.png"))
    end
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
using OpenCV
using Base.Test

include(joinpath(Pkg.dir("OpenCV"), "./test/jl/tests.jl"))


text:
using PooksoftOptionsKit
using PyPlot
using Dates
using Distributions

# simulation logic -
function simulate(pathToSimulationFile::String, currentClose::Float64)

    # initialize -
    lattice_model = nothing
    asset_set = nothing
    asset_price_array = nothing
    data_array_at_exp = nothing
    prob_array = Array{Float64,1}()

    # load the lattice model -
    result = build_simulation_lattice_data_structure(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        lattice_model = result.value
    end

    # load the asset set -
    result = build_simulation_asset_set(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        asset_set = result.value
    end

    # load asset price array -
    result = build_simulation_price_array(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        asset_price_array = result.value
    end

    # run simulation at expiration -
    # compute -
    result = compute_option_profit_and_loss_at_expiration(asset_set, asset_price_array)
    if (typeof(result.value) == PSError)
        return result.value
    else
        data_array_at_exp = result.value
    end

    # check - are we in the money?
    iv = lattice_model.volatility
    d = Laplace(current_close,iv)
    for value in asset_price_array
        prob = 1 - cdf(d,value)
        push!(prob_array,prob)
    end

    # return -
    return (data_array_at_exp, prob_array)
end

# setup -
current_close = 7.23
path_to_design_file = "$(pwd())/tests/config/IronCondor-F.json"
(data_array_at_exp,pa) = simulate(path_to_design_file, current_close)


text:
using Test
using Distributions

@testset "Truncated normal, mean and variance" begin
    d = Distributions.TruncatedNormal(0,1,100,115); @test mean(d) ≈ 100.00999800099926070518490239457545847490332879043
    d = Distributions.TruncatedNormal(0,1,50,70); @test var(d) ≈ 0.00039904318680389954790992722653605933053648912703600
    d = Distributions.TruncatedNormal(-2,3,50,70); @test mean(d) ≈ 50.171943499898757645751683644632860837133138152489
    d = Distributions.TruncatedNormal(-2,3,50,70); @test var(d) ≈ 0.029373438107168350377591231295634273607812172191712
end


text:
module adapters_postgresql_transactions_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo Schema Raw

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

Repo.debug_sql()

include("options.jl")

# @test_throws Repo.NeedsConnectError Repo.get(Employee, 2)

Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    Options.for_postgresql...
)

Repo.execute([DROP TABLE IF EXISTS Employee])
Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID BIGSERIAL,
        Name VARCHAR(255),
        Salary FLOAT(8),
        PRIMARY KEY (ID)
    )"""))

multiple_changes = [
    (Name="Tom", Salary=20000.25),
    (Name="Jim", Salary=30000.00),
]
Repo.insert!(Employee, multiple_changes)

df = Repo.query(Employee)
@test size(df) == (2,)

Repo.delete!(Employee, (Name="Tom",))

df = Repo.query(Employee)
@test size(df) == (1,)

Repo.execute([BEGIN])
Repo.insert!(Employee, (Name="Tom", Salary=20000.25))
Repo.execute([SAVEPOINT :my_savepoint])
try
    Repo.insert!(Employee, (Name2="Tom", Salary=20000.25))
catch
    Repo.execute([ROLLBACK TO SAVEPOINT :my_savepoint])
end
Repo.execute([COMMIT])

df = Repo.query(Employee)
@test size(df) == (2,)

Repo.disconnect()

end # module adapters_postgresql_transactions_test


text:
using OMEinsum, OMEinsumContractionOrders
using Test, Random

@testset "slicer" begin
    log2_sizes = [1, 2,3 ,4.0]
    s = Slicer(log2_sizes, 3)
    push!(s, 1)
    @test_throws AssertionError push!(s, 1)
    push!(s, 2)
    push!(s, 3)
    @test_throws AssertionError push!(s, 4)
    replace!(s, 1=>4)
    @test s.log2_sizes == [1, 0.0, 0.0, 0.0]
    @test s.legs == Dict(2=>2.0, 3=>3.0, 4=>4.0)
    @test_throws AssertionError replace!(s, 1=>4)
end

@testset "slice iterator" begin
    se = SlicedEinsum(Slicing(['i', 'l']), ein"(ij,jk),(kl,lm)->im")
    it = OMEinsumContractionOrders.SliceIterator(se, uniformsize(se, 2))
    @test length(it) == 4
    for (i, v) in enumerate(it)
        @test v == it[i]
    end
end

@testset "SlicedEinsum" begin
    se = SlicedEinsum(Slicing(['i', 'l']), ein"(ij,jk),(kl,lm)->im")
    @test OMEinsum.flatten(se) == OMEinsum.flatten(se.eins)
    @test OMEinsum.labeltype(se) == Char
    xs = (randn(2,3), randn(3,4), randn(4,5), randn(5,6))
    size_info = Dict{Char,Int}()
    @test OMEinsum.get_size_dict!(se, xs, size_info) == Dict('i'=>2, 'j'=>3, 'k'=>4, 'l'=>5, 'm'=>6)
    @test OMEinsum.getixsv(se) == [['i','j'],['j','k'],['k','l'],['l','m']]
    @test OMEinsum.getiyv(se) == ['i','m']
    @test OMEinsum.label_elimination_order(se) == ['j','l', 'k']
    @test se(xs...) ≈ se.eins(xs...)
    @test uniquelabels(se) == ['i', 'j', 'k', 'l', 'm']
    @test uniformsize(se, 2) == Dict(zip(['i', 'j', 'k', 'l', 'm'], ones(Int, 5).*2))
end