text:
using Graphs

function create_graph(start_node, end_node)
    @assert length(start_node)==length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]

        link_idx = link_dic[parent,current]

        if link_idx != 0
            x[link_idx] = 1
        end

        current = parent
    end

    return x
end

# preparing a graph
graph = create_graph(start_node, end_node);
link_dic = sparse(start_node, end_node, 1:numLinks);

function BPR(flowVec, fcoeffs)
    bpr = similar(flowVec)
    for a = 1:length(bpr)
        bpr[a] = free_flow_time[a] * sum([fcoeffs[i] * (flowVec[a]/capacity[a])^(i-1) for i = 1:length(fcoeffs)])
        if bpr[a] < free_flow_time[a]
            bpr[a] = free_flow_time[a]
        end
    end
    return bpr
end

function all_or_nothing(travel_time, travel_demand_dict)
    state = []
    path = []
    x = zeros(size(start_node))

    for r=1:numNodes
        # for each origin node r, find shortest paths to all destination nodes
        state = dijkstra_shortest_paths(graph, travel_time, r)

        for s=1:numNodes
            # for each destination node s, find the shortest-path vector
            # load travel demand
            x = x + travel_demand_dict[(r,s)] * get_vector(state, r, s, link_dic)
        end
    end

    return x
end

function tapMSA_Multi(travel_demand_car, travel_demand_truck, fcoeffs, numIter=1000, tol=1e-6)

    # Finding a starting feasible solution
    travel_time = BPR(zeros(numLinks), fcoeffs);

    travel_time_car = 1.0 * BPR(zeros(numLinks), fcoeffs);
    travel_time_truck = 1.1 * BPR(zeros(numLinks), fcoeffs);

    xl_car = all_or_nothing(travel_time_car, travel_demand_car);
    xl_truck = all_or_nothing(travel_time_truck, travel_demand_truck);

    l = 1

    while l < numIter
        l += 1

        # Finding yl_car
        xl_car_old = xl_car
        travel_time_car = 1.0 * BPR(xl_car + 2.0 * xl_truck, fcoeffs)
        yl_car = all_or_nothing(travel_time_car, travel_demand_car);
        xl_car = xl_car + (yl_car - xl_car)/l
        xl_car_new = xl_car

        # Finding yl_truck
        xl_truck_old = xl_truck
        travel_time_truck = 1.1 * BPR(xl_car + 2.0 * xl_truck, fcoeffs)
        yl_truck = all_or_nothing(travel_time_truck, travel_demand_truck);
        xl_truck = xl_truck + (yl_truck - xl_truck)/l
        xl_truck_new = xl_truck

        relative_gap_car = norm(xl_car_new - xl_car_old, 1) / norm(xl_car_new, 1)
        relative_gap_truck = norm(xl_truck_new - xl_truck_old, 1) / norm(xl_truck_new, 1)

        if relative_gap_car < tol && relative_gap_truck < tol
            break
        end
    end

    tapFlows = Dict{}()
    tapFlowsCar = Dict{}()
    tapFlowsTruck = Dict{}()
    tapFlowVect = zeros(2,length(xl_car))

    for i = 1:length(ta_data.start_node)
        key = (ta_data.start_node[i], ta_data.end_node[i])
        tapFlowsCar[key] = xl_car[i]
        tapFlowsTruck[key] = xl_truck[i]
    end
    tapFlows["car"] = tapFlowsCar
    tapFlows["truck"] = tapFlowsTruck

    for j = 1:length(xl_car)
        tapFlowVect[1,j] = xl_car[j]
        tapFlowVect[2,j] = xl_truck[j]
    end

    return tapFlows, tapFlowVect

end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
module JsonGrinder
using Mill, JSON, Printf, Flux
using Mill: paddedprint, COLORS
include("schema.jl")
include("html_show_tools.jl")

using Mill: ArrayNode, BagNode, TreeNode, catobs
include("extractors/extractarray.jl")
include("extractors/extractbranch.jl")
include("extractors/extractcategorical.jl")
include("extractors/extractscalar.jl")
include("extractors/extractstring.jl")
include("extractors/extractvector.jl")
include("extractors/extractonehot.jl")
include("extractors/multirepresentation.jl")

export ExtractScalar, ExtractCategorical, ExtractArray, ExtractBranch, ExtractOneHot, ExtractVector, MultipleRepresentation, ExtractString
export suggestextractor, schema, extractbatch, generate_html
end # module

ext: jl


text:
"Smallest allowed graph weight/distance (instead of 0, which does not play well with sparse graphs.)"
const ϵ = 1e-6

"Lower threshold for Euclidean distance calculations. Distances less than this are treated as 0."
const L2THR = 1e-7

"""
    EMD

Earth Mover's Distance, a.k.a. 1-p Monge-Wasserstein distance. The algorithm as
implemented here assumes *balanced* transport where both vectors have the same statistical mass
(they both sum to 1). If they are not already, input vectors will be normalized to sum to 1.

"""
struct EMD <: SemiMetric
    grids::Union{Nothing,Tuple}
end

"EMD distance with with no grid provided. Grids default to axes(u,1) and axes(v,1)."
EMD() = EMD(nothing)

"""

    EMD(u::AbstractVector{T}, v::AbstractVector{T}) where {T <: Real}

Calculate the Earth Mover Distance (EMD) a.k.a the 1-Wasserstein distance
between the two given vectors, accepting a default grid. `u` and `v` are treated as
weights on the grid. The default grid is equal to the first axis of `u` and `v`.

```@example 1
u = rand(10);
u .= u / sum(u);
sum(u)
```

`EMD` implements the Distances package convention of runnable types:

```@example 1
v = rand(10);
v .= v / sum(v);
d = EMD()(u,v)
```

"""
function EMD(u::AbstractVector{T}, v::AbstractVector{T}) where {T <: Real}
    gridu = collect(first(axes(u)))
    gridv = collect(first(axes(v)))
    EMD(gridu, gridv, u, v)
end

"""

    EMD(u,v,uw,vw)

Calculate the Earth Mover's Distance using an explicit grid. This method is not intended to be called directly.
Instead, specify a grid in the call to sequence, with the WASS1D constant.

```julia

julia> A = rand(50, 100)
50×100 Array{Float64,2}:
[...]
julia> g = collect(0.5:0.5:div(size(A,1),2))
50-element Array{Float64,1}:
   0.5
   1.0
   1.5
  [...]
  24.0
  24.5
  25.0
julia> sequence(A; grid = g)
┌ Info: Sequencing data with
│     shape: (50, 100)
│     metric(s): (Euclidean(0.0), EMD(nothing), KLDivergence(), Energy(nothing))
└     scale(s): (1, 2, 4)
[ Info: `Euclidean(0.0) at scale 1: η = 5.214 (3.4s)
[...]
```
"""
function EMD(u,v,uw,vw)
    ndims(u) == 1 && ndims(u) == 1 || error("Only 1-d data vectors are supported.")
    length(u) == length(uw) || error("u and u weights must have same length. Got u $(length(u)) and uw $(length(uw))")
    length(v) == length(vw) || error("v and v weights must have same length. Got v $(length(v)) and vw $(length(vw))")
    _cdf_distance(ecdf(float(u); weights = float(uw)), ecdf(float(v); weights = float(vw)), 1)
end

"Same as `EMD(u,v)` but using Distances-style runnable type syntax."
(m::EMD)(u,v=u) = (m.grids === nothing) ? EMD(u,v) : EMD(m.grids...,u,v)

"Same as `EMD(u,v,uw,vw)` but using Distances-style runnable type syntax."
(::EMD)(u,v,uw,vw) = EMD(u,v,uw,vw)

"""
Convenience method for `EMD(u,v)`.

See [`EMD(u,v)`](@ref)
"""
emd(u::AbstractArray, v = u) = EMD(u, v)

"""
Convenience method for `EMD(u,v,uw,vw)`.

See [`EMD(u,v,uw,vw)`](@ref)
"""
emd(u,v,uw,vw) = EMD(u,v,uw,vw)

# ******** ENERGY *********

"""
    Energy

Energy distance as defined by Székely. An explicit grid may be provided. Default is
`axes(u,1)` and `axes(v,1)`.

[Wikipedia page on Energy distance](https://en.wikipedia.org/wiki/Energy_distance)

"""
struct Energy <: SemiMetric
    grids::Union{Nothing,Tuple}
end

"Default constructor with no specified grid."
Energy() = Energy(nothing)

"Enables dispatching on the type directly"
(m::Energy)(u,v) = (m.grids === nothing) ? Energy(u,v) : Energy(m.grids...,u,v)

"Direct dispatch on an object of type Energy."
(::Energy)(u,v,uw,vw) = Energy(u,v,uw,vw)

"""

    Energy(u,v,uw,vw)

Calculate Székely's energy distance between the two given vectors, `u` and `v`,
whose weights `uw`, `vw` are treated as a empirical cumulative distribution function (CDF).
`u` and `v` must have the same length, respectively, as `uw` and `vw`.
"""
function Energy(u,v,uw,vw)
    length(u) ==  length(uw) || error("u and u weights must have same length. Got u $(length(u)) and uw $(length(uw))")
    length(v) ==  length(vw) || error("v and v weights must have same length. Got v $(length(v)) and vw $(length(vw))")
    ndims(u) == 1 && ndims(v) == 1 || error("Only 1-d data vectors are supported.")
    sqrt(2) * _cdf_distance(ecdf(float(u); weights = float(uw)), ecdf(float(v); weights = float(vw)), 2)
end

"""

    Energy(u,v)

Calculate Székely's energy distance between the two given vectors, accepting a default grid.
`u` and `v` are treated as weights on the grid. The default grid is equal to the first axis of `u` and `v`.
"""
function Energy(u::AbstractVector{T}, v::AbstractVector{T}) where {T <: Real}
    gridu = collect(first(axes(u)))
    gridv = collect(first(axes(v)))
    Energy(gridu,gridv,u,v)
end

"""
Convenience function for Energy with a supplied grid.
"""
energy(u, v, uw, vw) = Energy(u, v, uw, vw)

"""
Convenience function for Energy with a default unit grid.
"""
energy(u::AbstractVector, v=u) = Energy(u,v)

"""
Returns the statistical distance between two cumulative
distributions u et v. With dimension parameter p = 1,
the distance is the 1 Wasserstein (or Earth Mover's) Distance.
With p=2, the distance is Energy.

References (from the stats.py source code)
>    [1] "Energy distance", https://en.wikipedia.org/wiki/Energy_distance
>    [2] Szekely "E-statistics: The energy of statistical samples." Bowling
>        Green State University, Department of Mathematics and Statistics,
>        Technical Report 02-16 (2002).
>    [3] Rizzo, Szekely "Energy distance." Wiley Interdisciplinary Reviews:
>        Computational Statistics, 8(1):27-38 (2015).
>    [4] Bellemare, Danihelka, Dabney, Mohamed, Lakshminarayanan, Hoyer,
>        Munos "The Cramer Distance as a Solution to Biased Wasserstein
>        Gradients" (2017). :arXiv:`1705.10743`.
"""
function _cdf_distance(u::ECDF, v::ECDF, p::Int=1)

# values are pre-sorted ascending
    uv = u.sorted_values
    vv = v.sorted_values
    # merge and sort all grid values
    # We insert a zero to ensure that the first element's
    # grid distance survives the diff operation
    uplusv = vcat(0., uv, vv)
    sort!(uplusv)
    @debug "uplusv" uplusv
    Δx = diff(uplusv, dims = 1)
    @debug "Delta Force!" Δx

# Shortcut return in case the grids are identical
    if uv == vv # we can safely remove zeros from the delta
        Δx = filter(v->!iszero(v), Δx)
        UA = u(uv) # calculate u,v CDFs on their separate grids
        VA = v(vv)
    else
        A = uplusv[1:end-1]
        UA = u(A) # Calculate CDFs on whole domain
        VA = v(A)
    end

    if p == 1
        return sum(abs.(UA .- VA) .* Δx)
    elseif p == 2
        return sqrt(sum( (UA .- VA) .^ 2 .* Δx ))
    else
        return sum(abs.((UA .- VA)) .^ p .* Δx) ^ 1/p
    end
end

# *********** CONSTANTS ************

"""
Constant to use when specifying the Euclidean or L2 distance metric for a sequencing run. This
sets the default threshold to the value of [`L2THR`](@ref) = 1e-7.

To set a different threshold level, use the type constructor `Distances.Euclidean(<your value>)` directly
(i.e. instead of the `L2` constant).
"""
const L2 = Distances.Euclidean(L2THR)

"""
Constant for specifiying the Kullbach-Leibler Divergence metric.
"""
const KLD = Distances.KLDivergence()

"""
Constant for specifiying the Monge-Wasserstein a.k.a. 1-p Wasserstein a.k.a. Earth Mover's Distance (EMD) metric.
This metric is sensitive to the underlying grid. Default is unit grid taken from the axes of the data vector.

See [EMD]@ref
"""
const WASS1D = EMD()

"""
Constant for specifiying the Energy distance as defined by Székely.

See [Energy](@ref)
"""
const ENERGY = Energy()

"Convenience constant to represent all of L2, WASS1D, KLD, and ENERGY."
const ALL_METRICS = (L2, WASS1D, KLD, ENERGY)

"Less verbose output for PeriodicEuclidean."
Base.show(io::IO, s::PeriodicEuclidean) = write(io, "PeriodicEuclidean($(length(s.periods))) [$(minimum(s.periods)),$(maximum(s.periods))]")


text:
include("Linear.jl")
include("jacobian.jl")
using Base.Test

X = map(Float32, randn(5, 3))
l = Linear{Float32}(3, 4)

# Test grad. wrt. parameter W
Dfwd = jacobian_fwd(l, X, l.w)
Dbwd = jacobian_bwd(l, X, l.w, l.∂w)

@test_approx_eq Dfwd Dbwd

# Test grad. wrt. parameter b
Dfwd = jacobian_fwd(l, X, l.b)
Dbwd = jacobian_bwd(l, X, l.b, l.∂b)

@test_approx_eq Dfwd Dbwd

# Test grad. wrt. input X
Dfwd = jacobian_fwd(l, X)
Dbwd = jacobian_bwd(l, X)

@test_approx_eq Dfwd Dbwd


text:
using Meganet
using Base.Test

for TYPE=[Float64,Float32]
    S = getSoftMaxLoss(TYPE)
    display(S)

    @testset "softmax (with bias) $TYPE" begin
        testLossFunction(S,out=false)
    end

    S = getSoftMaxLoss(TYPE)
    S.addBias = false
    @testset "softmax (no bias) $TYPE" begin
        testLossFunction(S,out=false)
    end
end


text:
function bridge_sampling(
      samples::Array{Float64,2},
      log_density::Function;
      verbose = false,
)
      nr = size(samples, 1) ÷ 2
      samples_4_fit = permutedims(Array(samples[1:nr, :]))
      samples_4_iter = permutedims(Array(samples[(nr+1):end, :]))
      r0 = 0.5
      tol1 = 1e-10
      n_post = size(samples_4_iter, 2)
      m = vec(mapslices(mean, samples_4_fit, dims = 2))::Vector{Float64}
      V = diagm(ones(size(samples, 2)))
      q11 = Vector{Float64}(undef, n_post) # unnormalized posterior
      @threads for i = 1:n_post
            q11[i] = log_density(Array(samples_4_iter[:, i]))
      end
      q12 = Vector{Float64}(undef, n_post) # g
      @threads for i = 1:n_post
            q12[i] = logpdf(MvNormal(m, V), samples_4_iter[:, i])
      end
      gen_samples = rand(MvNormal(m, V), n_post)
      q21 = Vector{Float64}(undef, n_post) # unnormalized posterior
      @threads for i = 1:n_post
            q21[i] = log_density(gen_samples[:, i])
      end
      q22 = Vector{Float64}(undef, n_post) # g
      @threads for i = 1:n_post
            q22[i] = logpdf(MvNormal(m, V), gen_samples[:, i])
      end
      l1 = q11 - q12 # samples_4_iter, unnormalized posterior / g
      l2 = q21 - q22# gen_samples, unnormalized posterior / g
      lstar = median(l1)
      n_1 = length(l1)
      n_2 = length(l2)
      neff = (1.0 * n_1)::Float64
      try
            neff = minimum(summarystats(sim).value[:, end, 1])
      catch e
      end
      s1 = neff ./ (neff .+ n_2)
      s2 = n_2 ./ (neff .+ n_2)
      logr = log(r0)
      log_r_vals = [logr]::Vector{Float64}
      logml = logr + lstar

      function upd(logr)
            lognumi = logsumexp([
                  l2[i] - lstar - logaddexp(log(s1) + l2[i] - lstar, log(s2) + logr) for
                  i = 1:size(l2, 1)
            ])
            logdeni = logsumexp([
                  -logaddexp(log(s1) + l1[i] - lstar, log(s2) .+ logr) for
                  i = 1:size(l1, 1)
            ])
            log(n_1) - log(n_2) + lognumi - logdeni
      end

      ft = optimize(y -> (upd(y)-y)^2, -10000., 0.)
      ft.minimizer + lstar

end

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Distributions, StatsBase, IntervalSets, ValueShapes, ArraysOfArrays

@testset "hierarchial_distribution" begin
    let
        primary_dist = NamedTupleDist(
            foo = LogNormal(1, 0.3),
            bar = Normal(2.0, 1.0)
        )

        f = v -> NamedTupleDist(baz = fill(Normal(v.bar, v.foo), 3))

        @test @inferred(HierarchicalDistribution(f, primary_dist)) isa HierarchicalDistribution
        hd = HierarchicalDistribution(f, primary_dist)

        @test @inferred(unshaped(hd)) isa BAT.UnshapedHDist
        ud = unshaped(hd)

        @test @inferred(rand(hd)) isa NamedTuple
        @test @inferred(rand(hd, ())) isa ShapedAsNT
        @test @inferred(rand(ud)) isa AbstractVector{<:Real}
        @test @inferred(varshape(hd)) == NamedTupleShape(foo = ScalarShape{Real}(), bar = ScalarShape{Real}(), baz = ArrayShape{Real}(3))
        @test @inferred(varshape(ud)) == ArrayShape{Real}(5)

        ux = [2.7, 4.3, 8.7, 8.7, 8.7]
        @test @inferred(logpdf(ud, ux)) ≈ logpdf(primary_dist.foo, 2.7) + logpdf(primary_dist.bar, 4.3) + 3 * logpdf(Normal(4.3, 2.7), 8.7)
        @test @inferred(logpdf(hd, varshape(hd)(ux))) == logpdf(ud, ux)
        @test @inferred(logpdf(hd, BAT.strip_shapedasnt(varshape(hd)(ux)))) == logpdf(ud, ux)

        samples = bat_sample(hd, MCMCSampling(mcalg = HamiltonianMC(), trafo = NoDensityTransform(), nsteps = 10^4)).result
        @test isapprox(cov(unshaped.(samples)), cov(ud), rtol = 0.25)
    end

    let
        hd = HierarchicalDistribution(
            v -> NamedTupleDist(b = Normal(v.a, 1.2)),
            NamedTupleDist(a = Normal(2.3, 1.9))
        )

        cov_expected = [1.9^2 1.9^2; 1.9^2 1.9^2 + 1.2^2]

        @test isapprox(cov(unshaped(hd)), cov_expected, rtol = 0.05)
        @test isapprox(mean(unshaped.(rand(sampler(hd), 10^5))), [2.3, 2.3], rtol = 0.05)
    end
end


text:
using FastAnonymous
using Base.Test

immutable RGB
    r::Float64
    g::Float64
    b::Float64
end

offset = 1.2
f = @anon x->sqrt(x+offset)
@test_approx_eq f(7.8) 3.0

to_gray(x::Real) = x
to_gray(x::RGB) = 0.299*x.r + 0.587*x.g + 0.114*x.b
gray = @anon x->to_gray(x)
@test gray(0.8) == 0.8
@test gray(RGB(1,0,0)) == 0.299

h = @anon x->sin(x^2)
@test h(2.1) == sin((2.1)^2)

# Check that we are well optimized
# (allocation == 0 means no boxing)
x = rand(10)
y = similar(x)
map!(f, y, x)
@test (@allocated map!(f, y, x)) == 0

c = reinterpret(RGB, rand(30))
map!(gray, y, c)
@test (@allocated map!(gray, y, c)) == 0
map!(gray, y, x)
@test (@allocated map!(gray, y, x)) == 0

# Does it work in a function?
function testwrap(dest, A, offset)
    ff = @anon x->(x+offset)^2
    map!(ff, dest, A)  # compile map! with ff
    @allocated map!(ff, dest, A)
end

@test testwrap(y, x, 1.2) == 0
@test y[1] == (x[1]+1.2)^2
@test testwrap(y, x, 3.4) == 0
@test y[1] == (x[1]+3.4)^2

# Does it work across modules?
module TestMod
mysquare(x) = x^2
end

f = @anon x->TestMod.mysquare(x)
@test f(3) == 9

module TestMod2
export mycube
mycube(x) = x^3
end

using TestMod2
f = @anon x->mycube(x)
@test f(-2) == -8


text:
module PCBs

using Circuits

function translate!(sc::SubCircuit; dx=0, dy=0)
    for e in sc.c.elements
        for (d,k) in [(dx,:x), (dx,:centerx), (dy,:y), (dy,:centery)]
            if k ∈ keys(e.meta)
                e.meta[k] += d
            end
        end
    end
    sc
end

translate(obj; kwargs...) = translate!(copy(obj); kwargs...)

module Lisp
include("kicad_lisp.jl")
end
using .Lisp

module Settings
general() = Dict{Symbol,Number}(
    :thickness => 1.6,
    :drawings => 0,
    :tracks => 0,
    :zones => 0,
    :modules => 1,
    :nets => 1)

layers() = Dict{Int,Vector{Symbol}}(
    0 => [Symbol("F.Cu"), :signal],
    31 => [Symbol("B.Cu"), :signal],
    32 => [Symbol("B.Adhes"), :user],
    33 => [Symbol("F.Adhes"), :user],
    34 => [Symbol("B.Paste"), :user],
    35 => [Symbol("F.Paste"), :user],
    36 => [Symbol("B.SilkS"), :user],
    37 => [Symbol("F.SilkS"), :user],
    38 => [Symbol("B.Mask"), :user],
    39 => [Symbol("F.Mask"), :user],
    40 => [Symbol("Dwgs.User"), :user],
    41 => [Symbol("Cmts.User"), :user],
    42 => [Symbol("Eco1.User"), :user],
    43 => [Symbol("Eco2.User"), :user],
    44 => [Symbol("Edge.Cuts"), :user],
    45 => [Symbol("Margin"), :user],
    46 => [Symbol("B.CrtYd"), :user],
    47 => [Symbol("F.CrtYd"), :user],
    48 => [Symbol("B.Fab"), :user],
    49 => [Symbol("F.Fab"), :user])

setup() = Dict{Symbol,Any}(
    :last_trace_width => 0.25,
    :trace_clearance => 0.2,
    :zone_clearance => 0.508,
    :zone_45_only => :no,
    :trace_min => 0.2,
    :segment_width => 0.2,
    :edge_width => 0.15,
    :via_size => 0.8,
    :via_drill => 0.4,
    :via_min_size => 0.4,
    :via_min_drill => 0.3,
    :uvia_size => 0.3,
    :uvia_drill => 0.1,
    :uvias_allowed => :no,
    :uvia_min_size => 0.2,
    :uvia_min_drill => 0.1,
    :pcb_text_width => 0.3,
    :pcb_text_size => [1.5, 1.5],
    :mod_edge_width => 0.15,
    :mod_text_size => [1, 1],
    :mod_text_width => 0.15,
    :pad_size => [1.524, 1.524],
    :pad_drill => 0.762,
    :pad_to_mask_clearance => 0.051,
    :solder_mask_min_width => 0.25,
    :aux_axis_origin => [0, 0],
    :visible_elements => 0xFFFFFF7F,
    :pcbplotparams => Dict{Symbol, Any}(
        :layerselection => Symbol("0x010fc_ffffffff"),
        :usegerberextensions => false,
        :usegerberattributes => false,
        :usegerberadvancedattributes => false,
        :creategerberjobfile => false,
        :excludeedgelayer => true,
        :linewidth => 0.100000,
        :plotframeref => false,
        :viasonmask => false,
        :mode => 1,
        :useauxorigin => false,
        :hpglpennumber => 1,
        :hpglpenspeed => 20,
        :hpglpendiameter => 15.000000,
        :psnegative => false,
        :psa4output => false,
        :plotreference => true,
        :plotvalue => true,
        :plotinvisibletext => false,
        :padsonsilk => false,
        :subtractmaskfromsilk => false,
        :outputformat => 1,
        :mirror => false,
        :drillshape => 1,
        :scaleselection => 1,
        :outputdirectory => ""
    )
)

net_classes() = Dict{Tuple{Symbol,String},Dict{Symbol,Float64}}(
    (:Default,"This is the default net class.") => Dict{Symbol,Float64}(
        :clearance => 0.2,
        :trace_width => 0.25,
        :via_dia => 0.8,
        :via_drill => 0.4,
        :uvia_dia => 0.3,
        :uvia_drill => 0.1
    )
)

zone_settings() = Dict(:tstamp => 0,
                       :hatch => [:edge, 0.508],
                       :connect_pads => [:yes, Dict(:clearance => 0.508)],
                       :min_thickness => 0.254,
                       :fill => Dict(:arc_segments => 16,
                                    :thermal_gap => 0.508,
                                    :thermal_bridge_width => 0.508))

end

mutable struct Zone
    name::Union{Symbol,String}
    layer::Symbol
    polygon::Matrix{Float64}
    settings::Dict
end
function Zone(name, layer, polygon; settings = Settings.zone_settings())
    Zone(name, Symbol(layer), polygon, settings)
end

rectangular_zone(name, layer, x, y, w, h; kwargs...) =
    Zone(name, layer, [x y; x y+h; x+w y+h; x+w y]; kwargs...)

mutable struct Segment
    layer::Symbol
    net::Union{Symbol,String}
    coordinates::Matrix{Float64}
    width::Float64
end
Segment(layer::Union{String,Symbol}, net::Union{String,Symbol},
        x, y, dx, dy, width) =
            Segment(Symbol(layer), net, [x y; x+dx y+dy], width)

Base.copy(segment::Segment) =
    Segment(segment.layer, segment.net, copy(segment.coordinates), segment.width)

function translate!(segment::PCBs.Segment; dx=0, dy=0)
    segment.coordinates[:,1] .+= dx
    segment.coordinates[:,2] .+= dy
    segment
end

struct PCB
    circuit::Circuit
    filename::String
    directories::Vector{String}
    version::String
    host::String
    general::Dict{Symbol,Number}
    page::String
    layers::Dict{Int,Vector{Symbol}}
    setup::Dict{Symbol,Any}
    net_classes::Dict{Tuple{Symbol,String},Dict{Symbol,Float64}}
    zones::Vector{Zone}
    segments::Vector{Segment}
end

PCB(circuit::Circuit, filename::String;
    directories=[],
    version="20171130",
    host="(5.0.1-3-g963ef8bb5)",
    general=Settings.general(),
    page="A4",
    layers=Settings.layers(),
    setup=Settings.setup(),
    net_classes=Settings.net_classes(),
    zones=Zone[], segments=Segment[]) =
        PCB(circuit, filename, directories,
            version, host,
            general, page, layers, setup, net_classes,
            zones, segments)

function settings_dict_to_lisp(name::String, settings::Dict, sorted::Bool=false;
                               ks=collect(keys(settings)))
    node = Lisp.KNode(Lisp.KSym(name))
    sorted && sort!(ks)
    children = map(ks) do k
        v = settings[k]
        child = Lisp.KNode(Lisp.obj(k))
        if v isa Vector
            for e in v
                if e isa Dict
                    s = settings_dict_to_lisp("tmp", e)
                    append!(child.children, s.children[2:end])
                else
                    push!(child.children, Lisp.obj(e))
                end
            end
        elseif v isa Dict
            child = settings_dict_to_lisp(string(k), v)
        else
            push!(child.children, Lisp.obj(v))
        end
        child
    end
    append!(node.children, children)
    node
end

function define_nets(pcb::PCB)
    net_nodes = Lisp.KNode[]
    nets = node_names(pcb.circuit)
    unique_nets = unique(nets)

    pcb.general[:nets] += length(unique_nets)

    push!(net_nodes, settings_dict_to_lisp("general", pcb.general))
    append!(net_nodes, Lisp.Read("(page $(pcb.page))"))
    push!(net_nodes, settings_dict_to_lisp("layers", pcb.layers, true))
    push!(net_nodes, settings_dict_to_lisp("setup", pcb.setup))

    append!(net_nodes, Lisp.Read("(net 0 \"\")"))
    for (i,net) in enumerate(unique_nets)
        net_node = Lisp.KNode(ksym"net", Lisp.KInt(i), Lisp.obj(net))
        push!(net_nodes, net_node)
    end

    for ((name,descr),v) in pcb.net_classes
        nc = Lisp.KNode(ksym"net_class", Lisp.KSym(name), Lisp.KStr(descr))
        append!(nc.children, settings_dict_to_lisp("net_class", v).children[2:end])
        if name == :Default
            for (i,net) in enumerate(unique_nets)
                append!(nc.children, Lisp.Read("(add_net \"$(net)\")"))
            end
        end
        push!(net_nodes, nc)
    end
    net_nodes,nets,unique_nets
end

function read_kicad_lisp_file(filename)
    isfile(filename) || throw(ArgumentError("$filename missing"))
    @debug "Reading $filename"
    open(Lisp.Read, filename)[1]
end

function find_footprint(footprint, directories)
    directories = copy(directories)
    push!(directories, ENV["KISYSMOD"])

    footprint_path = split(footprint, ":")
    footprint_path[1] = "$(footprint_path[1]).pretty"
    footprint_path[2] = "$(footprint_path[2]).kicad_mod"
    footprint_path = joinpath(footprint_path...)

    @debug "Finding footprint $footprint"
    for dir in directories
        candidate = joinpath(dir, footprint_path)
        @debug "Trying $candidate"
        isfile(candidate) && return read_kicad_lisp_file(candidate)
    end
    throw(ArgumentError("Footprint $footprint not found"))
end

xy(x, y) = Lisp.KNode(ksym"xy", Lisp.obj(x), Lisp.obj(y))
xyz(x, y, z) = Lisp.KNode(ksym"xyz", Lisp.obj(x), Lisp.obj(y), Lisp.obj(z))

function get_model(filename::String; at=[0,0,0], scale=[1,1,1], rotate=[0,0,0])
    [ksym"model", Lisp.KSym(filename),
     Lisp.KNode(ksym"at", xyz(at...)),
     Lisp.KNode(ksym"scale", xyz(scale...)),
     Lisp.KNode(ksym"rotate", xyz(rotate...))]
end

get_model(model::Pair{String,<:Dict{Symbol}}) = get_model(model[1];model[2]...)

function replace_layers!(footprint::Lisp.KNode, layers::Vector{Pair{Regex,String}})
    mirror = Lisp.Read("(justify mirror)")
    do_mirror = any(map(isequal("B"), last.(layers)))
    Lisp.traverse_nodes(footprint) do node
        fc = first(node.children)
        if fc isa Lisp.KSym
            if fc.name ∈ ["layer","layers"]
                node.children = copy(node.children)
                for i = 2:length(node.children)
                    for (pat,subs) in layers
                        m = match(pat, node.children[i].name)
                        if m != nothing
                            node.children[i] = Lisp.KSym("$(subs).$(m[1])")
                            break # We don't want the reverse pattern to apply as well
                        end
                    end
                end
            elseif fc.name == "fp_text"
                node.children = copy(node.children)
                Lisp.traverse_nodes(node, "effects") do text_prop
                    text_prop.children = copy(text_prop.children)
                    append!(text_prop.children, mirror)
                end
            end
        end
    end
end

function replace_layers!(footprint::Lisp.KNode, layers::Vector{Pair{String,String}})
    layers = vcat(map(l -> Regex("$(l[1])\\.(.*)") => l[2], layers),
                  map(l -> Regex("$(l[2])\\.(.*)") => l[1], layers))
    replace_layers!(footprint, layers)
end

replace_layers!(footprint::Lisp.KNode, layers::Pair{String,String}) =
    replace_layers!(footprint, [layers])

function modify_text!(footprint::Lisp.KNode, text::Dict{String,<:Dict})
    Lisp.traverse_nodes(footprint, "fp_text") do node
        text_name = node.children[2].name
        for t in keys(text)
            if text_name == t
                node.children = copy(node.children)
                for (i,n) in enumerate(node.children)
                    if n isa Lisp.KNode
                        s = Symbol(n.children[1].name)
                        s ∈ keys(text[t]) || continue
                        node.children[i] = Lisp.KNode(vcat(Lisp.KSym(s), Lisp.obj(text[t][s])))
                    end
                end
            end
        end
    end
end

function elements_to_lisp(pcb::PCB, nets, unique_nets)
    footprints = Dict{String, Lisp.KNode}()
    element_nodes = Lisp.KNode[]

    for e in pcb.circuit.elements
        f = try
            e.meta[:footprint]
        catch
            throw(ArgumentError("Footprint missing from $e"))
        end
        footprint = copy(get(footprints, f, find_footprint(f, pcb.directories)))
        ks = keys(e.meta)

        x = (:x ∈ ks ? e.meta[:x]
             : :centerx ∈ ks && :offsetx ∈ ks ? e.meta[:centerx] + e.meta[:offsetx]
             : 0)
        y = (:y ∈ ks ? e.meta[:y]
             : :centery ∈ ks && :offsety ∈ ks ? e.meta[:centery] + e.meta[:offsety]
             : 0)
        rot = :rot ∈ ks ? e.meta[:rot] : 0

        at = Lisp.KNode(ksym"at")
        for c in [x,y,rot]
            push!(at.children, Lisp.obj(c))
        end
        push!(footprint.children, at)

        pads = map(enumerate(pins(pcb.circuit, e))) do (i,pi)
            net = nets[pi]
            net_index = findfirst(isequal(net), unique_nets)
            @debug "Connecting pin $i of $e (global pin $pi) to net $(net) (#$(net_index))"
            net,net_index
        end
        for (i,c) in enumerate(footprint.children)
            if c isa Lisp.KNode && length(c.children) > 2
                if c.children[1] == ksym"pad" && c.children[2] isa Lisp.KInt
                    pad = pads[c.children[2].i]
                    c.children = copy(c.children)
                    push!(c.children,
                          Lisp.KNode(ksym"net", Lisp.KInt(pad[2]), Lisp.KStr(pad[1])))
                elseif c.children[1] == ksym"model" && :model ∈ ks
                    c.children = get_model(e.meta[:model])
                elseif c.children[1] == ksym"fp_text" &&
                    c.children[2] == ksym"reference"
                    c.children = copy(c.children)
                    c.children[3] = Lisp.KStr(string(e))
                end
            end
        end
        :layers ∈ ks && replace_layers!(footprint, e.meta[:layers])
        :text ∈ ks && modify_text!(footprint, e.meta[:text])

        push!(element_nodes, footprint)
    end
    element_nodes
end

function segments_to_lisp(pcb::PCB, nets, unique_nets)
    segment_nodes = Lisp.KNode[]
    for segment in pcb.segments
        cos = segment.coordinates
        for i = 1:size(cos,1)-1
            push!(segment_nodes,
                  Lisp.KNode(ksym"segment",
                             Lisp.KNode(ksym"start", Lisp.obj.(cos[i,:])...), Lisp.KNode(ksym"end", Lisp.obj.(cos[i+1,:])...),
                             Lisp.KNode(ksym"width", Lisp.obj(segment.width)),
                             Lisp.KNode(ksym"layer", Lisp.obj(segment.layer)),
                             Lisp.KNode(ksym"net", Lisp.obj(findfirst(isequal(segment.net), unique_nets)))))
        end
    end
    segment_nodes
end

function zones_to_lisp(pcb::PCB, nets, unique_nets)
    zone_nodes = Lisp.KNode[]
    for zone in pcb.zones
        zone_node = settings_dict_to_lisp("zone", zone.settings)
        insert!(zone_node.children, 2, Lisp.KNode(ksym"layer", Lisp.KSym(zone.layer)))
        insert!(zone_node.children, 2, Lisp.KNode(ksym"net_name", Lisp.KSym(zone.name)))
        insert!(zone_node.children, 2, Lisp.KNode(ksym"net", Lisp.KInt(findfirst(isequal(zone.name), unique_nets))))
        polygon = Lisp.KNode(ksym"polygon", Lisp.KNode(ksym"pts"))
        for i = 1:size(zone.polygon,1)
            push!(polygon.children[2].children, xy(zone.polygon[i,:]...))
        end
        push!(zone_node.children, polygon)
        push!(zone_nodes, zone_node)
    end
    zone_nodes
end

function Base.show(io::IO, pcb::PCB)
    node = Lisp.Read("(kicad_pcb (version $(pcb.version)) (host pcbnew \"$(pcb.host)\"))")[1]

    net_nodes,nets,unique_nets = define_nets(pcb)
    append!(node.children, net_nodes)
    append!(node.children, elements_to_lisp(pcb, nets, unique_nets))
    append!(node.children, segments_to_lisp(pcb, nets, unique_nets))
    append!(node.children, zones_to_lisp(pcb, nets, unique_nets))

    show(io, node)
end

function save(pcb::PCB,overwrite::Bool=false)
    filename = "$(pcb.filename).kicad_pcb"
    if isfile(filename)
        if overwrite
            @warn "Overwriting existing $(filename)"
        else
            error("$(filename) already exists")
        end
    else
        @info "Saving $(pcb.circuit) to $filename"
    end
    open(filename, "w") do file
        show(file, pcb)
    end
end

export PCB, translate

end # module


text:
# Add tests to the quote for functions with  FastMath varients.
const FASTABLE_AST = quote
    @testset "Trig" begin
        @testset "Basics" for x = (Float64(π)-0.01, Complex(π, π/2))
            test_scalar(sin, x)
            test_scalar(cos, x)
            test_scalar(tan, x)
        end
        @testset "Hyperbolic" for x = (Float64(π)-0.01, Complex(π-0.01, π/2))
            test_scalar(sinh, x)
            test_scalar(cosh, x)
            test_scalar(tanh, x)
        end
        @testset "Inverses" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asin, x)
            test_scalar(acos, x)
            test_scalar(atan, x)
        end
        @testset "Multivariate" begin
            @testset "sincos" begin
                x, Δx, x̄ = randn(3)
                Δz = (randn(), randn())

                frule_test(sincos, (x, Δx))
                rrule_test(sincos, Δz, (x, x̄))
            end
        end
    end

    @testset "exponents" begin
        for x in (-0.1, 6.4)
            test_scalar(inv, x)

            test_scalar(exp, x)
            test_scalar(exp2, x)
            test_scalar(exp10, x)
            test_scalar(expm1, x)

            test_scalar(cbrt, x)

            if x >= 0
                test_scalar(sqrt, x)
                test_scalar(log, x)
                test_scalar(log2, x)
                test_scalar(log10, x)
                test_scalar(log1p, x)
            end
        end
    end

    @testset "Unary complex functions" begin
        for x in (-4.1, 6.4)
            test_scalar(abs, x)
            test_scalar(angle, x)
            test_scalar(abs2, x)
            test_scalar(conj, x)
        end
    end

    @testset "Unary functions" begin
        for x in (-4.1, 6.4)
            test_scalar(+, x)
            test_scalar(-, x)
        end
    end

    @testset "binary function ($f)" for f in (/, +, -, hypot, atan, rem, ^, max, min)
        x, Δx, x̄ = 10rand(3)
        y, Δy, ȳ = rand(3)
        Δz = rand()

        frule_test(f, (x, Δx), (y, Δy))
        rrule_test(f, Δz, (x, x̄), (y, ȳ))
    end

    @testset "sign" begin
        @testset "at points" for x in (-1.1, -1.1, 0.5, 100)
            test_scalar(sign, x)
        end

        @testset "Zero over the point discontinuity" begin
            # Can't do finite differencing because we are lying
            # following the subgradient convention.

            _, pb = rrule(sign, 0.0)
            _, x̄ = pb(10.5)
            @test extern(x̄) == 0

            _, ẏ = frule((Zero(), 10.5), sign, 0.0)
            @test extern(ẏ) == 0
        end
    end
end

# Now we generate tests for fast and nonfast versions
@eval @testset "fastmath_able Base functions" begin
    $FASTABLE_AST
end

@eval @testset "fastmath_able FastMath functions" begin
    $(Base.FastMath.make_fastmath(FASTABLE_AST))
end


text:
using Documenter
using Random
using Soss
using StableRNGs

DocMeta.setdocmeta!(Soss, :DocTestSetup,
    quote
        using Random
        using Soss
        using StableRNGs
        using MeasureTheory
        Random.seed!(3)
    end; recursive=true)

@testset "Doctests" begin
    doctest(Soss)
end


text:
using Debugger

function foo(n)
    x = n+1
    ((BigInt[1 1; 1 0])^x)[2,1]
end

foo(20)


text:
using BinDeps
using Compat

@BinDeps.setup

const librsync = library_dependency("librsync", aliases=["librsync"])

const librsyncName = "librsync-2.0.1"
provides(Sources,
         URI("https://github.com/librsync/librsync/archive/v2.0.1.zip"),
         librsync,
         unpacked_dir = librsyncName)

const prefix = joinpath(BinDeps.depsdir(librsync), "usr")
const srcdir = joinpath(BinDeps.depsdir(librsync), "src", librsyncName)
const builddir = joinpath(srcdir, "build")
const doDebug = true
const debugArg = doDebug ? "-DCMAKE_BUILD_TYPE=Debug " : "-DCMAKE_BUILD_TYPE=Release "
const buildRdiff = true
const RdiffArg = buildRdiff ? "-DBUILD_RDIFF=ON " : "-DBUILD_RDIFF=OFF "

provides(SimpleBuild,
  (@build_steps begin
    GetSources(librsync)
    CreateDirectory(builddir)
    @build_steps begin
      ChangeDirectory(builddir)
      `cmake -DCMAKE_INSTALL_PREFIX=$(prefix) $(debugArg)$(RdiffArg) ..`
      `make`
      `make check`
      `make install`
    end
  end),
  librsync,
  os = :Unix)

@BinDeps.install @compat Dict(:librsync => :librsync)


text:
# sample_target_mismatch-multivariate_uniform_vs_beta-kernel_vs_graph
#
# This script compares the values of the multivariate graph Stein
# discrepancy and kernel discrepancy obtained when samples are drawn i.i.d.
# from their targets.

include("experiment_utils.jl")

# we see first if we want to parallelize
@parseintcli parallel "p" "parallel" "false"
if parallel == "true"
    addprocs(CPU_CORES - 1)
end

using Distributions: Beta
using SteinDistributions:
    SteinUniform,
    gradlogdensity
using SteinDiscrepancy:
    SteinChampionLenardMillsKernel,
    ksd,
    gsd

# random number generator seed
# we allow this to be set on the commandline
@parseintcli seed "s" "seed" 7
# set the target distribution
@parsestringcli distname "t" "target" "uniform"
# set the discrepancy
@parsestringcli discrepancytype "d" "discrepancy" "championlenardmills"
# Specify dimension of sampled vectors
d = 2
# Specify a maximum number of samples
maxn = 200000
# uniform target
uniform_target = SteinUniform(d)
# beta target
beta_target = Beta(3, 3)
# the graphsolver
# setup kernel (if necessary)
graphsolver = nothing
kernel = nothing
if discrepancytype == "championlenardmills"
    kernel = SteinChampionLenardMillsKernel([0.0,1.0], 2)
elseif discrepancytype == "graph-gurobi"
    graphsolver = "gurobi"
elseif discrepancytype == "graph-clp"
    graphsolver = "clp"
end
# Draw samples from target distribution
if distname == "beta"
    target = beta_target
    X = @setseed reshape(rand(beta_target, d*maxn), maxn, d)
elseif distname == "uniform"
    target = uniform_target
    X = @setseed rand(uniform_target, maxn)
end
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# Sample sizes at which optimization problem will be solved
ns = vcat(10:1:min(10,maxn),
          100:100:min(1000,maxn),
          2000:1000:min(10000,maxn),
          20000:10000:min(100000,maxn),
          200000:100000:maxn)

# Solve optimization problem at each sample size
for i in ns
    @printf("[Beginning n=%d]:\n", i)
    # Create program with data subset
    Xi = X[1:i,:]
    discrepancy = nothing
    solvetime = nothing
    if kernel != nothing
        @printf("[Beginning kernel computations n=%d]\n", i)
        result = ksd(points=Xi, gradlogdensity=gradlogp, kernel=kernel)
        discrepancy = sqrt(result.discrepancy2)
        solvetime = result.solvetime
    elseif discrepancytype in ["graph-clp", "graph-gurobi"]
        @printf("[Beginning graph solver n=%d]\n", i)
        result = gsd(points=Xi,
                     gradlogdensity=gradlogp,
                     supportlowerbounds=zeros(d),
                     supportupperbounds=ones(d),
                     solver=graphsolver)
        discrepancy = sum(result.objectivevalue)
        solvetime = sum(result.solvetime)
    end
    # save data
    instance_data = Dict{Any, Any}(
        "target" => distname,
        "n" => i,
        "d" => d,
        "seed" => seed,
        "discrepancytype" => discrepancytype,
        "discrepancy" => discrepancy,
        "nprocs" => nprocs(),
        "solvetime" => solvetime,
    )
    save_json(
        instance_data;
        dir="sample_target_mismatch-multivariate_uniform_vs_beta-kernel_vs_graph",
        target=distname,
        discrepancytype=discrepancytype,
        n=i,
        d=d,
        ncores=nprocs(),
        seed=seed
    )
end

@printf("COMPLETE!")


text:
using DandelionWebSockets
using DandelionWebSockets: STATE_OPEN, STATE_CONNECTING, STATE_CLOSING, STATE_CLOSED
using DandelionWebSockets: SocketState, AbstractPonger, SendTextFrame, FrameFromServer
using DandelionWebSockets: handle
import DandelionWebSockets: write
import DandelionWebSockets: on_text, on_binary
import DandelionWebSockets: state_connecting, state_open, state_closing, state_closed
using Test

include("stubs.jl")
include("test_factories.jl")
include("issues_test.jl")
include("client_protocol_from_server_test.jl")
include("client_protocol_from_client_test.jl")
include("client_protocol_connection_test.jl")
include("masking_test.jl")
include("core_test.jl")
include("handshake_short_test.jl")
include("handshake_logic_test.jl")
include("security_considerations_test.jl")
include("data_frames_test.jl")
include("frame_sender_test.jl")
include("frame_writer_test.jl")
include("closing_behaviour_test.jl")
include("network_small_test.jl")
include("WebSocketUpgrade/WebSocketUpgrade_test.jl")
include("UniformResourceIdentifiers/UniformResourceIdentifiers_test.jl")


text:
include("./imports.jl")

@testset ExtendedTestSet "FluxTraining.jl" begin
    include("./testfunctional.jl")
    include("./testoptim.jl")
    include("./testprotected.jl")
    include("./testtraining.jl")
    include("./testcallbacks.jl")
end


text:
@testset "Core" begin
    e2 = Edge(1, 3)
    e3 = Edge(1, 4)
    @test @inferred(is_ordered(e2))
    @test @inferred(!is_ordered(reverse(e3)))

    gx = Graph(10)
    for g in testgraphs(gx)
        add_vertices!(g, 5)
        @test @inferred(nv(g)) == 15
    end

    g5w = WheelGraph(5)
    for g in testgraphs(g5w)
        @test @inferred(indegree(g, 1)) == @inferred(outdegree(g, 1)) == @inferred(degree(g, 1)) == 4
        @test @inferred(indegree(g)) == @inferred(outdegree(g)) == @inferred(degree(g)) == [4, 3, 3, 3, 3]

        @test @inferred(Δout(g)) == @inferred(Δin(g)) == @inferred(Δ(g)) == 4
        @test @inferred(δout(g)) == @inferred(δin(g)) == @inferred(δ(g)) == 3
        @test @inferred(weights(g)) == LightGraphs.DefaultDistance(nv(g))
        z1 = @inferred(degree_histogram(g))
        z2 = @inferred(degree_histogram(g, indegree))
        z3 = @inferred(degree_histogram(g, outdegree))

        @test z1 == z2 == z3 == Dict(4 => 1, 3 => 4)

        @test @inferred(neighbors(g, 2)) == @inferred(all_neighbors(g, 2)) == [1, 3, 5]
        @test @inferred(common_neighbors(g, 1, 5)) == [2, 4]

        gsl = copy(g)
        add_edge!(gsl, 3, 3)
        add_edge!(gsl, 2, 2)

        @test @inferred(!has_self_loops(g))
        @test @inferred(has_self_loops(gsl))
        @test @inferred(num_self_loops(g)) == 0
        @test @inferred(num_self_loops(gsl)) == 2

        @test @inferred(density(g)) == 0.8

        @test eltype(squash(g)) == UInt8
    end

    g5wd = WheelDiGraph(5)
    for g in testdigraphs(g5wd)
        @test @inferred(indegree(g, 2)) == 2
        @test @inferred(outdegree(g, 2)) == 1
        @test @inferred(degree(g, 2)) == 3
        @test @inferred(indegree(g)) == [0, 2, 2, 2, 2]
        @test @inferred(outdegree(g)) == [4, 1, 1, 1, 1]
        @test @inferred(degree(g)) == [4, 3, 3, 3, 3]
        @test @inferred(weights(g)) == LightGraphs.DefaultDistance(nv(g))

        @test @inferred(Δout(g)) == @inferred(Δ(g)) == 4
        @test @inferred(Δin(g)) == 2
        @test @inferred(δout(g)) == 1
        @test @inferred(δin(g)) == 0
        @test @inferred(δ(g)) == 3

        z1 = @inferred(degree_histogram(g))
        z2 = @inferred(degree_histogram(g, indegree))
        z3 = @inferred(degree_histogram(g, outdegree))
        @test z1 == Dict(4 => 1, 3 => 4)
        @test z2 == Dict(0 => 1, 2 => 4)
        @test z3 == Dict(4 => 1, 1 => 4)

        @test @inferred(neighbors(g, 2)) == [3]
        @test Set(@inferred(all_neighbors(g, 2))) == Set([1, 3, 5])
        @test @inferred(common_neighbors(g, 1, 5)) == [2]

        gsl = copy(g)
        add_edge!(gsl, 3, 3)
        add_edge!(gsl, 2, 2)

        @test @inferred(!has_self_loops(g))
        @test @inferred(has_self_loops(gsl))
        @test @inferred(num_self_loops(g)) == 0
        @test @inferred(num_self_loops(gsl)) == 2

        @test @inferred(density(g)) == 0.4
        @test eltype(squash(g)) == UInt8
        end
    end


text:
export plotsDir
export dataDir
export Beamsplitter
export classicalFreqData
export classicalSinglePlot
export classicalSumPlot

"""
    plotsDir(name::String,dirname::String = "plots")

Concatenates string with name of plots directory
"""
function plotsDir(name::String,dirname::String = "plots")
    return joinpath(dirname,name)
end

"""
    dataDir(name::String,dirname::String = "data")

Concatenates string with name of data directory
"""
function dataDir(name::String,dirname::String = "data")
    return joinpath(dirname,name)
end

"""
    γIntensityPlot(times::Array,γint::Array,prefix::String)

Plots average photon counts and saves to file. Returns relative path to resulting plot.
"""
function γIntensityPlot(times::Array,γint::Array,prefix::String)
    γplotName = string(prefix,"avg-photons-vs-time.svg")
    γplot = plot(times,γint,label=false)
    xlabel!("time (ns)")
    ylabel!("Avg photon counts")
    title!("Average Photon Counts vs Time")
    savefig(γplot,γplotName)
    @info "Saved $γplotName"
    return γplotName
end

export γIntensityPlot

"""
    γCountPlot(times::Array,γint::Array,prefix::String)

Plots 'measured' photon counts and saves to file. Returns relative path to resulting plot.
"""
function γCountPlot(times::Array,γcount::Array, γint::Array,prefix::String)
    γplotName = string(prefix,"photon-counts-vs-time.svg")
    γplot = bar(times,γcount,label="Simulated Counts")
    sig1 = sqrt.(γint)
    sig1low = map(x -> x < 1 ? x : sqrt(x),γint)
    plot!(γplot,times,γint,label = "Average Counts"; ribbon = (sig1low,sig1))
    xlabel!("time (ns)")
    ylabel!("photon counts")
    title!("Photon Counts vs Time")
    savefig(γplot,γplotName)
    @info "Saved $γplotName"
    return γplotName
end

export γCountPlot
"""
    γCorrTimePlot(timeVec::Vector, g2Vec::Vector, params::Dict, prefix::String)

Saves plot of photon correlation time series
"""
function γCorrTimePlot(timeDF::DataFrame, params::Dict, prefix::String)
    out = []
    γplotName = string(prefix,"time-domain-photon-correlation.svg")
    inzcorr  = timeDF[!,:corr1] .> 0
    nzcorr = timeDF[!,:corr1][inzcorr]
    nztime = timeDF[!,:time][inzcorr]
	γplot = plot(nztime,nzcorr,label = false)
    xlabel!(L"\tau \textrm{ (ns)}")
	ylabel!("\$g^{(2)}(\\tau)\$")
    title!("Photon Correlations vs Time Offset")
    savefig(γplot,γplotName)
    push!(out,γplotName)
    @info "Saved $(γplotName)"
    if params[:repeat] > 1
        γplotName = string(prefix,"time-domain-photon-correlation-sum.svg")
        nzcorr = timeDF[!,:sum][inzcorr]
        γplot = plot(nztime,nzcorr,label = false)
        xlabel!(L"\tau \textrm{ (ns)}")
        ylabel!("\$\\sum_{i=1}^{$(params[:repeat])}g_i^{(2)}(\\tau)\$")
        title!("Sum of Photon Correlations vs Time Offset")
        savefig(γplot,γplotName)
        push!(out,γplotName)
        @info "Saved $(γplotName)"
    end
    return out
end

export γCorrTimePlot

"""
    γCorrFreqPlot(freqVec::Vector, fftVec::Vector, params::Dict, prefix::String)

Saves plot of photon correlation Fourier transform
"""
function γCorrFreqPlot(freqDF::DataFrame, params::Dict, prefix::String)
    out = []
    γplotName = string(prefix,"frequency-domain-photon-correlation.svg")
	γplot = plot(freqDF[!,:freq],freqDF[!,:corr1],label = false)
	xlabel!("frequency (GHz)")
	ylabel!("\$\\hat{g}^{(2)}(\\nu)\$")
	title!("Photon Correlations vs Frequency")
    if length(params[:νm]) > 1
        shifts = map(x->abs(x[2]-x[1]),subsets(params[:νm],2))
        shiftname = length(params[:νm]) > 2 ? "frequency shifts" : "frequency shift"
        vline!(shifts,label=shiftname,ls=:dash)
    end
    savefig(γplot,γplotName)
    push!(out,γplotName)
    @info "Saved $(γplotName)"
    if params[:repeat] > 1
        γplotName = string(prefix,"frequency-domain-photon-correlation-sum.svg")
        γplot = plot(freqDF[!,:freq],freqDF[!,:sum],label = false)
        xlabel!("frequency (GHz)")
        ylabel!("\$\\sum_{i=1}^{$(params[:repeat])}\\hat{g}^{(2)}(\\nu)\$")
        title!("Sum of Photon Correlations vs Frequency")
        if length(params[:νm]) > 1
            shifts = map(x->abs(x[2]-x[1]),subsets(params[:νm],2))
            shiftname = length(params[:νm]) > 2 ? "frequency shifts" : "frequency shift"
            vline!(shifts,label=shiftname,ls=:dash)
        end
        savefig(γplot,γplotName)
        push!(out,γplotName)
        @info "Saved $(γplotName)"
    end
    return out
end

export γCorrFreqPlot


text:
using Test
using CubicSplines
using Random
using PyPlot # Plot the curve

@testset "cubic_splines.jl" begin

    # Create data that is symmetric about zero
    xdata = -2.0:1.0:2.0
    ydata = mod.(xdata, 2)
    @assert reverse(xdata) == -xdata
    @assert reverse(ydata) == +ydata

    # Test spline is symmetric about zero
    xs = -2.0:0.5:2.0
    ys = CubicSpline(xdata, ydata)[xs]
    @assert reverse(xs) == -xs
    @test reverse(ys) == +ys

    # Interpolate a sine function with white noise
    xdata = range(0,stop=4π,length=20) .+ 0.5*rand(MersenneTwister(1234),20)
    ydata = sin.(xdata)

    # Create the spline
    spline = CubicSpline(xdata, ydata)
    @test spline[π] ≈ 0.00652422 rtol = 1e-5
    @test spline[0.35] ≈ 0.346190 rtol = 1e-5

    # Non-allocating vector filling
    x = [1.0, 2.0]
    y = zeros(2)
    spline(y,x)
    for i in 1:length(y)
        @test y[i] == spline(x[i])
    end

    # Plot the spline
    # xs = range(xdata[1], stop=xdata[end], length=1000)
    # ys = spline[xs]
    # plt.plot(xdata, ydata, "o")
    # plt.plot(xs, ys)
    # plt.xlabel("x")
    # plt.ylabel("y")
    # plt.close()

    # Extrapolate the spline with a linear curve to the left
    spline = CubicSpline(xdata, ydata; extrapl=[1.0]) # Linear extrapolation with an incline of 1.0

    # Data from outside the interpolation range
    left_value = -1.0
    right_value = 14.0
    @test spline[left_value] ≈ -1.004278 rtol = 1e-5 # Extrapolate linearly to the left
    @test_throws String spline[right_value] # No extrapolation (throw an error)

    # Plot the spline
    # xs = range(left_value, stop=xdata[end], length=1000)
    # ys = spline[xs]
    # plt.plot(xdata, ydata, "o")
    # plt.plot(xs, ys)
    # plt.xlabel("x")
    # plt.ylabel("y")
    # plt.title("With linear extrapolation to the left")
    # plt.close()

    # Extrapolate the spline with a quadratic curve to the left and a cubic curve to the right
    spline = CubicSpline(xdata, ydata; extrapl=[0.0; 0.5], extrapr=[1.0; 0.0; 0.5])
    left_value = -1.0
    right_value = 14.0
    @test spline[left_value] ≈ 1.13020 rtol = 1e-5 # Extrapolate linearly to the left
    @test spline[right_value] ≈ 1.93991 rtol = 1e-5 # Extrapolate linearly to the left

    # Plot the spline
    # xs = range(left_value, stop=right_value, length=2000)
    # ys = spline[xs]
    # plt.plot(xdata, ydata, "o")
    # plt.plot(xs, ys)
    # plt.xlabel("x")
    # plt.ylabel("y")
    # plt.title("With extrapolation on both sides")
    # plt.close()

    # Gradient calculation
    xdata = range(0,stop=4pi,length=20)
    ydata = sin.(xdata)
    spline = CubicSpline(xdata, ydata; extrapl=[2.0; 1.0], extrapr=[1.0; 0.0; 0.5])

    # Plot the spline
    # xs = range(left_value, stop=right_value, length=2000)
    # ys = spline[xs]
    # plt.plot(xdata, ydata, "o")
    # plt.plot(xs, ys)
    # plt.xlabel("x")
    # plt.ylabel("y")
    # plt.title("With extrapolation on both sides")
    # plt.close()

    # Calculate the spline gradient at different points and evaluate it
    @test gradient(spline, -2.0) ≈ -2.0
    @test gradient(spline, -2.0, 2) ≈ 2.0
    @test gradient(spline, -2.0, 3) ≈ 0.0
    @test gradient(spline, pi/2, 1) ≈ 0.01851987 rtol = 1e-5
end


text:
struct Stardate
  sd::Float64
  canonical::String
  short::String
  originaltz::TimeZones.TimeZone
end

"""
    defaulttimezone()

Returns a TimeZone object corresponding to the
IPFS file `/t/env/TZ`.
"""
function defaulttimezone()
    tzstring = read(`ipfs files read /t/env/TZ`, String)
    TimeZones.TimeZone(tzstring)
end

function Stardate(x::Float64, tz1::TimeZone = tz"UTC")
  Stardate(
    x,
    Printf.@sprintf("%0.15f", x),
    Printf.@sprintf("%0.3f", x),
    tz1
  )
end

function Stardate(zdt::ZonedDateTime)
  tx = astimezone(zdt, tz"UTC")
  y0 = year(tx)
  t0 = Stardates.getstartofyear(y0).zoneddatetime
  t1 = Stardates.getstartofyear(y0 + 1).zoneddatetime
  Stardate(y0 + (tx - t0) / (t1 - t0), timezone(zdt))
end

function Stardate(dt::Dates.DateTime)
  # Assume dt is in UTC
  y0 = Dates.year(dt)
  t0 = Stardates.getstartofyear(y0).unix
  t1 = Stardates.getstartofyear(y0 + 1).unix
  Stardate(y0 + (Dates.datetime2unix(dt) - t0) / (t1 - t0))
end

function Stardate(d::Dates.Date, hr, mi, ss, tz)
  Stardate(ZonedDateTime(
    Dates.year(d),
    Dates.month(d),
    Dates.day(d),
    hr,
    mi,
    ss,
    tz,
  ))
end

function Stardate(;
    year,
    month,
    day,
    hour = 12,
    minute = 0,
    second = 0,
    tz = nothing,
    style = nothing,
    copy = false
)
  if tz == nothing
    tz1 = defaulttimezone()
  else
    tz1 = TimeZone(tz)
  end
  Stardate(ZonedDateTime(year, month, day, hour, minute, second, tz1))
end

function stardate(args...)
  Stardate(args...).sd
end

"Return the stardate of the closing bell of the date"
function nyse_closing_stardate(d::Dates.Date)
  stardate(d, 16, 0, 0, tz"America/New_York")
end

"Return the stardate of the closing bell of the date represented by iso8601"
function nyse_closing_stardate(iso8601)
  nyse_closing_stardate(parse(Dates.Date, iso8601))
end

function Stardate(st::Base.StatStruct)
  Stardate(Dates.unix2datetime(st.mtime))
end

function mstardate(fn::AbstractString)
  Stardate(stat(fn))
end

function Stardate(style::Symbol)
  @assert(style == :now)
  Stardate(now(defaulttimezone()))
end

function mediawiki(x::Stardate)
    "<stardate digits=\"3\" tz=\"" *
    TimeZones.name(x.originaltz) * "\">" *
    x.canonical *
    "</stardate>"
end


text:
import Base.Iterators:partition

function load_data(filename)

    open(filename) do f
        bags = Dict{AbstractString,Vector{AbstractString}}()
        for line in readlines(f)
            line = strip(replace(line, [',','.'] => ""))
            line = strip(replace(line, "bags" => ""))
            line = strip(replace(line, "bag" => ""))
            line = strip(replace(line, "no other" => ""))
            color, contain = split(strip(line),"contain")
            loads = String[]
            for x in partition(split(contain), 3)
                n , c1, c2 = x
                for i in 1:parse(Int, n)
                    push!(loads,string(c1," ",c2))
                end
            end
            bags[strip(color)] = loads
        end
        return bags
    end

end

function hold_shiny_gold( color, bags )

    if "shiny gold" in bags[color]
       return true
    elseif length(bags[color]) > 0
       return any([hold_shiny_gold( c, bags ) for c in unique(bags[color])])
    else
	   return false
    end

end

function contained_bags( color, bags, n)

    if length(bags[color]) == 0
        return n
    else
        s = length(bags[color])
        for c in bags[color]
           s += contained_bags(c, bags, n)
        end
        return s
    end

end

bags = load_data("input07.txt")
println(count([hold_shiny_gold(c, bags) for c in keys(bags)]))
println(contained_bags( "shiny gold", bags, 0))

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
using UnicodePlots: heatmap

function read_dots_and_folds(file)
    dots = Tuple{Int, Int}[]
    folds = Tuple{String, Int}[]
    for line in eachline(file)
        isempty(line) && continue
        if startswith(line, "fold")
            axis, at = split(replace(line, "fold along " => ""), "=")
            push!(folds, (axis, parse(Int, at)))
        else
            push!(dots, Tuple(parse.(Int, split(line, ","))))
        end
    end
    dots, folds
end

fold_x!(dots, x) = unique!(map!(dot -> (mirror(dot[1], x), dot[2]), dots, dots))

fold_y!(dots, y) = unique!(map!(dot -> (dot[1], mirror(dot[2], y)), dots, dots))

mirror(a, b) = a < b ? a : 2b - a

function fold!(dots, axis, at)
    axis == "x" && fold_x!(dots, at)
    axis == "y" && fold_y!(dots, at)
end

function plot(dots)
    x_max = maximum(x for x in first.(dots))
    y_max = maximum(y for y in last.(dots))
    paper = falses(y_max + 1, x_max + 1)
    for (x, y) in dots
        paper[y_max - y + 1, x + 1] = true
    end
    heatmap(paper)
end

dots, folds = read_dots_and_folds("day13/dots.txt")

first_fold = popfirst!(folds)
fold!(dots, first_fold...)
p1 = length(dots)

@show p1

for fold in folds
    fold!(dots, fold...)
end

plot(dots)
p2 = "PFKLKCFP"

@show p2


text:
using Base.Test
using Bebop

const sampRate = 44100
const duration = 5.0

initialize() do
    openStream(sampRate, 2, Int16) do stream
        record = readStream(stream, duration)
        sleep(1)
        writeStream(stream, record)
    end
end


text:
# function df_to_matrix(df::DataFrame) #Converts a data frame with missing values to a matrix with NaN
#     taxa = Vector{String}(df[!, :taxon])
#     n, p = size(df)
#     p = p - 1
#     nms = names(df)

#     data = fill(NaN, n, p)
#     for i = 1:p
#         for j = 1:n
#             x = df[j, nms[i + 1]]
#             if !ismissing(x)
#                 data[j, i] = x
#             end
#         end
#     end
#     return taxa, data
# end

function use_dates!(bx::BEASTXMLElement)
    data = get_data(bx)
    newick = get_newick(bx)
    data.use_times = true
    newick.fix_tree = false
end

function add_trait!(bx::BEASTXMLElement, data::Matrix{Float64}, trait::String)
    p = size(data, 2)
    data_el = get_data(bx)
    tm_el = get_treeModel(bx)

    add_trait!(data_el, data, trait)
    add_leaf_param!(tm_el, trait, p)
end

function get_id(el::XMLElement)
    return attribute(el, bn.ID)
end

function set_id!(el::XMLElement, id::String)
    set_attribute(el, bn.ID, id)
end

function add_matrix_parameter(pel::XMLElement, M::AbstractArray{T, 2};
        id::String = "") where T <: Number
    mat_el = new_element(bn.MATRIX_PARAMETER)
    if id != ""
        set_attribute(mat_el, bn.ID, id)
    end
    n, p = size(M)
    for i = 1:n
        param_el = new_element(bn.PARAMETER)
        set_attribute(param_el, bn.VALUE, join(M[i, :], ' '))
        add_child(mat_el, param_el)
    end
    add_child(pel, mat_el)
    return mat_el
end

function add_diagonal_matrix(pel::XMLElement, M::Vector{T};
        id::String = "", lower::String = "") where T <: Number

    mat_el = new_child(pel, bn.DIAGONAL_MATRIX)
    p_el = add_parameter(mat_el, value = M, id=id)
    if lower != ""
        set_attribute(p_el, bn.LOWER, lower)
    end
    return mat_el
end

# function add_parameter(pel::XMLElement, val::Vector{<:Number};
#                        lower::Vector{<:Number} = Float64[])
#     el = new_element(bn.PARAMETER)
#     set_attribute(el, bn.VALUE, join(val, ' '))
#     if length(lower) > 0
#         if length(lower) == 1 || length(lower) == length(val)
#             set_attribute(el, bn.LOWER, join(lower, ' '))
#         else
#             error("incompatible dimension between the parameter values " *
#                 "($(length(val))) and lower bounds ($(length(lower)))")
#         end
#     end

#     add_child(pel, el)
# end

function add_ref_el(pel::XMLElement, el::XMLElement;
            new_name::String = name(el))
    ref_el = reference_element(el, new_name)
    add_child(pel, ref_el)
    return ref_el
end

function add_ref_el(pel::XMLElement, name::String, id::String)
    ref_el = new_child(pel, name)
    set_attribute(ref_el, bn.IDREF, id)
end

function add_ref_els(pel::XMLElement, els::Array{XMLElement})
    for el in els
        add_ref_el(pel, el)
    end
end

function reference_element(el::XMLElement)
    nm = name(el)
    return reference_element(el, nm)
end

function reference_element(el::XMLElement, nm::String)
    id = attribute(el, bn.ID)
    if isnothing(id)
        id = attribute(el, bn.IDREF)
    end

    if isnothing(id)
        @warn "Element does not have id. Adding whole element."
        return el
    else
        ref_el = new_element(nm)
        set_attribute(ref_el, bn.IDREF, id)
        return ref_el
    end
end

function xml_vec(n::Int)
    v = Vector{XMLOrNothing}(undef, n)
    fill!(v, nothing)
    return v
end

function check_posdef(x::AbstractArray{Float64, 2})
    if !isposdef(x)
        error("Matrix must be positive definite.")
    end
end


text:
export states, batches

struct StatesIterator{Safe,T,S}
    states::T
    sz::S
end

"""
    states(v::StateCollection)

Returns an iterator to iterate all states in the batch/array of batches
`v`.

By default, the iteration is done using unsafe views for increased performance.
"""
function states(v::Union{AbstractArray,AbstractDoubled}, safe=Val{false})
    sz = state_size(v)
    sz = CartesianIndices(sz)

    do_safe = safe isa Val{true} ? true : false
    return StatesIterator{do_safe, typeof(v), typeof(sz)}(v,sz)
end

Base.length(iter::StatesIterator) = length(iter.sz)
Base.eltype(iter::StatesIterator) = state_eltype(iter.states)
Base.lastindex(iter::StatesIterator) = length(iter)

function Base.iterate(it::StatesIterator{false}, state=1)
    state > length(it.sz) && return nothing
    @inbounds v = state_uview(it.states, it.sz[state].I...)

    return (v, state+1)
end

function Base.getindex(it::StatesIterator{false}, i::Integer)
    @boundscheck checkbounds(it.sz, i)
    @inbounds v = state_uview(it.states, it.sz[i].I...)

    return v
end

####
struct BatchesIterator{Safe,T,S}
    states::T
    sz::S
    b_sz::Int  # size of batches
    b_num::Int # number of batches per entry
end

function batches(v::Union{AbstractArray,AbstractDoubled}, batch_size=batch_size(v), safe=Val{false})
    sz = state_size(v)
    batch_num, remainder = divrem(sz[1],batch_size)

    @assert remainder == 0

    sz = CartesianIndices(sz[2:end])
    do_safe = safe isa Val{true} ? true : false
    return BatchesIterator{do_safe, typeof(v), typeof(sz)}(v, sz, batch_size, batch_num)
end

Base.length(iter::BatchesIterator) = length(iter.sz)*iter.b_num
Base.eltype(iter::BatchesIterator) = state_eltype(iter.states)
Base.lastindex(iter::BatchesIterator) = length(iter)

function Base.iterate(it::BatchesIterator{false}, state=1)
    state > length(it.sz) && return nothing
    @inbounds v = state_uview(it.states, it.sz[state].I...)

    return (v, state+1)
end

function Base.getindex(it::BatchesIterator{false}, i::Integer)
    @boundscheck checkbounds(it.sz, i)
    @inbounds v = state_uview(it.states, it.sz[i].I...)

    return v
end


text:
import Base.Broadcast.broadcasted
export QolElemAtom, qol_elementwise, square, sumsquares, invpos, /
export sign, monotonicity, curvature, conic_form!

struct QolElemAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr}
    size::Tuple{Int, Int}

    function QolElemAtom(x::AbstractExpr, y::AbstractExpr)
        if x.size != y.size
            error("elementwise quad over lin must take two arguments of the same size")
        end
        children = (x, y)
        return new(:qol_elem, hash(children), children, x.size)
    end
end

function sign(q::QolElemAtom)
    return Positive()
end

function monotonicity(q::QolElemAtom)
    return (sign(q.children[1]) * Nondecreasing(), Nonincreasing())
end

function curvature(q::QolElemAtom)
    return ConvexVexity()
end

function evaluate(q::QolElemAtom)
    return (evaluate(q.children[1]).^2) ./ evaluate(q.children[2])
end

function conic_form!(q::QolElemAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, q)
        sz = q.children[1].size
        t = Variable(sz[1], sz[2])
        qol_objective = conic_form!(t, unique_conic_forms)
        x, y = q.children
        conic_form!(SOCElemConstraint(y + t, y - t, 2 * x), unique_conic_forms)
        # add implicit constraint y >= 0
        conic_form!(y >= 0, unique_conic_forms)
        cache_conic_form!(unique_conic_forms, q, qol_objective)
    end
    return get_conic_form(unique_conic_forms, q)
end

qol_elementwise(x::AbstractExpr, y::AbstractExpr) = QolElemAtom(x, y)

broadcasted(::typeof(^),x::AbstractExpr,k::Int) = k==2 ? QolElemAtom(x, Constant(ones(x.size[1], x.size[2]))) : error("raising variables to powers other than 2 is not implemented")

invpos(x::AbstractExpr) = QolElemAtom(Constant(ones(x.size[1], x.size[2])), x)
broadcasted(::typeof(/), x::Value, y::AbstractExpr) = DotMultiplyAtom(Constant(x), invpos(y))
/(x::Value, y::AbstractExpr) = size(y) == (1,1) ? MultiplyAtom(Constant(x), invpos(y)) : error("cannot divide by a variable of size $(size(y))")
sumsquares(x::AbstractExpr) = square(norm2(x))

function square(x::AbstractExpr)
    if sign(x) == ComplexSign()
        error("Square of complex number is not DCP. Did you mean square_modulus?")
    else
        QolElemAtom(x, Constant(ones(x.size[1], x.size[2])))
    end
end


text:
pids = addprocs(10)

@everywhere begin
    include("voiConsAuc.jl")
    f(n, p1, p2) = evsi_sim([1, .5, .5], [.8, 1.25, .75], [n * p1, n * p2, n * (1 - (p1 + p2))], nsims = 1000000)
    function g(n)
        s = 100
        x  = zeros(s, s)
        p1 = linspace(1, 0, s)
        p2 = linspace(0, 1, s)
        for i = 1:s
            for j = 1:i
                x[i, j] = f(n, p1[i], p2[j])
            end
        end
        y = ind2sub(size(x), indmax(x))
        [p1[y[1]], p2[y[2]]]
    end
end

y = @parallel hcat for n = linspace(0.5, 10, 50)
  g(n)
end
rmprocs(pids)

writedlm("evsi3sim.csv", y, ',')


text:
function descr(T, d)
    ex = quote
        MMI.docstring(::Type{<:$T}) = MMI.docstring_ext($T; descr=$d)
    end
    eval(ex)
    return
end

function meta(T; input=Unknown, target=Unknown, output=Unknown,
              weights::Bool=false, descr::String="")
    ex = quote
        MMI.input_scitype(::Type{<:$T})    = $input
        MMI.output_scitype(::Type{<:$T})   = $output
        MMI.target_scitype(::Type{<:$T})   = $target
        MMI.supports_weights(::Type{<:$T}) = $weights
        MMI.docstring(::Type{<:$T}) = MMI.docstring_ext($T; descr=$descr)
    end
    eval(ex)
    return
end