text:
"""
    plot(net::HybridNetwork; useEdgeLength=false, mainTree=false, showTipLabel=true,
         showNodeNumber=false, showEdgeLength=false, showGamma=false, edgeColor=colorant"black",
         majorHybridEdgeColor=colorant"deepskyblue4", minorHybridEdgeColor=colorant"deepskyblue",
         showEdgeNumber=false, showIntNodeLabel=false, edgeLabel=[], nodeLabel=[])

Plots a network, from left to right.

- useEdgeLength: if true, the tree edges and major hybrid edges are
  drawn proportionally to their length. Minor hybrid edges are not, however.
  Note that edge lengths in coalescent units may scale very poorly with time.
- mainTree: if true, the minor hybrid edges are ommitted.
- showTipLabel: if true, taxon labels are shown. You may need to zoom out to see them.
- showNodeNumber: if true, nodes are labelled with the number used internally.
- showEdgeLength: if true, edges are labelled with their length (above)
- showGamma: if true, hybrid edges are labelled with their heritability (below)
- edgeColor: color for tree edges. black by default.
- majorHybridEdgeColor: color for major hybrid edges
- minorHybridEdgeColor: color for minor hybrid edges
- showEdgeNumber: if true, edges are labelled with the number used internally.
- showIntNodeLabel: if true, internal nodes are labelled with their names.
  Useful for hybrid nodes, which do have tags like 'H1'.
- edgeLabel: dataframe with two columns: the first with edge numbers, the second with labels
  (like bootstrap values) to annotate edges. empty by default.
- nodeLabel: dataframe with two columns: the first with node numbers, the second with labels
  (like bootstrap values for hybrid relationships) to annotate nodes. empty by default.

Note that `plot` actually modifies some (minor) attributes of the network,
as it calls `directEdges!`, `preorder!` and `cladewiseorder!`.

If hybrid edges cross tree and major edges, you may choose to rotate some tree
edges to eliminate crossing edges, using `rotate!`.
"""
function Gadfly.plot(net::HybridNetwork; useEdgeLength=false::Bool,
        mainTree=false::Bool, showTipLabel=true::Bool, showNodeNumber=false::Bool,
        showEdgeLength=false::Bool, showGamma=false::Bool,
        edgeColor=colorant"black"::ColorTypes.Colorant,
        majorHybridEdgeColor=colorant"deepskyblue4"::ColorTypes.Colorant,
        minorHybridEdgeColor=colorant"deepskyblue"::ColorTypes.Colorant,
        showEdgeNumber=false::Bool, showIntNodeLabel=false::Bool,
        edgeLabel=DataFrame()::DataFrame, nodeLabel=DataFrame()::DataFrame)

    (edge_xB, edge_xE, edge_yB, edge_yE, node_x, node_y, node_yB, node_yE,
     xmin, xmax, ymin, ymax) = getEdgeNodeCoordinates(net, useEdgeLength)

    !net.node[net.root].leaf ||
        @warn "the root is leaf $(net.node[net.root].name): the plot will look weird..."

    mylayers = Layer[] # gadfly layers
    # one layers for each edge
    for i=1:net.numEdges
        if (!mainTree || net.edge[i].isMajor)
            col = edgeColor
            if net.edge[i].hybrid
              if (net.edge[i].isMajor) col = majorHybridEdgeColor;
              else col = minorHybridEdgeColor; end
            end
            push!(mylayers,
              layer(x = [edge_xB[i],edge_xE[i]],
                    y = [edge_yB[i],edge_yE[i]], Geom.line,
                    Theme(default_color=col))[1])
        end
    end
    # one layer for each (vertical) clade
    for i=1:net.numNodes
        if (net.node[i].leaf) continue; end
        push!(mylayers,
              layer(y = [node_yB[i],node_yE[i]],
                    x = [node_x[i], node_x[i]], Geom.line)[1])
    end
    labelnodes, nodeLabel = checkNodeDataFrame(net, nodeLabel)
    if (showTipLabel || showNodeNumber || showIntNodeLabel || labelnodes)
      # white dot beyond tip labels/name and root node label to force enough zoom out
      expfac = 0.1
      expfacy = 0.5 # additive expansion for y axis
      push!(mylayers, layer(x=[xmin-(xmax-xmin)*expfac,xmax+(xmax-xmin)*expfac],
                            y=[ymin- expfacy, ymax+ expfacy],
               Geom.point, Theme(default_color=colorant"white"))[1])
      # data frame to place tip names and node annotations (labels)
      ndf = prepareNodeDataFrame(net, nodeLabel, showNodeNumber,
              showIntNodeLabel, labelnodes, node_x, node_y)
      if (showTipLabel)
        push!(mylayers, layer(ndf[ndf[!,:lea], [:x,:y,:name]], y="y", x="x", label="name",
            Geom.label(position=:right ;hide_overlaps=true))[1])
      end
      if (showIntNodeLabel)
        push!(mylayers, layer(ndf[.!ndf[!,:lea], [:x,:y,:name]], y="y", x="x", label="name",
            Geom.label(position=:above ;hide_overlaps=true))[1])
      end
      if (showNodeNumber)
        push!(mylayers, layer(ndf, y="y", x="x", label="num",
            Geom.label(position=:dynamic ;hide_overlaps=true))[1])
      end
      if labelnodes
        push!(mylayers, layer(ndf[!,[:x,:y,:lab]], y="y", x="x", label="lab",
            Geom.label(position=:left ;hide_overlaps=false))[1])
      end
    end
    # data frame for edge annotations.
    labeledges, edf = prepareEdgeDataFrame(net, edgeLabel, mainTree,
                        edge_xB, edge_xE, edge_yB, edge_yE)
    if labeledges
      push!(mylayers, layer(edf[!,[:x,:y,:lab]], y="y", x="x", label="lab",
            Geom.label(position=:above ;hide_overlaps=false))[1])
    end
    if (showEdgeLength)
      push!(mylayers, layer(edf[!,[:x,:y,:len]], y="y", x="x", label="len",
            Geom.label(position=:below ;hide_overlaps=false))[1])
    end
    if (showGamma && net.numHybrids>0)
      if !mainTree
        push!(mylayers, layer(edf[edf[!,:hyb] .& edf[!,:min], [:x,:y,:gam]], y="y", x="x",label="gam",
              Geom.label(position=:below ;hide_overlaps=true),
              Theme(point_label_color=minorHybridEdgeColor))[1])
      end
      push!(mylayers, layer(edf[edf[!,:hyb] .& .!edf[!,:min], [:x,:y,:gam]], y="y", x="x",label="gam",
            Geom.label(position=:below ;hide_overlaps=true),
            Theme(point_label_color=majorHybridEdgeColor))[1])
    end
    if (showEdgeNumber)
      push!(mylayers, layer(edf[!,[:x,:y,:num]], y="y", x="x", label="num",
            Geom.label(position=:dynamic ;hide_overlaps=false))[1])
    end

    Gadfly.plot(mylayers, Guide.ylabel(nothing), Guide.xlabel(nothing), #("time"),
         Guide.xticks(ticks=:auto, label=false), # ticks=[xmin,xmin,xmax,xmax*1.1],
         Guide.yticks(ticks=:auto, label=false), # ticks=[ymin,ymax],
         Theme(default_color=edgeColor,grid_color=colorant"white",grid_line_width=0pt))
end


text:
export differentiate

"""
    differentiate(p::AbstractPolynomialLike, v::AbstractVariable, deg::Union{Int, Val}=1)

Differentiate `deg` times the polynomial `p` by the variable `v`.

    differentiate(p::AbstractPolynomialLike, vs, deg::Union{Int, Val}=1)

Differentiate `deg` times the polynomial `p` by the variables of the vector or
tuple of variable `vs` and return an array of dimension `deg`. It is recommended
to pass `deg` as a `Val` instance when the degree is known at compile time, e.g.
`differentiate(p, v, Val{2}())` instead of `differentiate(p, x, 2)`, as this
will help the compiler infer the return type.

    differentiate(p::AbstractArray{<:AbstractPolynomialLike, N}, vs, deg::Union{Int, Val}=1) where N

Differentiate the polynomials in `p` by the variables of the vector or tuple of variable `vs` and return an array of dimension `N+deg`.
If `p` is an `AbstractVector` this returns the Jacobian of `p` where the i-th row containts the partial
derivaties of `p[i]`.

### Examples

```julia
p = 3x^2*y + x + 2y + 1
differentiate(p, x) # should return 6xy + 1
differentiate(p, x, Val{1}()) # equivalent to the above
differentiate(p, (x, y)) # should return [6xy+1, 3x^2+1]
differentiate( [x^2+y, z^2+4x], [x, y, z]) # should return [2x 1 0; 4 0 2z]
```
"""
function differentiate end

# Fallback for everything else
differentiate(α::T, v::AbstractVariable) where T = zero(T)
differentiate(v1::AbstractVariable, v2::AbstractVariable) = v1 == v2 ? 1 : 0
differentiate(t::AbstractTermLike, v::AbstractVariable) = coefficient(t) * differentiate(monomial(t), v)
# The polynomial function will take care of removing the zeros
differentiate(p::APL, v::AbstractVariable) = polynomial(differentiate.(terms(p), v), SortedState())
differentiate(p::RationalPoly, v::AbstractVariable) = (differentiate(p.num, v) * p.den - p.num * differentiate(p.den, v)) / p.den^2

const ARPL = Union{APL, RationalPoly}

function differentiate(ps::AbstractArray{PT}, xs::AbstractArray) where {PT <: ARPL}
    differentiate.(reshape(ps, (size(ps)..., 1)), reshape(xs, 1, :))
end

function differentiate(ps::AbstractArray{PT}, xs::Tuple) where {PT <: ARPL}
    differentiate(ps, collect(xs))
end

# TODO: this signature is probably too wide and creates the potential
# for stack overflows
differentiate(p::ARPL, xs) = [differentiate(p, x) for x in xs]

# differentiate(p, [x, y]) with TypedPolynomials promote x to a Monomial
differentiate(p::ARPL, m::AbstractMonomial) = differentiate(p, variable(m))

# The `R` argument indicates a desired result type. We use this in order
# to attempt to preserve type-stability even though the value of `deg` cannot
# be known at compile time. For scalar `p` and `x`, we set R to be the type
# of differentiate(p, x) to give a stable result type regardless of `deg`. For
# vectors p and/or x this is impossible (since differentiate may return an array),
# so we just set `R` to `Any`
function (_differentiate_recursive(p, x, deg::Int, ::Type{R})::R) where {R}
    if deg < 0
        throw(DomainError(deg, "degree is negative"))
    elseif deg == 0
        return p
    else
        return differentiate(differentiate(p, x), x, deg-1)
    end
end

differentiate(p, x, deg::Int) = _differentiate_recursive(p, x, deg, Base.promote_op(differentiate, typeof(p), typeof(x)))
differentiate(p::AbstractArray, x,                              deg::Int) = _differentiate_recursive(p, x, deg, Any)
differentiate(p,                x::Union{AbstractArray, Tuple}, deg::Int) = _differentiate_recursive(p, x, deg, Any)
differentiate(p::AbstractArray, x::Union{AbstractArray, Tuple}, deg::Int) = _differentiate_recursive(p, x, deg, Any)

# This is alternative, Val-based interface for nested differentiation.
# It has the advantage of not requiring an conversion or calls to
# Base.promote_op, while maintaining type stability for any argument
# type.
differentiate(p, x, ::Val{0}) = p
differentiate(p, x, ::Val{1}) = differentiate(p, x)

@static if VERSION < v"v0.7.0-"
    # Marking this @pure helps julia v0.6 figure this out
    Base.@pure _reduce_degree(::Val{N}) where {N} = Val{N - 1}()
    function differentiate(p, x, deg::Val{N}) where N
        if N < 0
            throw(DomainError(deg))
        else
            differentiate(differentiate(p, x), x, _reduce_degree(deg))
        end
    end
else
    # In Julia v0.7 and above, we can remove the _reduce_degree trick
    function differentiate(p, x, deg::Val{N}) where N
        if N < 0
            throw(DomainError(deg))
        else
            differentiate(differentiate(p, x), x, Val{N - 1}())
        end
    end
end

ext: jl


text:
using RandomMatrices
using Base.Test

srand(1)
dx = 0.001

let
p = WhiteNoiseProcess(dx)
@test start(p) == nothing
S = next(p)
@test done(p, S) == false
end

let
p = BrownianProcess(dx)
@test start(p) == 0
S = next(p, 0.0)
@test done(p, S) == false
end

let
p = AiryProcess(dx, 2.0)
S = start(p)
@test S == fill(-2/dx^2, 1, 1)
S = next(p, S)
@test done(p, S) == false
end


text:
function lgd15()
#
#  test polynomial suggested by Goedecker
#   Legendre polynomial of degree 15
#
    p = lgd(15);
    z = [ 0
   0.98799251802063
   0.93727339240040
   0.84820658341064
   0.72441773136009
   0.57097217260857
   0.39415134707756
  -0.98799251802067
  -0.93727339240030
  -0.84820658341075
  -0.72441773136004
   0.20119409399743
  -0.57097217260858
  -0.39415134707756
  -0.20119409399743]
    z = [z ones(15)]

    p, PolyZeros(z)
end


text:
using GOF3R
using Test

@testset "GOF3R.jl" begin
    GOF3R.s3gof3r_jll.gof3r() do exe
        @test success(`$exe --version`)
    end
end


text:
@testset "Confidence intervals ................. " begin

  ss = SampleSpace(5:20, 50)

  params = MESS(
    ss, .2, .4, alpha = .05, beta = .2, pess = .4
  )

  design = optimaldesign(10, params, solver)
  supp   = support(design)

  cpe    = OCEstimator(design, qsolver)

  ci     = CPInterval(design, confidence = .9)
  @test minimum(coverage.(ci, linspace(0, 1))) > .85

  lim  = hcat(limits.(ci, supp[:, 1], supp[:, 2])...)'
  @test all((0 .<= lim[:, 1]) .& (1 .>= lim[:, 2]))

  ci     = ECPInterval(cpe, confidence = .9)
  @test minimum(coverage.(ci, linspace(0, 1))) >= .9
  lim  = hcat(limits.(ci, supp[:, 1], supp[:, 2])...)'
  @test all((0 .<= lim[:, 1]) .& (1 .>= lim[:, 2]))

  # do not test mmwinterval, not recommended for use

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
module Day21
__revise_mode__ = :eval

using Advent: save_input
using DataStructures: counter, inc!

function recurse(dict::Dict{S, Set{T}}, tocheck::Vector{S}) where S where T
    isempty(tocheck) && return dict
    sort!(tocheck; by=x -> length(dict[x]))

    if any(isempty, values(dict)) || length(union(values(dict))) < length(dict)
        return empty(dict)
    end

    k = popfirst!(tocheck)
    for j in dict[k]
        save = deepcopy(dict)
        dict[k] = Set([j])
        for (k2, v2) in dict
            k == k2 || setdiff!(v2, dict[k])
        end
        res = recurse(dict, tocheck)
        isempty(res) || return res
        merge!(dict, save)
    end
    return empty(dict)
end

function parse_input(filename)
    possible = Dict{String,Set{String}}()
    ingred_counts = counter(String)
    for line in eachline(filename)
        _ingreds, _allers = split(line, "(contains ")
        ingreds = split(_ingreds)
        map(Base.Fix1(inc!, ingred_counts), ingreds)

        allers = Set(split(_allers[1:end-1], ", "))
        for aller in allers
            intersect!(get!(possible, aller, Set(ingreds)), ingreds)
        end
    end

    possible, ingred_counts
end

function solve_a(filename)
    possible, ingred_counts = parse_input(filename)
    bad = reduce(union, values(possible))
    sum(v for (k, v) in ingred_counts if !(k in bad))
end

function solve_b(filename)
    possible, _ = parse_input(filename)
    tocheck = possible |> keys |> collect
    bad = Dict(k => only(v) for (k, v) in recurse(possible, tocheck))
    join((bad[k] for k in sort(collect(keys(possible)))), ",")
end

solve(filename = save_input(2020, 21)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using PkgTemplates
using PkgTemplates: FilePlugin, Template
import PkgTemplates: source, destination, view

struct Dockerfile <: FilePlugin
    tag::String
    file::String
    destination::String
    function Dockerfile(; tag = string(VERSION), with_jupyter = false)
        if with_jupyter
            new(tag, joinpath("templates", "with_jupyter", "Dockerfile"), "Dockerfile")
        else
            new(tag, joinpath("templates", "Dockerfile"), "Dockerfile")
        end
    end
end

source(p::Dockerfile) = p.file
destination(p::Dockerfile) = p.destination
view(f::Dockerfile, ::Template, pkg::AbstractString) = Dict("PKG" => pkg, "tag" => f.tag)

# ---

Base.@kwdef struct Jupytext <: FilePlugin
    file::String = "templates/jupytext.toml"
    destination::String = "jupytext.toml"
end

source(p::Jupytext) = p.file
destination(p::Jupytext) = p.destination

# ---

Base.@kwdef struct PlaygroundNotebook <: FilePlugin
    file::String = "templates/playground/notebook/.gitkeep"
    destination::String = "playground/notebook/.gitkeep"
end

source(p::PlaygroundNotebook) = p.file
destination(p::PlaygroundNotebook) = p.destination

# ---

Base.@kwdef struct PlaygroundPluto <: FilePlugin
    file::String = "templates/playground/pluto/.gitkeep"
    destination::String = "playground/pluto/.gitkeep"
end

source(p::PlaygroundPluto) = p.file
destination(p::PlaygroundPluto) = p.destination

# ---

Base.@kwdef struct Jupytext <: FilePlugin
    file::String = "templates/jupytext.toml"
    destination::String = "jupytext.toml"
end

source(p::Jupytext) = p.file
destination(p::Jupytext) = p.destination

# ---

Base.@kwdef struct DevContainer <: FilePlugin
    file::String = "templates/dot_devcontainer/devcontainer.json"
    destination::String = ".devcontainer/devcontainer.json"
end

source(p::DevContainer) = p.file
destination(p::DevContainer) = p.destination

# ---

struct DockerCompose <: FilePlugin
    file::String
    destination::String
    function DockerCompose(; with_jupyter = false)
        if with_jupyter
            new(joinpath("templates", "with_jupyter", "docker-compose.yml"), "docker-compose.yml")
        else
            new(joinpath("templates", "docker-compose.yml"), "docker-compose.yml")
        end
    end
end

source(p::DockerCompose) = p.file
destination(p::DockerCompose) = p.destination
view(::DockerCompose, ::Template, pkg::AbstractString) = Dict("DOCKER_IMAGE" => lowercase(pkg) * "jl", "PKG" => pkg)

# ---

Base.@kwdef struct Makefile <: FilePlugin
    file::String = "templates/Makefile"
    destination::String = "Makefile"
end

source(p::Makefile) = p.file
destination(p::Makefile) = p.destination
view(::Makefile, ::Template, pkg::AbstractString) = Dict("DOCKER_IMAGE" => lowercase(pkg) * "jl")

# ---

Base.@kwdef struct VSCodeExtensions <: FilePlugin
    file::String = "templates/dot_vscode/extensions.json"
    destination::String = ".vscode/extensions.json"
end

source(p::VSCodeExtensions) = p.file
destination(p::VSCodeExtensions) = p.destination


text:
"""
    PersistencePair{T, U}

fields:

* `birth::T`
* `death::T`
* `data::U`
"""
struct PersistencePair{T, U}
    birth ::T
    death ::T
    data  ::U
end

function PersistencePair(birth, death, data::U=nothing) where U
    b, d = promote(birth, death)
    T    = typeof(b)
    PersistencePair{T, U}(b, d, data)
end

function Base.show(io::IO, pair::PersistencePair{T, U}) where {T, U}
    compact = get(io, :compact, false)
    if !compact
        print(io, "PersistencePair{$T, $U}")
    end
    print(io, "($(birth(pair)), $(death(pair))")
    if U != Void
        print(io, "; ")
        showcompact(io, data(pair))
    end
    print(io, ")")
end

"""
    birth(pair)

Get the birth time of a `PersistencePair`.
"""
birth(pair::PersistencePair) = pair.birth
"""
    death(pair)

Get the death time of a `PersistencePair`.
"""
death(pair::PersistencePair) = pair.death
"""
    data(pair)

Get the data associated with a `PersistencePair`.
"""
data(pair::PersistencePair) = pair.data

Base.start(pair::PersistencePair) = 1
Base.next(pair::PersistencePair, i::Int) =
    if i == 1
        (birth(pair), 2)
    else
        (death(pair), 3)
    end
Base.done(pair::PersistencePair, i::Int) = i == 3

# ============================================================================ #
"""
    PersistenceBarcode{T, U, D, A<:AbstractVector{PersistencePair{T, U}}}

fields:

* `barcodes::SVector{D, A}`

The `barcodes` contains a vector of `PersistencePair`s for each dimension.

constructors:

`PersistenceBarcode(arr::AbstractVector{A})`
`PersistenceBarcode(arrs::Vararg{A})`
"""
struct PersistenceBarcode{T, U, D, A<:AbstractVector{PersistencePair{T, U}}}
    barcodes::SVector{D, A}
end

function PersistenceBarcode(arr::AbstractVector{A}) where
        A<:AbstractVector{PersistencePair{T, U}} where {T, U}
    D = length(arr)
    PersistenceBarcode{T, U, D, A}(SVector{D, A}(arr))
end

function PersistenceBarcode(arrs::Vararg{A, D}) where
        {A<:AbstractVector{PersistencePair{T, U}}, D} where {T, U}
    PersistenceBarcode{T, U, D, A}(SVector{D}([arrs...]))
end

function Base.show(io::IO, barcode::PersistenceBarcode{T, U, D}) where {T, U, D}
    compact = get(io, :compact, false)
    if compact
        print(io, "$(D-1)-d PersistenceBarcode{$T, $U}")
    else
        print(io, "$(D-1)-d PersistenceBarcode{$T, $U}:")
        for i in 0:D-1
            print(io, "\n dim $i:")
            for p in barcode[i]
                print(io, "\n  ")
                showcompact(io, p)
            end
        end
    end
end

function Base.:(==)(barcode1::PersistenceBarcode, barcode2::PersistenceBarcode)
    dim(barcode1) == dim(barcode2) || return false
    for i in eachindex(barcode1)
        barcode1[i] == barcode2[i] || return false
    end
    true
end

"""
    dim(barcode)

Get the dimension of a `PersistenceBarcode`.
"""
dim(::PersistenceBarcode{T, U, D}) where {T, U, D} = D - 1

# Indexing (zero-based!)
Base.getindex(barcode::PersistenceBarcode, i) = barcode.barcodes[i + 1]
Base.length(barcode::PersistenceBarcode{T, U, D}) where {T, U, D} = D
Base.endof(barcode::PersistenceBarcode) = endof(barcode.barcodes) - 1
Base.eachindex(barcode::PersistenceBarcode) = 0:dim(barcode)

# Iteration
Base.start(barcode::PersistenceBarcode) = start(barcode.barcodes)
Base.next(barcode::PersistenceBarcode, st) = next(barcode.barcodes, st)
Base.done(barcode::PersistenceBarcode, st) = done(barcode.barcodes, st)
Base.iteratorsize(::PersistenceBarcode) = Base.HasLength()
Base.iteratoreltype(barcode::PersistenceBarcode) = Base.HasEltype()
Base.eltype(barcode::PersistenceBarcode{T, U, D, A}) where {T, U, D, A} = A

Base.filter(f, barcode::PersistenceBarcode) =
    PersistenceBarcode(map(b -> filter(f, b), barcode.barcodes))

Base.map(f, barcode::PersistenceBarcode) =
    PersistenceBarcode(map(b -> map(f, b), barcode.barcodes))


text:
abstract type Flag end
struct UpdateBetaWithSkewT <: Flag end
struct UpdateLambdaWithSkewT <: Flag end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module Tutorial5
using FEniCS

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "P", 1)
#in the original tutorial, these are calculated via SymPy
#although it has no actual FEniCS code in calculating those, so I have simply
#stated the result
u_code = "x[0]+2*x[1]+1"
f_code = "-10*x[0]-20*x[1]-10"

calc(var) = 1+var*var

u_D = Expression(u_code, degree=2)
bc = DirichletBC(V, u_D, "on_boundary")
u=FeFunction(V)
v = TestFunction(V)
f = Expression(f_code, degree=2)
F = calc(u)*dot(grad(u), grad(v))*dx - f*v*dx
nlvsolve(F,u,bc)
end#module


text:
push!(LOAD_PATH, pwd())

using QuIPS

quip = [
    (:X, 3),
    (:X, 2),
    (:CCX, (3, 2, 1)),
]

QC = QCircuit(quip, 3)

@time run!(QC)

println(abs2.(QC.wfn))

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
module GuessworkQuantumSideInfo

using LinearAlgebra, SparseArrays # standard libraries
using Random: Random, randperm!, AbstractRNG # standard library
using Convex # Julia SDP solvers / interfaces
using Combinatorics: multiset_permutations, combinations # used for POVM outcomes
using UnPack # helper

# Basic functions for quantum states and examples
export ket, bra, ⊗, I, dm, randdm, randprobvec, iid_copies, BB84_states
include("quantum_states.jl")

# The basic SDP routine
export guesswork
include("SDP_formulation.jl")

export guesswork_lower_bound
include("lower_bound.jl")

# Re-run the SDP adding constraints one at a time
export guesswork_upper_bound
include("upper_bound.jl")

# reformulate as a mixed-integer SDP
export guesswork_MISDP
include("MISDP_formulation.jl")

# Utilities for analyzing the POVMs resulting from the SDPs
export pmfN
include("analyze_measurements.jl")

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Storage

Currently only electric storage is modeled but thermal storage will be added soon.

TODO update this doc string and perhaps the way we define arbitrary storage.

```julia
Base.@kwdef struct ElecStorage <: AbstractStorage
    min_kw::Float64 = 0.0
    max_kw::Float64 = 1.0e4
    min_kwh::Float64 = 0.0
    max_kwh::Float64 = 1.0e6
    internal_efficiency_pct::Float64 = 0.975
    inverter_efficiency_pct::Float64 = 0.96
    rectifier_efficiency_pct::Float64 = 0.96
    soc_min_pct::Float64 = 0.2
    soc_init_pct::Float64 = 0.5
    can_grid_charge::Bool = true
    installed_cost_per_kw::Float64 = 840.0
    installed_cost_per_kwh::Float64 = 420.0
    replace_cost_per_kw::Float64 = 410.0
    replace_cost_per_kwh::Float64 = 200.0
    inverter_replacement_year::Int = 10
    battery_replacement_year::Int = 10
    macrs_option_years::Int = 7
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    total_itc_pct::Float64 = 0.0
    total_rebate_per_kw::Float64 = 0.0
    total_rebate_per_kwh::Float64 = 0.0
end
```
"""
Base.@kwdef struct ElecStorage <: AbstractStorage
    min_kw::Float64 = 0.0
    max_kw::Float64 = 1.0e4
    min_kwh::Float64 = 0.0
    max_kwh::Float64 = 1.0e6
    internal_efficiency_pct::Float64 = 0.975
    inverter_efficiency_pct::Float64 = 0.96
    rectifier_efficiency_pct::Float64 = 0.96
    soc_min_pct::Float64 = 0.2
    soc_init_pct::Float64 = 0.5
    can_grid_charge::Bool = true
    installed_cost_per_kw::Float64 = 840.0
    installed_cost_per_kwh::Float64 = 420.0
    replace_cost_per_kw::Float64 = 410.0
    replace_cost_per_kwh::Float64 = 200.0
    inverter_replacement_year::Int = 10
    battery_replacement_year::Int = 10
    macrs_option_years::Int = 7
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    total_itc_pct::Float64 = 0.0
    total_rebate_per_kw::Float64 = 0.0
    total_rebate_per_kwh::Float64 = 0.0
end

struct Storage <: AbstractStorage
    types::Array{Symbol,1}
    raw_inputs::Dict{Symbol, AbstractStorage}
    min_kw::Dict{Symbol, Float64}
    max_kw::Dict{Symbol, Float64}
    min_kwh::Dict{Symbol, Float64}
    max_kwh::Dict{Symbol, Float64}
    charge_efficiency::Dict{Symbol, Float64}
    discharge_efficiency::Dict{Symbol, Float64}
    soc_min_pct::Dict{Symbol, Float64}
    soc_init_pct::Dict{Symbol, Float64}
    installed_cost_per_kw::Dict{Symbol, Float64}
    installed_cost_per_kwh::Dict{Symbol, Float64}
    can_grid_charge::Array{Symbol,1}
    grid_charge_efficiency::Dict{Symbol, Float64}
end

"""
    # function Storage(d::Dict{Symbol,Dict}, f::Financial)

Construct Storage struct from Dict with keys for each storage type (eg. :elec) and values with
    input dicts for each storage type. Note that the REopt inputs are indexed on the storage type.
"""
function Storage(d::Dict, f::Financial)  # nested dict
    types = Symbol[]
    can_grid_charge = Symbol[]
    raw_vals = Dict(zip(fieldnames(Storage), [Any[] for _ in range(1, stop=fieldcount(Storage))]))
    # raw_vals has keys = fieldnames(Storage) and values = arrays of values for each type in types

    for (storage_type, input_dict) in d

        push!(types, storage_type)
        struct_name = string(titlecase(string(storage_type)) * "Storage")  # eg. ElecStorage
        storage_instance = eval(Meta.parse(struct_name * "(;$input_dict...)"))

        if storage_instance.can_grid_charge
            push!(can_grid_charge, storage_type)
        end
        fill_storage_vals!(raw_vals, storage_instance, f)
    end

    storage_args = Dict(
        :types => types,
        :can_grid_charge => can_grid_charge,
    )
    d2 = Dict()  # Julia won't let me use storage_args: "unable to check bounds for indices of type Symbol"
    for k in keys(raw_vals)
        d2[k] = Dict(zip(types, raw_vals[k]))
    end

    grid_charge_efficiency = Dict(:elec => convert(Float64, d2[:charge_efficiency][:elec]))

    return Storage(
        storage_args[:types],
        d2[:raw_inputs],
        d2[:min_kw],
        d2[:max_kw],
        d2[:min_kwh],
        d2[:max_kwh],
        d2[:charge_efficiency],
        d2[:discharge_efficiency],
        d2[:soc_min_pct],
        d2[:soc_init_pct],
        d2[:installed_cost_per_kw],
        d2[:installed_cost_per_kwh],
        storage_args[:can_grid_charge],
        grid_charge_efficiency
    )
    # TODO expand for smart thermostat
end

"""
    fill_storage_vals!(d::Dict, s::AbstractStorage, f::Financial)

Fill `d`'s values using `s` and `f`. The dict `d` must have keys equivalent to the fieldnames(Storage)
"""
function fill_storage_vals!(d::Dict, s::AbstractStorage, f::Financial)
    push!(d[:raw_inputs], s)
    push!(d[:min_kw], s.min_kw)
    push!(d[:max_kw], s.max_kw)
    push!(d[:min_kwh], s.min_kwh)
    push!(d[:max_kwh], s.max_kwh)
    push!(d[:soc_min_pct], s.soc_min_pct)
    push!(d[:soc_init_pct], s.soc_init_pct)

    push!(d[:charge_efficiency], s.rectifier_efficiency_pct * s.internal_efficiency_pct^0.5)
    push!(d[:discharge_efficiency], s.inverter_efficiency_pct * s.internal_efficiency_pct^0.5)

    push!(d[:installed_cost_per_kw], effective_cost(;
        itc_basis=s.installed_cost_per_kw,
        replacement_cost=s.replace_cost_per_kw,
        replacement_year=s.inverter_replacement_year,
        discount_rate=f.owner_discount_pct,
        tax_rate=f.owner_tax_pct,
        itc=s.total_itc_pct,
        macrs_schedule = s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
        macrs_bonus_pct=s.macrs_bonus_pct,
        macrs_itc_reduction = s.macrs_itc_reduction,
        rebate_per_kw = s.total_rebate_per_kw
    ))
    push!(d[:installed_cost_per_kwh], effective_cost(;
        itc_basis=s.installed_cost_per_kwh,
        replacement_cost=s.replace_cost_per_kwh,
        replacement_year=s.inverter_replacement_year,
        discount_rate=f.owner_discount_pct,
        tax_rate=f.owner_tax_pct,
        itc=s.total_itc_pct,
        macrs_schedule=s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
        macrs_bonus_pct=s.macrs_bonus_pct,
        macrs_itc_reduction = s.macrs_itc_reduction
    ))
    d[:installed_cost_per_kwh][end] -= s.total_rebate_per_kwh
end


text:
#################################################################################################
# AIM 1: nanopolish's LRT vs Accuracy
#################################################################################################
## Deps
using Dierckx
using CodecZlib
using StatsPlots
using Distributions
using DelimitedFiles
using Plots.PlotMeasures

## Constants
const noise_levels = ["$(s)" for s in collect(2.0:0.5:3.5)]
const data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-1"

# Colors
const blnd_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]
const s_color_cde = Dict("0.5" => blnd_col[1], "1.0" => blnd_col[2], "1.5" => blnd_col[3], "2.0" => blnd_col[4],
                         "2.5" => blnd_col[5], "3.0" => blnd_col[6], "3.5" => blnd_col[7], "4.0" => blnd_col[8])

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Functions
#################################################################################################

# Function to compute accuracy
function comp_accu(truth, pred)

    # Return accuracy
    return sum(truth .== pred) / length(truth)

end

function get_perc_calls(s, lr_thrs)

    # Read in
    stream = GzipDecompressorStream(open("$(data_dir)/nanopolish/gm12878_chr22_sigma_$(s)_nanopolish_tuples_sample.tsv.gz"))
    in_data = readdlm(stream)
    close(stream)

    # Remove NaN
    kp_ind = .! map(isequal("n/a"), in_data[:,2])
    in_data = in_data[kp_ind,:]

    # Get data
    perc_calls = []
    n_pairs = length(in_data[:,2]) - 1
    for lr in lr_thrs
        println("$(lr)")
        # Get valid pairs
        val_pairs = [all(abs.(in_data[i:(i + 1),2]) .>= lr) for i = 1:n_pairs]
        # Push
        push!(perc_calls, sum(val_pairs) / n_pairs)
    end

    # Return
    return perc_calls

end

function get_acc_calls(s, lr_thrs)

    # Read in
    stream = GzipDecompressorStream(open("$(data_dir)/nanopolish/gm12878_chr22_sigma_$(s)_nanopolish_tuples_sample.tsv.gz"))
    in_data = readdlm(stream)
    close(stream)

    # Remove NaN
    kp_ind = .! map(isequal("n/a"), in_data[:,2])
    in_data = in_data[kp_ind,:]

    # Get all XX data
    n_pairs = length(in_data[:,2]) - 1
    true_xx = [in_data[i,1] * in_data[i + 1,1] for i = 1:n_pairs]
    pred_xx = [in_data[i,3] * in_data[i + 1,3] for i = 1:n_pairs]

    # Get data
    accs = []
    for lr in lr_thrs
        println("$(lr)")
        # Get valid pairs
        val_pairs = [all(abs.(in_data[i:(i + 1),2]) .>= lr) for i = 1:n_pairs]
        if sum(val_pairs) > 0
            # Keep valid pairs
            true_xx_thrs = true_xx[val_pairs]
            pred_xx_thrs = pred_xx[val_pairs]
            # Push
            push!(accs, comp_accu(true_xx_thrs, pred_xx_thrs))
        else
            # Push
            push!(accs, 1.0)
        end
    end

    # Return
    return accs

end

#################################################################################################
# Calls
#################################################################################################
max_lrt_plt = 20.0
lr_thrs = vcat(collect(0.0:0.1:10.0), collect(10.0:0.5:25.0), collect(25.0:25.0:125.0))

# Plot percentage of calls
plt_perc = plot(xlabel="LRT", ylabel="percentage of calls", xlim=(0, max_lrt_plt), ylim=(0, 1.0));
for s in noise_levels
    println("Sigma $(s)")

    # Get data
    perc_calls = get_perc_calls(s, lr_thrs)

    # Plot spline
    plot!(plt_perc, lr_thrs, perc_calls, label="\\sigma=$(s)", color=s_color_cde["$(s)"])

end

# Add vertical line
vline!(plt_perc,[2.5],linestyle=:dash,color="black",label="LRT=2.5")

# Save
plot(plt_perc)
# savefig("$(data_dir)/Nanopolish-LRT-vs-Percentage-Calls-XX.pdf")

# Plot accuracy of calls per noise level
plt_err = plot(xlabel="LRT", ylabel="percentage call error", xlim=(0, max_lrt_plt), ylim=(0, 0.31));
for s in noise_levels
    println("Sigma $(s)")

    # Get data
    acc_calls = get_acc_calls(s, lr_thrs)

    # Fit spline
    spl = Spline1D(lr_thrs, 1.0 .- acc_calls, k=3, s=1e-2)

    # Plot spline
    plot!(plt_err, lr_thrs, spl(lr_thrs), label="\\sigma=$(s)", color=s_color_cde["$(s)"])

end

# Add vertical line
vline!(plt_err,[2.5],linestyle=:dash,color="black",label="LRT=2.5")

# Save individual plot
plot(plt_err)
# savefig("$(data_dir)/Nanopolish-LRT-vs-Error-Calls-Per-Sigma-XX.pdf")

# Plot accuracy of calls per noise level
plt_err_acc = plot(xlabel="percentage calls", ylabel="percentage call error", xlim=(0, 1), ylim=(0, 0.31), legend=:topleft);
for s in noise_levels
    println("Sigma $(s)")

    # Get data
    perc_calls = get_perc_calls(s, lr_thrs)
    acc_calls = get_acc_calls(s, lr_thrs)

    # Add (0,0) point
    push!(perc_calls, 0.0)
    push!(acc_calls, 1.0)

    # Fit spline
    spl = Spline1D(1.0 .- perc_calls, 1.0 .- acc_calls, k=4, s=1e-1)

    # Plot spline
    plot!(plt_err_acc, 1.0 .- perc_calls, spl(perc_calls), label="\\sigma=$(s)", color=s_color_cde["$(s)"])

end

# Save
plot(plt_err_acc)
# savefig("$(data_dir)/Nanopolish-Perc-vs-Error-Calls-XX.pdf")

plot(plt_err,plt_perc,plt_err_acc,layout=(3, 1),size=(700, 1000))
savefig("$(data_dir)/Nanopolish-LRT-All-XX.pdf")


text:
import FillArrays: getindex_value, fillsimilar, fillzero

"""
    AbstractQuasiFill{T, N, Axes} <: AbstractQuasiArray{T, N}

Supertype for lazy array types whose entries are all equal to constant.
"""
abstract type AbstractQuasiFill{T, N, Axes} <: AbstractQuasiArray{T, N} end

Base.@propagate_inbounds @inline function _fill_getindex(::Type{IND}, F::AbstractQuasiFill, kj::IND) where IND
    @boundscheck checkbounds(F, kj...)
    getindex_value(F)
end

Base.@propagate_inbounds @inline function _fill_getindex(::Type{IND}, A::AbstractQuasiFill{<:Any,N}, I) where {N,IND}
    @boundscheck checkbounds(A, I...)
    shape = Base.index_shape(I...)
    fillsimilar(A, shape...)
end

==(a::AbstractQuasiFill, b::AbstractQuasiFill) = axes(a) == axes(b) && getindex_value(a) == getindex_value(b)

_getindex(::Type{IND}, F::AbstractQuasiFill, k::IND) where IND = _fill_getindex(IND, F, k)

@inline function _setindex!(::Type{IND}, F::AbstractQuasiFill, v, k::IND) where IND
    @boundscheck checkbounds(F, k...)
    v == getindex_value(F) || throw(ArgumentError("Cannot setindex! to $v for an AbstractQuasiFill with value $(getindex_value(F))."))
    F
end

@inline function fill!(F::AbstractQuasiFill, v)
    v == getindex_value(F) || throw(ArgumentError("Cannot fill! with $v an AbstractQuasiFill with value $(getindex_value(F))."))
    F
end

rank(F::AbstractQuasiFill) = iszero(getindex_value(F)) ? 0 : 1

"""
    QuasiFill{T, N, Axes}

A lazy representation of a quasi-array of dimension `N`
whose entries are all equal to a constant of type `T`,
with axes of type `Axes`.
Typically created by `QuasiFill`
```
"""
struct QuasiFill{T, N, Axes} <: AbstractQuasiFill{T, N, Axes}
    value::T
    axes::Axes

    QuasiFill{T,N,Axes}(x::T, sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
        new{T,N,Axes}(x,sz)
end

QuasiFill{T,N,Axes}(x, sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
    QuasiFill{T,N,Axes}(convert(T, x)::T, sz)

@inline QuasiFill{T, N}(x::T, sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
    QuasiFill{T,N,Axes}(x, sz)
@inline QuasiFill{T, N}(x, sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
    QuasiFill{T,N}(convert(T, x)::T, sz)
QuasiFill{T}(x, sz::NTuple{N,Any}) where {T,N} = QuasiFill{T,N}(x, sz)
QuasiFill(x::T, sz::NTuple{N,Any}) where {T,N} = QuasiFill{T,N}(x, sz)

QuasiFill{T,N}(x, sz...) where {T,N} = QuasiFill{T,N}(x, map(_inclusion, sz))
QuasiFill{T}(x, sz...) where T = QuasiFill{T}(x, map(_inclusion, sz))
QuasiFill(x, sz...) = QuasiFill(x, map(_inclusion, sz))

# We restrict to  when T is specified to avoid ambiguity with a QuasiFill of a QuasiFill
@inline QuasiFill{T}(F::QuasiFill{T}) where T = F
@inline QuasiFill{T,N}(F::QuasiFill{T,N}) where {T,N} = F
@inline QuasiFill{T,N,Axes}(F::QuasiFill{T,N,Axes}) where {T,N,Axes} = F

@inline axes(F::QuasiFill) = F.axes
@inline size(F::QuasiFill) = length.(F.axes)
MemoryLayout(::Type{<:QuasiFill}) = FillLayout()

@inline getindex_value(F::QuasiFill) = F.value

AbstractQuasiArray{T}(F::QuasiFill{T}) where T = F
AbstractQuasiArray{T,N}(F::QuasiFill{T,N}) where {T,N} = F
AbstractQuasiArray{T}(F::QuasiFill{V,N}) where {T,V,N} = QuasiFill{T}(convert(T, F.value)::T, F.axes)
AbstractQuasiArray{T,N}(F::QuasiFill{V,N}) where {T,V,N} = QuasiFill{T}(convert(T, F.value)::T, F.axes)

convert(::Type{AbstractQuasiArray{T}}, F::QuasiFill{T}) where T = F
convert(::Type{AbstractQuasiArray{T,N}}, F::QuasiFill{T,N}) where {T,N} = F
convert(::Type{AbstractQuasiArray{T}}, F::QuasiFill) where {T} = AbstractQuasiArray{T}(F)
convert(::Type{AbstractQuasiArray{T,N}}, F::QuasiFill) where {T,N} = AbstractQuasiArray{T,N}(F)
convert(::Type{AbstractQuasiFill}, F::AbstractQuasiFill) = F
convert(::Type{AbstractQuasiFill{T}}, F::AbstractQuasiFill) where T = convert(AbstractQuasiArray{T}, F)
convert(::Type{AbstractQuasiFill{T,N}}, F::AbstractQuasiFill) where {T,N} = convert(AbstractQuasiArray{T,N}, F)

copy(F::QuasiFill) = QuasiFill(F.value, F.axes)

# ambiguity fix
convert(::Type{T}, F::T) where T<:QuasiFill = F

getindex(F::QuasiFill{<:Any,0}) = getindex_value(F)

sort(a::AbstractQuasiFill; kwds...) = a
sort!(a::AbstractQuasiFill; kwds...) = a

for (Typ, funcs, func) in ((:QuasiZeros, :zeros, :zero), (:QuasiOnes, :ones, :one))
    @eval begin
        """ `$($Typ){T, N, Axes} <: AbstractQuasiFill{T, N, Axes}` (lazy `$($funcs)` with axes)"""
        struct $Typ{T, N, Axes} <: AbstractQuasiFill{T, N, Axes}
            axes::Axes
            @inline $Typ{T,N,Axes}(sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
                new{T,N,Axes}(sz)
            @inline $Typ{T,N}(sz::Axes) where Axes<:Tuple{Vararg{Any,N}} where {T, N} =
                new{T,N,Axes}(sz)
        end
        @inline $Typ{T}(sz::SZ) where SZ<:Tuple{Vararg{Any,N}} where {T, N} = $Typ{T, N}(sz)
        @inline $Typ(sz::SZ) where SZ<:Tuple{Vararg{Any,N}} where N = $Typ{Float64,N}(sz)
        @inline $Typ(sz...) = $Typ(map(_inclusion,sz))
        @inline $Typ{T}(sz...) where T = $Typ{T}(map(_inclusion,sz))

        @inline axes(Z::$Typ) = Z.axes
        @inline size(Z::$Typ) = length.(Z.axes)
        @inline getindex_value(Z::$Typ{T}) where T = $func(T)

        AbstractQuasiArray{T}(F::$Typ{T}) where T = F
        AbstractQuasiArray{T,N}(F::$Typ{T,N}) where {T,N} = F
        AbstractQuasiArray{T}(F::$Typ) where T = $Typ{T}(F.axes)
        AbstractQuasiArray{T,N}(F::$Typ{V,N}) where {T,V,N} = $Typ{T}(F.axes)
        convert(::Type{AbstractQuasiArray{T}}, F::$Typ{T}) where T = AbstractQuasiArray{T}(F)
        convert(::Type{AbstractQuasiArray{T,N}}, F::$Typ{T,N}) where {T,N} = AbstractQuasiArray{T,N}(F)
        convert(::Type{AbstractQuasiArray{T}}, F::$Typ) where T = AbstractQuasiArray{T}(F)
        convert(::Type{AbstractQuasiArray{T,N}}, F::$Typ) where {T,N} = AbstractQuasiArray{T,N}(F)

        copy(F::$Typ) = F

        getindex(F::$Typ{T,0}) where T = getindex_value(F)
    end
end

"""
    fillsimilar(a::AbstractQuasiFill, axes)

creates a fill object that has the same fill value as `a` but
with the specified axes.
For example, if `a isa QuasiZeros` then so is the returned object.
"""
fillsimilar(a::QuasiOnes{T}, axes...) where T = QuasiOnes{T}(axes)
fillsimilar(a::QuasiZeros{T}, axes...) where T = QuasiZeros{T}(axes)
fillsimilar(a::AbstractQuasiFill, axes...) = QuasiFill(getindex_value(a), axes)

fillsimilar(a::QuasiOnes{T}, axes::AbstractUnitRange{Int}...) where T = Ones{T}(axes)
fillsimilar(a::QuasiZeros{T}, axes::AbstractUnitRange{Int}...) where T = Zeros{T}(axes)
fillsimilar(a::AbstractQuasiFill, axes::AbstractUnitRange{Int}...) = Fill(getindex_value(a), axes)

rank(F::QuasiZeros) = 0
rank(F::QuasiOnes) = 1

diag(F::QuasiZeros{T,2}) where T = QuasiZeros{T}((axes(F,1),))

const QuasiEye{T,Axes} = QuasiDiagonal{T,QuasiOnes{T,1,Tuple{Axes}}}

@inline QuasiEye{T}(n) where T = QuasiDiagonal(QuasiOnes{T}(n))
@inline QuasiEye(n) = QuasiDiagonal(QuasiOnes(n))

# function iterate(iter::Eye, istate = (1, 1))
#     (i::Int, j::Int) = istate
#     m = size(iter, 1)
#     return i > m ? nothing :
#         ((@inbounds getindex(iter, i, j)),
#          j == m ? (i + 1, 1) : (i, j + 1))
# end

isone(::QuasiEye) = true

@inline QuasiEye{T}(A::AbstractQuasiMatrix) where T = QuasiEye{T}(axes(A))
@inline QuasiEye(A::AbstractQuasiMatrix) = QuasiEye{eltype(A)}(axes(A))

quasiscaling(J, ax) = QuasiDiagonal(QuasiFill(J.λ, ax))

for op in (:+, :-, :*, :/, :\)
    @eval begin
        $op(A::AbstractQuasiMatrix, J::UniformScaling) = $op(A, quasiscaling(J, axes(A,2)))
        $op(J::UniformScaling, A::AbstractQuasiMatrix) = $op(quasiscaling(J, axes(A,1)), A)
    end
end

########
# show
#######

function summary(io::IO, F::QuasiOnes)
    print(io, "ones(")
    summary(io, F.axes[1])
    for a in tail(F.axes)
        print(io, ", ")
        summary(io, a)
    end
    print(io, ")")
end

function summary(io::IO, F::QuasiZeros)
    print(io, "zeros(")
    summary(io, F.axes[1])
    for a in tail(F.axes)
        print(io, ", ")
        summary(io, a)
    end
    print(io, ")")
end

function summary(io::IO, F::QuasiFill)
    print(io, "fill($(F.value), ")
    summary(io, F.axes[1])
    for a in tail(F.axes)
        print(io, ", ")
        summary(io, a)
    end
    print(io, ")")
end

#########
#  Special matrix types
#########

function convert(::Type{QuasiDiagonal}, Z::QuasiZeros{T,2}) where T
    n,m = axes(Z)
    n ≠ m && throw(BoundsError(Z))
    QuasiDiagonal(QuasiZeros{T}((n,)))
end

function convert(::Type{QuasiDiagonal{T}}, Z::QuasiZeros{V,2}) where {T,V}
    n,m = axes(Z)
    n ≠ m && throw(BoundsError(Z))
    QuasiDiagonal(QuasiZeros{T}((n,)))
end

#########
# maximum/minimum
#########

for op in (:maximum, :minimum)
    @eval $op(x::AbstractQuasiFill) = getindex_value(x)
end

#########
# Cumsum
#########

sum(x::AbstractQuasiFill) = getindex_value(x)*measure(axes(x,1))
sum(x::QuasiZeros) = getindex_value(x)

# define `sum(::Callable, ::AbstractQuasiFill)` to avoid method ambiguity errors on Julia 1.0
sum(f, x::AbstractQuasiFill) = _sum(f, x)
sum(f::Base.Callable, x::AbstractQuasiFill) = _sum(f, x)
_sum(f, x::AbstractQuasiFill) = measure(x) * f(getindex_value(x))

#########
# zero
#########

zero(r::QuasiZeros{T,N}) where {T,N} = r
zero(r::QuasiOnes{T,N}) where {T,N} = QuasiZeros{T,N}(r.axes)
zero(r::QuasiFill{T,N}) where {T,N} = QuasiZeros{T,N}(r.axes)

#########
# in
#########
in(x, A::AbstractQuasiFill) = x == getindex_value(A)

####
# Algebra
####

+(a::AbstractQuasiFill) = a
-(a::AbstractQuasiFill) = QuasiFill(-getindex_value(a), axes(a))

# Fill +/- Fill
function +(a::AbstractQuasiFill{T, N}, b::AbstractQuasiFill{V, N}) where {T, V, N}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return QuasiFill(getindex_value(a) + getindex_value(b), axes(a))
end
-(a::AbstractQuasiFill, b::AbstractQuasiFill) = a + (-b)

## Transpose/Adjoint
# cannot do this for vectors since that would destroy scalar dot product

transpose(a::QuasiOnes{T,2}) where T = QuasiOnes{T}(reverse(a.axes))
adjoint(a::QuasiOnes{T,2}) where T = QuasiOnes{T}(reverse(a.axes))
transpose(a::QuasiZeros{T,2}) where T = QuasiZeros{T}(reverse(a.axes))
adjoint(a::QuasiZeros{T,2}) where T = QuasiZeros{T}(reverse(a.axes))
transpose(a::QuasiFill{T,2}) where T = QuasiFill{T}(transpose(a.value), reverse(a.axes))
adjoint(a::QuasiFill{T,2}) where T = QuasiFill{T}(adjoint(a.value), reverse(a.axes))

permutedims(a::AbstractQuasiFill{<:Any,1}) = fillsimilar(a, Inclusion(1), axes(a,1))
permutedims(a::AbstractQuasiFill{<:Any,2}) = fillsimilar(a, reverse(a.axes)...)

## Algebraic identities

+(a::QuasiZeros) = a
-(a::QuasiZeros) = a

# Zeros +/- Zeros
function +(a::QuasiZeros{T}, b::QuasiZeros{V}) where {T, V}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return QuasiZeros{promote_type(T,V)}(axes(a))
end
-(a::QuasiZeros, b::QuasiZeros) = -(a + b)
-(a::QuasiOnes, b::QuasiOnes) = QuasiZeros(a)+QuasiZeros(b)

# Zeros +/- Fill and Fill +/- Zeros
function +(a::AbstractQuasiFill{T}, b::QuasiZeros{V}) where {T, V}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return convert(AbstractQuasiFill{promote_type(T, V)}, a)
end
+(a::QuasiZeros, b::AbstractQuasiFill) = b + a
-(a::AbstractQuasiFill, b::QuasiZeros) = a + b
-(a::QuasiZeros, b::AbstractQuasiFill) = a + (-b)

# QuasiZeros +/- Array and Array +/- QuasiZeros
function +(a::QuasiZeros{T, N}, b::AbstractQuasiArray{V, N}) where {T, V, N}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractQuasiArray{promote_type(T,V),N}(b)
end
function +(a::QuasiArray{T, N}, b::QuasiZeros{V, N}) where {T, V, N}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractQuasiArray{promote_type(T,V),N}(a)
end

function -(a::QuasiZeros{T, N}, b::AbstractQuasiArray{V, N}) where {T, V, N}
    axes(a) ≠ axes(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::QuasiArray{T, N}, b::QuasiZeros{V, N}) where {T, V, N} = a + b

##
# view
##

Base.@propagate_inbounds _unsafe_view(::Type{IND}, A::AbstractQuasiFill{<:Any,N}, I::Tuple) where {IND,N} =
    _fill_getindex(IND, A, Base.to_indices(A,I))

# not getindex since we need array-like indexing
Base.@propagate_inbounds function _unsafe_view(::Type{IND}, A::AbstractQuasiFill{<:Any,N}, I::Vararg{IND, N}) where {N,IND}
    @boundscheck checkbounds(A, I...)
    fillsimilar(A)
end

####
# broadcast
####

### map

map(f::Function, r::AbstractQuasiFill) = QuasiFill(f(getindex_value(r)), axes(r))

### Unary broadcasting

function broadcasted(::AbstractQuasiArrayStyle{N}, op, r::AbstractQuasiFill{T,N}) where {T,N}
    return QuasiFill(op(getindex_value(r)), axes(r))
end

broadcasted(::AbstractQuasiArrayStyle, ::typeof(+), r::QuasiZeros) = r
broadcasted(::AbstractQuasiArrayStyle, ::typeof(-), r::QuasiZeros) = r
broadcasted(::AbstractQuasiArrayStyle, ::typeof(+), r::QuasiOnes) = r

broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(conj), r::QuasiZeros{T,N}) where {T,N} = r
broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(conj), r::QuasiOnes{T,N}) where {T,N} = r
broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(real), r::QuasiZeros{T,N}) where {T,N} = QuasiZeros{real(T)}(r.axes)
broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(real), r::QuasiOnes{T,N}) where {T,N} = QuasiOnes{real(T)}(r.axes)
broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(imag), r::QuasiZeros{T,N}) where {T,N} = QuasiZeros{real(T)}(r.axes)
broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof(imag), r::QuasiOnes{T,N}) where {T,N} = QuasiZeros{real(T)}(r.axes)

### Binary broadcasting

function broadcasted(::AbstractQuasiArrayStyle, op, a::AbstractQuasiFill, b::AbstractQuasiFill)
    val = op(getindex_value(a), getindex_value(b))
    return QuasiFill(val, broadcast_shape(axes(a), axes(b)))
end

_broadcasted_zeros(f, a, b) = QuasiZeros{Base.Broadcast.combine_eltypes(f, (a, b))}(broadcast_shape(axes(a), axes(b)))
_broadcasted_ones(f, a, b) = QuasiOnes{Base.Broadcast.combine_eltypes(f, (a, b))}(broadcast_shape(axes(a), axes(b)))
_broadcasted_nan(f, a, b) = QuasiFill(convert(Base.Broadcast.combine_eltypes(f, (a, b)), NaN), broadcast_shape(axes(a), axes(b)))

broadcasted(::AbstractQuasiArrayStyle, ::typeof(+), a::QuasiZeros, b::QuasiZeros) = _broadcasted_zeros(+, a, b)
broadcasted(::AbstractQuasiArrayStyle, ::typeof(+), a::QuasiOnes, b::QuasiZeros) = _broadcasted_ones(+, a, b)
broadcasted(::AbstractQuasiArrayStyle, ::typeof(+), a::QuasiZeros, b::QuasiOnes) = _broadcasted_ones(+, a, b)

broadcasted(::AbstractQuasiArrayStyle, ::typeof(-), a::QuasiZeros, b::QuasiZeros) = _broadcasted_zeros(-, a, b)
broadcasted(::AbstractQuasiArrayStyle, ::typeof(-), a::QuasiOnes, b::QuasiZeros) = _broadcasted_ones(-, a, b)
broadcasted(::AbstractQuasiArrayStyle, ::typeof(-), a::QuasiOnes, b::QuasiOnes) = _broadcasted_zeros(-, a, b)

broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(+), a::QuasiZeros{<:Any,1}, b::QuasiZeros{<:Any,1}) = _broadcasted_zeros(+, a, b)
broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(+), a::QuasiOnes{<:Any,1}, b::QuasiZeros{<:Any,1}) = _broadcasted_ones(+, a, b)
broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(+), a::QuasiZeros{<:Any,1}, b::QuasiOnes{<:Any,1}) = _broadcasted_ones(+, a, b)

broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(-), a::QuasiZeros{<:Any,1}, b::QuasiZeros{<:Any,1}) = _broadcasted_zeros(-, a, b)
broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(-), a::QuasiOnes{<:Any,1}, b::QuasiZeros{<:Any,1}) = _broadcasted_ones(-, a, b)

broadcasted(::AbstractQuasiArrayStyle, ::typeof(*), a::QuasiZeros, b::QuasiZeros) = _broadcasted_zeros(*, a, b)

# In following, need to restrict to <: Number as otherwise we cannot infer zero from type
# TODO: generalise to things like SVector
for op in (:*, :/)
    @eval begin
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros, b::QuasiOnes) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros, b::QuasiFill{<:Number}) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros, b::Number) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros, b::Base.Broadcast.Broadcasted) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros, b::AbstractQuasiArray{<:Number}) = _broadcasted_zeros($op, a, b)
    end
end

for op in (:*, :\)
    @eval begin
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiOnes, b::QuasiZeros) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiFill{<:Number}, b::QuasiZeros) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::Number, b::QuasiZeros) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::Base.Broadcast.Broadcasted, b::QuasiZeros) = _broadcasted_zeros($op, a, b)
        broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::AbstractQuasiArray{<:Number}, b::QuasiZeros) = _broadcasted_zeros($op, a, b)
    end
end

for op in (:*, :/, :\)
    @eval broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiOnes, b::QuasiOnes) = _broadcasted_ones($op, a, b)
end

for op in (:/, :\)
    @eval broadcasted(::AbstractQuasiArrayStyle, ::typeof($op), a::QuasiZeros{<:Number}, b::QuasiZeros{<:Number}) = _broadcasted_nan($op, a, b)
end

for op in (:+, -)
    @eval begin
        function broadcasted(::AbstractQuasiArrayStyle{N}, ::typeof($op), a::AbstractQuasiArray{T,N}, b::QuasiZeros{V,N}) where {T,V,N}
            broadcast_shape(axes(a), axes(b)) == axes(a) || throw(ArgumentError("Cannot broadcast $a and $b. Convert $b to a Vector first."))
            LinearAlgebra.copy_oftype(a, promote_type(T,V))
        end

        broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof($op), a::AbstractQuasiFill{T,1}, b::QuasiZeros{V,1}) where {T,V} =
            Base.invoke(broadcasted, Tuple{DefaultQuasiArrayStyle, typeof($op), AbstractQuasiFill, AbstractQuasiFill}, DefaultQuasiArrayStyle{1}(), $op, a, b)
    end
end

function broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(+), a::QuasiZeros{T,N}, b::AbstractQuasiArray{V,N}) where {T,V,N}
    broadcast_shape(axes(a), axes(b)) == axes(a) || throw(ArgumentError("Cannot broadcast $a and $b. Convert $a to a Vector first."))
    LinearAlgebra.copy_oftype(b, promote_type(T,V))
end

broadcasted(::AbstractQuasiArrayStyle{1}, ::typeof(+), a::QuasiZeros{V,1}, b::AbstractQuasiFill{T,1}) where {T,V} =
            Base.invoke(broadcasted, Tuple{DefaultQuasiArrayStyle, typeof(+), AbstractQuasiFill, AbstractQuasiFill}, DefaultQuasiArrayStyle{1}(), +, a, b)

# Need to prevent array-valued fills from broadcasting over entry
_broadcast_getindex_value(a::AbstractQuasiFill{<:Number}) = getindex_value(a)
_broadcast_getindex_value(a::AbstractQuasiFill) = Ref(getindex_value(a))

broadcasted(::AbstractQuasiArrayStyle{N}, op, r::AbstractQuasiFill{T,N}, x::Number) where {T,N} = QuasiFill(op(getindex_value(r),x), axes(r))
broadcasted(::AbstractQuasiArrayStyle{N}, op, x::Number, r::AbstractQuasiFill{T,N}) where {T,N} = QuasiFill(op(x, getindex_value(r)), axes(r))
broadcasted(::AbstractQuasiArrayStyle{N}, op, r::AbstractQuasiFill{T,N}, x::Ref) where {T,N} = QuasiFill(op(getindex_value(r),x[]), axes(r))
broadcasted(::AbstractQuasiArrayStyle{N}, op, x::Ref, r::AbstractQuasiFill{T,N}) where {T,N} = QuasiFill(op(x[], getindex_value(r)), axes(r))

###
# Mul
###

MemoryLayout(::Type{<:QuasiZeros}) = ZerosLayout()
MemoryLayout(::Type{<:QuasiOnes}) = OnesLayout()

_quasi_mul(M::Mul{ZerosLayout}, _) = QuasiZeros{eltype(M)}(axes(M))
_quasi_mul(M::Mul{QuasiArrayLayout,ZerosLayout}, _) = QuasiZeros{eltype(M)}(axes(M))
_quasi_mul(M::Mul{QuasiArrayLayout,ZerosLayout}, ::NTuple{N,OneTo{Int}}) where N = Zeros{eltype(M)}(axes(M))
fillzeros(::Type{T}, a::Tuple{AbstractQuasiVector,Vararg{Any}}) where T<:Number = QuasiZeros{T}(a)
fillzeros(::Type{T}, a::Tuple{Any,AbstractQuasiVector,Vararg{Any}}) where T<:Number = QuasiZeros{T}(a)

copy(M::MulAdd{<:AbstractFillLayout,<:AbstractFillLayout,<:AbstractFillLayout,<:Any,<:AbstractQuasiArray}) =
    QuasiFill(measure(axes(M.A,2))*M.α*getindex_value(M.A)*getindex_value(M.B) + M.β*getindex_value(M.C), axes(M))
copy(M::MulAdd{ZerosLayout,ZerosLayout,ZerosLayout,<:Any,<:AbstractQuasiArray}) =
    QuasiZeros{eltype(M)}(axes(M))
copy(M::MulAdd{ZerosLayout,<:AbstractFillLayout,ZerosLayout,<:Any,<:AbstractQuasiArray}) =
    QuasiZeros{eltype(M)}(axes(M))
copy(M::MulAdd{<:AbstractFillLayout,ZerosLayout,ZerosLayout,<:Any,<:AbstractQuasiArray}) =
    QuasiZeros{eltype(M)}(axes(M))

inv(D::QuasiEye) = D

zero(x::AbstractQuasiArray{T}) where {T} = QuasiZeros{T}(axes(x))
one(x::AbstractQuasiArray{T}) where {T} = QuasiOnes{T}(axes(x))

zeros(::Type{T}, x::Inclusion, y::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiZeros{T}((x, y...))
zeros(x::Inclusion, y::Union{OneTo,IdentityUnitRange,Inclusion}...) = zeros(Float64, x, y...)
zeros(::Type{T}, x::Union{OneTo,IdentityUnitRange}, y::Inclusion, z::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiZeros{T}((x, y, z...))
zeros(x::Union{OneTo,IdentityUnitRange}, y::Inclusion, z::Union{OneTo,IdentityUnitRange,Inclusion}...) = zeros(Float64, x, y, z...)
ones(::Type{T}, x::Inclusion, y::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiOnes{T}((x, y...))
ones(x::Inclusion, y::Union{OneTo,IdentityUnitRange,Inclusion}...) = ones(Float64, x, y...)
ones(::Type{T}, x::Union{OneTo,IdentityUnitRange}, y::Inclusion, z::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiOnes{T}((x, y, z...))
ones(x::Union{OneTo,IdentityUnitRange}, y::Inclusion, z::Union{OneTo,IdentityUnitRange,Inclusion}...) = ones(Float64, x, y, z...)
fill(c, x::Inclusion, y::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiFill(c, (x, y...))
fill(c, x::Union{OneTo,IdentityUnitRange}, y::Inclusion, z::Union{OneTo,IdentityUnitRange,Inclusion}...) where T = QuasiFill(c, (x, y, z...))

iszero(x::AbstractQuasiFill) = iszero(getindex_value(x))
isone(x::AbstractQuasiFill) = isone(getindex_value(x))


text:
using MacroTools: postwalk, @q

_padtuple(x, n) = length(x) < n ? (x..., ones(n - length(x))...) : x
isconstant(x) = occursin(r"^-?\d+$", string(x)) || occursin(r"^-?\d+\.\d+$", string(x))

function extractrtinfo!(m::Module, netlist::Netlist, innames, outname, f, isbroadcast, submod, args...)
    # evaluate function
    retval = isbroadcast ? f.(args...) : f(args...)

    # get types
    intypes = Symbol.(typeof.(args))
    outtype = Symbol(typeof(retval))

    # get sizes
    insizes = _padtuple.(size.(args), 2)
    outsize = _padtuple(size(retval), 2)

    # push types onto module and netlist
    inputs = Variable[]
    outputs = ismissing(outname) ? Variable[] : [Variable(outname, outtype)]
    for (inname, intype, insize) in zip(innames, intypes, insizes)
        push!(inputs, Variable(inname, intype))
        if !contains(netlist, string(inname)) && all(insize .!= 0)
            if isconstant(inname)
                update!(netlist, Net(name = string(inname), class = :constant, size = insize))
            elseif haskey(m.parameters, inname)
                update!(netlist, Net(name = string(inname), class = :parameter, size = insize))
            else
                update!(netlist, Net(name = string(inname), class = :internal, size = insize))
            end
        end
    end

    # update operator if operation is submodule
    if !isnothing(submod)
        v = findnode(m.dfg, string.(innames), [string(outname)], submod)
        operator = Symbol(typeof(f))
        !isempty(v) && set_prop!(m.dfg, v[1], :operator, operator)
    else
        operator = Symbol(f)
    end
    updatetype!(m, inputs, outputs, operator)

    # update output in netlist
    if !ismissing(outname) && !contains(netlist, string(outname)) && all(outsize .!= 0)
        update!(netlist, Net(name = string(outname), class = :internal, size = outsize))
    end

    # allow handler to access runtime information
    op = Operation([intypes...], ismissing(outname) ? Symbol[] : [outtype], operator)
    if haskey(m.handlers, op)
        handler = m.handlers[op]
    else
        handler = gethandler(op)()
        m.handlers[op] = handler
    end
    extractrtinfo!(handler, innames, outname, submod, args, retval, f)

    return retval
end

function setinputs!(netlist::Netlist, innames)
    for inname in innames
        i = find(netlist, inname)
        if isnothing(i)
            error("Cannot set $inname as input because it does not exist in netlist.")
        else
            update!(netlist, Net(name = inname,
                                 type = netlist[i].type,
                                 class = :input,
                                 signed = netlist[i].signed,
                                 size = netlist[i].size))
        end
    end

    return netlist
end

function setoutputs!(netlist::Netlist, outnames)
    for outname in outnames
        i = find(netlist, outname)
        if isnothing(i)
            error("Cannot set $outname as output because it does not exist in netlist.")
        else
            update!(netlist, Net(name = inname,
                                 type = netlist[i].type,
                                 class = :output,
                                 signed = netlist[i].signed,
                                 size = netlist[i].size))
        end
    end

    return netlist
end

function getparameters!(parametermap, params)
    fields = Expr[]
    for param in params
        @capture(param, var_::T_ => val_) ||
            error("Cannot parse parameter: $param (use syntax 'paramname::paramtype => paramvalue')")
        parametermap[var] = val
        push!(fields, @q($var::$T = $val))
    end

    return parametermap, fields
end

function getsubmodules!(submodulemap, submodules)
    fields = Expr[]
    for submodule in submodules
        @capture(submodule, var_::T_) ||
            error("Cannot parse submodule: $submodule (use syntax 'modname::modtype')")
        submodulemap[var] = Symbol(T)
        push!(fields, @q($var::$T))
    end

    return submodulemap, fields
end

createinnerconstructor(name, initdef, fields) = isnothing(initdef) ? nothing :
    @q begin
        function $name($(fields...))
            $initdef
        end
    end

createstruct(name, pfields, sfields, constructor) = isnothing(constructor) ?
    esc(@q begin
        Base.@kwdef struct $name
            $(pfields...)
            $(sfields...)
        end
    end) :
    esc(@q begin
        Base.@kwdef struct $name
            $(pfields...)
            $(sfields...)

            $constructor
        end
    end)

function stripbroadcast(x)
    m = match(r"^\.(.*)$", string(x))

    return isnothing(m) ? (false, x) : (true, Symbol(m.captures[1]))
end
stripparameters(x; prefix) = @capture(x, p_.s_) && p == prefix ? s : x
stripargument(x) = @capture(x, name_::T_) ? name : error("Cannot strip argument expression $x")

issymbol(s) = false
issymbol(s::Symbol) = true

function parsenode!(m::Module, opd, op, args, prefix, modcalls, modsym, netsym)
    inputs = map(x -> Variable(x, :Any), args)
    outputs = ismissing(opd) ? Variable[] : map(x -> Variable(x, :Any), [opd])
    isbroadcast, op = stripbroadcast(op)
    op = opalias(op)

    opname = stripparameters(op; prefix = prefix)
    addnode!(m, inputs, outputs, opname)

    submod = (opname == op) ? nothing : opname

    return @q BitSAD.HW.extractrtinfo!($modsym, $netsym, $args, $(QuoteNode(:($opd))), $op, $isbroadcast, $(QuoteNode(:($submod))), $(modcalls...))
end

function parseexpr!(m::Module, expr; prefix, modsym, netsym, opd = nothing, level = 1, depth = 1, counter = 1)
    opd = isnothing(opd) ? Symbol("net_$(level)_$(depth)_$counter") : opd

    if issymbol(stripparameters(expr; prefix = prefix))
        return stripparameters(expr; prefix = prefix), expr
    else
        @capture(expr, op_(args__)) || error("Cannot parse statement $expr.")
    end

    wires = Symbol[]
    modcalls = Union{Expr, Symbol}[]
    for (i, arg) in enumerate(args)
        wire, modcall = parseexpr!(m, arg; prefix = prefix, modsym = modsym, netsym = netsym,
                                           level = level, depth = depth + 1, counter = i)
        push!(wires, wire)
        push!(modcalls, modcall)
    end

    return opd, parsenode!(m, opd, op, wires, prefix, modcalls, modsym, netsym)
end

function parsecircuit!(m::Module, f)
    (@capture(f, (dut_::T_)(args__) -> body_) && T == m.name) ||
        error("Cannot parse circuit description.")

    # original function definition
    funcdef = @q begin
        function ($dut::$T)($(args...))
            $body
        end
    end

    # separate body into statements
    # body = postwalk(x -> stripparameters(x; prefix = dut), body)
    @capture(body, begin statements__ end) || error("Cannot parse circuit description.")

    # gensym for modified function definition
    modsym = gensym("m")
    netsym = gensym("netlist")

    # step through body of function and extract DFG
    modstatements = Expr[]
    rargs = Symbol[]
    for (i, statement) in enumerate(statements)
        @capture(statement, opd_ = rhs_) ?
            push!(modstatements,
                @q $opd = $(parseexpr!(m, rhs; prefix = dut, modsym = modsym, netsym = netsym, opd = opd, level = i)[2])) :
        @capture(statement, op_(arg__)) ?
            push!(modstatements,
                parseexpr!(m, statement; prefix = dut, modsym = modsym, netsym = netsym, opd = missing, level = i)[2]) :
        @capture(statement, return (rvals__)) ? vcat(rargs, rvals) :
        @capture(statement, return rval_) ? vcat(rargs, [rval]) :
        error("Cannot parse statement $statement.")
    end

    # create modified function definition (for type extraction at runtime)
    innames = string.(stripargument.(args))
    outnames = string.(Symbol.(rargs))
    modfdef = @q begin
        function extractrtinfo!($dut::$T, $modsym::BitSAD.HW.Module, $netsym::BitSAD.HW.Netlist, $(args...))
            $(modstatements...)
            BitSAD.HW.setinputs!($netsym, $innames)
            BitSAD.HW.setoutputs!($netsym, $outnames)
        end
    end

    return esc(funcdef), esc(modfdef)
end

"""
    @circuit

Create a new BitSAD circuit by specifying the name, parameters,
submodules, initialization, and implementation algorithm.

Returns a `Tuple{Module, Function}`. See [`Module`](@ref).
The function returned extracts runtime information to store in the `Module`.

The function has the signature where `T` is the name of the circuit generated:
```julia
extractrtinfo!(dut::T, m::BitSAD.HW.Module, netlist::BitSAD.HW.Netlist, args...)
```

# Examples
```julia
m = @circuit Foo begin
    parameters : [
        a => 0.125
        b => 10
    ]

    submodules : [
        m::SDM,
        n::CircularBuffer{DBit}
    ]

    initialize : begin
        dut = new(a, b, m, n)
        fill!(dut.n, zero(DBit))

        return dut
    end

    circuit : (dut::Foo)(x::DBit) -> begin
        y = dut.a * x
        z = y + dut.b

        return z
    end
end
```
"""
macro circuit(name, body)
    parametermap = Dict{Symbol, Number}()
    submodulemap = Dict{Symbol, Symbol}()
    pfields = Expr[]
    sfields = Expr[]
    funcdef = body
    modfdef = body
    initdef = nothing
    m = Module(name = Symbol("$name"))

    # walk through definition
    postwalk(body) do statement
        # extract parameters
        if @capture(statement, parameters : [params__])
            m.parameters, pfields = getparameters!(parametermap, params)
        end

        # extract submodules
        if @capture(statement, submodules : [submods__])
            m.submodules, sfields = getsubmodules!(submodulemap, submods)
        end

        # extract initialization
        if @capture(statement, initialize : initroutine_)
            initdef = rmlines(initroutine)
        end

        # extract circuit
        if @capture(statement, circuit : f_)
            funcdef, modfdef = parsecircuit!(m, f)
        end

        return statement
    end

    innerconstructor = createinnerconstructor(name, initdef, vcat(collect(keys(parametermap)), collect(keys(submodulemap))))
    structdef = createstruct(name, pfields, sfields, innerconstructor)

    return @q begin
        $structdef

        $funcdef

        $m, $modfdef
    end
end


text:
@testset "chainrules" begin
    x, Δx, x̄ = randn(3)
    y, Δy, ȳ = randn(3)
    z, Δz, z̄ = randn(3)
    Δu = randn()

    ỹ = expm1(y)
    frule_test(Bijectors.find_alpha, (x, Δx), (ỹ, Δy), (z, Δz); rtol=1e-3, atol=1e-3)
    rrule_test(Bijectors.find_alpha, Δu, (x, x̄), (ỹ, ȳ), (z, z̄); rtol=1e-3, atol=1e-3)
end


text:
# Create positive definite matrix
to_posdef(A::AbstractMatrix) = A * A' + I
to_posdef_diagonal(a::AbstractVector) = Diagonal(a .^ 2 .+ 1)

# Create vectors in probability simplex.
function to_simplex(x::AbstractArray)
    max = maximum(x; dims = 1)
    y = exp.(x .- max)
    y ./= sum(y; dims = 1)
    return y
end
to_simplex(x::AbstractArray{<:AbstractArray}) = to_simplex.(x)

# Utility for testing `adapt_randn`
function test_adapt_randn(rng, x::AbstractVector, ::Type{T}, dims::Int...) where {T}
    Random.seed!(rng, 100)
    y = DistributionsAD.adapt_randn(rng, x, dims...)
    @test y isa Array{T}
    @test size(y) == dims

    Random.seed!(rng, 100)
    @test y == randn(rng, T, dims...)
end


text:
using Dates
printstyled("  sync!\n", color=:light_green)

fs = 100.0
ns = 10
t = floor(Int64, time()-60.0)*sμ
δ3 = 1.0
dtμ = round(Int, 1.0e6/fs)

# for k = 1:10
# Test three seismic and three irregularly sampled channels
S = randSeisData(4, s=1.0)[2:4] + randSeisData(4, c=1.0)[2:4]

# seismic channels will be ns seconds of data sampled at 100 Hz with preset gaps
nx = round(Int, ns * fs)
for i = 1:3
  S.x[i] = rand(nx)       # Set S.x to be 1000 samples each
  S.fs[i] = fs            # Set S.fs to be uniformly 100 Hz
end
S.t[1] = [1 t; nx 0]
S.t[2] = [1 t-ns*1000000; nx 0]
S.t[3] = [1 t; 2 round(Int, δ3*1.0e6); nx 0]

# irregular channels will be 100 points each, randomly sampled between the start and the end
n_irr = div(nx,10)
for i = 4:6
  S.x[i] = rand(n_irr)
  S.fs[i] = 0.0
  r = UnitRange{Int64}(0:ns*1000000)
  S.t[i] = hcat(zeros(Int64, n_irr), t.+sort(rand(r, n_irr)))

  # Ensure the last sample is always <1s from the end
  τ = view(S.t[i], :, 2)
  if last(τ)-first(τ) < (ns-1)*1.0e6
    τ[end] = first(τ) + (ns-0.5)*1.0e6
  end
end
sx = Lx(S)

# Set derived variables
ts, te = get_edge_times(S)
ts₀ = minimum(ts)
ts₁ = maximum(ts)

# TEST 1 =====================================================================
# Sync to last start; don't sync ends
printstyled(stdout,"    sync to last start; don't sync ends\n", color=:light_green)
T = sync(S)
basic_checks(T)
ts_new, te_new = get_edge_times(T)

# Expectation: T.t[2] is deleted
@test T.n == 5
@test S.id[1] == T.id[1]
@test findid(S.id[2], T.id) == 0
for i = 3:6
  @test(S.id[i] == T.id[i-1])
end

# Expectation: the latest time in S.t is now the earliest start time
@test minimum(ts_new .≥ ts₁) == true
@test maximum(ts_new .== ts₁) == true

# TEST 2 =====================================================================
# Sync to last start; don't sync ends

# change trace 2 to begin only 1s earlier; resync
S.t[2] = [1 t-1000000; nx 0]
ts, te = get_edge_times(S)
ts₀ = minimum(ts)
ts₁ = maximum(ts)

T = sync(S)
ts_new, te_new = get_edge_times(T)

# Expectation: T.x[2] is at least 100 pts shorter and exactly 100 shorter than T.x[1]
@test T.n == 6
tx = Lx(T)
@test tx[2] ≤ nx-100
@test tx[2] == tx[1] - 100
@test tx[2] ≤ tx[3]
@test minimum(ts_new .≥ ts₁) == true
@test maximum(ts_new .== ts₁) == true
basic_checks(T)

# TEST 3 =====================================================================
# Sync to last start; sync to first end
printstyled(stdout,"    sync to last start; sync to first end\n", color=:light_green)
te₀ = minimum(te)
te₁ = maximum(te)

T = sync(S, t="first")
basic_checks(T)
ts_new, te_new = get_edge_times(T)

# Expectation: T.x[1] and T.x[2] are same length, T.x[3] is shorter
tx = Lx(T)
@test tx[1] == tx[2]
@test tx[2] > tx[3]
@test minimum(te_new .≤ te₀) == true
@test maximum(te_new .== te₀) == true
@test minimum(sx.-tx.≥0) == true                # vx is always longer

# TEST 4 =====================================================================
# Sync to first start; sync to first end
printstyled(stdout,"    sync to first start; sync to last end\n", color=:light_green)
U = sync(S, s="first", t="first")
basic_checks(U)
ts_new, te_new = get_edge_times(U)

# Expectation:
ux = Lx(U)
@test minimum(ux.-tx.≥0) == true                # vx is always longer
@test tx[1]+100 ≤ ux[1]                        # U.x[1] gains >100 samples
@test minimum(ts_new .≥ ts₀) == true           # Should still be true
@test 2 ∈ findall((ts_new .== ts₀).== true)   # Defined earliest start time
@test ux[2]==ux[1]                            # Due to how we set the gap
@test ux[1]-ux[3] == 100
@test minimum(te_new .≤ te₀) == true
@test maximum(te_new .== te₀) == true
for i = 4:6
  @test ux[i] ≥ tx[i]
end

# TEST 5 =====================================================================
# Sync to first start; sync to last end
V = sync(S, s="first", t="last")
basic_checks(V)
ts_new, te_new = get_edge_times(V)
vx = Lx(V)

# Expectation:
@test minimum(vx.-sx.≥0) == true                # vx is always longer
@test vx[1]-sx[1] == 199
@test vx[1]-ux[1] == 199
@test minimum(ts_new .≥ ts₀) == true           # Should still be true
@test 2 ∈ findall((ts_new .== ts₀).== true)   # Defined earliest start time
@test vx[2] == vx[1]                           # Due to how we set the gap
@test vx[2] - vx[3] == 99
@test minimum(te_new .≤ te₁) == true
@test maximum(te_new .== te₁) == true
for i = 4:6
  @test vx[i] ≥ tx[i]
  @test vx[i] == sx[i]
end

# TEST 6 =====================================================================
# Sync to s = DateTime (first trace start), t = "none"
# trace 3 should be 100 samples shorter
# so should trace 2
printstyled(stdout,"    sync to DateTime; don't sync ends\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + round(Int, 1.0e6*nx/fs)
ds₆ = u2d(ts₆*1.0e-6)
de₆ =  u2d(te₆*1.0e-6)
W = sync(S, s=ds₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                      # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

# Repeat with an end time
W = sync(S, s=ds₆, t=de₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[3]-wx[3] == 99                       # Trace 3 is 99 samples shorter
for i in [1,2,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths; 2 gets padded
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

# TEST 7 =====================================================================
# Sync to DateTime 1s before first trace, 0.01s after; is it 101 pts longer?

# Repeat with an increased window range
ds₆ -= Second(1)
de₆ += Millisecond(10)
X = sync(S, s=ds₆, t=de₆)
basic_checks(X)
ts_new, te_new = get_edge_times(X)
xx = Lx(X)

# Expectations:
@test xx[1]-sx[1] == 101    # Trace 1 is 101 samples longer
@test xx[2]-sx[2] == 101    # Trace 2 is also 101 samples longer
@test xx[3]-sx[3] == 2      # Trace 3 loses 99 samples at the end, but gains 100 at end and 1 at start ... net gain +2
for i = 4:6
  @test xx[i] == sx[i]
end

# Sync to DateTime 2s after first trace, t="first"
ds₆ += Second(3)
X = sync(S, s=ds₆, t=de₆)
basic_checks(X)
ts_new, te_new = get_edge_times(X)
xx = Lx(X)

# Expectations:
@test xx[1] == 801          # Trace 1 should be 801 samples
@test xx[1] == xx[2]        # Should be the same due to padding
# In fact trace 2 should have 101 points appended
found = findall([occursin("appended 101", i) for i in X.notes[2]])
@test isempty(found) == false

# TEST 8 =====================================================================
# A few simple combinations; do these work?
printstyled(stdout,"    sync start to DateTime; sync first end\n", color=:light_green)
X = sync(S, s=ds₆, t="first"); basic_checks(X)

printstyled(stdout,"    sync start to DateTime; sync to last end\n", color=:light_green)
X = sync(S, s=ds₆, t="last"); basic_checks(X)

printstyled(stdout,"    sync start to string time; sync to last end\n", color=:light_green)
ss = string(ds₆)
Y = sync(S, s=ds₆, t="last", v=3); basic_checks(Y)

# Expect: X != Y due to notes, but all other fields equal
for f in datafields
  if f != :notes
    @test isequal(getfield(X,f), getfield(Y,f))
  end
end

# TEST 9 =====================================================================
# Do we actually prune campaign data when all times are out of range?
printstyled(stdout,"    prune all irregular data when all times are out of range\n", color=:light_green)

ss = string(ds₆)
Z = deepcopy(S)
t = deepcopy(Z.t[5])
t = hcat(t[:,1:1], vcat(0, diff(t[:,2:2], dims=1)))
Z.t[5] = deepcopy(t)

redirect_stdout(out) do
  sync!(Z, v=3); basic_checks(Z)
end

# Expect: Z[5] is gone
for i in [1,2,3,4,6]
  @test any(Z.id.==S.id[i])
end
@test (any(Z.id.==S.id[5]) == false)

# ===========================================================================
# method extenson to SeisChannel
printstyled(stdout,"    SeisChannel method extension\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + round(Int, 1.0e6*nx/fs)
ds₆ = u2d(ts₆*1.0e-6)
de₆ =  u2d(te₆*1.0e-6)
C = deepcopy(S[1])
sync!(C, s=ds₆)
W = SeisData(C)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Repeat with an end time
C = deepcopy(S[1])
W = SeisData(sync(C, s=ds₆, t=de₆))
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# ===========================================================================
# method extenson to SeisEvent
printstyled(stdout,"    SeisEvent method extension\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + round(Int, 1.0e6*nx/fs)
ds₆ = u2d(ts₆*1.0e-6)
de₆ =  u2d(te₆*1.0e-6)
Ev = SeisEvent(hdr = randSeisHdr(), data = deepcopy(S))
sync!(Ev.data, s=ds₆)
W = Ev.data
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                      # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

Ev = SeisEvent(hdr = randSeisHdr(), data = deepcopy(S))
W = sync(Ev.data, s=ds₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                      # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time


text:
module BenchParticle

using Random
using BenchmarkTools
using HOODESolver

suite = BenchmarkGroup()

function fct(u, p, t)
    s1, c1 = sincos(u[1]/2)
    s2, c2 = sincos(u[2])
    s3, c3 = sincos(u[3])
    return [0, 0, u[6], c1*s2*s3/2, s1*c2*s3, s1*s2*c3]
end

Random.seed!(561909)
u0 = 2rand(BigFloat,6)-ones(BigFloat,6)
epsilon=big"1e-7"
ordmax=9
A = [0 0 0  1 0 0;
     0 0 0  0 1 0;
     0 0 0  0 0 0;
     0 0 0  0 1 0;
     0 0 0 -1 0 0;
     0 0 0  0 0 0]

t_0 = big"0.0"
t_max = big"1.0"
prob = HOODEProblem(fct, u0, (t_0, t_max), missing, A, epsilon)
nb = 10^3
n_tau = 32
suite["solve"] = @benchmarkable solve(prob, nb_t=nb, order=ordmax,
                                      getprecision=false, nb_tau=n_tau, dense=false)

end

BenchParticle.suite


text:
using Test, MPI

using CLIMA
using CLIMA.MPIStateArrays

MPI.Initialized() || MPI.Init()

const mpicomm = MPI.COMM_WORLD

@static if haspkg("CuArrays")
  using CuArrays
  # make sure that broadcasting is not being done by scalar indexing into CuArrays
  CuArrays.allowscalar(false)
  ArrayType = CuArray
else
  ArrayType = Array
end

@testset "MPIStateArray basics" begin
  Q = MPIStateArray{Float32}(mpicomm, ArrayType, 4, 6, 8)

  @test eltype(Q) == Float32
  @test size(Q) == (4, 6, 8)

  fillval = 0.5f0
  fill!(Q, fillval)

  haspkg("CuArrays") && CuArrays.allowscalar(true)
  @test Q[1] ==  fillval
  @test Q[2, 3, 4] == fillval
  @test Q[end] == fillval

  @test Array(Q) == fill(fillval, 4, 6, 8)

  Q[2, 3, 4] = 2fillval
  @test Q[2, 3, 4] != fillval
  haspkg("CuArrays") && CuArrays.allowscalar(false)

  Qp = copy(Q)

  @test typeof(Qp) == typeof(Q)
  @test eltype(Qp) == eltype(Q)
  @test size(Qp) == size(Q)
  @test Array(Qp) == Array(Q)

  Qp = similar(Q)

  @test typeof(Qp) == typeof(Q)
  @test eltype(Qp) == eltype(Q)
  @test size(Qp) == size(Q)

  copyto!(Qp, Q)
  @test Array(Qp) == Array(Q)
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using CompilerTools
using Base.Test

## Tests for CompilerTools.LivenessAnalysis

function test_liveness_1(x::Int64, y::Int64, s::Int64)
    C1 = (x - 256.0) / s + 0.6
    C2 = (y - 256.0) / s + 0.8
    acc1 = 0.0
    acc2 = 0.0
    I2 = 0.0
    R2 = 0.0
    n = 0
    while (R2+I2 < 1.0) && (n < 512)
       acc1 = (acc2+acc2)*acc1+C1;
       acc2 = R2-I2+C2;
       R2 = acc2*acc2;
       I2 = acc1*acc1;
       n += 1;
    end
    return n;
end

ast_lv_1 = CompilerTools.Helper.LambdaInfo(test_liveness_1, (Int64,Int64,Int64), code_typed(test_liveness_1, (Int64,Int64,Int64))[1])
# cfg_lv_1 = CompilerTools.CFGs.from_ast(ast_lv_1) :: CompilerTools.CFGs.CFG

linfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast_lv_1)

#CompilerTools.LivenessAnalysis.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)
lives_1 = CompilerTools.LivenessAnalysis.from_lambda(ast_lv_1)

# The AST changes too frequently to be testing anything but perhaps the starting block and even that is suspect.
@test length(lives_1.basic_blocks[lives_1.cfg.basic_blocks[-1]].live_in) == 3
@test length(lives_1.basic_blocks[lives_1.cfg.basic_blocks[-1]].live_out) == 7
@test length(lives_1.basic_blocks[lives_1.cfg.basic_blocks[-1]].def) == 7
@test length(lives_1.basic_blocks[lives_1.cfg.basic_blocks[-1]].use) == 3

deps = CompilerTools.TransitiveDependence.computeDependencies(lives_1)

#for entry in deps
#    vd = CompilerTools.LambdaHandling.getVarDef(entry[1], linfo)
#    println(entry[1], " ", entry[2], " ", vd.name)
#end


text:
import TPT

using Base.Test

using Plots; pyplot()

println("--- Unit WCA ---")

#
# Single component: Liquid Na
# Ref: I. H. Umar et al.: J. Phys. F: Metal Phys., Vol. 4 (1974), 1691-1706.
#
σ_Na = 6.235

ahs = TPT.AHS(σ = 6.235, ρ = 3.564e-3)
wca = TPT.WCA(ahs, 373)
pseudo = TPT.Ashcroft(1.0, 1.67)
nfe = TPT.NFE(wca, pseudo)
sys = TPT.TPTSystem(wca, nfe, m = 4.191e+4)

residue = sys.ref.residue
σ_estim = TPT.hsdiameter_estimate(nfe, 373.)[1,1]
σ_wca = TPT.hsdiameter(sys.ref)[1,1]

B = TPT.blipfunction(sys.ref)[1,1]
plot(r -> B(r)*r^2, 4, 8, ylims=:auto)

g_ahs = TPT.paircorrelation(ahs)[1,1]
g_wca = TPT.paircorrelation(sys)[1,1]

plot([g_ahs, g_wca], 2, 20, labels = ["AHS" "WCA"], ylims=:auto)

S = TPT.entropy(sys)
U = TPT.internal(sys)
F = TPT.helmholtz(sys)

@testset "Unit WCA" begin
  @test isapprox(residue, 0, atol=1e-4)
  @test isapprox(σ_wca, 6.17, atol=1e-2)
  @test isapprox(S, 7.56, atol=1e-2)
  @test isapprox(U, -0.237, atol=1e-3)
  @test isapprox(F, -0.244, atol=1e-3)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## from base/boot.jl:
#
# immutable String <: AbstractString
#     data::Vector{UInt8}
# end
#

## basic UTF-8 decoding & iteration ##

is_surrogate_lead(c::Unsigned) = ((c & ~0x003ff) == 0xd800)
is_surrogate_trail(c::Unsigned) = ((c & ~0x003ff) == 0xdc00)
is_surrogate_codeunit(c::Unsigned) = ((c & ~0x007ff) == 0xd800)
is_valid_continuation(c) = ((c & 0xc0) == 0x80)

const utf8_offset = [
    0x00000000, 0x00003080,
    0x000e2080, 0x03c82080,
    0xfa082080, 0x82082080,
]

const utf8_trailing = [
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,
]

## required core functionality ##

function endof(s::String)
    d = s.data
    i = length(d)
    @inbounds while i > 0 && is_valid_continuation(d[i])
        i -= 1
    end
    i
end

function length(s::String)
    d = s.data
    cnum = 0
    for i = 1:length(d)
        @inbounds cnum += !is_valid_continuation(d[i])
    end
    cnum
end

@noinline function slow_utf8_next(d::Vector{UInt8}, b::UInt8, i::Int)
    # potentially faster version
    # d = s.data
    # a::UInt32 = d[i]
    # if a < 0x80; return Char(a); end
    # #if a&0xc0==0x80; return '\ufffd'; end
    # b::UInt32 = a<<6 + d[i+1]
    # if a < 0xe0; return Char(b - 0x00003080); end
    # c::UInt32 = b<<6 + d[i+2]
    # if a < 0xf0; return Char(c - 0x000e2080); end
    # return Char(c<<6 + d[i+3] - 0x03c82080)

    if is_valid_continuation(b)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    trailing = utf8_trailing[b + 1]
    if length(d) < i + trailing
        return '\ufffd', i+1
    end
    c::UInt32 = 0
    for j = 1:(trailing + 1)
        c <<= 6
        c += d[i]
        i += 1
    end
    c -= utf8_offset[trailing + 1]
    return Char(c), i
end

@inline function next(s::String, i::Int)
    # function is split into this critical fast-path
    # for pure ascii data, such as parsing numbers,
    # and a longer function that can handle any utf8 data
    d = s.data
    b = d[i]
    if b < 0x80
        return Char(b), i + 1
    end
    return slow_utf8_next(d, b, i)
end

function first_utf8_byte(ch::Char)
    c = UInt32(ch)
    b = c < 0x80    ? c%UInt8 :
        c < 0x800   ? ((c>>6)  | 0xc0)%UInt8 :
        c < 0x10000 ? ((c>>12) | 0xe0)%UInt8 :
                      ((c>>18) | 0xf0)%UInt8
    return b
end

function reverseind(s::String, i::Integer)
    j = length(s.data) + 1 - i
    d = s.data
    while is_valid_continuation(d[j])
        j -= 1
    end
    return j
end

## overload methods for efficiency ##

sizeof(s::String) = sizeof(s.data)

isvalid(s::String, i::Integer) =
    (1 <= i <= endof(s.data)) && !is_valid_continuation(s.data[i])

const empty_utf8 = String(UInt8[])

function getindex(s::String, r::UnitRange{Int})
    isempty(r) && return empty_utf8
    i, j = first(r), last(r)
    d = s.data
    if i < 1 || i > length(s.data)
        throw(BoundsError(s, i))
    end
    if is_valid_continuation(d[i])
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    if j > length(d)
        throw(BoundsError())
    end
    j = nextind(s,j)-1
    String(d[i:j])
end

function search(s::String, c::Char, i::Integer)
    if i < 1 || i > sizeof(s)
        i == sizeof(s) + 1 && return 0
        throw(BoundsError(s, i))
    end
    d = s.data
    if is_valid_continuation(d[i])
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, d[i]))
    end
    c < Char(0x80) && return search(d, c%UInt8, i)
    while true
        i = search(d, first_utf8_byte(c), i)
        (i==0 || s[i] == c) && return i
        i = next(s,i)[2]
    end
end

function rsearch(s::String, c::Char, i::Integer)
    c < Char(0x80) && return rsearch(s.data, c%UInt8, i)
    b = first_utf8_byte(c)
    while true
        i = rsearch(s.data, b, i)
        (i==0 || s[i] == c) && return i
        i = prevind(s,i)
    end
end

function string(a::String...)
    if length(a) == 1
        return a[1]::String
    end
    # ^^ at least one must be UTF-8 or the ASCII-only method would get called
    data = Array{UInt8}(0)
    for d in a
        append!(data,d.data)
    end
    String(data)
end

function string(a::Union{String,Char}...)
    s = Array{UInt8}(0)
    for d in a
        if isa(d,Char)
            c = UInt32(d::Char)
            if c < 0x80
                push!(s, c%UInt8)
            elseif c < 0x800
                push!(s, (( c >> 6          ) | 0xC0)%UInt8)
                push!(s, (( c        & 0x3F ) | 0x80)%UInt8)
            elseif c < 0x10000
                push!(s, (( c >> 12         ) | 0xE0)%UInt8)
                push!(s, (((c >> 6)  & 0x3F ) | 0x80)%UInt8)
                push!(s, (( c        & 0x3F ) | 0x80)%UInt8)
            elseif c < 0x110000
                push!(s, (( c >> 18         ) | 0xF0)%UInt8)
                push!(s, (((c >> 12) & 0x3F ) | 0x80)%UInt8)
                push!(s, (((c >> 6)  & 0x3F ) | 0x80)%UInt8)
                push!(s, (( c        & 0x3F ) | 0x80)%UInt8)
            else
                # '\ufffd'
                push!(s, 0xef); push!(s, 0xbf); push!(s, 0xbd)
            end
        else
            append!(s,(d::String).data)
        end
    end
    String(s)
end

function reverse(s::String)
    dat = s.data
    n = length(dat)
    n <= 1 && return s
    buf = Vector{UInt8}(n)
    out = n
    pos = 1
    @inbounds while out > 0
        ch = dat[pos]
        if ch > 0xdf
            if ch < 0xf0
                (out -= 3) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out + 1], buf[out + 2], buf[out + 3] = ch, dat[pos + 1], dat[pos + 2]
                pos += 3
            else
                (out -= 4) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out+1], buf[out+2], buf[out+3], buf[out+4] = ch, dat[pos+1], dat[pos+2], dat[pos+3]
                pos += 4
            end
        elseif ch > 0x7f
            (out -= 2) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
            buf[out + 1], buf[out + 2] = ch, dat[pos + 1]
            pos += 2
        else
            buf[out] = ch
            out -= 1
            pos += 1
        end
    end
    String(buf)
end

write(io::IO, s::String) = write(io, s.data)

pointer(x::String) = pointer(x.data)
pointer(x::String, i::Integer) = pointer(x.data)+(i-1)

convert(::Type{String}, s::String) = s
convert(::Type{String}, v::Vector{UInt8}) = String(v)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia