text:
@doc raw"""
    Line2

An object `l` of the data type [`Line2`](@ref) is a directed straight line in
the two-dimensional Euclidean plane ``Ε^2``.

It is defined by the set of points with Cartesian coordinates ``(x, y)`` that
satisfy the equation

```math
    l:\; a\, x + b\, y + c = 0.
```

The line splits ``Ε^2`` in a *positive* and a *negative* side. A point p with
Cartesian coordinates ``(px, py)`` is on the positive side of `l`, iff
``a\, px + b\, py + c > 0``, it is on the negative side of `l`, iff
``a\, px + b\, py + c < 0``. The positive side is to the left of `l`.

# Example

Let us first define two Cartesian two-dimensional points in the Euclidean place
``Ε^2``. Their dimension and the fact that they are Cartesian is expressed by
the suffix `2` and the `C2` suffix in their representation.

```jldoctest
julia> p, q = Point2(1.0,1.0), Point2(4.0,7.0)
(PointC2(1, 1), PointC2(4, 7))
```

To define a line `l` we write:

```jldoctest
julia> l = Line2(p,q)
Line_2(-6, 3, 3)
```
"""
Line2

"""
    Line2(a::Real, b::Real, c::Real)

Introduces a line `l` with the line equation in Cartesian coordinates
``ax + by + c = 0``.
"""
Line2(a::Real, b::Real, c::Real) =
    Line2(convert(FT, a), convert(FT, b), convert(FT, c))

"""
    Line2(p::Point2, q::Point2)

Introduces a line `l` passing through the points `p` and `q`.

Line `l` is directed from `p` to `q`.
"""
Line2(p::Point2, q::Point2)

"""
    Line2(p::Point2, d::Direction2)

Introduces a line `l` passing through point `p` with direction `d`.
"""
Line2(p::Point2, d::Direction2)

"""
    Line2(p::Point2, d::Vector2)

Introduces a line `l` passing through point `p` and oriented by `v`.
"""
Line2(p::Point2, v::Vector2)

"""
    Line2(s::Segment2)

Introduces a line `l` supporting the segment `s`, oriented from source to
target.
"""
Line2(s::Segment2)

"""
    Line2(r::Ray2)

Introduces a line `l` supporting the ray `r`, with same orientation.
"""
Line2(r::Ray2)

"""
    ==(l::Line2, h::Line2)

Test for equality: two lines are equal, iff they have a non empty intersection
and the same direction.
"""
==(l::Line2, h::Line2)

"""
    a(l::Line2)

Returns the first coefficient of `l`.
"""
a(l::Line2)

"""
    b(l::Line2)

Returns the second coefficient of `l`.
"""
b(l::Line2)

"""
    c(l::Line2)

Returns the third coefficient of `l`.
"""
c(l::Line2)

"""
    point(l::Line2, i::Real)

Returns an arbitrary point on `l`.

It holds `point(l, i) == point(l, j)`, iff `i == j`. Furthermore, `l` is
directed from `point(l, i)` to `point(l, j)`, for all `i` ``<`` `j`.
"""
point(l::Line2, i::Real)

@cxxdereference point(l::Line2, i::Real) = point(l, convert(FT, i))

"""
    projection(l::Line2, p::Point2)

Returns the orthogonal projection of `p` onto `l`.
"""
projection(l::Line2, p::Point2)

"""
    x_at_y(l::Line2, y::Real)

Returns the ``x``-coordinate of the point at `l` with the given
``y``-coordinate.

!!! info "Precondition"

    `l` is not horizontal.
"""
x_at_y(l::Line2, y::Real)

@cxxdereference x_at_y(l::Line2, y::Real) = x_at_y(l, convert(FT, y))

"""
    y_at_x(l::Line2, x::Real)

Returns the ``y``-coordinate of the point at `l` with the given
``x``-coordinate.

!!! info "Precondition"

    `l` is not vertical.
"""
y_at_x(l::Line2, x)

@cxxdereference y_at_x(l::Line2, x::Real) = y_at_x(l, convert(FT, x))

"""
    is_degenerate(l::Line2)

Line `l` is degenerate, if the coefficients `a` and `b` of the line equation are
zero.
"""
is_degenerate(l::Line2)

"""
    is_horizontal(l::Line2)

Line `l` is horizontal, if coefficient `a` of the line equation is zero.
"""
is_horizontal(l::Line2)

"""
    is_vertical(l::Line2)

Line `l` is horizontal, if coefficient `b` of the line equation is zero.
"""
is_vertical(l::Line2)

"""
    oriented_side(l::Line2, p::Point2)

Returns [`ON_ORIENTED_BOUNDARY`](@ref), [`ON_NEGATIVE_SIDE`](@ref), or the
constant [`ON_POSITIVE_SIDE`](@ref), depending on the position of `p` relative
to the oriented line `l`.
"""
oriented_side(l::Line2, p::Point2)

"""
    has_on(l::Line2, p::Point2)

Returns `true`, iff `p` lies properly on `l`.
"""
has_on(l::Line2, p::Point2)

"""
    has_on_positive_side(l::Line2, p::Point2)

Returns `true`, iff `p` lies on the positive side of `l`.
"""
has_on_positive_side(l::Line2, p::Point2)

"""
    has_on_negative_side(l::Line2, p::Point2)

Returns `true`, iff `p` lies on the negative side of `l`.
"""
has_on_negative_side(c::Line2, p::Point2)

"""
    to_vector(l::Line2)

Returns a vector having the direction of `l`.
"""
to_vector(l::Line2)

"""
    direction(l::Line2)

Returns the direction of `l`.
"""
direction(l::Line2)

"""
    opposite(l::Line2)

Returns the line with opposite direction.
"""
opposite(l::Line2)

"""
    perpendicular(l::Line2, p::Point2)

Returns the line perpendicular to `l` and passing through `p`, where the
direction is the direction of `l` rotate counterclockwise by 90 degrees.
"""
perpendicular(l::Line2, p::Point2)

"""
    transform(l::Line2, t::AffTransformation2)

Returns the line obtained by applying `t` on a point on `l` and the direction of
`l`.
"""
transform(l::Line2, t::AffTransformation2)


text:
using MacroTools

"""
Expand:

    f(x, ...) = ...

Into:

    f(x, ...) = ...
    f(x::AbstractShim, ...) = f(super(x), ...)
"""
function inject_shim(ex)
    ex = shortdef(ex)
    if ex.head == :block
        exprs = ex.args
    else
        exprs = [ex]
    end
    return Expr(:block, Iterators.flatten(inject_shim1.(exprs))...)
end

function inject_shim1(ex)
    if @capture(ex, ((f_(all_args__; kwargs__))|f_(all_args__)) = body__)
        all_args = Expr(:tuple, all_args...)
        if @capture(all_args, (((self_::T_)|self_), args__))
        elseif @capture(all_args, ((self_::T_)|self_))
            args = []
        else
            error("Unexpected: $all_args")
        end
        if self isa Vector && length(self) == 1
            self, = self
        end
        if kwargs === nothing
            kwargs = []
        end
        return [
            ex,
            :($f($self::AbstractShim, $(args...); $(kwargs...)) =
              $f(super($self), $(args...); $(kwargs...))),
        ]
    else
        return [ex]
    end
end

macro shimmed(ex)
    esc(inject_shim(ex))
end

ext: jl


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

"""
The CMake module provides cross-platform installation of a
recent version of `cmake`, and exports a constant `cmake`
giving the path of this `cmake`` executable.
"""
module CMake
export cmake

depsjl = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsjl)
    error("CMake not properly installed. Please run\nPkg.build(\"CMake\")")
else
    include(depsjl)
end

@doc "`cmake` is a `String` giving the path of the `cmake` executable." cmake

end # module


text:
@symbol_func function tau_E(cur_reactor::AbstractReactor)
  cur_scaling = cur_reactor.mode_scaling

  cur_tau_E = cur_scaling[:constant]

  cur_tau_E *= cur_reactor.H

  cur_tau_E *= cur_reactor.I_P ^ cur_scaling[:I]

  cur_tau_E *= cur_reactor.R_0 ^ cur_scaling[:R]

  cur_tau_E *= a(cur_reactor) ^ cur_scaling[:a]

  cur_tau_E *= kappa_tau(cur_reactor) ^ cur_scaling[:kappa]

  cur_tau_E *= cur_reactor.n_bar ^ cur_scaling[:n]

  cur_tau_E *= cur_reactor.B_0 ^ cur_scaling[:B]

  cur_tau_E *= cur_reactor.A ^ cur_scaling[:A]

  cur_tau_E /= P_L(cur_reactor) ^ cur_scaling[:P]

  cur_tau_E
end


text:
include("solver_preamble.jl")
import SDPT3
factory = with_optimizer(SDPT3.Optimizer, printlevel = 0)
config = MOI.Test.Config(atol=1e-4, rtol=1e-4)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config, String[
        #  Expression: termination_status(model) == MOI.OPTIMAL
        #   Evaluated: MathOptInterface.SLOW_PROGRESS == MathOptInterface.OPTIMAL
        #  Expression: primal_status(model) == MOI.FEASIBLE_POINT
        #   Evaluated: MathOptInterface.NO_SOLUTION == MathOptInterface.FEASIBLE_POINT
        #  Expression: ≈(value(p), x ^ 3, atol=atol, rtol=rtol)
        #   Evaluated: 1.0173770753701414x³ + 1.5745867283645243e-11x² - 5.33707366495717e-12x + 7.1247443708677416e-12 ≈ x³ (atol=0.0001, rtol=0.0001)
        "sdsos_concave_then_convex_cubic"
    ])
end
@testset "SDP" begin
    Tests.sd_test(factory, config, String[
        #  Expression: JuMP.termination_status(model) == MOI.OPTIMAL
        #   Evaluated: MathOptInterface.OTHER_ERROR == MathOptInterface.OPTIMAL
        #  Expression: JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        #   Evaluated: MathOptInterface.NO_SOLUTION == MathOptInterface.FEASIBLE_POINT
        "chebyshev"
    ])
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/calc_sum_stats_pr.jl
## 2017 Keir A. Ashby
#  Calculates summary statistics that can be used to generate period-radius distribution plots
#  Based primarily off of calc_summary_stats_obs_binned_rates from christiansen_func.jl

function calc_sum_stats_pr(cat_obs::KeplerObsCatalog, param::SimParam, trueobs_cat::Bool = false, cat_phys::KeplerPhysicalCatalog=KeplerPhysicalCatalog([]))
  ssd = Dict{ASCIIString,Any}()
  cache = Dict{ASCIIString,Any}()

  if !trueobs_cat
    num_targets = get_int(param,"num_targets_sim_pass_one")
    ssd["num targets"] = num_targets
  else
    ssd["num targets"] = get_int(param,"num_kepler_targets")
  end

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet
# Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  if ( length( find(x::KeplerTargetObs-> length(x.obs) > max_tranets_in_sys, cat_obs.target[idx_tranets] ) ) > 0)   # Make sure max_tranets_in_sys is at least big enough for observed systems
    warn("Observational data has more transiting planets in one systems than max_tranets_in_sys allows.")
  end
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  num_sys_tranets = zeros(max_tranets_in_sys)                           # Since observed data, don't need to calculate probabilities.
  for n in 1:max_tranets_in_sys                                         # Make histogram of N-tranet systems
    num_sys_tranets[n] = length(idx_n_tranets[n])
  end
  ssd["num_sys_tranets"] = num_sys_tranets
  ssd["planets detected"] = num_tranets

    period_obs_list = zeros(num_tranets)
    weight_list = zeros(num_tranets)
    radius_obs_list = zeros(num_tranets)

  if trueobs_cat==true
    mult_true   = zeros(max_tranets_in_sys)

    n = 1    # tranet id
    for i in idx_tranets
      for j in 1:num_planets(cat_obs.target[i])
	if length(cat_obs.target[i].obs)==0 continue end
        period_obs_list[n] = cat_obs.target[i].obs[j].period
        weight_list[n] = 1.0
        radius_obs_list[n] = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius
        num_planets_in_sys = length(cat_obs.target[i].obs)
        mult_true[num_planets_in_sys]+=1
        n = n+1
      end
    end

    ssd["period_list"] = period_obs_list
    ssd["radius_list"] = radius_obs_list
    ssd["multiplicity_true"] = mult_true

  else  #Here we'll calculate periods and radii for the physical and observed catalogs as well as their planet multiplicities.
    period_phys_list = Float64[] # zeros(length(idx_phys))      #Holds the periods for every planet in the physical catalog
    radius_phys_list = Float64[] # zeros(length(idx_phys))	#Holds the radii for every planet in the physical catalog
    mult_phys        = zeros(max_tranets_in_sys) 		#Holds the number of systems for each variation of planet multiplicity in the physical catalog
    mult_obs         = zeros(max_tranets_in_sys) 		#Holds the number of systems for each variation of planet multiplicity in the observed catalog
    a = 1    # tranet id
    b = 1    # tranet id
    for i in 1:length(cat_phys.target) 				#get from physical catalog
       for s in 1:length(cat_phys.target[i].sys)
         if length(cat_phys.target[i].sys)==0 continue end
	 for p in 1:length(cat_phys.target[i].sys[s].planet)
	   if length(cat_phys.target[i].sys[s].planet)==0 continue end
	   #period_phys_list[a] = cat_phys.target[i].sys[s].orbit[p].P
           #radius_phys_list[a] = cat_phys.target[i].sys[s].planet[p].radius
	   push!(period_phys_list,cat_phys.target[i].sys[s].orbit[p].P)
           push!(radius_phys_list,cat_phys.target[i].sys[s].planet[p].radius)
	   num_planets_in_sys = length(cat_phys.target[i].sys[s].planet)
           mult_phys[num_planets_in_sys]+=1
           a+=1
         end
	end
    end
    for j in idx_tranets
      for h in 1:num_planets(cat_obs.target[j])
	if length(cat_obs.target[j].obs)==0 continue end
        period_obs_list[b] = cat_obs.target[j].obs[h].period
	weight_list[b] = 1.0
        radius_obs_list[b] = sqrt(cat_obs.target[j].obs[h].depth*cat_obs.target[j].star.radius)
        num_planets_in_sys = length(cat_obs.target[j].obs)
        mult_obs[num_planets_in_sys]+=1
        b+=1
      end
    end

    ssd["period_phys_list"]   = period_phys_list
    ssd["radius_phys_list"]   = radius_phys_list
    ssd["period_obs_list"]    = period_obs_list
    ssd["radius_obs_list"]    = radius_obs_list
    ssd["max_tranets_in_sys"] = max_tranets_in_sys
    ssd["multiplicity_phys"]  = mult_phys
    ssd["multiplicity_obs"]   = mult_obs
  end

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_obs_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_obs_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #ssd["planets detected"] = sum(np_bin)
  ssd["planets table"] = np_bin

  return CatalogSummaryStatistics(ssd, cache)
end


text:
"""

```
mddf_frame_self!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)
```

Computes the MDDF for a single frame, modifies the data in the `R` (`Result`) structure.

"""
function mddf_frame_self!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)

  # Simplify code by assigning some shortened names
  trajectory = framedata.trajectory
  volume_frame = framedata.volume_frame
  rdf_count_random_frame = framedata.rdf_count_random_frame
  md_count_random_frame = framedata.md_count_random_frame
  dc = framedata.dc
  dmin_mol = framedata.dmin_mol
  dref_mol = framedata.dref_mol
  x_solvent_random = framedata.x_solvent_random
  lc_solvent = framedata.lc_solvent
  solute = trajectory.solute
  solvent = trajectory.solvent
  x_solute = trajectory.x_solute
  x_solvent = trajectory.x_solvent

  # Reset counters for this frame
  reset!(volume_frame)
  @. rdf_count_random_frame = 0.
  @. md_count_random_frame = 0.

  # get pbc sides in this frame
  sides = getsides(trajectory,iframe)

  volume_frame.total = sides[1]*sides[2]*sides[3]
  R.volume.total = R.volume.total + volume_frame.total

  R.density.solute = R.density.solute + (solute.nmols / volume_frame.total)
  R.density.solvent = R.density.solvent + (solvent.nmols / volume_frame.total)

  # Add the box side information to the box structure, in this frame
  box = Box(options.lcell,sides,R.cutoff)

  # Will wrap everything relative to the reference atom of the first molecule
  # and move everything such that that center is in the origin. This is important
  # to simplify the computation of cell indexes, as the minimum coordinates are
  # automatically -side/2 at each direction
  solute_center = x_solute[R.irefatom]
  wrap!(x_solvent,sides,solute_center)
  center_to_origin!(x_solvent,solute_center)

  # Initialize linked cells
  initcells!(x_solvent,box,lc_solvent)

  # Check if the cutoff is not too large considering the periodic cell size
  if R.cutoff > sides[1]/2. || R.cutoff > sides[2]/2. || R.cutoff > sides[3]/2.
    error("""
      in MDDF: cutoff or dbulk > periodic_dimension/2 in frame: $iframe
               max(cutoff,dbulk) = $(R.cutoff)
               sides read from file = $(sides)
               If sides are zero it is likely that the PBC information is not available
               in the trajectory file.
    """)
  end

  n_solvent_in_bulk = 0.
  local n_dmin_in_bulk
  for isolvent in 1:solvent.nmols-1
    # We need to do this one solute molecule at a time to avoid exploding the memory
    # requirements
    x_this_solute = viewmol(isolvent,x_solvent,solvent)

    # Compute all distances between solute and solvent atoms which are smaller than the
    # cutoff (this is the most computationally expensive part), the distances are returned
    # in the dc structure
    cutoffdistances_self!(R.cutoff,x_this_solute,x_solvent,lc_solvent,box,dc,
                          solvent,isolvent)

    # For each solute molecule, update the counters (this is highly suboptimal, because
    # within updatecounters there are loops over solvent molecules, in such a way that
    # this will loop with cost nsolute*nsolvent. However, I cannot see an easy solution
    # at this point with acceptable memory requirements
    n_dmin_in_bulk, n_dref_in_bulk = updatecounters!(R,solute,solvent,dc,dmin_mol,dref_mol)
    n_solvent_in_bulk += n_dref_in_bulk
  end
  # The normalization below is tricky. The number that comes out from updatecounters is the
  # sum, for every solvent molecule (minus one) of the distances that were not found to be in
  # the solute domain. The total number of distances is n^2, because the sum inside updatecounter
  # is made for all molecules (the same function is used for cross-distribution), but we called updatecounters
  # only (n-1) times, so the actual sum of the distances considered is n(n-1). From this set
  # the number of distances returned in n_dref_in_bulk is r=n(n-1)-nd/2, where nd is the
  # number of distances in the domain. Thus, the number of distances in the domain, considering
  # symmetric terms, is nd=2n(n-1)-2r. The average number of distances in the domain, per
  # molecule, is thus nd/n=2(n-1)-2r/n. Finally, the number of solvent molecules in
  # bulk, for each molecule, is the total number of other molecules, (n-1), minus the
  # number of molecules in the domain, that is (n-1)-nd/n=(n-1)-2(n-1)+2r/n, which
  # finally simplifies to 2r/n-(n-1), which is the equation below.
  n_solvent_in_bulk = 2*n_solvent_in_bulk/solvent.nmols - (solvent.nmols-1)

  #
  # Computing the random-solvent distribution to compute the random minimum-distance count
  #
  for i in 1:options.n_random_samples
    # generate random solvent box, and store it in x_solvent_random
    for j in 1:solvent.nmols
      # Choose randomly one molecule from the bulk, if there are actually bulk molecules
      if n_dmin_in_bulk != 0
        jmol = dmin_mol[random(RNG,solvent.nmols-n_dmin_in_bulk+1:solvent.nmols)].jmol
      else
        jmol = random(RNG,1:solvent.nmols)
      end
      # Indexes of this molecule in the x_solvent array
      x_ref = viewmol(jmol,x_solvent,solvent)
      # Indexes of the random molecule in random array
      x_rnd = viewmol(j,x_solvent_random,solvent)
      # Generate new random coordinates (translation and rotation) for this molecule
      random_move!(x_ref,R.irefatom,sides,x_rnd,RNG)
    end

    # wrap random solvent coordinates to origin
    wrap!(x_solvent_random,sides)

    # Initialize linked cells
    initcells!(x_solvent_random,box,lc_solvent)

    # Choose randomly one solute molecule to be the solute in this sample
    i_rand_mol = random(RNG,1:solvent.nmols)
    x_this_solute = viewmol(i_rand_mol,x_solvent,solvent)

    # Compute all distances between solute and solvent atoms which are smaller than the
    # cutoff (this is the most computationally expensive part), the distances are returned
    # in the dc structure (here we do not use the self function)
    cutoffdistances!(R.cutoff,x_this_solute,x_solvent_random,lc_solvent,box,dc)

    # Update the counters and get the number of solvent molecules in bulk
    updatecounters!(R,rdf_count_random_frame,md_count_random_frame,
                    solvent,dc,dmin_mol,dref_mol)

  end # random solvent sampling

  # Update global counters with the data of this frame
  update_counters_frame!(R,rdf_count_random_frame,md_count_random_frame,volume_frame,
                         solute,solvent,n_solvent_in_bulk)

  nothing
end


text:
include("helperfunctions.jl")
include("randncm.jl")
include("ProjPSD.jl")
include("calllbfgsb.jl")
include("ncm.jl")
include("dualobj.jl")

function genprob(
    n,
    γ;
    seed = 0,
    gaussian_noise = false,
    memlim = 10,
    maxfgcalls = 100_000,
)
    U, G, H = randncm(
        n,
        γ = γ,
        seed = seed,
        gaussian_noise = gaussian_noise,
    )

    ncm = NCM(
        n,
        memlim = memlim,
        maxfgcalls = maxfgcalls,
    )

    return U, G, H, ncm
end

function runall(
    G,
    H,
    ncm;
    maxfgcalls = 100_000,
    tol = 1e-1,
    printlevel = 1,
    useXold = true,
)
    if useXold
        X, y = CorNewton3(G)
        ncm.Xold .= X
    end

    @time ncm(
        G,
        H,
        method = :IAPG,
        maxfgcalls = maxfgcalls,
        tol = tol,
        printlevel = printlevel,
        useXold = useXold,
    )

    if useXold
        ncm.Xold .= X
    end

    @time ncm(
        G,
        H,
        method = :IR,
        τ = 0.95,
        maxfgcalls = maxfgcalls,
        tol = tol,
        printlevel = printlevel,
        useXold = useXold,
    )

    if useXold
        ncm.Xold .= X
    end

    H2 = ncm.H2
    H2.data .= H .^ 2
    L = fronorm(H2, ncm.proj.work)

    @time ncm(
        G,
        H,
        method = :IER,
        α = 19/L,
        σ = 1.0,
        maxfgcalls = maxfgcalls,
        tol = tol,
        printlevel = printlevel,
        useXold = useXold,
    )

    return nothing
end

function tester(
    n,
    γ;
    maxfgcalls = 100_000,
    seed = 0,
    gaussian_noise = false,
    tol = 1e-1,
    printlevel = 1,
    useXold = true,
)
    U, G, H, ncm = genprob(
        n,
        γ,
        seed = seed,
        gaussian_noise = gaussian_noise,
        maxfgcalls = maxfgcalls,
    )
    runall(
        G,
        H,
        ncm,
        maxfgcalls = maxfgcalls,
        tol = tol,
        printlevel = printlevel,
        useXold = useXold,
    )
end

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
module NodeJS

export nodejs_cmd, npm_cmd

depsjl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl)
    error("NodeJS not properly installed. Please run\nPkg.build(\"NodeJS\")")
else
    include(depsjl)
end

"""
Return the full path of the node command.
"""
function nodejs_cmd()
    return `$node_executable_path`
end

"""
Return the full path of the npm script.
"""
function npm_script()
    return npm_script_path
end

"""
Return the full path of the npm command.
"""
function npm_cmd()
    return Sys.iswindows() ? Cmd(`$npm_executable_path`) : Cmd(`$(nodejs_cmd()) $(npm_script())`)
end

end # module


text:
module StrategyBase

  export filter_elements!, Strategy
  abstract type Strategy end
  function filter_elements!(s::Strategy)end

end


text:
# p38.jl - solve u_xxxx = exp(x), u(-1)=u(1)=u'(-1)=u'(1)=0
#         (compare p13.jl)

# Construct discrete biharmonic operator:
N = 15; (D,x) = cheb(N);
S = diagm([0; 1 ./(1-x[2:N].^2); 0]);
D4 = (diagm(1-x.^2)*D^4 - 8*diagm(x)*D^3 - 12*D^2)*S;
D4 = D4[2:N,2:N];

# Solve boundary-value problem and plot result:
f = @. exp(x[2:N]); u = D4\f; u = [0;u;0];
clf(); axes([.1,.4,.8,.5]);
plot(x,u,".",markersize=10);
axis([-1,1,-.01,.06]); grid(true);
xx = (-1:.01:1);
uu = (1-xx.^2).*polyval(polyfit(x,S*u),xx);
plot(xx,uu);

# Determine exact solution and print maximum error:
A = [1 -1 1 -1; 0 1 -2 3; 1 1 1 1; 0 1 2 3];
V = xx.^(0:3)';
c = A\exp.([-1,-1,1,1]); exact = exp.(xx) - V*c;
title("max err = $(signif(norm(uu-exact,Inf),5))",fontsize=12);


text:
module Zd_Arithmetics

using LinearAlgebra

import Base: +, -, *, convert, promote_rule, show, trailing_zeros, Int, AbstractFloat, inv

export

    ℤd,
    +,
    -,
    *,
    convert,
    promote_rule,
    show,
    trailing_zeros,
    Int,
    AbstractFloat,
    inv

include("datatype.jl")
include("operations.jl")

end # module


text:
# likelihood

logl(w, Ldx) = sum(log(Ldx*w))
dlogl(w, Ldx) = sum(Ldx./(Ldx*w), 1) |> vec

# entropy

repeatweights(w, zmult::Int) = repmat(w/zmult, zmult)
repeatweights(w, L::Matrix)  = repeatweights(w, Int(size(L,1) / length(w)))

# Monte Carlo approximation to the z-entropy
function hz(w, L, wz=repeatweights(w, L))
  rhoz = L * w           # \Int L(z|x) * pi(x) dx
  l = 0.
  for (r,w) in zip(rhoz, wz)
    r == 0 && continue
    l -= log(r)*w
  end
  l
end

# derivative of the MC approximation hz
function dhz(w, L, wz=repeatweights(w, L))
  nz, nx  = size(L)
  zmult   = Int(nz/nx)

  rhoz    = L * w
  factors = wz ./ rhoz

  d = zeros(w)
  @Threads.threads for k = 1:nx
    @simd for i = 1:nz
      @inbounds d[k] -= L[i,k] * factors[i]
    end
    @simd for m = 0:zmult-1
      @inbounds d[k] -= log(rhoz[k+nx*m]) / zmult
    end
  end
  d
end


text:
struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache
  k::uType
  half::uEltypeNoUnits
  c1::uEltypeNoUnits
  c2::uEltypeNoUnits
end
@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType} <: StochasticDiffEqMutableCache
  utmp::uType
  dutmp::uType
  k::uType
  gtmp::uType
  noise::rateNoiseType
  half::uEltypeNoUnits
  c1::uEltypeNoUnits
  c2::uEltypeNoUnits
end

function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype.x[1])
  c1 = exp(-alg.gamma*dt)
  c2 = sqrt(1 - c1^2)
  BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
end

function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  dutmp = zero(u.x[1])
  utmp = zero(u.x[2])
  k = zero(rate_prototype.x[1])

  gtmp = zero(rate_prototype.x[1])
  noise = zero(rate_prototype.x[1])

  half = uEltypeNoUnits(1//2)
  c1 = exp(-alg.gamma*dt)
  c2 = sqrt(1 - c1^2)

  BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2))
end


text:
# Mostly copy/pasted from StreamStats.jl

# v0.4 updates from
# https://github.com/PalladiumConsulting/StreamStats.jl/blob/master/src/hyper_log_log.jl

hash32(d::Any) = hash(d) % UInt32
maskadd32(x::UInt32, mask::UInt32, add::UInt32) = (x & mask) + add
ρ(s::UInt32) = UInt32(leading_zeros(s)) + 0x00000001
const toInt = Int
const toUInt = UInt

function α(m::UInt32)
    if m == 0x00000010
        return 0.673
    elseif m == 0x00000020
        return 0.697
    elseif m == 0x00000040
        return 0.709
    else # if m >= UInt32(128)
        return 0.7213 / (1 + 1.079 / m)
    end
end

"""
`HyperLogLog(b)`

Approximate count of distinct elements.  `HyperLogLog` differs from other OnlineStats
in that any input to `fit!(o::HyperLogLog, input)` is considered a singleton.  Thus,
a vector of inputs must be done by:

```julia
o = HyperLogLog(4)
for yi in y
    fit!(o, yi)
end
```
"""
type HyperLogLog <: OnlineStat
    m::UInt32
    M::Vector{UInt32}
    mask::UInt32
    altmask::UInt32
    n::Int
end

function HyperLogLog(b::Integer)
    if !(4 <= b <= 16)
        throw(ArgumentError("b must be an Integer between 4 and 16"))
    end

    m = 0x00000001 << b

    M = zeros(UInt32, m)

    mask = 0x00000000
    for i in 1:(b - 1)
        mask |= 0x00000001
        mask <<= 1
    end
    mask |= 0x00000001

    altmask = ~mask

    return HyperLogLog(m, M, mask, altmask, 0)
end

function Base.show(io::IO, counter::HyperLogLog)
    @printf(io, "A HyperLogLog counter w/ %d registers", Int(counter.m))
    return
end

function fit!(counter::HyperLogLog, v::Any)
    x = hash32(v)
    j = maskadd32(x, counter.mask, 0x00000001)
    w = x & counter.altmask
    counter.M[j] = max(counter.M[j], ρ(w))
    counter.n += 1
    return
end

nobs(o::HyperLogLog) = o.n

function value(counter::HyperLogLog)
    S = 0.0

    for j in 1:counter.m
        S += 1 / (2^counter.M[j])
    end

    Z = 1 / S

    E = α(counter.m) * toUInt(counter.m)^2 * Z

    if E <= 5//2 * counter.m
        V = 0
        for j in 1:counter.m
            V += toInt(counter.M[j] == 0x00000000)
        end
        if V != 0
            E_star = counter.m * log(counter.m / V)
        else
            E_star = E
        end
    elseif E <= 1//30 * 2^32
        E_star = E
    else
        E_star = -2^32 * log(1 - E / (2^32))
    end

    return E_star
end

# TODO: Figure out details here
# function confint(counter::HyperLogLog)
#     e = estimate(counter)
#     delta = e * 1.04 / sqrt(counter.m)
#     return e - delta, e + delta
# end


text:
abstract type AbstractILMProblem{DT,ST,DTP} end

function regenerate_problem() end

"""
$(TYPEDEF)

When defining a problem type with scalar data, make it a subtype of this.
"""
abstract type AbstractScalarILMProblem{DT,ST,DTP} <: AbstractILMProblem{DT,ST,DTP} end

"""
$(TYPEDEF)

When defining a problem type with vector data, make it a subtype of this.
"""
abstract type AbstractVectorILMProblem{DT,ST,DTP} <: AbstractILMProblem{DT,ST,DTP} end

"""
The `@ilmproblem` macro is used to automatically generate a type
particular to an immersed-layer problem, which can then be used for dispatch on
those types of problems. It takes two arguments: the name of the problem
(to which `Problem` will be automatically appended), and whether the problem
is of scalar or vector type. For example,

    @ilmproblem(StokesFlow,vector)

would generate a type `StokesFlowProblem` dealing with vector-valued data. The
resulting type then automatically has a constructor that allows one to pass
in the grid information and bodies, as well as optional choices for the
DDF function and the scaling type. For the example, this constructor would be

    StokesFlowProblem(grid,bodies[,ddftype=CartesianGrids.Yang3][,scaling=GridScaling])

Note that there is another constructor for problems with no surfaces that only
requires that the grid information be passed, e.g.,

    StokesFlowProblem(grid)

There are several keyword arguments for the problem constructor

- `ddftype = ` to set the DDF type. The default is `Yang3`.
- `scaling = ` to set the scaling type, `IndexScaling` (default) or `GridScaling`.
- `dtype = ` to set the element type to `Float64` (default) or `ComplexF64`.
- `phys_params = ` to pass in physical parameters
- `bc = ` to pass in boundary condition data or functions
- `forcing = ` to pass in forcing functions or data
- `motions = ` to provide function(s) that specify the velocity of the immersed surface(s). Note: if this keyword is used, it is assumed that surfaces will move.
- `timestep_func =` to pass in a function for time-dependent problems that provides the time-step size.
                  It is expected that this function takes in two arguments,
                  the `grid::PhysicalGrid` and `phys_params`, and returns the time step. It is up to the
                  user to decide how to determine this. It could also simply return a
                  constant value, regardless of the arguments.

"""
macro ilmproblem(name,vector_or_scalar)
    vs_string = string(vector_or_scalar)
    abtype = Symbol("Abstract"*uppercasefirst(vs_string)*"ILMProblem")
    typename = Symbol(string(name)*"Problem")
    return esc(quote

          @doc """
              $($typename)

          ILM problem type dealing with $($vs_string)-type data.
          """
          struct $typename{DT,ST,DTP,BLT,PHT,BCF,FF,DTF,MTF} <: $abtype{DT,ST,DTP}
             g :: PhysicalGrid
             bodies :: BLT
             phys_params :: PHT
             bc :: BCF
             forcing :: FF
             timestep_func :: DTF
             motions :: MTF
             $typename(g::PT,bodies::BodyList;ddftype=ImmersedLayers.DEFAULT_DDF,
                                              scaling=ImmersedLayers.DEFAULT_SCALING,
                                              dtype=ImmersedLayers.DEFAULT_DATA_TYPE,
                                              phys_params=nothing,
                                              bc=nothing,
                                              forcing=nothing,
                                              timestep_func=nothing,
                                              motions=nothing) where {PT <: PhysicalGrid} =
                    new{ddftype,scaling,dtype,typeof(bodies),typeof(phys_params),typeof(bc),typeof(forcing),typeof(timestep_func),
                                        typeof(ImmersedLayers._list(motions))}(
                                              g,bodies,phys_params,bc,forcing,timestep_func,ImmersedLayers._list(motions))
          end

          $typename(g::PhysicalGrid,body::Body;kwargs...) = $typename(g,BodyList([body]); kwargs...)
          $typename(g::PhysicalGrid;kwargs...) = $typename(g,BodyList(); kwargs...)

          ImmersedLayers.regenerate_problem(sys::ILMSystem{S,P},bl::BodyList) where {S,P<:$typename} =
                                      $typename(sys.base_cache.g,bl,
                                                ddftype=ImmersedLayers._ddf_type(P),
                                                scaling=ImmersedLayers._scaling_type(P),
                                                dtype=ImmersedLayers._element_type(P),
                                                phys_params=sys.phys_params,
                                                bc=sys.bc,
                                                forcing=sys.forcing,
                                                timestep_func=sys.timestep_func,
                                                motions=sys.motions)

           nothing
     end)

end

regenerate_problem(sys,body::Body) = regenerate_problem(sys,BodyList([body]))

_list(m::RigidBodyTools.AbstractMotion) = MotionList([m])
_list(m::MotionList) = m
_list(::Nothing) = nothing
_list(::Any) = nothing

_ddf_type(::AbstractILMProblem{DT}) where {DT} = DT
_ddf_type(::Type{P}) where P <: AbstractILMProblem{DT} where {DT} = DT
_scaling_type(::AbstractILMProblem{DT,ST}) where {DT,ST} = ST
_scaling_type(::Type{P}) where P <: AbstractILMProblem{DT,ST} where {DT,ST} = ST
_element_type(::AbstractILMProblem{DT,ST,DTP}) where {DT,ST,DTP} = DTP
_element_type(::Type{P}) where P <: AbstractILMProblem{DT,ST,DTP} where {DT,ST,DTP} = DTP

@ilmproblem BasicScalarILM scalar
@ilmproblem BasicVectorILM vector

# Extend this function for other problem types in order to create an extra cache
# of variables and operators for the problem
"""
    prob_cache(prob,base_cache::BasicILMCache)

This function is called by [`__init`](@ref) to generate a problem-specific extra
cache. Extend this function in order to generate an extra cache for a
user-defined problem type. The user must define the cache type itself.
""" prob_cache

prob_cache(prob::BasicScalarILMProblem,base_cache::BasicILMCache) = nothing
prob_cache(prob::BasicVectorILMProblem,base_cache::BasicILMCache) = nothing


text:
__precompile__()

module GrowableArrays
  using EllipsisNotation
  import Base: setindex!, getindex, push!, endof
  immutable GrowableArray{T,A,N} <: AbstractArray{T,N}
      data::Vector{A}
  end
  #=
  function GrowableArray(elem::Number)
    data = Vector{typeof(elem)}(0)
    push!(data,elem)
    GrowableArray{typeof(elem),typeof(elem),1}(data)
  end
  =#

  function GrowableArray(elem;initvalue=true)
    data = Vector{typeof(elem)}(0)
    if initvalue
      if typeof(elem) <: GrowableArray || typeof(elem) <: StackedArray #Copying GrowableArrays changes them!
        push!(data,elem)
      else
        push!(data,copy(elem))
      end
    end

    if typeof(elem) <: AbstractArray
      GrowableArray{eltype(elem),typeof(elem),ndims(elem)+1}(data)
    else
      GrowableArray{typeof(elem),typeof(elem),1}(data)
    end
  end
  Base.sizehint!(G::GrowableArray,i::Int) = sizehint!(G.data,i)
  Base.length(G::GrowableArray) = length(G.data)
  Base.push!(G::GrowableArray,Garr::GrowableArray) = push!(G.data,Garr) #Copying GrowableArrays changes them!
  Base.push!(G::GrowableArray,arr::AbstractArray)  = push!(G.data,copy(arr))
  Base.size(G::GrowableArray) = (length(G.data), size(G.data[1])...)
  Base.getindex(G::GrowableArray, i::Int) = G.data[i] # expand a linear index out
  Base.getindex(G::GrowableArray, i::Int, I::Int...) = G.data[i][I...]
  function Base.setindex!(G::GrowableArray, elem,i::Int) ##TODO: Add type checking on elem
    G.data[ind2sub(size(G), i)...] = elem
  end
  function Base.setindex!(G::GrowableArray, elem,i::Int,I::Int...)
    G.data[i][I...] = elem
  end
  function Base.push!(G::GrowableArray,elem)
    push!(G.data,elem)
  end

  immutable StackedArray{T,N,A} <: AbstractArray{T,N}
      data::A
      dims::NTuple{N,Int}
  end
  StackedArray(vec::AbstractVector) = StackedArray(vec, (length(vec), size(vec[1])...)) # TODO: ensure all elements are the same size
  StackedArray{A<:AbstractVector, N}(vec::A, dims::NTuple{N}) = StackedArray{eltype(eltype(A)),N,A}(vec, dims)
  Base.size(S::StackedArray) = S.dims
  Base.getindex(S::StackedArray, i::Int) = S.data[ind2sub(size(S), i)...] # expand a linear index out
  Base.getindex(S::StackedArray, i::Int, I::Int...) = S.data[i][I...]
  Base.push!(G::GrowableArray,Sarr::StackedArray)  = push!(G.data,Sarr.data)

  export StackedArray, GrowableArray, setindex!, getindex, push!, endof

end # module


text:
using CUDA
CUDA.versioninfo()
a = CuArray([1,2,3,4])
a .+= 1

task_local_storage(:ScalarIndexing, CUDA.GPUArrays.ScalarWarn)
for i in eachindex(a)
    a[i] += 1fend
a
CUDA.allowscalar(false)
a[1]
CUDA.@allowscalar a[1]
a = CuArray{Float32}(undef, 10)
fill!(a, 42)
c = CUDA.ones(length(a))
using LinearAlgebra
mul!(d, CUDA.rand(10,10), d)
out = Ref{Cint}()
CUBLAS.cublasIsamin_v2(CUBLAS.handle(), length(d), d, stride(d, 1), out)
out[]

map(a) do x
    x + 1
end

reduce(+, a)

accumulate(+, a)

broadcast(a, eachcol(a)) do x,Y
    x + sum(Y)
end

using Tullio
using CUDA, CUDAKernels, KernelAbstractions

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

path = "examples/examples_files/doublependulum.urdf"
mech = Mechanism(path)

storage = simulate!(mech, 10., record = true)
visualize(mech, storage)


text:
include("include.jl")
include("read_lp.jl")

#println("GUROBI NOT WORKING")
#using Gurobi

### SHOULD REALLY BE LOOKING AT DUAL
###

#for tol in [1e-8]

function lp_feasible(tol, problem_list)
  status_list = Array{String,1}()
  for problem_name in problem_list
      mod = read_lp_into_JuMP(problem_name, -tol, false)
      setsolver(mod, GurobiSolver(FeasibilityTol=1e-9,BarHomogeneous=1))
      JuMP.build(mod)
      status =solve(mod)
      string_status = (status == :Optimal) ? "true" : "false"
      push!(status_list, string_status)
  end

  return status_list
end
tol_set = [10.0, 1e-4, 1e-6, 1e-8]
status_list = Dict{Float64,Array{String,1}}()
problem_list = load_netlib(Int(1e5))

for tol in tol_set
  status_list[tol] = lp_feasible(tol, problem_list)
end

#end
#using CSV
#CSV.write("seq-ipm-figure/table.csv",results_table)
# Pkg.clone("https://github.com/farrellm/YTables.jl.git")
using DataFrames, YTables

@show interior, no_interior

status_list = results[tol]
divide = Int(length(status_list)/2)

#tol1 =

tol1 = 1e-4
tol2 = 1e-6
tol3 = 1e-8
status11 = status_list[tol1][1:divide]
status12 = status_list[tol2][1:divide]
names1 = problem_list[1:divide]
status21 = status_list[tol1][(divide+1):end]
status22 = status_list[tol2][(divide+1):end]
names2 = problem_list[(divide+1):end]

df = DataFrame(problem_name1 = names1, feasible11 = status11, problem_name2 = names2, feasible21 = status21 )

#writetable("seq-ipm-figures/table.csv", df)

latex(df)

interior = Dict()
no_interior = Dict()
for tol in tol_set
  interior[tol] = sum("true" .== status_list[tol])
  no_interior[tol] = sum("false" .== status_list[tol])
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
module Utils

import ACE

import ACE: PolyTransform, transformed_jacobi, Rn1pBasis,
            NaiveTotalDegree, init1pspec!, Ylm1pBasis,
            Product1pBasis

# import ACE.PairPotentials: PolyPairBasis

# - simple ways to construct a radial basis
# - construct a descriptor
# - simple wrappers to generate RPI basis functions (ACE + relatives)

@doc raw"""
`Rn_basis` construct a ``R_n`` basis; all arguments are keyword arguments with
defaults:
* `r0 = 1.0`
* `trans = PolyTransform(2, r0)`
* `D = NaiveTotalDegree()`
* `maxdeg = 6`
* `rcut = 2.5`
* `rin = 0.5 * r0`
* `pcut = 2`
* `pin = 0`
* `constants = false`
"""
function Rn_basis(;
      # transform parameters
      r0 = 1.0,
      trans = PolyTransform(2, r0),
      # degree parameters
      D = NaiveTotalDegree(),
      maxdeg = 6,
      # radial basis parameters
      rcut = 2.5,
      rin = 0.5 * r0,
      pcut = 2,
      pin = 0,
      constants = false)

   J = transformed_jacobi(maxdeg, trans, rcut, rin; pcut=pcut, pin=pin)
   return Rn1pBasis(J)
end

@doc raw"""
Construct a ``R_n * Y_l^m`` 1-particle basis.
All arguments are keyword arguments; see documentation of `ACE.Utils.Rn_basis`.
"""
function RnYlm_1pbasis(; maxdeg=6, kwargs...)
   Rn = Rn_basis(; maxdeg = maxdeg, kwargs...)
   Ylm = Ylm1pBasis(maxdeg)
   B1p = ACE.Product1pBasis((Rn, Ylm))
   init1pspec!(B1p)
   return B1p
end

invariant_basis(; kwargs...) =
      symm_basis(ACE.Invariant(); kwargs...)

symm_basis(φ; maxν = 3, maxdeg = 6, kwargs...) =
      ACE.SymmetricBasis(φ,
                         RnYlm_1pbasis(; maxdeg=maxdeg, kwargs...),
                         maxν,
                         maxdeg)

# function rpi_basis(; species = :X, N = 3,
#       # transform parameters
#       r0 = 2.5,
#       trans = PolyTransform(2, r0),
#       # degree parameters
#       D = SparsePSHDegree(),
#       maxdeg = 8,
#       # radial basis parameters
#       rcut = 5.0,
#       rin = 0.5 * r0,
#       pcut = 2,
#       pin = 0,
#       constants = false,
#       rbasis = transformed_jacobi(get_maxn(D, maxdeg, species), trans, rcut, rin;
#                                   pcut=pcut, pin=pin),
#       # one-particle basis
#       Basis1p = RnYlm1pBasis,
#       basis1p = Basis1p(rbasis; species = species, D = D) )
#
#    return RPIBasis(basis1p, N, D, maxdeg, constants)
# end
#
# descriptor = rpi_basis
# ace_basis = rpi_basis
#
# function pair_basis(; species = :X,
#       # transform parameters
#       r0 = 2.5,
#       trans = PolyTransform(2, r0),
#       # degree parameters
#       maxdeg = 8,
#       # radial basis parameters
#       rcut = 5.0,
#       rin = 0.5 * r0,
#       pcut = 2,
#       pin = 0,
#       rbasis = transformed_jacobi(maxdeg, trans, rcut, rin; pcut=pcut, pin=pin))
#
#    return PolyPairBasis(rbasis, species)
# end

end


text:
# Sample postprocessing script for virtual satellite trackings.
# Simultaneous virtual satellite trackings at uniformly sampled locations on a 2D plane.
# Outputs are stored in binary format for sharing within Julia.
#
# Hongyang Zhou, hyzhou@umich.edu

using Glob, Vlasiator
using JLD2: jldsave

"Select cells in 2D `meta` with uniform distance `dx`."
function sample(meta, dx)
   dcell = floor(Int, dx ÷ meta.dcoord[1])

   cellid = meta.cellid[meta.cellindex]
   cellid = reshape(cellid, meta.ncells[1], meta.ncells[2])
   cellid_select = cellid[1+dcell:dcell:end-dcell, 1+dcell:dcell:end-dcell]
end

function extract_vars(files, dx)
   nfiles = length(files)

   meta = load(files[1])
   ids = sample(meta, dx)
   locations = [getcellcoordinates(meta, id) for id in ids]

   nsize = size(ids)
   println("Number of virtual satellites: $(length(ids))")

   # variables to be extracted
   t   = zeros(Float32, nfiles)
   rho = zeros(Float32, nfiles, nsize[1], nsize[2])
   vx  = zeros(Float32, nfiles, nsize[1], nsize[2])
   vy  = zeros(Float32, nfiles, nsize[1], nsize[2])
   vz  = zeros(Float32, nfiles, nsize[1], nsize[2])
   p   = zeros(Float32, nfiles, nsize[1], nsize[2])
   bx  = zeros(Float32, nfiles, nsize[1], nsize[2])
   by  = zeros(Float32, nfiles, nsize[1], nsize[2])
   bz  = zeros(Float32, nfiles, nsize[1], nsize[2])
   ex  = zeros(Float32, nfiles, nsize[1], nsize[2])
   ey  = zeros(Float32, nfiles, nsize[1], nsize[2])
   ez  = zeros(Float32, nfiles, nsize[1], nsize[2])

   # Extract data from each frame
   for i = eachindex(files)
      meta = load(files[i])
      t[i] = meta.time
      rho[i,:,:] = readvariable(meta, "proton/vg_rho", ids)
      v = readvariable(meta, "proton/vg_v", ids)
      vx[i,:,:] = v[1,:]
      vy[i,:,:] = v[2,:]
      vz[i,:,:] = v[3,:]
      p[i,:,:] = readvariable(meta, "vg_pressure", ids)
      b = readvariable(meta, "vg_b_vol", ids)
      bx[i,:,:] = b[1,:]
      by[i,:,:] = b[2,:]
      bz[i,:,:] = b[3,:]
      e = readvariable(meta, "vg_e_vol", ids)
      ex[i,:,:] = e[1,:]
      ey[i,:,:] = e[2,:]
      ez[i,:,:] = e[3,:]
   end

   # Save into binary file
   jldsave("satellites_uniform_sampled.jld2";
      locations, t, rho, vx, vy, vz, p, bx, by, bz, ex, ey, ez)
end

#####

files = glob("bulk*.vlsv", "./")

dx = 5RE # uniform sampling distance [m]

println("Number of files: $(length(files))")
println("Running with $(Threads.nthreads()) threads...")

@time extract_vars(files, dx)

println("Virtual satellite extraction done!")


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: kinesis_analytics

using AWS.Compat
using AWS.UUIDs
"""
    AddApplicationCloudWatchLoggingOption()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Adds a CloudWatch log stream to monitor application configuration errors. For more information about using CloudWatch log streams with Amazon Kinesis Analytics applications, see Working with Amazon CloudWatch Logs.

# Required Parameters
- `ApplicationName`: The Kinesis Analytics application name.
- `CloudWatchLoggingOption`: Provides the CloudWatch log stream Amazon Resource Name (ARN) and the IAM role ARN. Note: To write application messages to CloudWatch, the IAM role that is used must have the PutLogEvents policy action enabled.
- `CurrentApplicationVersionId`: The version ID of the Kinesis Analytics application.

"""

add_application_cloud_watch_logging_option(ApplicationName, CloudWatchLoggingOption, CurrentApplicationVersionId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationCloudWatchLoggingOption", Dict{String, Any}("ApplicationName"=>ApplicationName, "CloudWatchLoggingOption"=>CloudWatchLoggingOption, "CurrentApplicationVersionId"=>CurrentApplicationVersionId); aws_config=aws_config)
add_application_cloud_watch_logging_option(ApplicationName, CloudWatchLoggingOption, CurrentApplicationVersionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationCloudWatchLoggingOption", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CloudWatchLoggingOption"=>CloudWatchLoggingOption, "CurrentApplicationVersionId"=>CurrentApplicationVersionId), args)); aws_config=aws_config)

"""
    AddApplicationInput()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.   Adds a streaming source to your Amazon Kinesis application. For conceptual information, see Configuring Application Input.  You can add a streaming source either when you create an application or you can use this operation to add a streaming source after you create an application. For more information, see CreateApplication. Any configuration update, including adding a streaming source using this operation, results in a new version of the application. You can use the DescribeApplication operation to find the current application version.  This operation requires permissions to perform the kinesisanalytics:AddApplicationInput action.

# Required Parameters
- `ApplicationName`: Name of your existing Amazon Kinesis Analytics application to which you want to add the streaming source.
- `CurrentApplicationVersionId`: Current version of your Amazon Kinesis Analytics application. You can use the DescribeApplication operation to find the current application version.
- `Input`: The Input to add.

"""

add_application_input(ApplicationName, CurrentApplicationVersionId, Input; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationInput", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "Input"=>Input); aws_config=aws_config)
add_application_input(ApplicationName, CurrentApplicationVersionId, Input, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationInput", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "Input"=>Input), args)); aws_config=aws_config)

"""
    AddApplicationInputProcessingConfiguration()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Adds an InputProcessingConfiguration to an application. An input processor preprocesses records on the input stream before the application's SQL code executes. Currently, the only input processor available is AWS Lambda.

# Required Parameters
- `ApplicationName`: Name of the application to which you want to add the input processing configuration.
- `CurrentApplicationVersionId`: Version of the application to which you want to add the input processing configuration. You can use the DescribeApplication operation to get the current application version. If the version specified is not the current version, the ConcurrentModificationException is returned.
- `InputId`: The ID of the input configuration to add the input processing configuration to. You can get a list of the input IDs for an application using the DescribeApplication operation.
- `InputProcessingConfiguration`: The InputProcessingConfiguration to add to the application.

"""

add_application_input_processing_configuration(ApplicationName, CurrentApplicationVersionId, InputId, InputProcessingConfiguration; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationInputProcessingConfiguration", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "InputId"=>InputId, "InputProcessingConfiguration"=>InputProcessingConfiguration); aws_config=aws_config)
add_application_input_processing_configuration(ApplicationName, CurrentApplicationVersionId, InputId, InputProcessingConfiguration, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationInputProcessingConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "InputId"=>InputId, "InputProcessingConfiguration"=>InputProcessingConfiguration), args)); aws_config=aws_config)

"""
    AddApplicationOutput()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Adds an external destination to your Amazon Kinesis Analytics application. If you want Amazon Kinesis Analytics to deliver data from an in-application stream within your application to an external destination (such as an Amazon Kinesis stream, an Amazon Kinesis Firehose delivery stream, or an AWS Lambda function), you add the relevant configuration to your application using this operation. You can configure one or more outputs for your application. Each output configuration maps an in-application stream and an external destination.  You can use one of the output configurations to deliver data from your in-application error stream to an external destination so that you can analyze the errors. For more information, see Understanding Application Output (Destination).   Any configuration update, including adding a streaming source using this operation, results in a new version of the application. You can use the DescribeApplication operation to find the current application version. For the limits on the number of application inputs and outputs you can configure, see Limits. This operation requires permissions to perform the kinesisanalytics:AddApplicationOutput action.

# Required Parameters
- `ApplicationName`: Name of the application to which you want to add the output configuration.
- `CurrentApplicationVersionId`: Version of the application to which you want to add the output configuration. You can use the DescribeApplication operation to get the current application version. If the version specified is not the current version, the ConcurrentModificationException is returned.
- `Output`: An array of objects, each describing one output configuration. In the output configuration, you specify the name of an in-application stream, a destination (that is, an Amazon Kinesis stream, an Amazon Kinesis Firehose delivery stream, or an AWS Lambda function), and record the formation to use when writing to the destination.

"""

add_application_output(ApplicationName, CurrentApplicationVersionId, Output; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationOutput", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "Output"=>Output); aws_config=aws_config)
add_application_output(ApplicationName, CurrentApplicationVersionId, Output, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationOutput", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "Output"=>Output), args)); aws_config=aws_config)

"""
    AddApplicationReferenceDataSource()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Adds a reference data source to an existing application. Amazon Kinesis Analytics reads reference data (that is, an Amazon S3 object) and creates an in-application table within your application. In the request, you provide the source (S3 bucket name and object key name), name of the in-application table to create, and the necessary mapping information that describes how data in Amazon S3 object maps to columns in the resulting in-application table.  For conceptual information, see Configuring Application Input. For the limits on data sources you can add to your application, see Limits.   This operation requires permissions to perform the kinesisanalytics:AddApplicationOutput action.

# Required Parameters
- `ApplicationName`: Name of an existing application.
- `CurrentApplicationVersionId`: Version of the application for which you are adding the reference data source. You can use the DescribeApplication operation to get the current application version. If the version specified is not the current version, the ConcurrentModificationException is returned.
- `ReferenceDataSource`: The reference data source can be an object in your Amazon S3 bucket. Amazon Kinesis Analytics reads the object and copies the data into the in-application table that is created. You provide an S3 bucket, object key name, and the resulting in-application table that is created. You must also provide an IAM role with the necessary permissions that Amazon Kinesis Analytics can assume to read the object from your S3 bucket on your behalf.

"""

add_application_reference_data_source(ApplicationName, CurrentApplicationVersionId, ReferenceDataSource; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationReferenceDataSource", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "ReferenceDataSource"=>ReferenceDataSource); aws_config=aws_config)
add_application_reference_data_source(ApplicationName, CurrentApplicationVersionId, ReferenceDataSource, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("AddApplicationReferenceDataSource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "ReferenceDataSource"=>ReferenceDataSource), args)); aws_config=aws_config)

"""
    CreateApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.   Creates an Amazon Kinesis Analytics application. You can configure each application with one streaming source as input, application code to process the input, and up to three destinations where you want Amazon Kinesis Analytics to write the output data from your application. For an overview, see How it Works.  In the input configuration, you map the streaming source to an in-application stream, which you can think of as a constantly updating table. In the mapping, you must provide a schema for the in-application stream and map each data column in the in-application stream to a data element in the streaming source. Your application code is one or more SQL statements that read input data, transform it, and generate output. Your application code can create one or more SQL artifacts like SQL streams or pumps. In the output configuration, you can configure the application to write data from in-application streams created in your applications to up to three destinations.  To read data from your source stream or write data to destination streams, Amazon Kinesis Analytics needs your permissions. You grant these permissions by creating IAM roles. This operation requires permissions to perform the kinesisanalytics:CreateApplication action.   For introductory exercises to create an Amazon Kinesis Analytics application, see Getting Started.

# Required Parameters
- `ApplicationName`: Name of your Amazon Kinesis Analytics application (for example, sample-app).

# Optional Parameters
- `ApplicationCode`: One or more SQL statements that read input data, transform it, and generate output. For example, you can write a SQL statement that reads data from one in-application stream, generates a running average of the number of advertisement clicks by vendor, and insert resulting rows in another in-application stream using pumps. For more information about the typical pattern, see Application Code.  You can provide such series of SQL statements, where output of one statement can be used as the input for the next statement. You store intermediate results by creating in-application streams and pumps. Note that the application code must create the streams with names specified in the Outputs. For example, if your Outputs defines output streams named ExampleOutputStream1 and ExampleOutputStream2, then your application code must create these streams.
- `ApplicationDescription`: Summary description of the application.
- `CloudWatchLoggingOptions`: Use this parameter to configure a CloudWatch log stream to monitor application configuration errors. For more information, see Working with Amazon CloudWatch Logs.
- `Inputs`: Use this parameter to configure the application input. You can configure your application to receive input from a single streaming source. In this configuration, you map this streaming source to an in-application stream that is created. Your application code can then query the in-application stream like a table (you can think of it as a constantly updating table). For the streaming source, you provide its Amazon Resource Name (ARN) and format of data on the stream (for example, JSON, CSV, etc.). You also must provide an IAM role that Amazon Kinesis Analytics can assume to read this stream on your behalf. To create the in-application stream, you need to specify a schema to transform your data into a schematized version used in SQL. In the schema, you provide the necessary mapping of the data elements in the streaming source to record columns in the in-app stream.
- `Outputs`: You can configure application output to write data from any of the in-application streams to up to three destinations. These destinations can be Amazon Kinesis streams, Amazon Kinesis Firehose delivery streams, AWS Lambda destinations, or any combination of the three. In the configuration, you specify the in-application stream name, the destination stream or Lambda function Amazon Resource Name (ARN), and the format to use when writing data. You must also provide an IAM role that Amazon Kinesis Analytics can assume to write to the destination stream or Lambda function on your behalf. In the output configuration, you also provide the output stream or Lambda function ARN. For stream destinations, you provide the format of data in the stream (for example, JSON, CSV). You also must provide an IAM role that Amazon Kinesis Analytics can assume to write to the stream or Lambda function on your behalf.
- `Tags`: A list of one or more tags to assign to the application. A tag is a key-value pair that identifies an application. Note that the maximum number of application tags includes system tags. The maximum number of user-defined application tags is 50. For more information, see Using Tagging.
"""

create_application(ApplicationName; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("CreateApplication", Dict{String, Any}("ApplicationName"=>ApplicationName); aws_config=aws_config)
create_application(ApplicationName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("CreateApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName), args)); aws_config=aws_config)

"""
    DeleteApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Deletes the specified application. Amazon Kinesis Analytics halts application execution and deletes the application, including any application artifacts (such as in-application streams, reference table, and application code). This operation requires permissions to perform the kinesisanalytics:DeleteApplication action.

# Required Parameters
- `ApplicationName`: Name of the Amazon Kinesis Analytics application to delete.
- `CreateTimestamp`:  You can use the DescribeApplication operation to get this value.

"""

delete_application(ApplicationName, CreateTimestamp; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplication", Dict{String, Any}("ApplicationName"=>ApplicationName, "CreateTimestamp"=>CreateTimestamp); aws_config=aws_config)
delete_application(ApplicationName, CreateTimestamp, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CreateTimestamp"=>CreateTimestamp), args)); aws_config=aws_config)

"""
    DeleteApplicationCloudWatchLoggingOption()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Deletes a CloudWatch log stream from an application. For more information about using CloudWatch log streams with Amazon Kinesis Analytics applications, see Working with Amazon CloudWatch Logs.

# Required Parameters
- `ApplicationName`: The Kinesis Analytics application name.
- `CloudWatchLoggingOptionId`: The CloudWatchLoggingOptionId of the CloudWatch logging option to delete. You can get the CloudWatchLoggingOptionId by using the DescribeApplication operation.
- `CurrentApplicationVersionId`: The version ID of the Kinesis Analytics application.

"""

delete_application_cloud_watch_logging_option(ApplicationName, CloudWatchLoggingOptionId, CurrentApplicationVersionId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationCloudWatchLoggingOption", Dict{String, Any}("ApplicationName"=>ApplicationName, "CloudWatchLoggingOptionId"=>CloudWatchLoggingOptionId, "CurrentApplicationVersionId"=>CurrentApplicationVersionId); aws_config=aws_config)
delete_application_cloud_watch_logging_option(ApplicationName, CloudWatchLoggingOptionId, CurrentApplicationVersionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationCloudWatchLoggingOption", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CloudWatchLoggingOptionId"=>CloudWatchLoggingOptionId, "CurrentApplicationVersionId"=>CurrentApplicationVersionId), args)); aws_config=aws_config)

"""
    DeleteApplicationInputProcessingConfiguration()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Deletes an InputProcessingConfiguration from an input.

# Required Parameters
- `ApplicationName`: The Kinesis Analytics application name.
- `CurrentApplicationVersionId`: The version ID of the Kinesis Analytics application.
- `InputId`: The ID of the input configuration from which to delete the input processing configuration. You can get a list of the input IDs for an application by using the DescribeApplication operation.

"""

delete_application_input_processing_configuration(ApplicationName, CurrentApplicationVersionId, InputId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationInputProcessingConfiguration", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "InputId"=>InputId); aws_config=aws_config)
delete_application_input_processing_configuration(ApplicationName, CurrentApplicationVersionId, InputId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationInputProcessingConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "InputId"=>InputId), args)); aws_config=aws_config)

"""
    DeleteApplicationOutput()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Deletes output destination configuration from your application configuration. Amazon Kinesis Analytics will no longer write data from the corresponding in-application stream to the external output destination. This operation requires permissions to perform the kinesisanalytics:DeleteApplicationOutput action.

# Required Parameters
- `ApplicationName`: Amazon Kinesis Analytics application name.
- `CurrentApplicationVersionId`: Amazon Kinesis Analytics application version. You can use the DescribeApplication operation to get the current application version. If the version specified is not the current version, the ConcurrentModificationException is returned.
- `OutputId`: The ID of the configuration to delete. Each output configuration that is added to the application, either when the application is created or later using the AddApplicationOutput operation, has a unique ID. You need to provide the ID to uniquely identify the output configuration that you want to delete from the application configuration. You can use the DescribeApplication operation to get the specific OutputId.

"""

delete_application_output(ApplicationName, CurrentApplicationVersionId, OutputId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationOutput", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "OutputId"=>OutputId); aws_config=aws_config)
delete_application_output(ApplicationName, CurrentApplicationVersionId, OutputId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationOutput", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "OutputId"=>OutputId), args)); aws_config=aws_config)

"""
    DeleteApplicationReferenceDataSource()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Deletes a reference data source configuration from the specified application configuration. If the application is running, Amazon Kinesis Analytics immediately removes the in-application table that you created using the AddApplicationReferenceDataSource operation.  This operation requires permissions to perform the kinesisanalytics.DeleteApplicationReferenceDataSource action.

# Required Parameters
- `ApplicationName`: Name of an existing application.
- `CurrentApplicationVersionId`: Version of the application. You can use the DescribeApplication operation to get the current application version. If the version specified is not the current version, the ConcurrentModificationException is returned.
- `ReferenceId`: ID of the reference data source. When you add a reference data source to your application using the AddApplicationReferenceDataSource, Amazon Kinesis Analytics assigns an ID. You can use the DescribeApplication operation to get the reference ID.

"""

delete_application_reference_data_source(ApplicationName, CurrentApplicationVersionId, ReferenceId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationReferenceDataSource", Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "ReferenceId"=>ReferenceId); aws_config=aws_config)
delete_application_reference_data_source(ApplicationName, CurrentApplicationVersionId, ReferenceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DeleteApplicationReferenceDataSource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "CurrentApplicationVersionId"=>CurrentApplicationVersionId, "ReferenceId"=>ReferenceId), args)); aws_config=aws_config)

"""
    DescribeApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Returns information about a specific Amazon Kinesis Analytics application. If you want to retrieve a list of all applications in your account, use the ListApplications operation. This operation requires permissions to perform the kinesisanalytics:DescribeApplication action. You can use DescribeApplication to get the current application versionId, which you need to call other operations such as Update.

# Required Parameters
- `ApplicationName`: Name of the application.

"""

describe_application(ApplicationName; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DescribeApplication", Dict{String, Any}("ApplicationName"=>ApplicationName); aws_config=aws_config)
describe_application(ApplicationName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DescribeApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName), args)); aws_config=aws_config)

"""
    DiscoverInputSchema()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Infers a schema by evaluating sample records on the specified streaming source (Amazon Kinesis stream or Amazon Kinesis Firehose delivery stream) or S3 object. In the response, the operation returns the inferred schema and also the sample records that the operation used to infer the schema.  You can use the inferred schema when configuring a streaming source for your application. For conceptual information, see Configuring Application Input. Note that when you create an application using the Amazon Kinesis Analytics console, the console uses this operation to infer a schema and show it in the console user interface.   This operation requires permissions to perform the kinesisanalytics:DiscoverInputSchema action.

# Optional Parameters
- `InputProcessingConfiguration`: The InputProcessingConfiguration to use to preprocess the records before discovering the schema of the records.
- `InputStartingPositionConfiguration`: Point at which you want Amazon Kinesis Analytics to start reading records from the specified streaming source discovery purposes.
- `ResourceARN`: Amazon Resource Name (ARN) of the streaming source.
- `RoleARN`: ARN of the IAM role that Amazon Kinesis Analytics can assume to access the stream on your behalf.
- `S3Configuration`: Specify this parameter to discover a schema from data in an Amazon S3 object.
"""

discover_input_schema(; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DiscoverInputSchema"; aws_config=aws_config)
discover_input_schema(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("DiscoverInputSchema", args; aws_config=aws_config)

"""
    ListApplications()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Returns a list of Amazon Kinesis Analytics applications in your account. For each application, the response includes the application name, Amazon Resource Name (ARN), and status. If the response returns the HasMoreApplications value as true, you can send another request by adding the ExclusiveStartApplicationName in the request body, and set the value of this to the last application name from the previous response.  If you want detailed information about a specific application, use DescribeApplication. This operation requires permissions to perform the kinesisanalytics:ListApplications action.

# Optional Parameters
- `ExclusiveStartApplicationName`: Name of the application to start the list with. When using pagination to retrieve the list, you don't need to specify this parameter in the first request. However, in subsequent requests, you add the last application name from the previous response to get the next page of applications.
- `Limit`: Maximum number of applications to list.
"""

list_applications(; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("ListApplications"; aws_config=aws_config)
list_applications(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("ListApplications", args; aws_config=aws_config)

"""
    ListTagsForResource()

Retrieves the list of key-value tags assigned to the application. For more information, see Using Tagging.

# Required Parameters
- `ResourceARN`: The ARN of the application for which to retrieve tags.

"""

list_tags_for_resource(ResourceARN; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("ListTagsForResource", Dict{String, Any}("ResourceARN"=>ResourceARN); aws_config=aws_config)
list_tags_for_resource(ResourceARN, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN), args)); aws_config=aws_config)

"""
    StartApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Starts the specified Amazon Kinesis Analytics application. After creating an application, you must exclusively call this operation to start your application. After the application starts, it begins consuming the input data, processes it, and writes the output to the configured destination.  The application status must be READY for you to start an application. You can get the application status in the console or using the DescribeApplication operation. After you start the application, you can stop the application from processing the input by calling the StopApplication operation. This operation requires permissions to perform the kinesisanalytics:StartApplication action.

# Required Parameters
- `ApplicationName`: Name of the application.
- `InputConfigurations`: Identifies the specific input, by ID, that the application starts consuming. Amazon Kinesis Analytics starts reading the streaming source associated with the input. You can also specify where in the streaming source you want Amazon Kinesis Analytics to start reading.

"""

start_application(ApplicationName, InputConfigurations; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("StartApplication", Dict{String, Any}("ApplicationName"=>ApplicationName, "InputConfigurations"=>InputConfigurations); aws_config=aws_config)
start_application(ApplicationName, InputConfigurations, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("StartApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "InputConfigurations"=>InputConfigurations), args)); aws_config=aws_config)

"""
    StopApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Stops the application from processing input data. You can stop an application only if it is in the running state. You can use the DescribeApplication operation to find the application state. After the application is stopped, Amazon Kinesis Analytics stops reading data from the input, the application stops processing data, and there is no output written to the destination.  This operation requires permissions to perform the kinesisanalytics:StopApplication action.

# Required Parameters
- `ApplicationName`: Name of the running application to stop.

"""

stop_application(ApplicationName; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("StopApplication", Dict{String, Any}("ApplicationName"=>ApplicationName); aws_config=aws_config)
stop_application(ApplicationName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("StopApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName), args)); aws_config=aws_config)

"""
    TagResource()

Adds one or more key-value tags to a Kinesis Analytics application. Note that the maximum number of application tags includes system tags. The maximum number of user-defined application tags is 50. For more information, see Using Tagging.

# Required Parameters
- `ResourceARN`: The ARN of the application to assign the tags.
- `Tags`: The key-value tags to assign to the application.

"""

tag_resource(ResourceARN, Tags; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("TagResource", Dict{String, Any}("ResourceARN"=>ResourceARN, "Tags"=>Tags); aws_config=aws_config)
tag_resource(ResourceARN, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes one or more tags from a Kinesis Analytics application. For more information, see Using Tagging.

# Required Parameters
- `ResourceARN`: The ARN of the Kinesis Analytics application from which to remove the tags.
- `TagKeys`: A list of keys of tags to remove from the specified application.

"""

untag_resource(ResourceARN, TagKeys; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("UntagResource", Dict{String, Any}("ResourceARN"=>ResourceARN, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(ResourceARN, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN, "TagKeys"=>TagKeys), args)); aws_config=aws_config)

"""
    UpdateApplication()

 This documentation is for version 1 of the Amazon Kinesis Data Analytics API, which only supports SQL applications. Version 2 of the API supports SQL and Java applications. For more information about version 2, see Amazon Kinesis Data Analytics API V2 Documentation.  Updates an existing Amazon Kinesis Analytics application. Using this API, you can update application code, input configuration, and output configuration.  Note that Amazon Kinesis Analytics updates the CurrentApplicationVersionId each time you update your application.  This operation requires permission for the kinesisanalytics:UpdateApplication action.

# Required Parameters
- `ApplicationName`: Name of the Amazon Kinesis Analytics application to update.
- `ApplicationUpdate`: Describes application updates.
- `CurrentApplicationVersionId`: The current application version ID. You can use the DescribeApplication operation to get this value.

"""

update_application(ApplicationName, ApplicationUpdate, CurrentApplicationVersionId; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("UpdateApplication", Dict{String, Any}("ApplicationName"=>ApplicationName, "ApplicationUpdate"=>ApplicationUpdate, "CurrentApplicationVersionId"=>CurrentApplicationVersionId); aws_config=aws_config)
update_application(ApplicationName, ApplicationUpdate, CurrentApplicationVersionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_analytics("UpdateApplication", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationName"=>ApplicationName, "ApplicationUpdate"=>ApplicationUpdate, "CurrentApplicationVersionId"=>CurrentApplicationVersionId), args)); aws_config=aws_config)


text:
module TestTimeResp
using CustomTest
using Base.Test
using ControlSystems

A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = eye(2)
R = eye(1)
L = lqr(sys,Q,R)

u(i,x) = -L*x # Form control law
t=0:0.1:50
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0)
@test sum(abs(x[end,:])) < eps()

#Do a manual simulation with uout
ym, tm, xm = lsim(sys, uout, t, x0)
@test_approx_eq y ym
@test_approx_eq x xm

# Now compare to closed loop
# Discretization is needed before feedback
sysd = c2d(sys, 0.1)[1]
# Create the closed loop system
sysdfb = ss(sysd.A-sysd.B*L, sysd.B, sysd.C, sysd.D, 0.1)
#Simulate without input
yd, td, xd = lsim(sysdfb, zeros(501), t, x0)

@test_approx_eq y yd
@test_approx_eq x xd

#Test step and impulse
t0 = 0:0.05:2
systf = [tf(1,[1,1]) 0; 0 tf([1,-1],[1,2,1])]
sysss = ss([-1 0 0; 0 -2 -1; 0 1 0], [1 0; 0 1; 0 0], [1 0 0; 0 1 -1], 0)
yreal = zeros(length(t0), 2, 2)
xreal = zeros(length(t0), 3, 2)

#Step tf
y, t, x = step(systf, t0)
yreal[:,1,1] = 1-exp(-t)
yreal[:,2,2] = -1+exp(-t)+2*exp(-t).*t
@test_approx_eq_eps y yreal 1e-14
#Step ss
y, t, x = step(sysss, t)
@test_approx_eq_eps y yreal 1e-14
xreal[:,1,1] = yreal[:,1,1]
xreal[:,2,2] = exp(-t).*t
xreal[:,3,2] = exp(-t).*(-t-1) + 1
@test_approx_eq_eps x xreal 1e-14

#Impulse tf
y, t, x = impulse(systf, t)
yreal[:,1,1] = exp(-t)
yreal[:,2,2] = exp(-t).*(1 - 2.*t)
@test_approx_eq_eps y yreal 1e-14
#Impulse ss
y, t, x = impulse(sysss, t)
@test_approx_eq_eps y yreal 1e-14
xreal[:,1,1] = yreal[:,1,1]
xreal[:,2,2] = -exp(-t).*t + exp(-t)
xreal[:,3,2] = exp(-t).*t
@test_approx_eq_eps x xreal 1e-14

#Step response of discrete system with specified final time
G = tf([1], [1; zeros(3)], 1)
y, t2, x = step(G, 10)
@test_approx_eq_eps y [zeros(3); ones(8)] 1e-14
@test_approx_eq_eps t2 0:1:10 1e-14

#Impulse response of discrete system to final time that is not mulitple of the sample time
G = tf([1], [1; zeros(3)], 0.3)
y, t2, x = step(G, 2)
@test_approx_eq_eps y [zeros(3); ones(4)] 1e-14
@test_approx_eq_eps t2 0:0.3:1.8 1e-14

#Make sure t was never changed
@test t0 == t
end


text:
require("Color")
require("DataFrames")
require("Iterators")
require("JSON")
require("Mustache")

module Compose

using Color
typealias ColorOrNothing Union(ColorValue, Nothing)

using DataFrames
using Iterators
using JSON

import JSON.json

import Base.+, Base.-, Base.*, Base./, Base.|, Base.convert,
       Base.length, Base.==, Base.<, Base.<=, Base.>=, Base.isempty, Base.insert,
       Base.start, Base.next, Base.done, Base.copy, Base.isless, Base.max,
       Base.<<, Base.>>, Base.show, Base.hex

export |, <<, >>, pad, hstack, vstack, compose, combine, contents, decompose

import Mustache
import Iterators

using Color

# Empty combine. This violates the rules a bit, since nothing is not the
# identity element in any of the monoids, but it's sometimes convenient if it
# behaves as such.
combine() = nothing

include("util.jl")
include("list.jl")
include("measure.jl")
include("backend.jl")
include("property.jl")
include("form.jl")
include("canvas.jl")

# If available, pango and fontconfig are used to compute text extents and match
# fonts. Otherwise a simplistic pure-julia fallback is used.
try
    # Trigger an exception if unavailable.
    dlopen("libfontconfig")
    dlopen("libpangocairo-1.0")
    dlopen("libpango-1.0")

    include("fontconfig.jl")
    include("pango.jl")
catch
    include("fontfallback.jl")
end

# Use cairo for the PNG, PS, PDF if its installed.
try
    require("Cairo")
    include("cairo.jl")
catch
    PNG(::String, ::MeasureOrNumber, ::MeasureOrNumber) =
        error("Cairo must be installed to use the PNG backend.")
    PS(::String, ::MeasureOrNumber, ::MeasureOrNumber) =
        error("Cairo must be installed to use the PS backend.")
    PDF(::String, ::MeasureOrNumber, ::MeasureOrNumber) =
        error("Cairo must be installed to use the PDF backend.")
end

include("svg.jl")
include("d3.jl")
include("dataform.jl")

typealias Composable Union(Form, Property, Canvas)

# Compose operator
<<(a::Form,   b::Property)     = compose(a, b)
<<(a::Form,   b::DataProperty) = compose(a, b)
<<(a::Canvas, b::Form)         = compose(a, b)
<<(a::Canvas, b::Property)     = compose(a, b)
<<(a::Canvas, b::DataProperty) = compose(a, b)
<<(a::Canvas, b::Canvas)       = compose(a, b)

>>(b::Property,     a::Form)   = compose(a, b)
>>(b::DataProperty, a::Form)   = compose(a, b)
>>(b::Form,         a::Canvas) = compose(a, b)
>>(b::Property,     a::Canvas) = compose(a, b)
>>(b::DataProperty, a::Canvas) = compose(a, b)
>>(b::Canvas,       a::Canvas) = compose(a, b)

# Combine operator
|(xs::Property...) = combine(xs...)
|(xs::Form...)     = combine(xs...)

# Compose over hetergenous lists of things.
compose(x) = x
compose(xs::Tuple) = compose(xs...)
compose(xs::Array) = compose(xs...)
compose(x, y, zs...) = compose(compose(compose(x), compose(y)), zs...)

# Make nothings go away.
compose(::Nothing, ::Nothing) = nothing
compose(::Nothing, x::Composable) = x
compose(x::Composable, ::Nothing) = x

# Helpful functions# Create a new canvas containing the given canvas with margins on all sides of
# size u.
function pad(c::Canvas, u::MeasureOrNumber)
    u = size_measure(u)
    compose(canvas(), (canvas(u, u, 1w - 2u, 1h - 2u), c))
end

# TODO: pad_top, pad_bottom, pad_left, pad_right

# Create a new canvas containing the given canvases stacked horizontally.
#
# Args:
#  x0: X-position of the new root canvas
#  y0: Y-position of the new root canvas
#  height: Height of the root canvas.
#  aligned_canvases: One or more canveses accompanied with a vertical alignment
#                    specifier, giving the vertical positioning of the canvas.
#
function hstack(x0::MeasureOrNumber, y0::MeasureOrNumber,height::MeasureOrNumber,
                aligned_canvases::(Canvas, VAlignment)...)

    # To get the expected results, we scale width units, so that everything
    # fits.
    total_width_units = 0
    for (canvas, _)  in aligned_canvases
        if typeof(canvas.box.width) == SimpleMeasure{WidthUnit}
            total_width_units += canvas.box.width.value
        elseif typeof(canvas.box.width) == CompoundMeasure &&
               has(canvas.box.width.values, WidthUnit)
            total_width_units += canvas.box.width.values[WidthUnit]
       end
    end

    width = CompoundMeasure() + 0w

    if length(aligned_canvases) > 0
        width += sum([canvas.box.width for (canvas, _) in aligned_canvases])
    end

    if width.values[WidthUnit] > 0
        width.values[WidthUnit] /= total_width_units
    end

    height = y_measure(height)

    root = canvas(x0, y0, width, height)
    x = 0cx
    for (canvas, aln) in aligned_canvases
        canvas = copy(canvas)
        canvas.box = copy(canvas.box)

        if typeof(canvas.box.width) == SimpleMeasure{WidthUnit}
            canvas.box.width.value /= total_width_units
        elseif typeof(canvas.box.width) == CompoundMeasure &&
               has(canvas.box.width.values, WidthUnit)
            canvas.box.width.values[WidthUnit] /= total_width_units
        end

        # Should we interpret vbottom to mean 0?
        canvas.box.x0 = x
        if aln == vtop
            canvas.box.y0 = 0cy
        elseif aln == vcenter
            canvas.box.y0 = (height / 2) - (canvas.box.height / 2)
        elseif aln == vbottom
            canvas.box.y0 = height - canvas.box.height
        end

        root <<= canvas
        x += canvas.box.width
    end

    root
end

hstack() = canvas()

# Create a new canvas containing the given canvases stacked horizontally.
#
# This is the simple version of hstack. The root canvas will be placed on 0cx,
# 0cy, and its height will be the maximum of the canvases it contains. All
# canvases will be centered vertically.
#
function hstack(canvases::Canvas...)
    height = max([canvas.box.height for canvas in canvases])
    hstack(0, 0, height, [(canvas, vcenter) for canvas in canvases]...)
end

# Create a new canvas containing the given canvases stacked vertically.
#
# Args:
#  x0: X-position of the new root canvas
#  y0: Y-position of the new root canvas
#  width: Height of the root canvas.
#  aligned_canvases: One or more canveses accompanied with a horizontal alignment
#                    specifier, giving the horizontal positioning of the canvas.
#
function vstack(x0::MeasureOrNumber, y0::MeasureOrNumber, width::MeasureOrNumber,
                aligned_canvases::(Canvas, HAlignment)...)

    # Scale height units
    total_height_units = 0
    for (canvas, _)  in aligned_canvases
        if typeof(canvas.box.height) == SimpleMeasure{HeightUnit}
            total_height_units += canvas.box.height.value
        elseif typeof(canvas.box.height) == CompoundMeasure &&
               has(canvas.box.height.values, HeightUnit)
            total_height_units += canvas.box.height.values[HeightUnit]
       end
    end

    width = x_measure(width)

    height = CompoundMeasure() + 0h

    if length(aligned_canvases) > 0
        height += sum([canvas.box.height for (canvas, _) in aligned_canvases])
    end

    if height.values[HeightUnit] > 0
        height.values[HeightUnit] /= total_height_units
    end

    root = canvas(x0, y0, width, height)
    y = 0cy
    for (canvas, aln) in aligned_canvases
        canvas = copy(canvas)
        canvas.box = copy(canvas.box)

        if typeof(canvas.box.height) == SimpleMeasure{HeightUnit}
            canvas.box.height.value /= total_height_units
        elseif typeof(canvas.box.height) == CompoundMeasure &&
               has(canvas.box.height.values, HeightUnit)
            canvas.box.height.values[HeightUnit] /= total_height_units
        end

        canvas.box.y0 = y
        if aln == hleft
            canvas.box.x0 = 0cx
        elseif aln == hcenter
            canvas.box.x0 = (width / 2) - (canvas.box.width / 2)
        elseif aln == hright
            canvas.box.x0 = width - canvas.box.width
        end

        root <<= canvas
        y += canvas.box.height
    end

    root
end

vstack() = canvas()

# Create a new canvas containing the given canvases stacked horizontally.
#
# The simple version of vstack. The root canvas will be placed on 0cx, 0cy, and
# its width will be the maximum of the canvases it contains. All canvases will
# be centered horizontally..
#
function vstack(canvases::Canvas...)
    width = max([canvas.box.width for canvas in canvases])
    vstack(0, 0, width, [(canvas, hcenter) for canvas in canvases]...)
end

# Turn the tree represented by a ComposeType into a nested array (S-expression)
# representation.
#
# This operation is fully invertable with `compose`. Thas is
# decompose(compose(x)) = x should always be true if x is either a Form or a
# Canvas.
#
function decompose(p::Property)
    v = Array(Any, length(p))
    i = 1
    while !is(p, empty_property)
        v[i] = PropertySeq(p.primitive)
        i += 1
        p = p.next
    end
    v
end

decompose(f::EmptyForm) = {}

function decompose(f::FormTree)
    pv = decompose(f.property)
    cv = {decompose(child) for child in f.children}

    f = copy(f)
    f.property = empty_property
    f.children = ListNil{FormTree}()

    v = Array(Any, length(pv) + length(cv) + 1)
    v[1] = f

    i = 2
    for c in cv
        if length(c) == 1
            v[i] = c[1]
        else
            v[i] = c
        end
        i += 1
    end

    for p in pv
        v[i] = p
        i += 1
    end

    v
end

decompose(c::EmptyCanvas) = {}

function decompose(c::CanvasTree)
    pv = decompose(c.property)
    fv = decompose(c.form)
    cv = {decompose(child) for child in c.children}

    c = copy(c)
    c.property = empty_property
    c.form = empty_form
    c.children = ListNil{Canvas}()

    v = Array(Any, length(pv) + length(fv) + length(cv) + 1)
    v[1] = c

    i = 2
    for c in cv
        if length(c) == 1
            v[i] = c[1]
        else
            v[i] = c
        end
        i += 1
    end

    for f in fv
        if length(f) == 1
            v[i] = f[1]
        else
            v[i] = f
        end
        i += 1
    end

    for p in pv
        v[i] = p
        i += 1
    end

    v
end

# Emitting graphics.
#
# Calling draw on a backend without an explicit output files causes the graphic
# to be "emitted". What the means depends on the context. By default, it writes
# the file to a temporary file and opens it in the appropriate viewer.
#
# But this behavior may be changed. In particular when weaving documents, we can
# change the function used to emit graphics to instead embed the graphic in the
# document.
#
# More explicitly, this is just dynamic multiple dispatch. That is multiple
# dispath with bindings that we can change on a whim.
#

type Emitable
    mime::String
    data
end

emitters = Dict{String, Function}()

function emit(emitable::Emitable)
    if has(emitters, emitable.mime)
        emitter = emitters[emitable.mime]
        emitter(emitable.data)
    else
        warn("Unable to emit data of type ", string(typeof(emitable)))
    end
end

# Default emitters

function emitsvg(data::String)
    templ = readall(joinpath(Pkg.dir("Compose"), "src", "show.html"))
    htmlout_path, htmlout = mktemp()
    write(htmlout, Mustache.render(templ, {"svgdata" => data}))
    close(htmlout)
    htm_path = htmlout_path * ".html"
    mv(htmlout_path, htm_path)
    open_browser(htm_path)
end

emitters["image/svg+xml"] = emitsvg

end # module Compose