text:
#' # Modeling data 1
#'
#' Machine learning and data science is about modeling data. **Modeling** is the representation of an idea with some parameters and  a mathematical representation which we will encode in software. All machine learning methods are about training a computer to fit a model to some data. Even the fanciest neural networks are simply choices for models. In this notebook, we will begin to start building our first computational model of data.

#' ## Modeling data is hard!
#'
#' Let's pick up where we left off in notebook 1 with fruit. We were left with a riddle: when we load images of apples and bananas,

import Images
import Statistics
S = Statistics
# -

#' Load apple
apple = Images.load("data/10_100.jpg")
#' Load banana
banana = Images.load("data/104_100.jpg")
#' and then compare their average value for the color red, we end up with something that is perhaps surprising:

apple_red_amount = S.mean(Float64.(Images.red.(apple)))
banana_red_amount = S.mean(Float64.(Images.red.(banana)));

"The average value of red in the apple is $apple_red_amount, " *
"while the average value of red in the banana is $banana_red_amount."

#' We see that the banana's mean red value is higher than the apple's, even though the apple looks much redder. Can you guess why?
#'
#' There are actually two reasons. One of the reasons is the background: the image of the banana has a lot more white background than the apple, and that white background has a red value of 1! In our minds we ignore the background and say "the banana is bright yellow, the apple is dark red", but a computer just has a bundle of numbers and does not know where it should be looking.
#'
#' The other issue is that "bright yellow" isn't a color that exists in a computer. The computer has three colors: red, green, and blue. "Bright yellow" in a computer is a mixture of red and green, and it just so happens that to get this color yellow, it needs more red than the apple!

"The amount of red in the apple at (60, 60) is $(Float64(Images.red(apple[60, 60]))), " *
"while the amount of red in the banana at (60, 60) is $(Float64(Images.red(banana[60, 60])))."
#'
apple[60,60]
#'
banana[60,60]
#' This is a clear example that modeling data is hard!

#' ### A note on string interpolation
#'
#' In the last two input cells, we *interpolated a string*. This means that when we write the string using quotation marks (`"  "`), we insert a placeholder for some **value** we want the string to include. When the string is evaluated, the value we want the string to include replaces the placeholder. For example, in the following string,
#'
#' ```julia
#' mystring = "The average value of red in the apple is $apple_red_amount"
#' ```
#'
#' `$apple_red_amount` is a placeholder for the value stored in the variable `apple_red_amount`. Julia knows that we want to use the value bound to the variable `apple_red_amount` and *not* the word "apple_red_amount" because of the dollar sign, `$`, that comes before `apple_red_amount`.

#' #### Exercise 1
#'
#' Execute the following code to see what the dollar sign does:
#'
mypi = 3.14159
println("I have a variable called mypi that has a value of $mypi.")

#' #### Exercise 2
#'
#' Alter and execute the code that creates `mystring` below
#'
apple_blue_amount = S.mean(Float64.(Images.blue.(apple)))
mystring = "The average amount of blue in the apple is $apple_blue_amount"

#'
#' so that `println(mystring)` prints a string that reports the mean value of blue coloration in our image of an apple.

#' ## Take some time to think about the data
#'
#' Apples and bananas are very different, but how could we use the array of RGB values (which is how the images are represented in the computer, as we saw in notebook 1) to tell the difference between the two? Here are some quick ideas:
#'
#' - We could use the shape of the object in the image. But how can we encode ideas about shape from an array?
#' - We could use the size of the object in the image. But how do we calculate that size?
#' - We could use another color, or combinations of colors, from the image. Which colors?
#'
#' Let's go with the last route. The banana is yellow, which is a combination of red and green, while the apple is red. This means that the color that clearly differentiates between the two is not red, but green!

apple_green_amount = S.mean(Float64.(Images.green.(apple)))
banana_green_amount = S.mean(Float64.(Images.green.(banana)));

"The average value of green in the apple is $apple_green_amount, " *
"while the average value of green in the banana is $banana_green_amount."

#' The processes that we just went through are assigned fancy names: feature selection and data munging.
#'
#' **Feature selection** is the process of subsetting the data to a more relevant and informative set. We took the full image data and decided to select out the green channel.
#'
#' **Data munging** is transforming the data into a format more suitable for modeling. Here, instead of keeping the full green channel, we transformed it down to a single data point: the average amount of green.

#' ## Building a model
#'
#' We want to model the connection between "the average amount of green" and "is an apple or banana".
#' ![Drawing](data/data_flow.png)
#'
#' This model is a mathematical function which takes in our data and spits out a number that we will interpret as "is an apple" or "is a banana".
#' ![Drawing](data/what_is_model.png)
#'
#'
#' We will interpret the output of the function as "is an apple" if the output is close to 0, and "is a banana" if it's close to 1. Anything in the middle is something we are unsure about. Here we're using a mathematical function to perform a **classification**.
#'
#' Knowing how to declare and work with functions will allow us to model our data in the coming sections, so this is the subject of the next notebook!


text:
using Omega
# using Omega: samplemeanᵣ
using UnicodePlots
using Distributions

# function test()
#   μ1 = bernoulli(0.5)
#   μ2 = bernoulli(0.5)
#   μ = μ1 + μ2

#   c1 = rcd(normal(μ, 1.0), μ1)
#   nsamples = 100
#   means1 = [rand(samplemeanᵣ(c1, nsamples)) for i = 1:1000]
#   println(UnicodePlots.histogram(means1))

#   c2 = rcd(normal(μ, 1.0), μ2)
#   means2 = [rand(samplemeanᵣ(c2, nsamples)) for i = 1:1000]
#   println(UnicodePlots.histogram(means2))

#   c3 = rcd(normal(μ, 1.0), μ)
#   means3 = [rand(samplemeanᵣ(c3, nsamples)) for i = 1:1000]
#   println(UnicodePlots.histogram(means3))
# end

# test()

# function testbeta()
#   plotbeta(beta) = lineplot(i->Distributions.pdf(beta, i), 0.0001, 0.999)
#   α = uniform(0.001, 5.0)
#   β = uniform(0.001, 5.0)
#   b = betarv(α, β)
#   brcd = b ∥ₛ (α, β)
#   samples = rand((α, β), samplemeanᵣ(brcd) ==ₛ 0.5, SSMH; n = 100)
#   s = Distributions.Beta(rand(samples)...); plotbeta(s)
#   s = Distributions.Beta(rand(samples)...); plotbeta(s)
#   samples2 = rand((α, β), mean(brcd) ==ₛ α, SSMH; n = 100)
# end

# testbeta()

# # RCD here is not correct
# function testbad()
#   α = rademacher(0.5)
#   β = rademacher(0.5) * α
#   γ = rademacher(0.5) + α + β
#   samples = [rand(mean(γ ∥ (α, β))) for i = 1:1000]
# end

# testbad()

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
	sphinterpolate(cmd0::String="", arg1=[], kwargs...)

Spherical gridding in tension of data on a sphere

Full option list at [`sphinterpolate`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate .html)

Parameters
----------

- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = sphinterpolate(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate.html#g)
- **I** : **inc** : -- Str or Number --

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate.html#i)
- **Q** : **tension** : -- Number or Str --     Flags = mode[/options]

    Specify one of four ways to calculate tension factors to preserve local shape properties or satisfy arc constraints.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate .html#q)
- **T** : **var_tension** : -- Bool or Str --

    Use variable tension (ignored with -Q0 [constant]
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate.html#t)
- **Z** : **scale** : -- Bool or Str --

    Before interpolation, scale data by the maximum data range [no scaling].
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/sphinterpolate.html#z)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function sphinterpolate(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("sphinterpolate ", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'Q', d, [:Q :tension])
	cmd = add_opt(cmd, 'T', d, [:T :nodetable])
	cmd = add_opt(cmd, 'Z', d, [:Z :scale])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "sphinterpolate", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
sphinterpolate(arg1=[], cmd0::String=""; kw...) = sphinterpolate(cmd0, arg1; kw...)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Meshes

using Compat
using GeometryTypes
using Meshing
using MeshIO
using FileIO

export HomogenousMesh
export Mesh
export isosurface
export load, save #FileIO

typealias Mesh HomogenousMesh

#include("simplification.jl")

end # module Meshes


text:
using Random
using LightGraphs: SimpleGraph, is_connected

# TODO: this should probably be a concrete type <: AbstractMatrix that operates
# as an indexer
"""
    build_perm_matrix(p)

Builds a sparse permutation matrix form permutation `p` such that `P*x == x[p]`
"""
function build_perm_matrix(p)
    isnothing(p) && return I
    n = length(p)
    P = spzeros(n, n)
    for i in 1:n
        P[i, p[i]] = 1.0
    end
    return P
end

# Generates dual form:
#   min c^Tx st ΣF_ix_i + G ∈ PSDCONE()
#   Build from KKT conditions:
#       Fx = sum(F[i]*xstar[i] for i in 1:n) + G
#       all(eigvals(Matrix(Fx)) .>= 0)
#       all(eigvals(D) .>= 0)
#       all([0 .== c[i] - tr(F[i]*D) for i in 1:n])
#       tr(Fx*D) .== 0
"""
    generate_random_sdp(n; rand_seed=0)

Generates a random dual form SDP with side dimension `n`:
`min c'*x s.t. sum(F[i]*x[i]) + G ⪰ 0`

Returns `c, F, G, xstar, D`, where `xstar` and `D` are optimal primal and dual
variables respectively
"""
function generate_random_sdp(n; rand_seed=0)
    Random.seed!(rand_seed)

    D = diagm(1 .+ rand(n))
    F = Vector{SparseMatrixCSC{Float64, Int}}(undef, n)
    c = Vector{Float64}(undef, n)
    for i in 1:n
        F[i] = spzeros(n, n)
        block_size = randn() < 1.5 ? 2 : 10
        F[i][i:min(i+block_size,n), i:min(i+block_size,n)] .= 1
        c[i] = tr(D*F[i])
    end
    xstar = rand(n)
    Fx = sum(F[i]*xstar[i] for i in 1:n)
    G = -Fx

    return c, F, G, xstar, D
end

"""
    unzip(a)

Unzips a list of tuples `a`.

# Example
```julia-repl
julia> unzip([(1,2), (3,4), (5,6)])
([1, 3, 5], [2, 4, 6])
```
"""
unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))

"""
    is_separable(sp::SparseMatrixCSC)

Returns true if the sparsity pattern is separable (i.e., block diagonal).
"""
is_separable(sp::SparseMatrixCSC) = !is_connected(SimpleGraph(sp))


text:
export initialise, generator, plan, planner, policy, scan, objective, objectives

"""
    initialise()

Returns a populated model initialised and ready to run.
"""
function initialise(;
    griddims = (25, 25),
    experiment = Experiment(),
    municipalities = Dict("main" => (Governance(), 100)),
    lake_setup = lake_initial_state(Clear, Martin),
)
    @assert prod(griddims) >=
            sum(last(v) for v in values(municipalities)) + length(municipalities) "Total agents (municipalities + households) cannot be greater than the grid size"

    lake_state, lake_parameters = lake_setup
    prob =
        OrdinaryDiffEq.ODEProblem(lake_dynamics!, lake_state, (0.0, Inf), lake_parameters)

    space = GridSpace(griddims, periodic = false)
    properties = type2dict(experiment)
    merge!(properties, type2dict(Outcomes(); prefix = "outcomes"))
    push!(
        properties,
        :lake => OrdinaryDiffEq.init(prob, OrdinaryDiffEq.Tsit5()),
        :year => 0,
        :init_nutrients => lake_parameters.nutrients,
        :init_pike_mortality => lake_parameters.mp,
    )
    if properties[:nutrient_series] isa Noise && !isempty(properties[:nutrient_series].process.S₁.data)
        # Store the year in which the nutrient_series must transition to stable in the
        # event of a brownian bridge
        properties[:nutrient_stabilise] = Int(properties[:nutrient_series].process.S₁.data[1][1])
    end
    model = ABM(
        Union{Household,Municipality},
        space;
        properties = properties,
        scheduler = Schedulers.by_type((Household, Municipality), true),
        warn = false,
        rng = Random.RandomDevice()
    )

    total_houses = sum(last(m) for m in values(municipalities))
    real_estate_x = 1
    for (name, (gov, houses)) in municipalities
        #TODO: Do a proper voronoi partition, not this primitive segregation
        juristiction_x = Int(round(first(griddims) * (houses / total_houses)) - 1) # Separation of municipalities in x
        # Place Municipality headquarters in the middle of its jurisdiction
        municipality_id = nextid(model)
        municipality_pos = (
            Int(round((juristiction_x / 2) + real_estate_x)),
            Int(round(last(griddims) / 2)),
        )
        municipality = Municipality(
            municipality_id,
            municipality_pos,
            name,
            gov.information,
            gov.legislation,
            gov.regulate,
            gov.respond_direct,
            gov.threshold_variable,
            gov.interventions,
            gov.policies,
            gov.agents_uniform,
            gov.houseowner_type,
            gov.willingness_to_upgrade,
            gov.tolerance_level_affectors,
            gov.neighbor_distance,
            houses,
            0,
        )
        #TODO: This check needs to be updated now that knowledge is global.
        #for val in municipality.knowledge
        #    if val == :vegetation_imbalance
        #        @assert haskey(municipality.policies, Trawling) "Municipality must have a Trawling policy set to use :vegetation_imbalance"
        #    end
        #end
        add_agent_pos!(municipality, model)
        # Place houses within municipality jurisdiction
        for _ in 1:houses
            pos = (
                rand(model.rng, real_estate_x:min(real_estate_x + juristiction_x, first(griddims))),
                rand(model.rng, 1:last(griddims)),
            )
            while !isempty(pos, model) || pos == municipality_pos
                pos = (
                    rand(model.rng, real_estate_x:(real_estate_x+juristiction_x)),
                    rand(model.rng, 1:last(griddims)),
                )
            end
            compliance = gov.agents_uniform ? rand(model.rng) : gov.willingness_to_upgrade
            house =
                Household(nextid(model), pos, compliance, false, false, 0, municipality_id)
            add_agent_pos!(house, model)
        end
        real_estate_x += juristiction_x
    end
    return model
end

"""
    generator(model, replicates)

Creates `replicates` number of models with a different random seed. Useful for `ensemblerun!` calls.
"""
function generator(model, replicates)
    is_device = typeof(model.rng) <: Random.RandomDevice
    models = [deepcopy(model) for _ in 1:replicates]
    if is_device
        for m in models
            @set m.rng = Random.RandomDevice()
        end
    else
        for m in models
            seed!(m)
        end
    end
    return models
end

## Intervention planner

"""
    planner(plan(Angling))
    planner(plan(Planting; rate=5e-3),
            plan(Trawling, 1:3))

Provides a complete schedule of interventions for a [Municipality](@ref). Must
be used in conjunction with [`plan`](@ref).
"""
function planner(plan::Dict{Int,Vector{Intervention}}...)
    @assert length(unique(typeof(first(first(values(p)))) for p in plan)) == length(plan) "Cannot parse more than one plan per intervention."
    merge(vcat, plan...)
end

# This method generates a Dict{Integer, Vector{Intervention}} even though the value is
# always singular. The resultant `Dict` is not expected to be used as-is, but rather
# merged with other interventions to create a master plan.
"""
    plan(Angling) # Assume always on
    plan(Angling; rate = 2.5e-3) # Always on with custom rate
    plan(Angling, 7; rate = 3.2e-3) # Only one year (custom rate)
    plan(Angling, 3:5) # Only years 3 to 5
    plan(Angling, [(period = 1:4, ),
                   (year = 5, rate = 7.2e-3),
                   (period = 7:9, )]) # Active in years 1-5, 7-9 with a custom rate
                                      # in year 5

Helper that provides complex scheduling for interventions with a simple interface.

The return type is `Dict{Int, Vector{Intervention}}`, were the key is each year
the collection of interventions will be active. As a convention, year `-1` denotes an
'always active' intervention.

The result can be provided to `municipality.interventions`, although
this function should almost always used in conjunction with [`planner`](@ref).
"""
plan(::Type{I}; kwargs...) where {I<:Intervention} = plan(I, -1; kwargs...)

function plan(::Type{I}, year::Int; kwargs...) where {I<:Intervention}
    schedule = Dict{Int,Vector{Intervention}}()
    schedule[year] = [I(; kwargs...)]
    schedule
end

function plan(::Type{I}, period::UnitRange{Int}; kwargs...) where {I<:Intervention}
    schedule = Dict{Int,Vector{Intervention}}()
    for year in period
        schedule[year] = [I(; kwargs...)]
    end
    schedule
end

function plan(::Type{I}, values::Vector{<:NamedTuple}) where {I<:Intervention}
    @assert all(haskey(v, :year) || haskey(v, :period) for v in values) "Each entry must contain a `year` or `period`."
    schedule = Dict{Int,Vector{Intervention}}()
    filter = (year = 5, period = 1:7) # Dummy tuple which we can use to dump our temporal data
    for entry in values
        params = Base.structdiff(entry, filter)
        if haskey(entry, :year)
            schedule[entry.year] = [I(; params...)]
        else
            for year in entry.period
                schedule[year] = [I(; params...)]
            end
        end
    end
    schedule
end

"""
    policy(scan(Trawling), scan(Planting; rate = (1e-5, 2.6e-3)))

Enables the decision module to alter suggested planner values, optimising the system
state when possible. Used to set the `policies` of each [Municipality](@ref).
"""
policy(scans::Dict{Type{<:Intervention},NamedTuple}...) = merge(vcat, scans...)

"""
    scan(Trawling) # Activate decisions on trawling using default ranges.
    scan(Trawling; rate = (5e-3, 2e-2)) # Activate trawling decisions with custom
                                        # search range for `rate`.

Use the decision making optimiser to fine tune values in the planner.
`scan` expects a lower and upper bound of a range to scan for each intervention
property.
For the moment, this strategy excludes `WastewaterTreatment`, and does not
adjust active years. These must still be set in the planner.

Should be used in conjunction with [`policy`](@ref).
"""
function scan(::Type{I}; kwargs...) where {I<:Intervention}
    return Dict{Type{<:Intervention},NamedTuple}(I => default_policy(I; kwargs...))
end

default_policy(::Type{Planting}; rate = (0.0, 1e-2)) = (rate = rate,)
default_policy(::Type{Trawling}; rate = (0.0, 1e-2)) = (rate = rate,)
default_policy(::Type{Angling}; rate = (2.25e-3, 2.7e-3)) = (rate = rate,)

"""
    objectives(objective(min_time), objective(min_cost, 0.5))

Constructs a set of objectives which will be used to optimise policy towards the chosen
target.
"""
function objectives(objs::Tuple{Function,<:Real}...)
    fns, weights = collect(zip(objs...))
    weights = Float64.(weights ./ sum(weights))
    return Tuple(zip(fns, weights))
end

"""
    objective(min_time)
    objective(min_acceleration, 2)

Provides a weighted objective to be used when optimising policy. Should be used in
conjunction with [`objectives`](@ref).

When providing a weight, this value can be any positive `Real` number which will
rank this objective higher or lower than other objectives. It is not imperative
to make sure all numbers entered in the policy to add up to 100%, but bear in mind
that the values will be normalised for the optimisation.
"""
objective(obj::Function, weight::R = 1.0) where {R<:Real} = (obj, weight)

## Helpers

"""
    type2dict(struct; prefix = "")

Converts a `struct` into a `Dict`. Borrowed from DrWatson, but extended to allow some
prefix to be attached. This is helpful when merging multiple structs into one `Dict`.
"""
function type2dict(dt; prefix = "")
    di = Dict{Symbol,Any}()
    for n in propertynames(dt)
        id = isempty(prefix) ? n : Symbol(prefix, "_", n)
        di[id] = getproperty(dt, n)
    end
    di
end


text:
using Documenter
using Literate

using CuArrays

# generate tutorials
OUTPUT = joinpath(@__DIR__, "src/tutorials/generated")
Literate.markdown(joinpath(@__DIR__, "src/tutorials/intro.jl"), OUTPUT;
    repo_root_url="https://github.com/JuliaGPU/CuArrays.jl/blob/master")

makedocs(
    modules = [CuArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CuArrays.jl",
    pages = [
        "Home" => "index.md",
        "Tutorials"  => [
            "tutorials/generated/intro.md"
        ],
    ],
    doctest = true
)


text:
# Minimum segment length for ray tracing
const minimum_segment_length = 1e-4

# QuadraticSegment_2D
# ---------------------------------------------------------------------------------------------------
# Threshold beneath which the segment is treated as linear.
# Equal to 2|d⃗|θ, where d⃗ is the displacement of x⃗₃ from the midpoint of the linear segment (x⃗₁, x⃗₂)
# and θ is the angle between the intersecting line and d⃗
const QuadraticSegment_2D_1_intersection_ϵ = 5e-6

# UnstructuredMesh_2D
# ---------------------------------------------------------------------------------------------------
# Cell types are the same as VTK
const UnstructuredMesh_2D_linear_cell_types = UInt32[5, # Triangle
                                                     9  # Quadrilateral
                                                    ]
const UnstructuredMesh_2D_quadratic_cell_types = UInt32[22, # Triangle6
                                                        23  # Quadrilateral8
                                                       ]
const UnstructuredMesh_2D_cell_types = vcat(UnstructuredMesh_2D_linear_cell_types,
                                            UnstructuredMesh_2D_quadratic_cell_types)

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
Clips a plane by a tetrahedron.
The plane takes the form of a static row matrix
The tet takes the form of a static matrix
"""
function clip_plane_tet(plane::SMatrix{1,4,T,4}, tet::SMatrix{4,4,T,16}) where {T}
    function tet_mat_to_tuple(tet::SMatrix{4,4,T,16}) where {T}
        v_1 = SVector{3,T}(tet[1],  tet[2],  tet[3])
        v_2 = SVector{3,T}(tet[5],  tet[6],  tet[7])
        v_3 = SVector{3,T}(tet[9],  tet[10], tet[11])
        v_4 = SVector{3,T}(tet[13], tet[14], tet[15])
        return (v_1, v_2, v_3, v_4)
    end

    proj = plane * tet

    bool_neg = proj .< 0.0
    n_neg = sum(bool_neg)

    bool_pos = 0.0 .< proj
    n_pos = sum(bool_pos)

    if (n_pos == 0) || (n_neg == 0)
        return poly_eight{3,T}()
    else
        v = tet_mat_to_tuple(tet)
        if n_pos == 1
            bool_pos[1] && (return clip_plane_tet_1(proj, v))
            bool_pos[2] && (return clip_plane_tet_2(proj, v))
            bool_pos[3] && (return clip_plane_tet_3(proj, v))
            bool_pos[4] && (return clip_plane_tet_4(proj, v))
        elseif n_neg == 1
            bool_neg[1] && (return clip_plane_tet_1(proj, v))
            bool_neg[2] && (return clip_plane_tet_2(proj, v))
            bool_neg[3] && (return clip_plane_tet_3(proj, v))
            bool_neg[4] && (return clip_plane_tet_4(proj, v))
        else  # if (n_pos == 2) && (n_neg == 2)
            (bool_pos[1] == bool_pos[2]) && (return clip_plane_tet_12(proj, v))
            (bool_pos[1] == bool_pos[3]) && (return clip_plane_tet_13(proj, v))
            (bool_pos[1] == bool_pos[4]) && (return clip_plane_tet_14(proj, v))
        # else
        #     error("something is wrong")
        end
    end
end

function proj_weight_poly(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}, i_1::Int64, i_2::Int64) where {T}
    return weightPoly(v[i_1], v[i_2], proj[i_1], proj[i_2])
end

function clip_plane_tet_1(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 2, 1)
    b = proj_weight_poly(proj, v, 4, 1)
    c = proj_weight_poly(proj, v, 3, 1)
    tup = 0.0 < proj[1] ? (a, b, c) : (c, b, a)
    return poly_eight(tup)
end
function clip_plane_tet_2(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 1, 2)
    b = proj_weight_poly(proj, v, 3, 2)
    c = proj_weight_poly(proj, v, 4, 2)
    tup = 0.0 < proj[2] ? (a, b, c) : (c, b, a)
    return poly_eight(tup)
end
function clip_plane_tet_3(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 1, 3)
    b = proj_weight_poly(proj, v, 4, 3)
    c = proj_weight_poly(proj, v, 2, 3)
    tup = 0.0 < proj[3] ? (a, b, c) : (c, b, a)
    return poly_eight(tup)
end
function clip_plane_tet_4(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 1, 4)
    b = proj_weight_poly(proj, v, 2, 4)
    c = proj_weight_poly(proj, v, 3, 4)
    tup = 0.0 < proj[4] ? (a, b, c) : (c, b, a)
    return poly_eight(tup)
end

function clip_plane_tet_12(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 2, 3)
    b = proj_weight_poly(proj, v, 2, 4)
    c = proj_weight_poly(proj, v, 1, 4)
    d = proj_weight_poly(proj, v, 1, 3)
    tup = 0.0 < proj[1] ? (a, b, c, d) : (d, c, b, a)
    return poly_eight(tup)
end

function clip_plane_tet_13(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 1, 2)
    b = proj_weight_poly(proj, v, 1, 4)
    c = proj_weight_poly(proj, v, 3, 4)
    d = proj_weight_poly(proj, v, 3, 2)
    tup = 0.0 < proj[1] ? (a, b, c, d) : (d, c, b, a)
    return poly_eight(tup)
end

function clip_plane_tet_14(proj::SMatrix{1,4,T,4}, v::NTuple{4,SVector{3,T}}) where {T}
    a = proj_weight_poly(proj, v, 1, 3)
    b = proj_weight_poly(proj, v, 1, 2)
    c = proj_weight_poly(proj, v, 4, 2)
    d = proj_weight_poly(proj, v, 4, 3)
    tup = 0.0 < proj[1] ? (a, b, c, d) : (d, c, b, a)
    return poly_eight(tup)
end


text:
"""
	planar biped
		model inspired by Spring Flamingo
			from "Exploiting Inherent Robustness and Natural Dynamics in the Control of Bipedal Walking Robots"
	q = (x, z, t, q1, ..., q6)
		x - lateral position
		z - vertical position
		t - body orientation
		q1 - leg 1 thigh angle (absolute)
		q2 - leg 1 calf angle  (absolute)
		q3 - leg 1 foot angle  (absolute)
		q4 - leg 2 thigh angle (absolute)
		q5 - leg 2 calf angle  (absolute)
		q6 - leg 2 foot angle  (absolute)
"""
mutable struct Biped{T} <: Model{T}
    # dimensions
    nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
	nc::Int # contact points

    # torso
    l_torso::T
    lc_torso::T
    m_torso::T
    J_torso::T

    # leg 1
        # thigh
    l_thigh1::T
    lc_thigh1::T
    m_thigh1::T
    J_thigh1::T

        # calf
    l_calf1::T
    lc_calf1::T
    m_calf1::T
    J_calf1::T

		# foot
	lt_foot1::T # toe length
	lh_foot1::T # heel length
	m_foot1::T
	J_foot1::T

    # leg 2
        # thigh
    l_thigh2::T
    lc_thigh2::T
    m_thigh2::T
    J_thigh2::T

        # calf
    l_calf2::T
    lc_calf2::T
    m_calf2::T
    J_calf2::T

		# foot
	lt_foot2::T # toe length
	lh_foot2::T # heel length
	m_foot2::T
	J_foot2::T

	# joint friction
	friction_joint::Vector{T}

	# environment
	friction_body_world::Vector{T}
	friction_foot_world::Vector{T}
	gravity::T
end

function kinematics_body(model::Biped, q; mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if mode == :ee
		l = model.l_torso
	elseif mode == :com
		l = model.lc_torso
	elseif mode == :hip
		l = 0.0
	else
		@error "incorrect mode specified"
	end

	return [
			x - l * sin(t);
			z + l * cos(t);
		   ]
end

function kinematics_jacobian_body(model::Biped, q; mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if mode == :ee
		l = model.l_torso
	elseif mode == :com
		l = model.lc_torso
	elseif mode == :hip
		l = 0.0
	else
		@error "incorrect mode specified"
	end

	jac = zeros(eltype(q), 2, model.nq)

	jac[1, 1] = 1.0
	jac[1, 3] = -l * cos(t)
	jac[2, 2] = 1.0
	jac[2, 3] = -l * sin(t)

	return jac
end

function kinematics_thigh(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		le = model.l_thigh1
		lc = model.lc_thigh1
		idx = 4
	elseif leg == :leg2
		le = model.l_thigh2
		lc = model.lc_thigh2
		idx = 7
	else
		@error "incorrect leg specified"
	end

	if mode == :ee
		l = le
	elseif mode == :com
		l = lc
	else
		@error "incorrect mode specified"
	end

	return [
			x + l * sin(q[idx]);
			z - l * cos(q[idx]);
		   ]
end

function kinematics_jacobian_thigh(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		le = model.l_thigh1
		lc = model.lc_thigh1
		idx = 4
	elseif leg == :leg2
		le = model.l_thigh2
		lc = model.lc_thigh2
		idx = 7
	else
		@error "incorrect leg specified"
	end

	if mode == :ee
		l = le
	elseif mode == :com
		l = lc
	else
		@error "incorrect mode specified"
	end

	jac = zeros(eltype(q), 2, model.nq)
	jac[1, 1] = 1.0
	jac[1, idx] = l * cos(q[idx])
	jac[2, 2] = 1.0
	jac[2, idx] = l * sin(q[idx])

	return jac
end

function kinematics_calf(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		le = model.l_calf1
		lc = model.lc_calf1
		idx = 5
	elseif leg == :leg2
		le = model.l_calf2
		lc = model.lc_calf2
		idx = 8
	else
		@error "incorrect leg specified"
	end

	if mode == :ee
		l = le
	elseif mode == :com
		l = lc
	else
		@error "incorrect mode specified"
	end

	k_thigh = kinematics_thigh(model, q, leg=leg, mode=:ee)
	return k_thigh + [l * sin(q[idx]); -l * cos(q[idx])]
end

function kinematics_jacobian_calf(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		le = model.l_calf1
		lc = model.lc_calf1
		idx = 5
	elseif leg == :leg2
		le = model.l_calf2
		lc = model.lc_calf2
		idx = 8
	else
		@error "incorrect leg specified"
	end

	if mode == :ee
		l = le
	elseif mode == :com
		l = lc
	else
		@error "incorrect mode specified"
	end

	jac = kinematics_jacobian_thigh(model, q; leg=leg, mode=:ee)

	jac[1, idx] += l * cos(q[idx])
	jac[2, idx] += l * sin(q[idx])

	return jac
end

function kinematics_foot(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		lt = model.lt_foot1
		lh = model.lh_foot1
		idx = 6
	elseif leg == :leg2
		lt = model.lt_foot2
		lh = model.lh_foot2
		idx = 9
	else
		@error "incorrect leg specified"
	end

	if mode == :toe
		l = lt
	elseif mode == :heel
		l = -lh
	elseif mode == :com
		l = 0.0
	else
		@error "incorrect mode specified"
	end

	k_calf = kinematics_calf(model, q, leg=leg, mode=:ee)
	return k_calf + [l * cos(q[idx]); l * sin(q[idx])]
end

function kinematics_jacobian_foot(model::Biped, q; leg=:none, mode=:none)
	x = q[1]
	z = q[2]
	t = q[3]

	if leg == :leg1
		lt = model.lt_foot1
		lh = model.lh_foot1
		idx = 6
	elseif leg == :leg2
		lt = model.lt_foot2
		lh = model.lh_foot2
		idx = 9
	else
		@error "incorrect leg specified"
	end

	if mode == :toe
		l = lt
	elseif mode == :heel
		l = -lh
	elseif mode == :com
		l = 0.0
	else
		@error "incorrect mode specified"
	end

	jac = kinematics_jacobian_calf(model, q, leg=leg, mode=:ee)
	jac[1, idx] += -l * sin(q[idx])
	jac[2, idx] += l * cos(q[idx])

	return jac
end

function kinematics(model::Biped, q)
	# foot
	p_toe_1 = kinematics_foot(model, q, leg=:leg1, mode=:toe)
	p_heel_1 = kinematics_foot(model, q, leg=:leg1, mode=:heel)
	p_toe_2 = kinematics_foot(model, q, leg=:leg2, mode=:toe)
	p_heel_2 = kinematics_foot(model, q, leg=:leg2, mode=:heel)

	# knee
	p_knee_1 = kinematics_thigh(model, q, leg=:leg1, mode=:ee)
	p_knee_2 = kinematics_thigh(model, q, leg=:leg2, mode=:ee)

	# body
	p_body_1 = q[1:2]
	p_body_2 = kinematics_body(model, q, mode=:ee)

	[
	 p_toe_1;
	 p_heel_1;
	 p_toe_2;
	 p_heel_2;
	 p_knee_1;
	 p_knee_2;
	 p_body_1;
	 p_body_2;
	]
end

# Lagrangian
function lagrangian(model::Biped, q, q̇)
	L = 0.0

	# torso
	p_torso = kinematics_body(model, q, mode=:com)
	J_torso = kinematics_jacobian_body(model, q, mode=:com)
	v_torso = J_torso * q̇

	L += 0.5 * model.m_torso * transpose(v_torso) * v_torso
	L += 0.5 * model.J_torso * q̇[3]^2.0
	L -= model.m_torso * model.gravity * p_torso[2]

	# thigh 1
	p_thigh_1 = kinematics_thigh(model, q, leg=:leg1, mode=:com)
	J_thigh_1 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:com)
	v_thigh_1 = J_thigh_1 * q̇

	L += 0.5 * model.m_thigh1 * transpose(v_thigh_1) * v_thigh_1
	L += 0.5 * model.J_thigh1 * q̇[4]^2.0
	L -= model.m_thigh1 * model.gravity * p_thigh_1[2]

	# calf 1
	p_calf_1 = kinematics_calf(model, q, leg=:leg1, mode=:com)
	J_calf_1 = kinematics_jacobian_calf(model, q, leg=:leg1, mode=:com)
	v_calf_1 = J_calf_1 * q̇

	L += 0.5 * model.m_calf1 * transpose(v_calf_1) * v_calf_1
	L += 0.5 * model.J_calf1 * q̇[5]^2.0
	L -= model.m_calf1 * model.gravity * p_calf_1[2]

	# foot 1
	p_foot_1 = kinematics_foot(model, q, leg=:leg1, mode=:com)
	J_foot_1 = kinematics_jacobian_foot(model, q, leg=:leg1, mode=:com)
	v_foot_1 = J_foot_1 * q̇

	L += 0.5 * model.m_foot1 * transpose(v_foot_1) * v_foot_1
	L += 0.5 * model.J_foot1 * q̇[6]^2.0
	L -= model.m_foot1 * model.gravity * p_foot_1[2]

	# thigh 2
	p_thigh_2 = kinematics_thigh(model, q, leg=:leg2, mode=:com)
	J_thigh_2 = kinematics_jacobian_thigh(model, q, leg=:leg2, mode=:com)
	v_thigh_2 = J_thigh_2 * q̇

	L += 0.5 * model.m_thigh2 * transpose(v_thigh_2) * v_thigh_2
	L += 0.5 * model.J_thigh2 * q̇[7]^2.0
	L -= model.m_thigh2 * model.gravity * p_thigh_2[2]

	# calf 2
	p_calf_2 = kinematics_calf(model, q, leg=:leg2, mode=:com)
	J_calf_2 = kinematics_jacobian_calf(model, q, leg=:leg2, mode=:com)
	v_calf_2 = J_calf_2 * q̇

	L += 0.5 * model.m_calf2 * transpose(v_calf_2) * v_calf_2
	L += 0.5 * model.J_calf2 * q̇[8]^2.0
	L -= model.m_calf2 * model.gravity * p_calf_2[2]

	# foot 2
	p_foot_2 = kinematics_foot(model, q, leg=:leg2, mode=:com)
	J_foot_2 = kinematics_jacobian_foot(model, q, leg=:leg2, mode=:com)
	v_foot_2 = J_foot_2 * q̇

	L += 0.5 * model.m_foot2 * transpose(v_foot_2) * v_foot_2
	L += 0.5 * model.J_foot2 * q̇[9]^2.0
	L -= model.m_foot2 * model.gravity * p_foot_2[2]

	return L
end

function signed_distance(model::Biped, q)
	k = kinematics(model, q)

	return [k[2], k[4], k[6], k[8], k[10], k[12], k[14], k[16]]
end

function input_jacobian(::Biped, q)
	[0.0  0.0 -1.0  1.0  0.0  0.0  0.0  0.0  0.0;
	 0.0  0.0  0.0 -1.0  1.0  0.0  0.0  0.0  0.0;
	 0.0  0.0  0.0  0.0 -1.0  1.0  0.0  0.0  0.0;
	 0.0  0.0 -1.0  0.0  0.0  0.0  1.0  0.0  0.0;
	 0.0  0.0  0.0  0.0  0.0  0.0 -1.0  1.0  0.0;
	 0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  1.0]
end

function contact_jacobian(model::Biped, q)
	# foot
	J_toe_1 = kinematics_jacobian_foot(model, q, leg=:leg1, mode=:toe)
	J_heel_1 = kinematics_jacobian_foot(model, q, leg=:leg1, mode=:heel)
	J_toe_2 = kinematics_jacobian_foot(model, q, leg=:leg2, mode=:toe)
	J_heel_2 = kinematics_jacobian_foot(model, q, leg=:leg2, mode=:heel)

	# knee
	J_knee_1 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:ee)
	J_knee_2 = kinematics_jacobian_thigh(model, q, leg=:leg2, mode=:ee)

	# body
	J_body_1 = kinematics_jacobian_body(model, q, mode=:hip)
	J_body_2 = kinematics_jacobian_body(model, q, mode=:ee)

	[
	 J_toe_1;
	 J_heel_1;
	 J_toe_2;
	 J_heel_2;
	 J_knee_1;
	 J_knee_2;
	 J_body_1;
	 J_body_2;
	]
end

# nominal configuration
function nominal_configuration(model::Biped)
	[0.0; 1.0; 0.01 * π; -0.5 * π; 0.0 * π; -0.01 * π; 0.5 * π; 0.0 * π; -0.0 * π]
end

# friction_coefficients
function friction_coefficients(model::Biped)
	return [model.friction_foot_world; model.friction_body_world]
end

# dimensions
nq = 3 + 2 * 3            # generalized coordinates dimension
nu = 2 * 3                # control dimension
nc = 8                    # number of contact points
nw = 0                    # parameter dimension

# model parameters
m_torso = 12.0
m_thigh = 0.4598
m_calf = 0.306
m_foot = 0.3466

l_torso = 0.385
l_thigh = 0.42
l_calf = 0.45
l_toe = 0.1725 # distance ankle-toe
l_heel = 0.0525 # distance ankle-heel

lc_torso = 0.20
lc_thigh = 0.5 * l_thigh
lc_calf = 0.5 * l_calf

J_torso = 0.10
J_thigh = 0.01256
J_calf = 0.00952
J_foot = 0.0015

friction_joint = [0.1; 0.1; 0.1; 0.1; 0.1; 0.1] # q1-q6

# world parameters
friction_body_world = [0.5; 0.5; 0.5; 0.5] # knee1, knee2, body1, body2
friction_foot_world = [0.5; 0.5; 0.5; 0.5] # foot1, foot2, foot3, foot4
gravity = 9.81

biped = Biped(nq, nu, nw, nc,
			  l_torso, lc_torso, m_torso, J_torso,
			  l_thigh, lc_thigh, m_thigh, J_thigh,
			  l_calf, lc_calf, m_calf, J_calf,
			  l_toe, l_heel, m_foot, J_foot,
			  l_thigh, lc_thigh, m_thigh, J_thigh,
			  l_calf, lc_calf, m_calf, J_calf,
			  l_toe, l_heel, m_foot, J_foot,
			  friction_joint,
			  friction_body_world,
			  friction_foot_world,
			  gravity)

biped_contact_kinematics = [q -> kinematics_foot(biped, q, leg=:leg1, mode=:toe),
	q -> kinematics_foot(biped, q, leg=:leg1, mode=:heel),
	q -> kinematics_foot(biped, q, leg=:leg2, mode=:toe),
	q -> kinematics_foot(biped, q, leg=:leg2, mode=:heel),
	q -> kinematics_thigh(biped, q, leg=:leg1, mode=:ee),
	q -> kinematics_thigh(biped, q, leg=:leg2, mode=:ee),
	q -> q[1:2],
	q -> kinematics_body(biped, q, mode=:ee)]

biped_contact_kinematics_jacobians = [q -> kinematics_jacobian_foot(biped, q, leg=:leg1, mode=:toe),
	q -> kinematics_jacobian_foot(biped, q, leg=:leg1, mode=:heel),
	q -> kinematics_jacobian_foot(biped, q, leg=:leg2, mode=:toe),
	q -> kinematics_jacobian_foot(biped, q, leg=:leg2, mode=:heel),
	q -> kinematics_jacobian_thigh(biped, q, leg=:leg1, mode=:ee),
	q -> kinematics_jacobian_thigh(biped, q, leg=:leg2, mode=:ee),
	q -> kinematics_jacobian_body(biped, q, mode=:hip),
	q -> kinematics_jacobian_body(biped, q, mode=:ee)]

name(::Biped) = :biped


text:
using Test
using Sockets

include("../../src/JuliaBolt.jl")

using .JuliaBolt

# these tests are completed using a Neo4j Community Edition server
# running on the local machine with the default user "neo4j" and the
# password set to "password"

@testset "Connection Tests" begin

    @testset "test_connection_open_close" begin
        connection = bolt_connect("localhost", 7687, auth=("neo4j", "password"))
        @test isa(connection, Connection)
        bolt_close(connection)
    end

    @testset "test_connection_simple_run" begin
        records = []
        metadata = Dict()
        connection = bolt_connect("localhost", 7687, auth=("neo4j", "password"))
        bolt_run(connection, "RETURN 1", Dict(), on_success=(response, result)->merge!(metadata, result))
        bolt_pull_all(connection, on_success=(response, result)->merge!(metadata, result), on_records=(response, result)->append!(records, result))
        bolt_sync(connection)

        @test records == [[1]]

        bolt_close(connection)
    end

    @testset "test_multiple_chunk_response" begin
        b = ones(UInt8, 16365)
        records = []
        metadata = Dict()
        connection = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_run(connection, "CREATE (a) SET a.foo = \$x RETURN a", Dict("x"=>b), on_success=(_, result)->merge!(metadata, result))
        bolt_pull_all(connection, on_success=(_, result)->merge!(metadata, result), on_records=(_, result)->append!(records, result))
        bolt_sync(connection)

        foo = records[1][1][3]["foo"]

        @test b == foo

        bolt_close(connection)
    end

    @testset "test_return_1" begin
        records = []
        metadata = Dict()
        connection = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_run(connection, "RETURN \$x", Dict("x"=>1), on_success=(_, result)->merge!(metadata, result))
        bolt_pull_all(connection, on_success=(_, result)->merge!(metadata, result), on_records=(_, result)->append!(records, result))
        bolt_sync(connection)

        @test records == [[1]]

        bolt_close(connection)
    end

    @testset "test_return_1_in_tx" begin
        records = []
        metadata = Dict()
        connection = bolt_connect("localhost", 7687, auth=("neo4j", "password"))
        bolt_begin(connection, on_success=(_, result)->merge!(metadata, result))
        bolt_run(connection, "RETURN \$x", Dict("x"=>1), on_success=(_, result)->merge!(metadata, result))
        bolt_pull_all(connection, on_success=(_, result)->merge!(metadata, result), on_records=(_, result)->append!(records, result))
        bolt_commit(connection, on_success=(_, result)->merge!(metadata, result))
        bolt_sync(connection)

        @test records == [[1]]
        @test startswith(metadata["bookmark"], "neo4j:bookmark:")
        @test metadata["fields"] == ["\$x"]
        @test isa(metadata["t_first"], Integer)
        @test isa(metadata["t_last"], Integer)
        @test metadata["type"] == "r"

        bolt_close(connection)
    end

    @testset "test_begin_with_timeout" begin
        records = []
        metadata = Dict()
        cx1 = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_run(cx1, "CREATE (a:Node)")
        bolt_discard_all(cx1)
        bolt_sync(cx1)

        cx2 = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_begin(cx1, on_failure=(_, f)->println(f))
        bolt_run(cx1, "MATCH (a:Node) SET a.property = 1", on_failure=(_, f)->println(f))
        bolt_sync(cx1)

        bolt_begin(cx2, timeout=1, on_failure=(_, f)->println(f))
        bolt_run(cx2, "MATCH (a:Node) SET a.property = 2", on_failure=(_, f)->println(f))

        @test_throws TransientError bolt_sync(cx2)

        bolt_close(cx1)
        bolt_close(cx2)
    end

    @testset "test_run_with_timeout" begin
        cx1 = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_run(cx1, "CREATE (a:Node)")
        bolt_discard_all(cx1)
        bolt_sync(cx1)

        cx2 = bolt_connect("localhost", 7687, auth=("neo4j", "password"))

        bolt_begin(cx1, on_failure=(_, f)->println(f))
        bolt_run(cx1, "MATCH (a:Node) SET a.property = 1", on_failure=(_, f)->println(f))
        bolt_sync(cx1)

        bolt_run(cx2, "MATCH (a:Node) SET a.property = 2", timeout=1, on_failure=(_, f)->println(f))

        @test_throws TransientError bolt_sync(cx2)

        bolt_close(cx1)
        bolt_close(cx2)
    end
end

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
"""
FinEtools (C) 2017, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics.
"""
module myFinEtools

###########################################################################
# General facilities
###########################################################################
using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
# Exported: basic numerical types, type of data dictionary
export FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

using FinEtools.BoxModule: inbox, initbox!, updatebox!, boundingbox, inflatebox!, boxesoverlap
# Exported: methods for manipulating and testing boxes
export inbox, initbox!, updatebox!, boundingbox, inflatebox!, boxesoverlap

using FinEtools.PhysicalUnitModule: phun
# Exported: function for evaluating physical units
export phun

using FinEtools.RotationUtilModule: rotmat3!, skewmat!, cross3!, cross2
# Exported: functions for 3D rotation matrix computation, skew symmetric matrix computation,  in-place  cross product of 3-vectors, and  cross product of 2-vectors
export rotmat3!, skewmat!, cross3!, cross2

using FinEtools.CSysModule: CSys, updatecsmat!
# Exported: type  for coordinate systems, methods to invoke the update callback
export CSys, updatecsmat!

using FinEtools.FESetModule: FESet,  FESet0Manifold,  FESet1Manifold,  FESet2Manifold,  FESet3Manifold, manifdim, nodesperelem, count, fromarray!, connasarray, setlabel!, subset, cat, updateconn!, bfun, bfundpar, map2parametric, inparametric, centroidparametric,  FESetP1, FESetL2, FESetL3, FESetT3, FESetQ4, FESetQ9, FESetQ8, FESetT6, FESetH8, FESetH20, FESetH27, FESetT4, FESetT10
# Exported: types of finite element sets, abstract and concrete
export FESet,  FESet0Manifold,  FESet1Manifold,  FESet2Manifold,  FESet3Manifold, FESetP1, FESetL2, FESetL3, FESetT3, FESetQ4, FESetQ9, FESetQ8, FESetT6, FESetH8, FESetH20, FESetH27, FESetT4, FESetT10
# Exported: methods for accessing dimensions and counts
export manifdim, nodesperelem, count
# Exported: methods for  manipulating connectivity  and labels
export fromarray!, connasarray, setlabel!, subset, cat, updateconn!
# Exported: methods for computing basis function values and derivatives of basis functions  with respect to the parametric coordinates, and methods for working with parametric coordinates
export bfun, bfundpar, map2parametric, inparametric, centroidparametric

using FinEtools.FENodeSetModule: FENodeSet, spacedim, xyz3, count
# Exported: type for FE node sets, methods for accessing dimensions  and counts
export FENodeSet, spacedim, xyz3, count

using FinEtools.MeshSelectionModule: connectednodes, selectnode, selectelem, findunconnnodes
# Exported: functions to select (find) nodes and elements
export connectednodes, selectnode, selectelem, findunconnnodes

using FinEtools.MeshExportModule: vtkexportmesh
# Exported: VTK export
export vtkexportmesh

using FinEtools.MeshModificationModule: meshboundary,  fusenodes,  compactnodes,  mergemeshes, mergenmeshes, mergenodes,  renumberconn!,  meshsmoothing, mirrormesh, nodepartitioning, interior2boundary
# Exported: extraction of boundary, fusing of nodes  and merging of meshes, mesh smoothing,  node partitioning
export  meshboundary,  fusenodes,  compactnodes,  mergemeshes, mergenmeshes, mergenodes,  renumberconn!,  meshsmoothing, mirrormesh, nodepartitioning, interior2boundary

###########################################################################
# Mesh-generation functionality for various shapes
###########################################################################
using FinEtools.MeshQuadrilateralModule: Q4annulus, Q4quadrilateral, Q4elliphole, Q4block, Q4blockx, Q4refine, Q8block, Q4toQ8, Q8annulus, Q8blockx
# Exported: mesh generation functions for quadrilateral elements
export Q4annulus, Q4quadrilateral, Q4elliphole, Q4block, Q4blockx, Q4refine, Q8block, Q4toQ8, Q8annulus, Q8blockx

# My new sphere-generating function
include("Q4spheren.jl")
export Q4spheren

end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_issues_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_forks_repo_name: PetrKryslUCSD/FinEtoolsUseCase


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    n = 15
    x, y, color = rand(n), rand(n), rand(n)
    cmaps=[:cool,:viridis,:plasma,:inferno,:thermal,:leonardo,:winter,:spring,:ice]

    function FigGridHeat()
        fig = Figure(resolution = (1200, 800))
        axes = []
        c = 1
        for i in 1:2, j in 1:2:5
            ax=Axis(fig[i, j],aspect = 1,xgridvisible = false,ygridvisible = false)
            hmap = heatmap!(rand(10,10), colormap=cmaps[c], colorrange=(0, 1))
            ax.xticks = [1,10]
            ax.yticks = [1,10]
            ax.xticklabelsize = 20
            ax.yticklabelsize = 20
            cbar = Colorbar(fig,hmap,height = Relative(0.5),tickwidth = 2,
                                ticklabelsize = 18)
            cbar.ticks = [0,0.5,1]
            fig[i, j+1] = cbar
            c+=1
            push!(axes, ax) # just in case you need them later.
            fig
        end
        fig
    end
    fig = FigGridHeat()
    fig
    save(joinpath(@__DIR__, "output", "heatmapGrid.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
m = PoolModel("ss1")
fp = dirname(@__FILE__)

m <= Setting(:dataroot, "$(fp)/../reference/")
m <= Setting(:data_vintage, "190822")
observables = OrderedDict{Symbol,Observable}()
identity_transform_805 = function(levels)
    return levels[!,:p805]
end
identity_transform_904 = function(levels)
    return levels[!,:p904]
end
observables[:Model904] = Observable(:Model904, [:p904__wrongorigmatlab],
                                    identity_transform_904,
                                    identity_transform_904,
                                    "Model 904 predictive density",
                                    "Model 904 conditional predictive density scores")
observables[:Model805] = Observable(:Model805, [:p805__wrongorigmatlab],
                                    identity_transform_805,
                                    identity_transform_805,
                                    "Model 805 predictive density",
                                    "Model 805 conditional predictive density scores")

@testset "Test various forms of data loading" begin
    m.observable_mappings = observables
    @info "The following warnings are expected"
    df1 = load_data(m; try_disk = false)
    df1 = load_data(m)
    df2 = CSV.read(joinpath(dataroot(m), "raw/wrongorigmatlab_190822.csv"))
    df2[!,:date] = Vector{Dates.Date}(df2[!,:date])
    df2[!,:p904] = Vector{Float64}(df2[!,:p904])
    df2[!,:p805] = Vector{Float64}(df2[!,:p805])

    @test df1[!,:Model904] == df2[!,:p904]
    @test df1[!,:Model805] == df2[!,:p805]
end

nothing


text:
using DifferentialEquations
using OrdinaryDiffEq
using Distributions
using MonteCarloMeasurements
using StatsBase
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

δt = 1.0
tmax = 40.0
tspan = (0.0,tmax);

N = 1000.0;

n_samples = 1000; # Number of samples

p = [Particles(n_samples,Uniform(0.01,0.1)),
      Particles(n_samples,Uniform(5,20.0)),
      Particles(n_samples,Uniform(0.1,1.0))]

I₀=Particles(n_samples,Uniform(1.0,50.0))
u0 = [N-I₀,I₀,0.0]

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

sol_ode = solve(prob_ode, Tsit5(), dt=δt);

s20 = sol_ode(20.0)

l = @layout [a b c]
binwidth = 50
pl1 = histogram(s20[1],bins=0:binwidth:N, title="S(20)", xlabel="S", ylabel="Frequency", color=:blue)
pl2 = histogram(s20[2],bins=0:binwidth:N, title="I(20)", xlabel="I", ylabel="Frequency", color=:red)
pl3 = histogram(s20[3],bins=0:binwidth:N, title="R(20)", xlabel="R", ylabel="Frequency", color=:green)
plot(pl1,pl2,pl3,layout=l,legend=false)

corkendall(hcat(Array(p),Array(I₀)),Array(s20))


text:
#=
henonHeiles:
- Julia version:
- Author: ymocquar
- Date: 2019-10-22
=#
"""
 henon_heiles(u, p) : Henon-Heiles differential equation with parameters such that
it can be solved by the Julia solver ODEProblem and solve

Parameters :
    - u : the vector
    - p : the parameter of this version of Henon-Heiles equation, that is epsilon
    - t : the time that is unused.
Ouput :
    - the derivative vector
"""
# function henon_heiles(u,p)
#     return [u[3]/p, u[4], -u[1]/p-2u[1]*u[2], -u[2]-u[1]^2+u[2]^2]
# end
# henon_heiles_jul(u,p,t)=henon_heiles(u,p)
function henon_heiles(u, p, t)
    return [0, u[4], -2u[1] * u[2], -u[2] - u[1]^2 + u[2]^2]
end
function henon_heiles_julia(u, p, t)
    return [u[3] / p, u[4], -u[1] / p - 2u[1] * u[2], -u[2] - u[1]^2 + u[2]^2]
end


text:
struct ServiceMap
    thermal::Dict{Base.UUID,PSY.ThermalGen}
    hydro::Dict{Base.UUID,PSY.HydroGen}
    re::Dict{Base.UUID,PSY.RenewableGen}
    storage::Dict{Base.UUID,PSY.Storage}
end

function ServiceMap()
    return ServiceMap1(Dict{Base.UUID,PSY.ThermalGen}(),
            Dict{Base.UUID,PSY.HydroGen}(),
            Dict{Base.UUID,PSY.RenewableGen}(),
            Dict{Base.UUID,PSY.Storage}())
end

add_device(smap::ServiceMap, d::PSY.ThermalGen) =smap.thermal[IS.getuuid(d)] = d
add_device(smap::ServiceMap, d::PSY.HydroGen) = smap.hydro[IS.get_uuid(d)] = d
add_device(smap::ServiceMap, d::PSY.RenewableGen) = smap.re[IS.get_uuid(d)] = d
add_device(smap::ServiceMap, d::PSY.Storage) = smap.storage[IS.get_uuid(d)] = d

get_fields(smap::ServiceMap) = [smap.thermal,smap.hydro,smap.re,smap.storage]


text:
using Test
using VPTrees
using Random
using StringDistances

function hamming(a::Integer, b::Integer)
    count_ones(xor(a, b))
end

@testset "VPTree.jl" begin
    @testset "find in radius" begin
        Random.seed!(1)
        data = [UInt(1), UInt(15)]
        metric = hamming
        vptree = VPTree(data, metric)
        @test [1] == find(vptree, UInt(3), 1)
        @test Set([1, 2]) == Set(find(vptree, UInt(3), 2))
        data = [2,3,6]
        vptree = VPTree(data, hamming);
        @test Set([d for d in data if (hamming(d, 2)) <= 1]) == Set(data[find(vptree, 2, 1)])
        data = collect(1:200)
        vptree = VPTree(data, hamming);
        for f in data[find(vptree, 2, 4)]
            @test hamming(2, f) <= 4
        end
    end

    @testset "find n neighbors" begin
        Random.seed!(1)
        data = [UInt(1), UInt(15)]
        metric = hamming
        vptree = VPTree(data, metric)
        @test [1] == find_nearest(vptree, UInt(3), 1)
        @test Set([1, 2]) == Set(find_nearest(vptree, UInt(3), 2))
        data = collect(1:7)
        vptree = VPTree(data, hamming);
        find_nearest(vptree, 2, 2)
        [string(i, base=2) for i in data[find_nearest(vptree, 2, 2)]]
        for f in data[find_nearest(vptree, 2, 3)]
            @test hamming(2, f) <= 1
        end
    end

    @testset "quickselect" begin
        Random.seed!(1)
        a = [5,7,6,90,7,-1,3]
        k = 4
        distances = [abs(d + 5) for d in a]
        VPTrees.select!(a, k, distances)
        @test a[4] == 6
    end

    @testset "euclidean distance" begin
        Random.seed!(1)
        data = [(1,2),(15,16)]
        metric = (a, b) -> sqrt(sum((a .- b).^2))
        vptree = VPTree(data, metric)
        @test [1] == find(vptree, (3,3), 3.)
        @test Set([1, 2]) == Set(find(vptree, (3,3), 100.))
        @test Set([(1,2)]) == Set(data[find_nearest(vptree, (2,2), 1)])
    end

    @testset "levenshtein distance" begin
        Random.seed!(1)
        data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
        metric = (a, b) -> evaluate(Levenshtein(),a,b)
        vptree = VPTree(data, metric)
        query="blau"
        @test Set(["bla", "blub"]) == Set(data[find(vptree, query, 2)])
        @test Set(["bla", "blub", "baube"]) == Set(data[find_nearest(vptree, query, 3)])
    end
end


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

mutable struct IoK8sApiAppsV1beta1ScaleStatus <: SwaggerModel
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    selector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: selector
    targetSelector::Any # spec type: Union{ Nothing, String } # spec name: targetSelector

    function IoK8sApiAppsV1beta1ScaleStatus(;replicas=nothing, selector=nothing, targetSelector=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1ScaleStatus, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1beta1ScaleStatus, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiAppsV1beta1ScaleStatus, Symbol("targetSelector"), targetSelector)
        setfield!(o, Symbol("targetSelector"), targetSelector)
        o
    end
end # type IoK8sApiAppsV1beta1ScaleStatus

const _property_map_IoK8sApiAppsV1beta1ScaleStatus = Dict{Symbol,Symbol}(Symbol("replicas")=>Symbol("replicas"), Symbol("selector")=>Symbol("selector"), Symbol("targetSelector")=>Symbol("targetSelector"))
const _property_types_IoK8sApiAppsV1beta1ScaleStatus = Dict{Symbol,String}(Symbol("replicas")=>"Int32", Symbol("selector")=>"Dict{String, String}", Symbol("targetSelector")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1ScaleStatus }) = collect(keys(_property_map_IoK8sApiAppsV1beta1ScaleStatus))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1ScaleStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta1ScaleStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1ScaleStatus }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1ScaleStatus[property_name]

function check_required(o::IoK8sApiAppsV1beta1ScaleStatus)
    (getproperty(o, Symbol("replicas")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta1ScaleStatus }, name::Symbol, val)
end


text:
function xdot(x,y,z; sigma=10)
    return sigma*(y-x)
end

function ydot(x,y,z; rho=48)
    return x*(rho-z)-y
end

function zdot(x,y,z; beta=8/3)
    return x*y - beta*z
end

function make_step(x,y,z; delta=0.01)
    xp = x + xdot(x,y,z)*delta
    yp = y + ydot(x,y,z)*delta
    zp = z + zdot(x,y,z)*delta

    xn = x + 0.5*(xdot(x,y,z)+xdot(xp,yp,zp))*delta
    yn = y + 0.5*(ydot(x,y,z)+ydot(xp,yp,zp))*delta
    zn = z + 0.5*(zdot(x,y,z)+zdot(xp,yp,zp))*delta

    return xn, yn, zn
end

function run_sim(;runs=3000, count_z = false, timing=false, thresh = 0.0)
    x, y, z = 0, 1, 0
    X, Y, Z = Float64[0,], Float64[1,], Float64[0,]
    #dist = fill(0,2)
    cnt = 0
    zless = true
    cnt_list = Int[]

    for i in 1:runs
        x,y,z = make_step(x,y,z)
        push!(X, x)
        push!(Y, y)
        push!(Z, z)

        if (z < thresh && zless) || (z > thresh && !zless)
            cnt += 1
        else
            push!(cnt_list, cnt)
            zless = !zless
            cnt = 0
        end

        # if count_z
        #     if z > thresh
        #         dist[1] += 1
        #     else
        #         dist[2] += 1
        #     end
        # end
    end
    return X, Y, Z, cnt_list #,dist
end

X, Y, Z, times = run_sim(;runs=30000, timing=true, count_z = false, thresh = 24)

using Plots
Plots.plot(X,Y)
Plots.plot(Z)
Plots.histogram(Z;bins=100)

Plots.histogram((times*0.01)[3:end],bins=20)

function get_maxima(ls; add_min=false)
    maxs = Float64[]
    if add_min
        mins = Float64[]
        if ls[1]<ls[2]
            push!(mins,ls[1])
        end
    end
    if ls[1]>ls[2]
        push!(maxs,ls[1])
    end
    for i=2:length(ls)-1
        if ls[i-1]<ls[i]>ls[i+1]
            push!(maxs,ls[i])
        end
        if add_min
            if ls[i-1]>ls[i]<ls[i+1]
                push!(mins,ls[i])
            end
        end
    end
    if ls[end]>ls[end-1]
        push!(maxs,ls[end])
    end
    if add_min
        if ls[end]<ls[end-1]
            push!(mins,ls[end])
        end
        return maxs, mins
    else
        return maxs
    end
end

# get subsequent maxima plot
z_max = get_maxima(Z)
z_plot = Tuple{Float64, Float64}[]
for i=1:length(z_max)-1
    push!(z_plot,(z_max[i], z_max[i+1]))
end
i1 = [x[1] for x in z_plot]
i2 = [x[2] for x in z_plot]
Plots.scatter(i1,i2)

# get sequential max/min plot
z_max, z_min = get_maxima(Z, add_min=true)
Plots.scatter(z_max, z_min)

#=
points = Observable(Point3f[])
colors = Observable(Int[])
set_theme!(theme_black())

fig, ax, l = lines(points, color = colors,
    colormap = :inferno, transparency = true,
    axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
        viewmode = :fit, limits = (-30, 30, -30, 30, 0, 50)))

frames = 1:120
for frame in frames
    for i in 1:50
        push!(points[], step!(attractor))
        push!(colors[], frame)
    end
    ax.azimuth[] = 1.7pi + 0.3 * sin(2pi * frame / 120)
    notify.((points, colors))
    l.colorrange = (0, frame)
    sleep(0.1)
    #println(frame)
end
=#

max_stars_repo_name: m-geo/lorenz

max_issues_repo_name: m-geo/lorenz

max_forks_repo_name: m-geo/lorenz


text:
# convex relaxation (envelope) of cos function
function cv_cos(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  r::T = zero(T)
  kL::Int64 = Base.ceil((-0.5-xL/(2.0*pi)))
  if (x<=(pi-2.0*pi*kL))
    xL1::T = convert(T,xL+2.0*pi*kL)
    if (xL1 >= pi/2.0)
      return cos(x),-sin(x)
    end
    xU1::T = convert(T,min(xU+2.0*pi*kL,pi))
    if ((xL1>=(-pi/2))&&(xU1<=(pi/2)))
      if (abs(xL-xU)<MC_param.env_tol)
        r = zero(T)
      else
        r = (cos(xU)-cos(xL))/(xU-xL)
      end
      return cos(xL)+r*(x-xL),r
    end
    return cv_cosin(convert(T,x+(2.0*pi)*kL),xL1,xU1)
  end
  kU::Int64 = Base.floor((0.5-xU/(2.0*pi)))
  if (x>=(-pi-2.0*pi*kU))
    xU2::T = convert(T,xU+2.0*pi*kU)
    if (xU2<=-pi/2.0)
      return cos(x),-sin(x)
    end
    return cv_cosin(convert(T,x+2.0*pi*kU),convert(T,max(xL+2.0*pi*kU,-pi)),xU2)
  end
  return -one(T),zero(T)
end
# function for computing convex relaxation over nonconvex and nonconcave regions
function cv_cosin(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  xj::T = -Inf
  if (abs(xL)<=abs(xU))
    left::Bool = false
    x0::T = xU
    xm::T = xL
  else
    left = true
    x0 = xL
    xm = xU
  end
  try
    xj = newton(x0,xL,xU,cv_cosenv,dcv_cosenv,xm,zero(x0))
  catch e
    if isa(e, ErrorException)
      xj = golden_section(xL,xU,cv_cosenv,xm,zero(x0))
    end
  end
  if ((left && x<=xj)||((~left) && x>=xj))
    return cos(x),-sin(x)
  else
    if abs(xm-xj)<MC_param.env_tol
      r::T = zero(T)
    else
      r = (cos(xm)-cos(xj))/(xm-xj)
    end
    return cos(xm)+r*(x-xm),r
  end
end
# pivot point calculation function for convex relaxation of cosine
function cv_cosenv(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (x-y)*sin(x)+cos(x)-cos(y)
end
function dcv_cosenv(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (x-y)*cos(x)
end
# concave relaxation (envelope) of cos function
function cc_cos(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  temp = cv_cos(convert(T,x-pi),convert(T,xL-pi),convert(T,xU-pi))
  return -temp[1],-temp[2]
end
function cos_arg(xL::T,xU::T) where {T<:AbstractFloat}
  kL::Int64 = Base.ceil(-0.5-xL/(2.0*pi))
  xL1::T = convert(T,xL+2.0*pi*kL)
  xU1::T = convert(T,xU+2.0*pi*kL)
  if ~((xL1>=-pi)&&(xL1<=pi))
    error("Cosine Argument Calculation: xL out of bounds.")
  end
  if (xL1<=zero(T))
    if (xU1<=zero(T))
      arg1::T = xL
      arg2::T = xU
    elseif (xU1>=pi)
      arg1 = convert(T,pi-2.0*pi*kL)
      arg2 = convert(T,-2.0*pi*kL)
    else
      arg1 = (cos(xL1)<=cos(xU1)) ? xL : xU
      arg2 = convert(T,-2.0*pi*kL)
    end
  end
  if (xU1<=pi)
    arg1 = xU
    arg2 = xL
  elseif (xU1>=(2.0*pi))
    arg1 = convert(T,pi-2.0*pi*kL)
    arg2 = convert(T,2.0*pi-2.0*pi*kL)
  else
    arg1 = convert(T,pi-2.0*pi*kL)
    arg2 = (cos(xL1)>=cos(xU1)) ? xL : xU
  end
  return [arg1,arg2]
end

function cos(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T,eps_min::T = cos_arg(x.Intv.lo,x.Intv.hi)
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
	midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_cos(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_cos(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu>=1)
    gcc1::T,gdcc1::T = cc_cos(x.cv,x.Intv.lo,x.Intv.hi)
    gcv1::T,gdcv1::T = cv_cos(x.cv,x.Intv.lo,x.Intv.hi)
    gcc2::T,gdcc2::T = cc_cos(x.cc,x.Intv.lo,x.Intv.hi)
    gcv2::T,gdcv2::T = cv_cos(x.cc,x.Intv.lo,x.Intv.hi)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, cos(x.Intv),x.cnst, x.IntvBox, x.xref)
end

function sin(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  cos(x-pi/2.0)
end

# pivot point calculation function for convex relaxation of complimentary error function
@inline function tan_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (x-y)-(tan(x)-tan(y))/(one(T)+tan(x)^2)
end
# derivative of pivot point calculation function for convex relaxation of complimentary error function
@inline function tan_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return 2*tan(x)/(one(T)+tan(x)^2)*(tan(x)-tan(y))
end
# convex relaxation (envelope) of tangent function
@inline function cv_tan(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return tan(x),sec(x)^2
  elseif (xU<=zero(T))
    return line_seg(x,xL,tan(xL),xU,tan(xU)),dline_seg(x,xL,tan(xL),xU,tan(xU),sec(x)^2)
  else
    try
      p = secant(zero(T),xU,zero(T),xU,tan_env,xL,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(T),xU,tan_env,xL,zero(T))
      end
    end
    if (x<=p)
      return line_seg(x,xL,tan(xL),p,tan(p)),dline_seg(x,xL,tan(xL),p,tan(p),sec(x)^2)
    else
      return tan(x),sec(x)^2
    end
  end
end
# concave relaxation (envelope) of tangent function
@inline function cc_tan(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,tan(xL),xU,tan(xU)),dline_seg(x,xL,tan(xL),xU,tan(xU),sec(x)^2)
  elseif (xU<=zero(T))
    return tan(x),sec(x)^2
  else
    try
      p = secant(zero(T),xL,xL,zero(T),tan_env,xU,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),tan_env,xU,zero(T))
      end
    end
    if (x<=p)
       return tan(x),sec(x)^2
    else
       return line_seg(x,p,tan(p),xU,tan(xU)),dline_seg(x,p,tan(p),xU,tan(xU),sec(x)^2)
     end
  end
end
@inline function tan(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  if ((x.Intv.lo==-Inf)||(x.Intv.hi==Inf))
    error("Function unbounded on domain")
  end
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_tan(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_tan(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_tan(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_tan(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_tan(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_tan(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, tan(x.Intv),x.cnst, x.IntvBox, x.xref)
end

@inline function acos(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  return asin(-x)+pi/2.0
end

# pivot point calculation function for convex relaxation of arcsine
@inline function asin_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (asin(x)-asin(y))/(x-y)-one(T)/sqrt(one(T)-x^2)
end
# derivative of pivot point calculation function for convex relaxation of arcsine
@inline function asin_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return one(T)/((x-y)*sqrt(one(T)-x^2))-x/(one(T)-x^2)^(3/2)-(asin(x)-asin(y))/(x-y)^2
end
# convex relaxation (envelope) of arcsine function
@inline function cv_asin(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return asin(x),one(T)/sqrt(one(T)-x^2)
  elseif (xU<=zero(T))
    return line_seg(x,xL,asin(xL),xU,asin(xU)),dline_seg(x,xL,asin(xL),xU,asin(xU),one(T)/sqrt(one(T)-x^2))
  else
    try
      p = newton(xU,xL,xU,asin_env,asin_envd,xL,xU)
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(x),xU,asin_env,xL,zero(T))
      end
    end
    if (x<=p)
      return line_seg(x,xL,asin(xL),p,asin(p)),dline_seg(x,xL,asin(xL),p,asin(p),one(T)/sqrt(one(T)-x^2))
    else
      return asin(x),one(T)/sqrt(one(T)-x^2)
    end
  end
end
# concave relaxation (envelope) of arcsine function
@inline function cc_asin(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,asin(xL),xU,asin(xU)),dline_seg(x,xL,asin(xL),xU,asin(xU),one(T)/sqrt(one(T)-x^2))
  elseif (xU<=zero(T))
    return asin(x),one(T)/sqrt(one(T)-x^2)
  else
    try
      p = secant(zero(T),xL,xL,zero(T),asin_env,xU,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),asin_env,xU,zero(T))
      end
    end
    if (x<=p)
      return asin(x),one(T)/sqrt(one(T)-x^2)
    else
      return line_seg(x,p,asin(p),xU,asin(xU)),dline_seg(x,p,asin(p),xU,asin(xU),one(T)/sqrt(one(T)-x^2))
    end
  end
end
@inline function asin(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  if ((x.Intv.lo==-Inf)||(x.Intv.hi==Inf))
    error("Function unbounded on domain")
  end
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_asin(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_asin(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_asin(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_asin(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_asin(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_asin(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, asin(x.Intv),x.cnst, x.IntvBox, x.xref)
end

# pivot point calculation function for convex relaxation of arctangent
@inline function atan_env(x::T,y::T,z::T) where {T<:AbstractFloat}
  return (x-y)-(one(T)+x^2)*(atan(x)-atan(y))
end
# derivative of pivot point calculation function for convex relaxation of arctangent
@inline function atan_envd(x::T,y::T,z::T) where {T<:AbstractFloat}
  return -2*x*(atan(x)-atan(y))
end
# convex relaxation (envelope) of arctan function
@inline function cv_atan(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return line_seg(x,xL,atan(xL),xU,atan(xU)),dline_seg(x,xL,atan(xL),xU,atan(xU),one(T)/(x^2+one(T)))
  elseif (xU<=zero(T))
    return atan(x),one(T)/(x^2+one(T))
  else
    try
      p = newton(xL,xL,zero(T),atan_env,atan_envd,xU,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(xL,zero(T),atan_env,xU,zero(T))
      end
    end
    if (x<=p)
      return atan(x),one(T)/(x^2+one(T))
    else
      return line_seg(x,p,atan(p),xU,atan(xU)),dline_seg(x,p,atan(p),xU,atan(xU),one(T)/(x^2+one(T)))
    end
  end
end
# concave relaxation (envelope) of arctan function
@inline function cc_atan(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  p::T = zero(T)
  if (xL>=zero(T))
    return atan(x),one(T)/(x^2+one(T))
  elseif (xU<=zero(T))
    return line_seg(x,xL,atan(xL),xU,atan(xU)),dline_seg(x,xL,atan(xL),xU,atan(xU),one(T)/(x^2+one(T)))
  else
    try
      p = newton(xU,zero(T),xU,atan_env,atan_envd,xL,zero(T))
    catch e
      if isa(e, ErrorException)
        p = golden_section(zero(T),xU,atan_env,xL,zero(T))
      end
    end
    if (x<=p)
      return line_seg(x,xL,atan(xL),p,atan(p)),dline_seg(x,xL,atan(xL),p,atan(p),one(T)/(x^2+one(T)))
    else
      return atan(x),one(T)/(x^2+one(T))
    end
  end
end
@inline function atan(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cc_atan(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cv_atan(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cc_atan(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcv1::T,gdcv1::T = cv_atan(x.cv,x.Intv.lo,x.Intv.hi,c)
    gcc2::T,gdcc2::T = cc_atan(x.cc,x.Intv.lo,x.Intv.hi,c)
    gcv2::T,gdcv2::T = cv_atan(x.cc,x.Intv.lo,x.Intv.hi,c)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, atan(x.Intv),x.cnst, x.IntvBox, x.xref)
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestConceptRDF
using Test

using Serd
using Serd.RDF: Triple, Literal, Resource
using Catlab

using SemanticFlowGraphs

const R = Resource

@present TestPres(Monocl) begin
  A::Ob
  B::Ob
  C::Ob

  A0::Ob
  B0::Ob
  ::SubOb(A0,A)
  ::SubOb(B0,B)

  f::Hom(A,B)
  g::Hom(A,otimes(B,C))

  f0::Hom(A0,B0)
  ::SubHom(f0,f,SubOb(A0,A),SubOb(B0,B))
end

prefix = RDF.Prefix("ex", "http://www.example.org/#")
stmts = presentation_to_rdf(TestPres, prefix)
#write_rdf(stdout, stmts)

@test Triple(R("ex","A"), R("rdf","type"), R("monocl","TypeConcept")) in stmts
@test Triple(R("ex","A"), R("rdf","type"), R("monocl","TypeConcept")) in stmts
@test Triple(R("ex","A0"), R("rdf","type"), R("monocl","TypeConcept")) in stmts
@test Triple(R("ex","f"), R("rdf","type"), R("monocl","FunctionConcept")) in stmts

@test Triple(R("ex","A"), R("monocl","id"), RDF.Literal("A")) in stmts
@test Triple(R("ex","f"), R("monocl","id"), RDF.Literal("f")) in stmts

@test Triple(R("ex","A0"), R("monocl","subtypeOf"), R("ex","A")) in stmts
@test Triple(R("ex","f0"), R("monocl","subfunctionOf"), R("ex","f")) in stmts

@test Triple(R("ex","g"), R("monocl","inputs"), R("ex","g:in1")) in stmts
@test Triple(R("ex","g"), R("monocl","hasInput"), R("ex","g:in1")) in stmts
@test Triple(R("ex","g:in1"), R("list","hasContents"), R("ex","A")) in stmts
@test Triple(R("ex","g:in1"), R("list","index"), Literal(1)) in stmts

@test Triple(R("ex","g"), R("monocl","outputs"), R("ex","g:out1")) in stmts
@test Triple(R("ex","g"), R("monocl","hasOutput"), R("ex","g:out1")) in stmts
@test Triple(R("ex","g"), R("monocl","hasOutput"), R("ex","g:out2")) in stmts
@test Triple(R("ex","g:out1"), R("list","hasContents"), R("ex","B")) in stmts
@test Triple(R("ex","g:out1"), R("list","hasNext"), R("ex","g:out2")) in stmts
@test Triple(R("ex","g:out1"), R("list","index"), Literal(1)) in stmts
@test Triple(R("ex","g:out2"), R("list","hasContents"), R("ex","C")) in stmts
@test Triple(R("ex","g:out2"), R("list","index"), Literal(2)) in stmts

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
abstract PropertyPrimitive

# Meaningless isless function used to sort in optimize_batching
function Base.isless{T <: PropertyPrimitive}(a::T, b::T)
    for field in fieldnames(T)
        x = getfield(a, field)
        y = getfield(b, field)
        if isa(x, Colorant)
            if color_isless(x, y)
                return true
            elseif color_isless(y, x)
                return false
            end
        else
            if x < y
                return true
            elseif x > y
                return false
            end
        end
    end
    return false
end

immutable Property{P <: PropertyPrimitive} <: ComposeNode
    primitives::Vector{P}
end

function isempty(p::Property)
    return isempty(p.primitives)
end

function isscalar(p::Property)
    return length(p.primitives) == 1
end

# Some properties can be applied multiple times, most cannot.
function isrepeatable(p::Property)
    return false
end

function resolve{T}(box::AbsoluteBox, units::UnitBox, t::Transform, p::Property{T})
    return Property{T}([resolve(box, units, t, primitive) for primitive in p.primitives])
end

# Property primitive catchall: most properties don't need measure transforms
function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::PropertyPrimitive)
    return primitive
end

# Stroke
# ------

immutable StrokePrimitive <: PropertyPrimitive
	color::RGBA{Float64}
end

typealias Stroke Property{StrokePrimitive}

function stroke(c::(@compat Void))
    return Stroke([StrokePrimitive(RGBA{Float64}(0, 0, 0, 0))])
end

function stroke(c::@compat(Union{Colorant, AbstractString}))
	return Stroke([StrokePrimitive(parse_colorant(c))])
end

function stroke(cs::AbstractArray)
	return Stroke([StrokePrimitive(c == nothing ? RGBA{Float64}(0, 0, 0, 0) : parse_colorant(c)) for c in cs])
end

prop_string(::Stroke) = "s"

# Fill
# ----

immutable FillPrimitive <: PropertyPrimitive
	color::RGBA{Float64}
end

typealias Fill Property{FillPrimitive}

function fill(c::(@compat Void))
    return Fill([FillPrimitive(RGBA{Float64}(0.0, 0.0, 0.0, 0.0))])
end

function fill(c::@compat(Union{Colorant, AbstractString}))
	return Fill([FillPrimitive(parse_colorant(c))])
end

function fill(cs::AbstractArray)
	return Fill([FillPrimitive(c == nothing ? RGBA{Float64}(0.0, 0.0, 0.0, 0.0) : parse_colorant(c)) for c in cs])
end

prop_string(::Fill) = "f"

# StrokeDash
# ----------

immutable StrokeDashPrimitive <: PropertyPrimitive
    value::Vector{Measure}
end

typealias StrokeDash Property{StrokeDashPrimitive}

function strokedash(values::AbstractArray)
    return StrokeDash([StrokeDashPrimitive(collect(Measure, values))])
end

function strokedash(values::AbstractArray{AbstractArray})
    return StrokeDash([StrokeDashPrimitive(collect(Measure, value)) for value in values])
end

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::StrokeDashPrimitive)
    return StrokeDashPrimitive([resolve(box, units, t, v)
                                for v in primitive.value])
end

prop_string(::StrokeDash) = "sd"

# StrokeLineCap
# -------------

abstract LineCap
immutable LineCapButt <: LineCap end
immutable LineCapSquare <: LineCap end
immutable LineCapRound <: LineCap end

immutable StrokeLineCapPrimitive <: PropertyPrimitive
    value::LineCap

    function StrokeLineCapPrimitive(value::LineCap)
        return new(value)
    end

    function StrokeLineCapPrimitive(value::Type{LineCap})
        return new(value())
    end
end

typealias StrokeLineCap Property{StrokeLineCapPrimitive}

function strokelinecap(value::@compat(Union{LineCap, Type{LineCap}}))
    return StrokeLineCap([StrokeLineCapPrimitive(value)])
end

function strokelinecap(values::AbstractArray)
    return StrokeLineCap([StrokeLineCapPrimitive(value) for value in values])
end

prop_string(::StrokeLineCap) = "slc"

# StrokeLineJoin
# --------------

abstract LineJoin
immutable LineJoinMiter <: LineJoin end
immutable LineJoinRound <: LineJoin end
immutable LineJoinBevel <: LineJoin end

immutable StrokeLineJoinPrimitive <: PropertyPrimitive
    value::LineJoin

    function StrokeLineJoinPrimitive(value::LineJoin)
        return new(value)
    end

    function StrokeLineCapPrimitive(value::Type{LineJoin})
        return new(value())
    end
end

typealias StrokeLineJoin Property{StrokeLineJoinPrimitive}

function strokelinejoin(value::@compat(Union{LineJoin, Type{LineJoin}}))
    return StrokeLineJoin([StrokeLineJoinPrimitive(value)])
end

function strokelinejoin(values::AbstractArray)
    return StrokeLineJoin([StrokeLineJoinPrimitive(value) for value in values])
end

prop_string(::StrokeLineJoin) = "slj"

# LineWidth
# ---------

immutable LineWidthPrimitive <: PropertyPrimitive
    value::Measure

    function LineWidthPrimitive(value)
        return new(size_measure(value))
    end
end

typealias LineWidth Property{LineWidthPrimitive}

function linewidth(value::@compat(Union{Measure, Number}))
    return LineWidth([LineWidthPrimitive(value)])
end

function linewidth(values::AbstractArray)
    return LineWidth([LineWidthPrimitive(value) for value in values])
end

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::LineWidthPrimitive)
    return LineWidthPrimitive(resolve(box, units, t, primitive.value))
end

prop_string(::LineWidth) = "lw"

# Visible
# -------

immutable VisiblePrimitive <: PropertyPrimitive
    value::Bool
end

typealias Visible Property{VisiblePrimitive}

function visible(value::Bool)
    return Visible([VisiblePrimitive(value)])
end

function visible(values::AbstractArray)
    return Visible([VisiblePrimitive(value) for value in values])
end

prop_string(::Visible) = "v"

# FillOpacity
# -----------

immutable FillOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function FillOpacityPrimitive(value_::Number)
        value = @compat Float64(value_)
        if value < 0.0 || value > 1.0
            error("Opacity must be between 0 and 1.")
        end
        return new(value)
    end
end

typealias FillOpacity Property{FillOpacityPrimitive}

function fillopacity(value::Float64)
    return FillOpacity([FillOpacityPrimitive(value)])
end

function fillopacity(values::AbstractArray)
    return FillOpacity([FillOpacityPrimitive(value) for value in values])
end

prop_string(::FillOpacity) = "fo"

# StrokeOpacity
# -------------

immutable StrokeOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function StrokeOpacityPrimitive(value_::Number)
        value = @compat Float64(value_)
        if value < 0.0 || value > 1.0
            error("Opacity must be between 0 and 1.")
        end
        return new(value)
    end
end

typealias StrokeOpacity Property{StrokeOpacityPrimitive}

function strokeopacity(value::Float64)
    return StrokeOpacity([StrokeOpacityPrimitive(value)])
end

function strokeopacity(values::AbstractArray)
    return StrokeOpacity([StrokeOpacityPrimitive(value) for value in values])
end

prop_string(::StrokeOpacity) = "so"

# Clip
# ----

immutable ClipPrimitive{P <: Vec} <: PropertyPrimitive
    points::Vector{P}
end

typealias Clip Property{ClipPrimitive}

function clip()
    return Clip([ClipPrimitive(Array(Vec, 0))])
end

function clip{T <: XYTupleOrVec}(points::AbstractArray{T})
    XM, YM = narrow_polygon_point_types(Vector[points])
    if XM == Any
        XM = Length{:cx, Float64}
    end
    if YM == Any
        YM = Length{:cy, Float64}
    end
    VecType = Tuple{XM, YM}

    prim = ClipPrimitive(VecType[(x_measure(point[1]), y_measure(point[2])) for point in points])
    return Clip(typeof(prim)[prim])
end

function clip(point_arrays::AbstractArray...)
    XM, YM = narrow_polygon_point_types(point_arrays)
    VecType = XM == YM == Any ? Vec : Vec{XM, YM}
    PrimType = XM == YM == Any ? ClipPrimitive : ClipPrimitive{VecType}

    clipprims = Array(PrimType, length(point_arrays))
    for (i, point_array) in enumerate(point_arrays)
        clipprims[i] = ClipPrimitive(VecType[(x_measure(point[1]), y_measure(point[2]))
                                             for point in point_array])
    end
    return Property{PrimType}(clipprims)
end

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::ClipPrimitive)
    return ClipPrimitive{AbsoluteVec2}(
        AbsoluteVec2[resolve(box, units, t, point) for point in primitive.points])
end

prop_string(::Clip) = "clp"

# Font
# ----

immutable FontPrimitive <: PropertyPrimitive
    family::AbstractString
end

typealias Font Property{FontPrimitive}

function font(family::AbstractString)
    return Font([FontPrimitive(family)])
end

function font(families::AbstractArray)
    return Font([FontPrimitive(family) for family in families])
end

prop_string(::Font) = "fnt"

function Base.hash(primitive::FontPrimitive, h::UInt64)
    return hash(primitive.family, h)
end

function ==(a::FontPrimitive, b::FontPrimitive)
    return a.family == b.family
end

# FontSize
# --------

immutable FontSizePrimitive <: PropertyPrimitive
    value::Measure

    function FontSizePrimitive(value)
        return new(size_measure(value))
    end
end

typealias FontSize Property{FontSizePrimitive}

function fontsize(value::@compat(Union{Number, Measure}))
    return FontSize([FontSizePrimitive(value)])
end

function fontsize(values::AbstractArray)
    return FontSize([FontSizePrimitive(value) for value in values])
end

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::FontSizePrimitive)
    return FontSizePrimitive(resolve(box, units, t, primitive.value))
end

prop_string(::FontSize) = "fsz"

# SVGID
# -----

immutable SVGIDPrimitive <: PropertyPrimitive
    value::AbstractString
end

typealias SVGID Property{SVGIDPrimitive}

function svgid(value::AbstractString)
    return SVGID([SVGIDPrimitive(value)])
end

function svgid(values::AbstractArray)
    return SVGID([SVGIDPrimitive(value) for value in values])
end

prop_string(::SVGID) = "svgid"

function Base.hash(primitive::SVGIDPrimitive, h::UInt64)
    return hash(primitive.value, h)
end

function ==(a::SVGIDPrimitive, b::SVGIDPrimitive)
    return a.value == b.value
end

# SVGClass
# --------

immutable SVGClassPrimitive <: PropertyPrimitive
    value::Compat.UTF8String
end

typealias SVGClass Property{SVGClassPrimitive}

function svgclass(value::AbstractString)
    return SVGClass([SVGClassPrimitive(value)])
end

function svgclass(values::AbstractArray)
    return SVGClass([SVGClassPrimitive(value) for value in values])
end

function prop_string(svgc::SVGClass)
    if isscalar(svgc)
        return string("svgc(", svgc.primitives[1].value, ")")
    else
        return string("svgc(", svgc.primitives[1].value, "...)")
    end
end

function Base.hash(primitive::SVGClassPrimitive, h::UInt64)
    return hash(primitive.value, h)
end

function ==(a::SVGClassPrimitive, b::SVGClassPrimitive)
    return a.value == b.value
end

# SVGAttribute
# ------------

immutable SVGAttributePrimitive <: PropertyPrimitive
    attribute::Compat.ASCIIString
    value::Compat.ASCIIString
end

typealias SVGAttribute Property{SVGAttributePrimitive}

function svgattribute(attribute::AbstractString, value)
    return SVGAttribute([SVGAttributePrimitive(attribute, string(value))])
end

function svgattribute(attribute::AbstractString, values::AbstractArray)
    return SVGAttribute([SVGAttributePrimitive(attribute, string(value))
                         for value in values])
end

function svgattribute(attributes::AbstractArray, values::AbstractArray)
    return SVGAttribute(
        @makeprimitives SVGAttributePrimitive,
            (attribute in attributes, value in values),
            SVGAttributePrimitive(attribute, string(value)))
end

prop_string(::SVGAttribute) = "svga"

function Base.hash(primitive::SVGAttributePrimitive, h::UInt64)
    h = hash(primitive.attribute, h)
    h = hash(primitive.value, h)
    return h
end

function ==(a::SVGAttributePrimitive, b::SVGAttributePrimitive)
    return a.attribute == b.attribute && a.value == b.value
end

# JSInclude
# ---------

immutable JSIncludePrimitive <: PropertyPrimitive
    value::AbstractString
    jsmodule::@compat(Union{(@compat Void), @compat Tuple{AbstractString, AbstractString}})
end

typealias JSInclude Property{JSIncludePrimitive}

function jsinclude(value::AbstractString, module_name=nothing)
    return JSInclude([JSIncludePrimitive(value, module_name)])
end

# Don't bother with a vectorized version of this. It wouldn't really make #
# sense.

prop_string(::JSInclude) = "jsip"

# JSCall
# ------

immutable JSCallPrimitive <: PropertyPrimitive
    code::AbstractString
    args::Vector{Measure}
end

typealias JSCall Property{JSCallPrimitive}

function jscall(code::AbstractString, arg::Vector{Measure}=Measure[])
    return JSCall([JSCallPrimitive(code, arg)])
end

function jscall(codes::AbstractArray,
                args::AbstractArray{Vector{Measure}}=Vector{Measure}[Measure[]])
    return JSCall(
        @makeprimitives JSCallPrimitive,
            (code in codes, arg in args),
            JSCallPrimitive(code, arg))
end

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 primitive::JSCallPrimitive)
    # we are going to build a new string by scanning across "code" and
    # replacing %x with translated x values, %y with translated y values
    # and %s with translated size values.
    newcode = IOBuffer()

    i = 1
    validx = 1
    while true
        j = search(primitive.code, '%', i)

        if j == 0
            write(newcode, primitive.code[i:end])
            break
        end

        write(newcode, primitive.code[i:j-1])
        if j == length(primitive.code)
            write(newcode, '%')
            break
        elseif primitive.code[j+1] == '%'
            write(newcode, '%')
        elseif primitive.code[j+1] == 'x'
            val = resolve(box, units, t, (primitive.args[validx], 0mm))
            write(newcode, svg_fmt_float(val[1].value))
            validx += 1
        elseif primitive.code[j+1] == 'y'
            val = resolve(box, units, t, (0mm, primitive.args[validx]))
            write(newcode, svg_fmt_float(val[2].value))
            validx += 1
        elseif primitive.code[j+1] == 's'
            val = resolve(box, units, t, primitive.args[validx])
            write(newcode, svg_fmt_float(val.value))
            validx += 1
        else
            write(newcode, '%', primitive.code[j+1])
        end

        i = j + 2
    end

    return JSCallPrimitive(takebuf_string(newcode), Measure[])
end

function isrepeatable(p::JSCall)
    return true
end

function Base.isless(a::FillPrimitive, b::FillPrimitive)
    return color_isless(a.color, b.color)
end

function Base.isless(a::StrokePrimitive, b::StrokePrimitive)
    return color_isless(a.color, b.color)
end

prop_string(::JSCall) = "jsc"


text:
struct HeatMapLike end
struct WireframeLike end
struct SeriesLike end
struct HistogramLike end
struct ViolinLike end

struct DimensionalPlot end

@recipe function f(A::AbstractDimArray)
    DimensionalPlot(), A
end

@recipe function f(::DimensionalPlot, A::AbstractArray)
    A_fwd = _forwardorder(A)
    sertype = get(plotattributes, :seriestype, :none)
    if !(sertype in [:marginalhist])
        :title --> refdims_title(A_fwd)
    end
    if ndims(A) > 2
        parent(A)
    elseif sertype in [:heatmap, :contour, :volume, :marginalhist,
                       :surface, :contour3d, :wireframe, :scatter3d]
        HeatMapLike(), A_fwd
    elseif sertype in [:histogram, :stephist, :density, :barhist, :scatterhist, :ea_histogram]
        HistogramLike(), A_fwd
    elseif sertype in [:hline]
        :yguide --> label(A_fwd)
        parent(A_fwd)
    elseif sertype in [:vline, :andrews]
        :xguide --> label(A_fwd)
        parent(A_fwd)
    elseif sertype in [:violin, :dotplot, :boxplot]
        ViolinLike(), A_fwd
    elseif sertype in [:plot, :histogram2d, :none, :line, :path, :shape, :steppre,
                       :steppost, :sticks, :scatter, :hexbin, :barbins, :scatterbins,
                       :stepbins, :bins2d, :bar]
        SeriesLike(), A_fwd
    else
        parent(A_fwd)
    end
end

@recipe function f(s::SeriesLike, A::AbstractArray{T,1}) where T
    dim = dims(A, 1)
    :xguide --> label(dim)
    :yguide --> label(A)
    _xticks!(plotattributes, s, dim)
    _withaxes(dim, A)
end
@recipe function f(s::SeriesLike, A::AbstractArray{T,2}) where T
    A = permutedims(A, commondims(>:, (TimeDim, XDim, IndependentDim, YDim, ZDim, DependentDim, Dimension, Dimension), dims(A)))
    ind, dep = dims(A)
    :xguide --> label(ind)
    :yguide --> label(A)
    :legendtitle --> label(dep)
    :label --> permutedims(val(dep))
    :tickfontalign --> :left
    _xticks!(plotattributes, s, ind)
    _withaxes(ind, A)
end
@recipe function f(s::SeriesLike, A::AbstractArray{T,3}) where T
    A = permutedims(A, commondims(>:, (TimeDim, XDim, IndependentDim, YDim, ZDim, DependentDim, DependentDim, Dimension, Dimension, Dimension), dims(A)))
    ind, dep1, dep2 = dims(A)
    :xguide --> label(ind)
    :yguide --> label(A)
    :legendtitle --> label(dep)
    :label --> permutedims(val(dep))
    :tickfontalign --> :left
    _xticks!(plotattributes, s, ind)
    _withaxes(ind, A)
end

@recipe function f(s::HistogramLike, A::AbstractArray{T,1}) where T
    dim = dims(A, 1)
    :xguide --> label(A)
    _withaxes(dim, A)
end
@recipe function f(s::HistogramLike, A::AbstractArray{T,2}) where T
    A = permutedims(A, commondims(>:, (ZDim, YDim, DependentDim, IndependentDim, XDim, TimeDim, IndependentDim, Dimension, Dimension), dims(A)))
    ind, dep = dims(A)
    :xguide --> label(A)
    :legendtitle --> label(dep)
    :label --> permutedims(index(dep))
    _withaxes(ind, A)
end

@recipe function f(::ViolinLike, A::AbstractArray{T,1}) where T
    dim = dims(A, 1)
    :yguide --> label(A)
    parent(A)
end
@recipe function f(s::ViolinLike, A::AbstractArray{T,2}) where T
    A = permutedims(A, commondims(>:, (ZDim, YDim, DependentDim, XDim, TimeDim, IndependentDim, Dimension, Dimension), dims(A)))
    dep, ind = dims(A)
    :xguide --> label(ind)
    :yguide --> label(A)
    :legendtitle --> label(ind)
    :label --> permutedims(index(ind))
    _xticks!(plotattributes, s, ind)
    parent(A)
end

@recipe function f(s::HeatMapLike, A::AbstractArray{T,1}) where T
    dim = dims(A, 1)
    :xguide --> label(dim)
    :yguide --> label(A)
    _xticks!(plotattributes, s, dim)
    parent(A)
end

@recipe function f(s::HeatMapLike, A::AbstractArray{T,2}) where T
    A = permutedims(A, commondims(>:, (ZDim, YDim, XDim, TimeDim, Dimension, Dimension), dims(A)))
    y, x = dims(A)
    :xguide --> label(x)
    :yguide --> label(y)
    :zguide --> label(A)
    :colorbar_title --> label(A)
    _xticks!(plotattributes, s, x)
    _yticks!(plotattributes, s, y)
    _withaxes(x, y, A)
end

_withaxes(dim::Dimension, A::AbstractDimArray) =
    _withaxes(mode(dim), index(dim), parent(A))
_withaxes(::NoIndex, index, A::AbstractArray) = A
_withaxes(::IndexMode, index, A::AbstractArray) = index, A
_withaxes(::Categorical, index, A::AbstractArray) = eachindex(index), A

_withaxes(dx::Dimension, dy::Dimension, A::AbstractDimArray) =
    _withaxes(mode(dx), mode(dy), index(dx), index(dy), parent(A))
_withaxes(::IndexMode, ::IndexMode, ix, iy, A) = ix, iy, A
_withaxes(::NoIndex, ::IndexMode, ix, iy, A) = axes(A, 2), iy, A
_withaxes(::IndexMode, ::NoIndex, ix, iy, A) = ix, axes(A, 1), A
_withaxes(::NoIndex, ::NoIndex, ix, iy, A) = axes(A, 2), axes(A, 1), A

_xticks!(attr, s, d::Dimension) = _xticks!(attr, s, mode(d), index(d))
_xticks!(attr, s, ::Categorical, index) =
    RecipesBase.is_explicit(attr, :xticks) || (attr[:xticks] = (eachindex(index), index))
_xticks!(attr, s, ::IndexMode, index) = nothing

_yticks!(attr, s, d::Dimension) = _yticks!(attr, s, mode(d), index(d))
_yticks!(attr, s, ::Categorical, index) =
    RecipesBase.is_explicit(attr, :yticks) || (attr[:yticks] = (eachindex(index), index))
_yticks!(attr, s, ::IndexMode, index) = nothing

_forwardorder(A::AbstractArray) =
    reorder(A, ForwardIndex) |> a -> reorder(a, ForwardRelation)

"""
    refdims_title(A::AbstractDimArray)
    refdims_title(refdims::Tuple)
    refdims_title(refdim::Dimension)

Generate a title string based on reference dimension values.
"""
function refdims_title end
refdims_title(A::AbstractArray; kw...) = refdims_title(refdims(A); kw...)
function refdims_title(refdims::Tuple; kw...)
    join(map(rd -> refdims_title(rd; kw...), refdims), ", ")
end
function refdims_title(refdim::Dimension; kw...)
    string(name(refdim), ": ", refdims_title(mode(refdim), refdim; kw...))
end
function refdims_title(mode::AbstractSampled, refdim::Dimension; kw...)
    start, stop = map(string, bounds(refdim))
    if start == stop
        start
    else
         "$start to $stop"
    end
end
refdims_title(mode::IndexMode, refdim::Dimension; kw...) = string(val(refdim))