text:
# file containing methods for requesting ancillary information relating to simulations from GEOGloWS

"""
    available_data()

returns the regions available and their corresponding available dates.
See also `available_regions()` and `available_dates()`
"""
function available_data()
    endpoint = "AvailableData"

    r = _make_request(baseurl,endpoint)

    JSON3.read(r.body,Dict)
end

"""
    available_regions()

returns returns the available regions.
"""
function available_regions()
    endpoint = "AvailableRegions"

    r = _make_request(baseurl,endpoint)

    JSON3.read(r.body,Dict)
end

"""
    available_dates(region::AbstractString)

returns the available forecast dates.
"""
function available_dates(region::AbstractString)
    endpoint = "AvailableDates"

    r = _make_request(baseurl,endpoint,Dict(:region => region))

    JSON3.read(r.body,Dict)
end


text:
# using PyPlot
# close("all")
using jInv.Mesh
using ForwardHelmholtz
using KrylovMethods
using Base.Test

function fictitiousSourceTest2D(u,m,w,rhs)
  refineMesh(M::RegularMesh) = (getRegularMesh(M.domain, M.n*2))

  M  = getRegularMesh([-1.0 1.0 -1.0 1.0],[32,48])
  N  = 4
  err = zeros(N,2)
  # k=3
  uk = 0.0; sk= 0.0; rk = 0.0;k=1;ut=0
   for k=1:N
      M  = refineMesh(M)
  	  xn = getNodalGrid(M)

	  gamma = getABL(M,true,5*ones(Int64,M.dim),0.1);

  	  uk = u(xn[:,1],xn[:,2])
      mk = m(xn[:,1],xn[:,2]);

      rk = rhs(xn[:,1],xn[:,2]) + 1im*(w^2).*gamma[:].*mk[:].*uk;

	  A = GetHelmholtzOperator(M,mk[:],w,gamma[:],true,false);
	  V  = getVolume(M)

  	  ut = A\(-rk);

  	  res = norm(A*ut + rk,Inf);

      err[k,1] = abs(V[1,1]*dot((ut-uk),A*(ut-uk)))
      err[k,2] = norm(ut-uk,Inf)

	  # figure()
	  # imshow(reshape(mk,M.n[1]+1,M.n[2]+1))
	  # imshow(gamma'); colorbar()
	  # subplot(2,2,1)
	  # imshow(reshape(real(ut),M.n[1]+1,M.n[2]+1));;colorbar()
	  # title("ut")
	  # subplot(2,2,2)
	  # imshow(reshape(real(uk),M.n[1]+1,M.n[2]+1));;colorbar()
	  # title("uk")

	  # subplot(2,2,3)
	  # imshow(reshape(abs(uk-ut),M.n[1]+1,M.n[2]+1));colorbar()
	  # title("error")

      @printf "k=%d, n=[%d,%d], l2_err=%1.3e, factor=%1.3f linf_err=%1.3e\n" k M.n[1] M.n[2] err[k,1] err[max(k-1,1),1]/err[k,1] err[k,2]
   end

  # @test countnz(diff(log(err[:,1])).<-1.8) > 4
end

# Constant slowness
f = 1e-0;
w = 2*pi*f;
slowsq(x,y) = ones(size(x));
fictitiousSourceTest2D((x,y)->cos.(pi*x).*cos.(pi*y), slowsq,w,
      (x,y) -> (- pi.^2.*cos.(pi*x).*cos.(pi*y) - pi.^2.*cos.(pi*x).*cos.(pi*y) + slowsq(x,y).*cos.(pi*x).*cos.(pi*y)*(w.^2)  ))

# # VARYING SLOWNESS TEST
slowsq(x,y) = exp.(-2.0*(x.^2 + y.^2));
fictitiousSourceTest2D((x,y)->cos.(pi*x).*cos.(pi*y), slowsq,w,
      (x,y) -> (- pi.^2.*cos.(pi*x).*cos.(pi*y) - pi.^2.*cos.(pi*x).*cos.(pi*y) + slowsq(x,y).*cos.(pi*x).*cos.(pi*y)*(w.^2)  ))

println("\t== passed ! ==")

ext: jl


text:
using ITensors
using ITensorMakie
using Graphs
using GLMakie

g = grid((5,))
tn = itensornetwork(g; linkspaces=10, sitespaces=2)
@visualize fig tn

fig


text:
"""
    ConfidenceInterval

Abstract base type for all confidence interval types.
"""
abstract type ConfidenceInterval end

# these two enable array broadcasting of all methods!
Base.size(::ConfidenceInterval) = ()

Base.getindex(ci::ConfidenceInterval, i) = ci

"""
    limits{T<:Integer}(ci::ConfidenceInterval, x1::T, x2::T)

Return the confidence interval limits for observed x1 and x2.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| ci           | any ConfidenceInterval object |
| x1           | stage-one responses |
| x2           | stage-two responses |

# Return Value

Two element Real vector of limits [lower, upper].

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> ci = ClopperPearsonConfidenceInterval(est, confidence = .9)
julia> limits(ci, 0, 0)
```
"""
limits(ci::ConfidenceInterval, x1::T, x2::T) where {T<:Integer} = error("not implemented!")

confidence(ci::ConfidenceInterval) = try return ci.confidence catch error("not implemented") end

design(ci::ConfidenceInterval) = try return ci.design catch error("not implemented") end

"""
    coverage(ci::ConfidenceInterval, p::T; orientation::String = "overall") where {T<:Real}

Return coverage of given confidence interval and response rate `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| ci           | any ConfidenceInterval object |
| p            | response rate |
| orientation  | string indicating the coverage type - "overall", "lower", or "upper" |

# Return Value

Coverage probability given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> ci = ClopperPearsonConfidenceInterval(est, confidence = .9)
julia> coverage(ci, .5)
```
"""
function coverage(ci::ConfidenceInterval, p::T; orientation::String = "overall") where {T<:Real} # or upper or lower

  @checkprob p
  supp   = support(design(ci))
  res    = 0.0
  for i in 1:size(supp, 1)
    x1, x2 = supp[i, :]
    lim = limits(ci, x1, x2)
    if (orientation == "overall") & (lim[1] <= p) & (lim[2] >= p)
        res += pdf(design(ci), x1, x2, p)
    elseif (orientation == "upper") & (lim[2] >= p)
        res += pdf(design(ci), x1, x2, p)
    elseif (orientation == "lower") & (lim[1] <= p)
        res += pdf(design(ci), x1, x2, p)
    end
  end
  return res

end

"""
    meanwidth(ci::ConfidenceInterval, p::T) where {T<:Real}

Mean width of given confidence interval and response rate `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| ci           | any ConfidenceInterval object |
| p            | response rate |

# Return Value

Mean width given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> ci = ClopperPearsonConfidenceInterval(est, confidence = .9)
julia> meanwidth(ci, .5)
```
"""
function meanwidth(ci::ConfidenceInterval, p::T) where {T<:Real}

  d = design(ci)
  supp   = support(d)
  res    = 0.0
  for i in 1:size(supp, 1)
    x1, x2 = supp[i, :]
    res = res + pdf(d, x1, x2, p)*(limits(ci, x1, x2)[2] - limits(ci, x1, x2)[1])
  end
  return res

end

"""
    meaninterval(ci::ConfidenceInterval, p::T) where {T<:Real}

Mean interval (average limits) of given confidence interval and response rate `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| ci           | any ConfidenceInterval object |
| p            | response rate |

# Return Value

Mean interval given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> ci = ClopperPearsonConfidenceInterval(est, confidence = .9)
julia> meaninterval(ci, .5)
```
"""
function meaninterval(ci::ConfidenceInterval, p::T) where {T<:Real}

  d = design(ci)
  supp   = support(d)
  mean_limits = [0.0; 0.0]
  for i in 1:size(supp, 1)
      x1, x2      = supp[i, :]
      mean_limits = mean_limits + pdf(d, x1, x2, p)*limits(ci, x1, x2)
  end
  return mean_limits

end

"""
    findinconsistencies(ci::ConfidenceInterval, p0::T) where {T<:Real}

Return outcomes where the confidence interval contradicts the designs test decision.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| ci           | any ConfidenceInterval object |
| p            | response rate |

# Return Value

Mean interval given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> ci = ClopperPearsonConfidenceInterval(est, confidence = .9)
julia> findinconsistencies(ci)
```
"""
function findinconsistencies(ci::ConfidenceInterval, p0::T) where {T<:Real}

  d = design(ci)
  supp   = support(d)
  res    = []
  for i in 1:1:size(supp, 1)
      x1, x2 = supp[i, :]
      if (limits(ci, x1, x2)[1] <= p0) & (x1 + x2 > d.c[x1 + 1])
          push!(res, [x1 x2 d.c[x1 + 1] limits(ci, x1, x2)[1] p0])
      end
      if (limits(ci, x1, x2)[1] > p0) & (x1 + x2 <= d.c[x1 + 1])
          push!(res, [x1 x2 d.c[x1 + 1] limits(ci, x1, x2)[1] p0])
      end
  end
  res = vcat(res...)
  return res

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
function technique(tuning_run::Run,
                   channel::RemoteChannel,
                   iterate::Function;
                   name = "Technique Base")
    iteration  = 1
    cost_calls = tuning_run.cost_evaluations

    stop       = RemoteChannel(()->Channel{Bool}(1))
    report     = RemoteChannel(()->Channel{Bool}(1))

    @spawn tuning_run.stopping_criterion(tuning_run.duration, stop)
    @spawn tuning_run.reporting_criterion(tuning_run.report_after / 2, report)

    while !take!(stop)
        result                     = iterate(tuning_run)
        iteration                 += 1

        cost_calls                += result.cost_calls
        result.cost_calls          = cost_calls
        result.start               = tuning_run.starting_point
        result.technique           = name
        result.iterations          = iteration
        result.current_iteration   = iteration

        tuning_run.starting_point  = result.minimum
        tuning_run.starting_cost   = result.cost_minimum

        if take!(report)
            put!(channel, result)
        end
    end

    close(stop)
end


text:
a1 = GeoLocation(-33.8308, 151.223, 0.0)
a2 = GeoLocation(-33.8293, 151.221, 0.0)
A = [a1, a2]
b1 = GeoLocation(-33.8294, 150.22, 0.0)
b2 = GeoLocation(-33.8301, 151.22, 0.0)
B = [b1, b2]
node_a1 = Node(1,a1,nothing)
node_a2 = Node(1,a2,nothing)
node_b1 = Node(1,b1,nothing)
node_b2 = Node(1,b2,nothing)

@testset "Distance tests" begin
    @test isapprox(LightOSM.euclidean(a1, b1), 92.6448020780204)
    @test all(isapprox.(LightOSM.euclidean(A, B), [92.6448020780204, 0.1282389369277561]))
    @test isapprox(LightOSM.haversine(a1, b1), 92.64561837286445)
    @test all(isapprox.(LightOSM.haversine(A, B), [92.64561837286445, 0.1282389369295829]))
    @test LightOSM.euclidean(a1, b1) == LightOSM.euclidean(node_a1, node_b1)
    @test LightOSM.euclidean(A, B) == LightOSM.euclidean([node_a1, node_a2], [node_b1, node_b2])
    @test LightOSM.haversine(a1, b1) == LightOSM.haversine(node_a1, node_b1)
    @test LightOSM.haversine(A, B) == LightOSM.haversine([node_a1, node_a2], [node_b1, node_b2])
    @test LightOSM.euclidean([a1.lat, a1.lon], [b1.lat, b1.lon]) == LightOSM.euclidean(node_a1, node_b1)
    @test LightOSM.haversine([a1.lat, a1.lon], [b1.lat, b1.lon]) == LightOSM.haversine(node_a1, node_b1)

    @test isapprox(distance(a1, b1, :euclidean), 92.6448020780204)
    @test isapprox(distance(a1, b1, :haversine), 92.64561837286445)
    @test all(isapprox.(distance(A, B, :euclidean), [92.6448020780204, 0.1282389369277561]))
    @test all(isapprox.(distance(A, B, :haversine), [92.64561837286445, 0.1282389369295829]))
    @test_throws ArgumentError distance(a1, b1, :unknown_method)
end

@testset "Heading tests" begin
    a = GeoLocation(-33.8308, 151.223, 0.0)
    b = GeoLocation(-33.8294, 151.22, 0.0)
    # Used http://instantglobe.com/CRANES/GeoCoordTool.html to manually calculate headings and distance
    # Note: heading function returns bearings in range of [-180, 180], to convert to [0, 360] scale we need to adjust by (θ + 360) % 360
    @test abs((heading(a, b, :degrees) + 360) % 360 - 299.32559505087795) < 0.01
    @test sum((abs.(heading(A, B, :degrees) .+ 360) .% 360 - [299.32559505087795, 226.07812206538435])) < 0.01

    # Node methods
    @test heading(A, B, :degrees) == heading([node_a1, node_a2], [node_b1, node_b2], :degrees)

    # Radians
    @test isapprox(heading(a, b, :radians), deg2rad(heading(a, b, :degrees)))

    # Error handling
    @test_throws ArgumentError heading(a, b, :unknown_units)
end

@testset "calculate_location tests" begin
    a = GeoLocation(-33.8308, 151.223, 0.0)
    b = GeoLocation(-33.8294, 151.22, 0.0)
    dist_a = 100.0
    dist_b = 100.0
    heading_a = 10.0
    heading_b = 10.0
    end_a = calculate_location(a, heading_a, dist_a)
    @test isapprox(end_a.lat, -32.945001009035984)
    @test isapprox(end_a.lon, 151.40908284685628)
    @test isapprox(end_a.alt, 0.0)
    ends = calculate_location([a, b], [heading_a, heading_b], [dist_a, dist_b])
    @test ends[1] == end_a

    @test calculate_location(a, heading_a, 0.0) == a

    # Node methods
    node_a = Node(1, a, Dict{String, Any}())
    node_b = Node(1, b, Dict{String, Any}())
    @test calculate_location([a, b], [heading_a, heading_b], [dist_a, dist_b]) == calculate_location([node_a, node_b], [heading_a, heading_b], [dist_a, dist_b])
end


text:
# Machinery for complex plane mappings

"zr = rotate_and_scale(za, zb, z)"
function rotate_and_scale(za, zb, z)
    # computes
    # zr = M(z)
    # where
    # M(za) = -1
    # M(zb) = 1
    c_map = ((za+zb)/(za-zb), 2/(zb-za))
    zr = c_map[1] + z*c_map[2]
    return zr
end

function z_edges(grid, i)
    za = grid.edges[1,i] + 1im*grid.edges[2,i]
    zb = grid.edges[3,i] + 1im*grid.edges[4,i]
    return za, zb
end

function map_panels(grid::DiscreteCurve)
    L = legendre_matrix(grid.panelorder)
    ffit(z) = L*z
    coeffs = Array{Complex{Float64}}(grid.panelorder, grid.numpanels)
    for i=1:grid.numpanels
        za, zb = z_edges(grid, i)
        idx = (1:grid.panelorder) + grid.panelorder*(i-1)
        zpanel = grid.points[1,idx] + 1im*grid.points[2,idx]
        # First rescale [za,zb] to [-1,1], which helps Newton search later on
        zscaled = rotate_and_scale(za, zb, zpanel)
        c = ffit(zscaled)
        coeffs[:,i] = c
    end
    return coeffs
end

function invert_map(grid, coeffs, panel_idx, z; fall_back_to_initial::Bool=true)
    # Step 1: Rotate and scale
    za, zb = z_edges(grid, panel_idx)
    zr = rotate_and_scale(za, zb, z)
    # Step 2: Get the right expansion and solve with Newton
    c = coeffs[:, panel_idx]
    c = c[1 : min(length(c), 16)]     # Limit to 16 modes
    t0 = zr
    maxiter = 30
    converged = true
    t, relres, iter = newton_legendre(c, zr, t0, maxiter, 1e-13) # Try to hit 1e-13
    if iter == maxiter && relres > 1e-10 # ...but consider things converged if only 1e-10
        info("Newton did not converge, relres=", @sprintf("%.2e", relres), ". t=$t, t0=$t0")
        converged = false
        if fall_back_to_initial
            t = t0 # Fall back to initial guess
        end
    end
    tlocal = t
    # Step 3: Rescale t back to global parametrization
    ta = grid.t_edges[1, panel_idx]
    tb = grid.t_edges[2, panel_idx]
    tglobal = ta + (t+1)/2*(tb-ta)
    return tlocal, tglobal, converged
end


text:
"""
    Executor

An executor is a realization of a symbolic architecture defined by a `SymbolicNode`.
The actual forward and backward computation specified by the network architecture can
be carried out with an executor.
"""
type Executor
  handle :: MX_ExecutorHandle
  symbol :: SymbolicNode
  arg_arrays  :: Vector{NDArray}
  grad_arrays :: Vector{Union{Void,NDArray}}
  aux_arrays  :: Vector{NDArray}
  outputs     :: Vector{NDArray}
  arg_dict    :: Dict{Base.Symbol, NDArray}
  aux_dict    :: Dict{Base.Symbol, NDArray}
end
function Executor(hdr :: MX_ExecutorHandle, symbol :: SymbolicNode,
                  arg_arrays :: Vector{NDArray}, grad_arrays :: Vector{Union{Void,NDArray}},
                  aux_arrays :: Vector{NDArray})
  # get output arrays
  ref_size = Ref{MX_uint}(0)
  ref_hdrs = Ref{Ptr{MX_handle}}(0)
  @mxcall(:MXExecutorOutputs, (MX_handle, Ref{MX_uint}, Ref{Ptr{MX_handle}}),
          hdr, ref_size, ref_hdrs)
  out_hdrs = unsafe_wrap(Array, ref_hdrs[], ref_size[])
  out_arrays = [NDArray(MX_NDArrayHandle(x)) for x in out_hdrs]

  arg_names = list_arguments(symbol)
  @assert(length(arg_names) == length(unique(arg_names)), "Duplicated names in arguments: $arg_names")
  arg_dict = Dict{Base.Symbol,NDArray}(zip(arg_names, arg_arrays))

  aux_names = list_auxiliary_states(symbol)
  @assert(length(aux_names) == length(unique(aux_names)), "Duplicated names in auxiliary states: $aux_names")
  aux_dict = Dict{Base.Symbol,NDArray}(zip(aux_names, aux_arrays))

  Executor(hdr, symbol, arg_arrays, grad_arrays, aux_arrays, out_arrays, arg_dict, aux_dict)
end

function Base.unsafe_convert(::Type{MX_handle}, obj::Executor)
  Base.unsafe_convert(MX_handle, obj.handle)
end
Base.convert(t::Type{MX_handle}, obj::Executor) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::Executor) = Base.unsafe_convert(t, obj)

function _get_ndarray_inputs(arg_key::AbstractString, args::Vector{NDArray}, arg_names::Vector{Base.Symbol}, allow_missing::Bool)
  @assert(length(args) == length(arg_names), "Length of $arg_key does not match number of arguments")
  return (MX_handle[args...], args)
end
function _get_ndarray_inputs(arg_key::AbstractString, args::Dict{Base.Symbol,NDArray}, arg_names::Vector{Base.Symbol}, allow_missing::Bool)
  args_vec = map(arg_names) do name
    arr = get(args, name, nothing)
    if !allow_missing
      @assert(!isa(arr, Void), "Must specify all arguments in $arg_key ($name is missing)")
    end
    arr
  end
  # help the type inference
  if allow_missing
    args_vec = Union{NDArray,Void}[args_vec...]
  else
    args_vec = NDArray[args_vec...]
  end
  args_hdr = MX_handle[(isa(x,Void) ? MX_handle(0) : x) for x in args_vec]
  return (args_hdr, args_vec)
end

"""
    bind(sym, ctx, args; args_grad=Dict(), aux_states=Dict(), grad_req=GRAD_WRITE)

Create an `Executor` by binding a `SymbolicNode` to concrete `NDArray`.

# Arguments
* `sym::SymbolicNode`: the network architecture describing the computation graph.
* `ctx::Context`: the context on which the computation should run.
* `args`: either a list of `NDArray` or a dictionary of name-array pairs. Concrete
          arrays for all the inputs in the network architecture. The inputs typically include
          network parameters (weights, bias, filters, etc.), data and labels. See [`list_arguments`](@ref)
          and [`infer_shape`](@ref).
* `args_grad`:
* `aux_states`:
* `grad_req`:
"""
@enum GRAD_REQ GRAD_NOP=0 GRAD_WRITE=1 GRAD_ADD=3
function bind(self :: SymbolicNode, ctx :: Context, args :: Union{Vector{NDArray},Dict{Base.Symbol,NDArray}};
              args_grad  :: Union{Vector{NDArray},Dict{Base.Symbol,NDArray}} = Dict{Base.Symbol,NDArray}(),
              aux_states :: Union{Vector{NDArray},Dict{Base.Symbol,NDArray}} = Dict{Base.Symbol,NDArray}(),
              grad_req   :: Union{GRAD_REQ,Vector{GRAD_REQ},Dict{Base.Symbol,GRAD_REQ}} = GRAD_WRITE)

  arg_names = list_arguments(self)

  args_hdr, args           = _get_ndarray_inputs("args", args, arg_names, false)
  args_grad_hdr, args_grad = _get_ndarray_inputs("args_grad", args_grad, arg_names, true)
  aux_args_hdr, aux_states = _get_ndarray_inputs("aux_states", aux_states, list_auxiliary_states(self), false)

  if isa(grad_req, GRAD_REQ)
    reqs = MX_uint[grad_req for i=1:length(args)]
  elseif isa(grad_req, Vector{GRAD_REQ})
    @assert(length(grad_req) == length(args))
    reqs = MX_uint[grad_req...]
  elseif isa(grad_req, Dict{Base.Symbol, GRAD_REQ})
    reqs = MX_uint[get(grad_req, name, GRAD_NOP) for name in arg_names]
  end

  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXExecutorBind,
          (MX_handle, Cint, Cint, MX_uint, Ptr{MX_handle}, Ptr{MX_handle}, Ptr{MX_uint},
           MX_uint, Ptr{MX_handle}, Ref{MX_handle}),
          self, ctx.device_type, ctx.device_id, length(args), args_hdr,
          args_grad_hdr, reqs, length(aux_states), aux_args_hdr, ref_hdr)
  args_grad = convert(Vector{Union{Void,NDArray}}, args_grad)
  executor = Executor(MX_ExecutorHandle(ref_hdr[]), self,
                      args, args_grad, aux_states)
end
function bind(self :: SymbolicNode; kwargs...)
  kwargs = Dict(kwargs)
  @assert(haskey(kwargs, :args), "Must specify args")
  args = pop!(kwargs, :args)
  if haskey(kwargs, :context)
    context = pop!(kwargs, :context)
  else
    context = cpu()
  end
  bind(self, context, args; kwargs...)
end

function simple_bind(self :: SymbolicNode, ctx :: Context;
                     grad_req :: Union{GRAD_REQ, Dict{Symbol, GRAD_REQ}}=GRAD_WRITE,
                     kwargs...)
  arg_shapes, out_shapes, aux_shapes = infer_shape(self; kwargs...)
  @assert(!isa(arg_shapes, Void), "Information not enough to perform complete shape inference")

  arg_arrays = NDArray[zeros(shape, ctx) for shape in arg_shapes]
  arg_names  = list_arguments(self)

  grad_arrays = Dict{Symbol,NDArray}()

  if grad_req != GRAD_NOP
    shapes = zip(arg_names, arg_shapes)

    # if not in provided data, should be parameters
    provided_data_names = [x[1] for x in kwargs]
    shapes = filter(x -> !in(x[1], provided_data_names), shapes)

    # Remove all gradients for nop params
    # if isa(grad_req, Dict{Symbol, GRAD_REQ})
    #  shapes = filter(x -> grad_req[x[1]] != GRAD_NOP,shapes)
    # end

    for (name, shape) in shapes
      grad_arrays[name] = zeros(shape, ctx)
    end
  end

  aux_arrays = [zeros(shape, ctx) for shape in aux_shapes]
  return bind(self, ctx, arg_arrays, args_grad=grad_arrays, grad_req=grad_req, aux_states=aux_arrays)
end

function forward(self :: Executor; is_train::Bool=false, kwargs...)
  for (k,v) in kwargs
    @assert(k ∈ self.arg_dict, "Unknown argument $k")
    @assert(isa(v, NDArray), "Keyword argument $k must be an NDArray")
    copy!(self.arg_dict[k], v)
  end

  @mxcall(:MXExecutorForward, (MX_handle, Cint), self, is_train)
end

function backward(self :: Executor)
  backward(self, NDArray[])
end
function backward(self :: Executor, out_grad :: NDArray)
  backward(self, [out_grad])
end
function backward(self :: Executor, out_grads :: Vector{NDArray})
  out_grads = MX_handle[out_grads...]
  @mxcall(:MXExecutorBackward, (MX_handle, MX_uint, Ptr{MX_handle}), self, length(out_grads), out_grads)
end

function copy_params_from(self::Executor, arg_params::Dict{Base.Symbol,NDArray},
                          aux_params::Union{Void,Dict{Base.Symbol,NDArray}}=nothing;
                          allow_extra_params::Bool=false)
  for (name, array) in arg_params
    if haskey(self.arg_dict, name)
      copy!(self.arg_dict[name], array)
    else
      @assert(allow_extra_params, "Extra params $name not in the arguments")
    end
  end

  if !isa(aux_params, Void)
    for (name, array) in aux_params
      if haskey(self.aux_dict, name)
        copy!(self.aux_dict[name], array)
      else
        @assert(allow_extra_params, "Extra auxiliary state $name not recognized")
      end
    end
  end
end

"""
Get a debug string about internal execution plan.

Can be used to get an estimated about the memory cost.
```julia
  net = ... # Symbol
  dProvider = ... # DataProvider
  exec = mx.simple_bind(net, mx.cpu(), data=size(dProvider.data_batch[1]))
  dbg_str = mx.debug_str(exec)
  println(split(ref, ['\\n'])[end-2])
```
"""
function debug_str(self :: Executor)
  s_ref = Ref{Cstring}()
  @mxcall(:MXExecutorPrint, (MX_handle, Ptr{Cstring}), self.handle, s_ref)
  unsafe_string(s_ref[])
end


text:
## Translate Expr to Mxpr

# Input lines are lexed and parsed by the Julia parser into an AST stored
# in the Julia Expr type. The code here reinterprets the Expr into an Symata AST,
# which is stored in the Julia type Mxpr.

## NB: Checking for help characters, ?,  is done in doc.jl

# At the command line h"topic" prints help
macro h_str(s)
    try
        check_autoload(Symbol(strip(s)))
    catch
    end
    reg = eval(Expr(:macrocall, Symbol("@r_str"), Base.LineNumberNode(1), strip(s), "i"))
    print_matching_topics(reg)
    :Null
end

## Julia already has some of these. We need to consolidate

macro BI_str(s) parse(BigInt,s) end
macro BF_str(s) parse(BigFloat,s) end

"""
    symparsestring(s::String)::Vector{Any}

Parse `s` into an array of Symata expressions, with one expression per array element.

The Symata expressions are not evaluated. For example
`symparsestring("b + b")` returns the Symata expression
`Plus(b, b)`, rather than  `Times(2, b)`.

Note that the phrase *Symata expression* includes numbers, symbols, etc.

See `symparseeval`.
"""
function symparsestring(s)
    mxprs = Vector{Any}(undef, 0)
    s = sjpreprocess_string(s)
    i = 1
    local expr
    while ! (i > ncodeunits(s))
        symata_syntax_deprecation_warnings(false) do
            expr, i = Meta.parse(s,i)
        end
        try
            mx = extomx(expr)
            push!(mxprs, mx)
        catch e
            symprintln("Reading string: got error ", e)
        end
    end
    return mxprs
end

## Complicated:
## 1. preprocess :> to .>, because :> cannot be parsed
##  on output write .> as :>, so that it can be read again. But, .> also works as input
## 2. We want -> for Rule. But, IIRC, the precedence required parenthesizing more than we want.
##   So, we use =>. We would then like to preprocess -> to => on input and the reverse on output.
##   But, we want to be able to write:
## symata > f = :( (x) -> x^2 ), which would be rewritten f = :( (x) => x^2 ).
## We could, when analyzing the parsed :( (x) => x^2 ), convert => back to ->, but I don't want to risk
## messing the precedence of ->, which is asymmetric and special. Mostly because I guess this example will
## be one of the most common uses of :( expr ), to wrap a Julia function.
## So... what to do. Let us try simply breaking with Mma and using => for Rule.
## Perhaps we can also allow -> for Rule as well.
## We would have to do more than simple replacement form string preprocessing. I don't want to do that.
const PREPROCESS_SYMBOL_TRANSLATION = Dict{Any,Any}(
#                                                    "->" => "=>",   #disable
                                                    ":>" => ".>",
                                                    "^:=" => "&=")

const REVERSE_PREPROCESS_SYMBOL_TRANSLATION = Dict{Symbol,Symbol}(
#                                                                  :Rule => :(->), # disabled
                                                                  :Rule => :(=>),
                                                                  :RuleDelayed => Symbol(":>"),
                                                                  :UpSetDelayed => Symbol("^:="))

# Symata language expressions are first processed here (at
# least interactvley now. The string is rewritten and then
# parse to an AST which is the input to symataevaluate. This is
# mostly done to rewrite the string into legal Julia syntax.
# Now, we only look for the help symbol "?"
# sjpreprocess_interactive is called in REPL_symata.jl
const MAGIC_HELP_QUERY_STRING = "HELPQUERY!!"
const MAGIC_HELP_QUERY_SYMBOL = Symbol(MAGIC_HELP_QUERY_STRING)
function sjpreprocess_interactive(line::AbstractString)
    if length(line) > 1 && first(line) == '?'  # User wants documentation
        line =  MAGIC_HELP_QUERY_STRING * ", " *
            SubString(line, 2, lastindex(line)) # We add a comma so that the julia parse will accept it.
    end
    return sjpreprocess_string(line)
end

function sjpreprocess_string(line::AbstractString)
    for (k,v) in PREPROCESS_SYMBOL_TRANSLATION
        line = replace(line, k => v)
    end
    return line
end

## Underscore is not allowed in symbols. Instead,
## they signify `Blank`s.
##
## TODO. This is the only way we construct patterns at the moment.
## But the name of the pattern, the first arg, can refer to not just a Blank,
## but an entire pattern expression. Mma does this with a colon.
## f:(_^_), or f:_^_  --> Pattern[x, Power[Blank[], Blank[]]].
## We are using colon for Span. I don't know what we can do.
## Nov 2016. Moved away from :( ) for Julia expressions. Will remove it soon.
## This frees :( ) and :s for other things
## Mma does  Fullform[a::b] --> MessageName[a, "b"]. We could take
## :: for Pattern
function parseblank(s::AbstractString)
    a = split(s,['_'], keepempty=true)
    length(a) > 4 && error("parseblank: Illegal Pattern expression '$s'")
    if length(a) == 2
        blanktype = :Blank
        (blankhead,blankname) = (a[1],a[2])
    elseif length(a) == 3
        a[2] != "" && error("parseblank: Illegal Pattern expression '$s'")
        blanktype = :BlankSequence
        (blankhead,blankname) = (a[1],a[3])
    else
        a[2] != "" && a[3] != "" && error("parseblank: Illegal Pattern expression '$s'")
        blanktype = :BlankNullSequence
        (blankhead,blankname) = (a[1],a[4])
    end
    blank = isempty(blankname) ? mxpr(blanktype) : mxpr(blanktype,Symbol(blankname))
    return isempty(blankhead) ? blank : mxpr(:Pattern,Symbol(blankhead),blank)
end

# Pattern::argrx: Pattern called with 3 arguments; 2 arguments are expected.
function parsepattern(ex)
    nargs = newargs(length(ex.args))
    mxpr(:Pattern, map!(extomx, nargs, ex.args))
end

extomxarr!(ain, aout) =  foreach(x -> push!(aout, extomx(x)), ain)

# We currently have two kinds of symbols.
# 1) Those tagged by Julia symbols, forming a single namespace
# 2) Those tagged by a.b or Qsym(a,b) in which a is a context or namespace.
# The latter is an experiment, not well integrated with the rest of the language.
# Interactively, the second argument is an Expr.
# Read from a file, it is a QuoteNode. Don't know why.
function parse_qualified_symbol(ex::Expr)
    args = ex.args
    length(args) != 2 && error("extomx: We can only handle context qualifications like this: a.b")
    typeof(args[1]) == Symbol || error("extomx: expecting symbol as first argument to context qualification")
    a2 = args[2]
    if typeof(a2) != Expr
        if typeof(a2) != QuoteNode
            dump(ex)
            error("extomx: error parsing second argument $a2 of ", ex, " , Expected an Expr or QuoteNode, got ", typeof(a2))
        end
        return Qsym(args[1],a2.value)
    end
    typeof(a2.head) != Symbol && error("extomx: error parsing second argument of ", ex, ". Expected a symbol.")
    a2.head != :quote && error("extomx: error parsing second argument of ", ex, ". Expected symbol 'quote'.")
    length(a2.args) == 1 || error("extomx: error parsing second argument of ", ex, ". Expected one arg in quote node.")
    typeof(a2.args[1]) == Symbol || error("extomx: second argument of context qualification must be a symbol, got ", typeof(a2.args[1]))
    return Qsym(args[1],a2.args[1])
end

function parse_quoted(ex::Expr,newa)
    # Quotes are wrapped in Jxpr which is evaluated by Julia eval()
    @symwarn(":( ) for Julia code is deprecated. Use J( ) instead")
    head = :Jxpr           # This allows running Julia code from within Symata.
    push!(newa,ex.args[1]) # We evaluate the expression only whenever the Jxpr is evaled
                           # But, this is the same effect as evaling ex
    return head
end

function rewrite_Jxpr!(a, newa)
    for i in 2:length(a)  ## TODO: use `rest` iterator instead of indices
        x = a[i]
        if isa(x,Expr) && x.head == :kw  # in  J( x =  1) ,  x=1 is interpeted as keyword, value pair. But the user means to do assignment
            x.head = :(=)                # will this break legit uses of keywords ?
        end
        push!(newa, a[i]) # will be interpreted as Julia code, so don't translate it.
    end
    return mxpr(:Jxpr, newa)
end

function parse_pattern(pattern_expr)
    if isa(pattern_expr, Expr)  # assume it is a Function
        if pattern_expr.head == :call && length(pattern_expr.args) > 1 && pattern_expr.args[1] == :J
            #pattern = Core.eval(Symata, pattern_expr.args[2])
            pattern = eval(pattern_expr.args[2])
        else
            pattern = eval(eval(pattern_expr)) # first eval gets Symbol from Expr, Second gives Function.
        end
    else
        pattern = pattern_expr
    end
    isa(pattern, Symbol) || isa(pattern, Function) || typeof(pattern) <: Function ||
        error("extomx: argument to PatternTest must be a Symbol or a Function")
    return pattern
end

# FIXME: refactor after testing
function parse_colon!(a, newa)
    local mxhead
    if length(a) == 2
        lhs = extomx(a[1])
        if typeof(lhs) <: Union{BlankXXX, Mxpr{:Pattern}}  # This is probably not good enough!
            # Maybe : represents too many things.
            return mxpr(:Optional, lhs, extomx(a[2]))
        end
        mxhead = :Span      # Span syntax like:  a(::10), a(1::2), etc. clash with use of colon above
        extomxarr!(a, newa) # We may have to change the syntax
    else
        mxhead = :Span
        extomxarr!(a, newa)
    end
    return mxpr(mxhead, newa)
end

## TODO: handle head :parameters here.
## Check if the first argument (ex.args[2]) is a :parameters expression.
## If so, rewrite ex. We probably want f(a,b;c,d;e,f) to be three compound expressions with two expressions each.
function parse_call!(ex, newa)
    a = ex.args
    operator = a[1]
    if operator == :(:)
        return parse_colon!(@view(a[2:end]), newa)
    end
    mxhead = extomx(operator)
    if mxhead == :J
        return rewrite_Jxpr!(a, newa)
    end
    @inbounds for i in 2:length(a) push!(newa, extomx(a[i])) end
    return mxpr(mxhead, newa)
end

## Translate this separately otherwise the body always inserts a CompoundExpression
"""
    parse_function(ex::Expr)

Create an object of type `Mxpr{:Function}` from `ex`.
(Add someting more here.)
"""
function parse_function(ex::Expr)
    mxpr(:Function, extomx(ex.args[1]), extomx(ex.args[2].args[2]))
end

"""
    parse_macrocall(ex::Expr, newa::MxprArgs)

If `ex` is a Julia `cmd` macrocall (i.e. entered with backticks)
return an object of type `Mxpr{:PatternTest}`. Otherwise Julia-evaluate
the macrocall.
"""
function parse_macrocall(ex::Expr, newa::MxprArgs)
    fullmacroname = string(ex.args[1])
    r = r"^@(.*)_?cmd$"
    m = match(r, fullmacroname) # check for x`foo` (cmd macro)
    m == nothing && return(eval(ex)) # Not a cmd macro
    macroname_with_underscore = m.captures[1]
    macroname = SubString(macroname_with_underscore, 1, lastindex(macroname_with_underscore)-1)
    cmd_string = ex.args[3]
    pattern_expr = Meta.parse(cmd_string)
    pattern = parse_pattern(pattern_expr)
    push!(newa, extomx(Symbol(macroname)), pattern)
    return mxpr(:PatternTest, newa)
end

"""
    iscall(ex::Expr)

true if `ex` has head `:call`

    iscall(ex::Expr, op::Symbol)

true if the first argument, the operator is `op`.

    iscall(ex::Expr, op::Symbol, len::Int)

true if number of args (including op) is `len`

    iscall(ex::Expr, len::Int)

true if number of args (including op) is `len` for any `op`.
"""
iscall(ex::Expr) = ex.head == :call
# is ex a call with operator op ?
iscall(ex::Expr, op::Symbol) = ex.head == :call && ex.args[1] == op
# is ex a call with operator op and len args (including the op) ?
iscall(ex::Expr, op::Symbol, len::Int) = ex.head == :call && ex.args[1] == op && length(ex.args) == len
# is ex a call with len args (including the op) ?
iscall(ex::Expr, len::Int) = iscall(ex) && length(ex.args) == len
# is of the form    a op b , where op is <,>, etc.  i.e. this is  not a chained comparison
is_single_comparison(ex::Expr, op::Symbol) = ex.head == :comparison && length(ex.args) == 3 && ex.args[2] == op
# We check for :call repeatedly. We can optimize this later.
is_binary_minus(ex::Expr) = iscall(ex, :-, 3)
is_unary_minus(ex::Expr) = iscall(ex, :-, 2)
# number of args != 3 will pass through. But probably can't be used
is_division(ex::Expr) = iscall(ex, :/,3)
is_power(ex::Expr) = iscall(ex, :^)
is_sqrt(ex::Expr) = iscall(ex, :Sqrt)
# This may need to be refined
is_range(ex::Expr) = iscall(ex, :(:), 3)
is_range(args...) = false

"""
    range_start(ex::Expr)

return the start value assume `ex` evaluates to a `Range`.
"""
range_start(ex::Expr) = ex.args[2]
range_stop(ex::Expr) = ex.args[3]

## Needed for Julia v0.6.0. Pairs are parsed differently in different
## versions of Julia (as are many expressions).
## The change for Pair is here: https://github.com/JuliaLang/julia/pull/20327
## Goal is to encapsulate interpretation in one place. Here.
# """
#     is_pair(ex)
# """
# function is_pair(ex::Expr)
#     if iscall(ex, :(=>), 3) # In Julia v0.6
#         return true
#     elseif ex.head == :(=>) && length(ex.args) == 2
#         return true
#     end
#     return false
# end

# """
#     rewrite_pair(ex::Expr)
# Input is v6.0 Expr for :( a => b ). Output
# is v5.0 Expr for :( a => b ).
# """
# function rewrite_pair(ex::Expr)
#     Expr(:(=>),ex.args[2],ex.args[3])
# end

## save time, don't make Symata meval convert //(a,b) to rational
## we should also detect sums/products of like numbers, etc. and
## combine them. No* Sometimes we want to Hold expressions involving numbers.
function is_rational(ex::Expr)
    iscall(ex, :(//), 3) && isa(ex.args[2],Number) &&
        isa(ex.args[3],Number)
end

function is_complex(ex::Expr)
    iscall(ex, :Complex, 3) && isa(ex.args[2],Number) &&
        isa(ex.args[3],Number)
end

# In extomx, we first rewrite some Math to canonical forms
# a - b  -->  a + (-b)
rewrite_binary_minus(ex::Expr) = Expr(:call, :+, ex.args[2], Expr(:call,:(-),ex.args[3]))
#  - b  -->  -1 * b
rewrite_unary_minus(ex::Expr) = ex.args[2] == :Inf ? -Inf : Expr(:call, :*, -1 , ex.args[2])
# a / b -->  a * b^(-1)
rewrite_division(ex::Expr) = Expr(:call, :*, ex.args[2], Expr(:call,:^,ex.args[3],-1))

# not need in julia v0.4, but the parser has changed
# ==(a,b) --> comparison  a,==,b
# also,  '<:' is the head in v0.5. We don't yet handle this
#rewrite_to_comparison(ex::Expr) = Expr(:comparison, ex.args[2], ex.args[1], ex.args[3])

function rewrite_to_comparison(ex::Expr)
    op = ex.args[1]
    if haskey(comparison_translation,op)
        return Expr(:call, get(comparison_translation,op,:nothing), ex.args[2], ex.args[3])
    end
    return  Expr(:comparison, ex.args[2], ex.args[1], ex.args[3])
end

rewrite_comparison(ex::Expr) = ex
## The following does what we want. e.g. we get Less(x,y,z)
## But, we have not yet implemented the logic to reduce the expression.
## So, we disable this function and fall back to the older Comparison, for which
## we do have logic
# function rewrite_comparison(ex::Expr)
#     args = ex.args
#     ops = @view args[4:2:end]
#     op1 = args[2]
#     if all( x -> x == op1,  ops)
#         Expr(:call,get(comparison_translation,op1,op1),args[1:2:end]...)
#     else
#         ex
#     end
# end

# There is no binary minus, no division, and no sqrt in Mxpr's.
# Concrete example: a - b --> a + -b.
# We definitely need to dispatch on a hash query, or types somehow
# Other rewrites needed, but not done.

"""
    rewrite_expr(ex::Expr)

Preprocess `ex` to prepare for conversion to type `Mxpr`.
"""
function rewrite_expr(inex::Expr)
    if inex.head == :block  # changed in Julia v0.7, need to rmlines
        ex = MacroTools.rmlines(inex)
    else
        ex = inex
    end
    for i in 1:length(ex.args)
        x = ex.args[i]
        if isa(x,Expr) && x.head == :macrocall
            ex.args[i] = parse_macrocall(x, newargs())
        end
        if haskey(INSYMTRANS,x)
            ex.args[i] = INSYMTRANS[x]
        end
    end
    if iscall(ex) && length(ex.args) > 0 && is_comparison_symbol(ex.args[1])
        ex = rewrite_to_comparison(ex)
    elseif is_single_comparison(ex, :(.>))    # julia 0.4 parser does this. In 0.5, this is already a call
        return Expr(:call, :(.>), ex.args[1], ex.args[3])
    elseif ex.head == :comparison
        ex = rewrite_comparison(ex)
    elseif is_unary_minus(ex)    #  - b --> -1 * b
        ex = rewrite_unary_minus(ex)
    elseif is_binary_minus(ex)  #  a - b --> a + -b.
        ex = rewrite_binary_minus(ex)
    elseif is_division(ex) # a / b --> a + b^(-1)
        ex = rewrite_division(ex)
    elseif iscall(ex, :Exp, 2)  # Exp(x) --> E^x
        ex = Expr(:call, :^, :E, ex.args[2])
    elseif iscall(ex,:Sqrt,2) # This should happen at Mxpr level, and be optimized
        ex = Expr(:call, :^, ex.args[2], Expr(:call,:(//), 1,2))
    elseif iscall(ex, :(=>), 3)  # Convert v0.6 Pair to v0.5 Pair.
        ex = Expr(:(=>), ex.args[2], ex.args[3])
    elseif is_rational(ex)
        return eval(ex)  # TODO: don't do eval, use //
    elseif is_complex(ex)
        (real,imag) = (ex.args[2],ex.args[3])
        return isa(real,Real) && isa(imag,Real) ? complex(real,imag) : ex
    end
    return ex
end

## extomx

## Main translation routine
# We use Julia for lexing/parsing. But we change the semantics:
# sometimes a little, sometimes a lot.
"""
    extomx(expr::Expr)

Translate a Julia expression `expr`, to a Symata expression of
type `Mxpr{T}`, where `T` is typically determined by the operator
in `expr`.
"""
function extomx(ex::Expr)
    newa = newargs()
    local mxhead::Any
    ex = rewrite_expr(ex)
    isa(ex, Expr) || return ex  # may be a rational
    ohead = ex.head
    a = ex.args
    # We usually set the head and args in the conditional and construct Mxpr at the end
    if iscall(ex) return parse_call!(ex, newa)
    elseif ohead == :block && typeof(a[1]) == LineNumberNode  # g(x_Integer) = "int". julia finds line number node in rhs.
        return extomx(a[2])
    elseif ohead == :line return nothing # Ignore line number. part of misinterpretation of g(x_Integer) = "int".
    elseif ohead == :(->) return parse_function(ex)
    elseif haskey(JTOMSYM,ohead)
        mxhead = JTOMSYM[ohead]
        check_autoload(mxhead)
        extomxarr!(a,newa)
    elseif ohead == :kw  # Interpret keword as Set, but Expr is different than when ohead == :(=)
        mxhead = :Set
        extomxarr!(a,newa)
    elseif ohead == :(::)
        return parsepattern(ex)
    elseif ohead == :(:) # Eg the colon here: g(x_Integer:?(EvenQ)) := x
        error("This should never happen. Put assert or something here")
        if length(a) == 2
            lhs = extomx(a[1])
            if typeof(lhs) <: Union{BlankXXX, Mxpr{:Pattern}}  # This is probably not good enough!
                # Maybe : represents too many things.
                return mxpr(:Optional,lhs,extomx(a[2]))
            end
            mxhead = :Span      # Span syntax like:  a(::10), a(1::2), etc. clash with use of colon above
            extomxarr!(a,newa) # We may have to change the syntax
        else   # FIXME: refactor after testing
            mxhead = :Span
            extomxarr!(a,newa)
        end
    elseif ohead == :(.)  # don't need parens, but this is readable
        return parse_qualified_symbol(ex)
    elseif ohead == :quote
        mxhead = parse_quoted(ex, newa)
    elseif ohead == :macrocall
#        return eval(ex)
        return parse_macrocall(ex, newa)
    elseif ohead == :string
        mxhead = :StringInterpolation
        extomxarr!(a,newa)
    else
        dump(ex)
        error("extomx: No translation for Expr head '$(ohead)' in $ex")
    end
    mx = mxpr(mxhead,newa)  # Create the Mxpr
end

"""
    extomx(x::Integer)

Return `x` or convert `x` to `BigInt` if automatic
conversion is enabled with `BigIntInput(True)`.
"""
extomx(x::Integer) = getkerneloptions(:bigint_input) ? BigInt(x) : x
"""
    extomx(x::AbstractFloat)

Return `x` or convert `x` to `BigFloat` if automatic
conversion is enabled with `BigFloatInput(True)`.
"""
extomx(x::AbstractFloat) = getkerneloptions(:bigfloat_input) ? BigFloat(rationalize(x)) : x
extomx(x) = x
# This system needs to be rationalized.
# Comment out lines are no longer needed
# We can move the rest out of here too
# jtomsym is what picks these out.
function extomx(s::Symbol)
    s == :I && return complex(0,1)
    s == :∑ && return :Sum
    s == :True && return true
    s == :False && return false
    ss = string(s)
    return occursin("_", ss)  ? parseblank(ss) : getsym(jtomsym(s))
end

"""
    @extomx(ex)

inserts translation of Julia expression `ex` to Symata into code. This is used
in Symata code. For instance in `Collect` in `algebra.jl`.

`@extomx` differs from `@sym` in that it does not evaluate the resulting Symata expression.
"""
macro extomx(ex)
    :($(extomx(ex)))
end

## FIXME: Find a good name for this macro and document.
## eg: f = @jul x -> x^2
macro jul(ex) extomx(Expr(:call, :J, ex)) end


text:
# Fast sampling from truncated normal distributions

using Distributions

include("randtn_zigtables.jl")

const log2π = 1.83787706641
const tnzig_kmin = 5

# returns sample x ~ N(0,1) truncated at [a,b], using rejection sampling with
# a truncated exponential proposal
function randtn_truncexp(a::Real, b::Real)
    twoasq = 2abs2(a)
    c = exp(-a*(b - a)) - 1
    while true
        z = log(1 + c * rand())
        if -twoasq * log(rand()) > abs2(z)
            return a - z/a
        end
    end
end

# returns sample x ~ N(0,1) truncated at [a,b], using rejection sampling with
# a standard normal proposal
function randtn_norm(a::Real, b::Real)
    while true
        z = randn()
        if z >= a && z <= b
            return z
        end
    end
end

# returns sample x ~ N(0,1) truncated at [a,b], assuming b > a and |a| < |b|.
function randtnstdord(a::Real, b::Real, tp::Real)
    # a in left tail, or large mass, implies b > -a, use normal proposal
    if a < tnzig_xmin || tp > 0.3
        return randtn_norm(a, b)
    end

    # a in right tail implies b > a, use truncated exponential proposal
    if a > tnzig_xmax
        return randtn_truncexp(a, b)
    end

    # tnzig_xmin <= a <= tnzig_xmax, use Chopin's algorithm
    N = length(tnzig_x) - 1
    # a in interval ka or ka+1
    ka = tnzig_ncell[tnzig_i0 + ifloor(a * tnzig_invh)]
    kb = b < tnzig_xmax ?               # b in interval kb-1 or kb
        min(tnzig_ncell[tnzig_i0 + floor(b * tnzig_invh)] + 1, N) :
        N + 1                                     # b in right tail

    # for small kb - ka use rejection sampling with truncated exponential prop.
    if abs(kb - ka) < tnzig_kmin
        return randtn_truncexp(a, b)
    end

    while true
        # sample integer k from range [ka, kb]
        k = ka + ifloor((kb - ka + 1) * rand())

        if k == N + 1
            # right tail, rejection sampling with truncated exponential proposal
            xN1 = tnzig_x[end]
            z = -log(rand()) / xN1
            if abs2(z) < -2log(rand()) && z < b - xN1
                return xN1 + z
            end

        elseif k > ka + 1 && (k < kb - 1 || b >= tnzig_xmax)
            # all other intervals, away from boundaries
            u = rand()
            simy = tnzig_yu[k] * u;
            if simy < tnzig_yl[k]
                # happens most of the time
                return tnzig_x[k] + u * tnzig_delta[k]
            end
            # check if we are below the pdf curve
            sim = tnzig_x[k] + tnzig_d[k] * rand()
            if abs2(sim) + 2log(simy) + log2π < 0
                return sim
            end

        else
            # boundary intervals ka, ka+1, kb-1, kb
            sim = tnzig_x[k] + tnzig_d[k] * rand()
            if sim >= a && sim <= b
                simy = tnzig_yu[k] * rand()
                if simy < tnzig_yl[k] || abs2(sim) + 2log(simy) + log2π < 0
                    return sim
                end
            end
        end
    end
end
# same as randtnstdord, but does not require |a| < |b|
randtnstd(a::Real, b::Real, tp::Real) = abs(a) > abs(b) ?
                                        -randtnstdord(-b, -a, tp) :
                                        randtnstdord(a, b, tp)

function randtn(d::Truncated{Normal})
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    z = abs(a) > abs(b) ? -randtnstdord(-b, -a, d.tp) : randtnstdord(a, b, d.tp)
    return μ + σ * z
end

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn


text:
# Autogenerated wrapper script for Ccache_jll for x86_64-linux-gnu-cxx03
export ccache

using Zlib_jll
using Zstd_jll
JLLWrappers.@generate_wrapper_header("Ccache")
JLLWrappers.@declare_executable_product(ccache)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Zstd_jll)
    JLLWrappers.@init_executable_product(
        ccache,
        "bin/ccache",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
immutable SVD
  U::Matrix{Float64}
  V::Matrix{Float64}
  s::Vector{Float64}
  transp::Bool
end

function svd_singular_val_creation!(::Type{Val{1}}, e::Vector{Float64}, s::Vector{Float64}, ::Matrix{Float64}, V::Matrix{Float64},
                                    p::Int, ::Int, iter::Int, k::Int, m::Int, n::Int)
  # Deflate negligible s(p)
  f = e[p-1]
  e[p-1] = 0.0
  @inbounds for j = p-1:-1:k
    t = hypot(s[j],f)
    cs = s[j]/t
    sn = f/t
    s[j] = t
    if j != k
      f = -sn*e[j-1]
      e[j-1] = cs*e[j-1]
    end
    @inbounds for i = 1:n
      t = cs*V[i,j] + sn*V[i,p]
      V[i,p] = -sn*V[i,j] + cs*V[i,p]
      V[i,j] = t
    end
  end

  return iter, p, k
end

function svd_singular_val_creation!(::Type{Val{2}}, e::Vector{Float64}, s::Vector{Float64}, U::Matrix{Float64}, ::Matrix{Float64},
                                    p::Int, ::Int, iter::Int, k::Int, m::Int, n::Int)
  # Split at negligible s(k)
  f = e[k-1]
  e[k-1] = 0.0
  @inbounds for j = k:p
    t = hypot(s[j],f)
    cs = s[j]/t
    sn = f/t
    s[j] = t
    f = -sn*e[j]
    e[j] = cs*e[j]
    @inbounds for i = 1:m
      t = cs*U[i,j] + sn*U[i,k-1]
      U[i,k-1] = -sn*U[i,j] + cs*U[i,k-1]
      U[i,j] = t
    end
  end
  return iter, p, k
end

function svd_singular_val_creation!(::Type{Val{3}}, e::Vector{Float64}, s::Vector{Float64}, U::Matrix{Float64}, V::Matrix{Float64},
                                    p::Int, ::Int, iter::Int, k::Int, m::Int, n::Int)
  # Perform one qr step

  # Calculate the shift
  scale_ = max(
            max(
              max(
                max(abs(s[p]), abs(s[p-1])),
                  abs(e[p-1])),
                abs(s[k])),
              abs(e[k]))

  sp = s[p]/scale_
  spm1 = s[p-1]/scale_
  epm1 = e[p-1]/scale_
  sk = s[k]/scale_
  ek = e[k]/scale_
  b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0
  c = (sp*epm1)*(sp*epm1)
  shift = 0.0
  if b != 0.0 || c != 0.0
    shift = sqrt(b*b + c)
    if b < 0.0
      shift = -shift
    end
    shift = c/(b + shift)
  end
  f = (sk + sp)*(sk - sp) + shift
  g = sk*ek

  # Chase zeros
  @inbounds for j = k:p-1
    t = hypot(f,g)
    cs = f/t
    sn = g/t
    if j != k
      e[j-1] = t
    end
    f = cs*s[j] + sn*e[j]
    e[j] = cs*e[j] - sn*s[j]
    g = sn*s[j+1]
    s[j+1] = cs*s[j+1]
    @inbounds for i = 1:n
      t = cs*V[i,j] + sn*V[i,j+1]
      V[i,j+1] = -sn*V[i,j] + cs*V[i,j+1]
      V[i,j] = t
    end
    t = hypot(f,g)
    cs = f/t
    sn = g/t
    s[j] = t
    f = cs*e[j] + sn*s[j+1]
    s[j+1] = -sn*e[j] + cs*s[j+1]
    g = sn*e[j+1]
    e[j+1] = cs*e[j+1]
    if j < m
      @inbounds for i = 1:m
        t = cs*U[i,j] + sn*U[i,j+1];
        U[i,j+1] = -sn*U[i,j] + cs*U[i,j+1]
        U[i,j] = t
      end
    end
  end
  e[p-1] = f
  iter = iter + 1

  return iter, p, k
end

function svd_singular_val_creation!(::Type{Val{4}}, e::Vector{Float64}, s::Vector{Float64}, U::Matrix{Float64},
                                    V::Matrix{Float64}, p::Int, pp::Int, iter::Int, k::Int, m::Int, n::Int)
  # Make the singular values positive
  if s[k] <= 0.0
    s[k] = s[k] < 0.0 ? -s[k] : 0.0
    for i = 1:pp+1
      V[i,k] = -V[i,k]
    end
  end

  # Order the singular values
  @inbounds while k <= pp
    if s[k] >= s[k+1]
      break
    end
    s[k], s[k+1] = s[k+1], s[k]
    if k < n
      @inbounds for i = 1:n
        V[i, k], V[i, k+1] = V[i, k+1], V[i, k]
      end
    end
    if k < m
      @inbounds for i = 1:m
        U[i, k], U[i, k+1] = U[i, k+1], U[i, k]
      end
    end
    k += 1
  end

  iter = 0
  p -= 1

  return iter, p, k
end

function SVD(M::Matrix{Float64})
  m, n = size(M)

  if m >= n
    A = copy(M)
    transp = false
  else
    A = M'
    transp = true
  end

  s = Vector{Float64}(n)
  U = zeros(m, n)
  V = Matrix{Float64}(n, n)
  e = Vector{Float64}(n)
  work = Vector{Float64}(m)
  i = j = k = 0

  # reduce A to its bidiagonal form, storing the diagonal elements
  # in s, and super-diagonal elements in e
  nct = min(m-1, n)
  nrt = max(0, n-2)

  @inbounds @simd for k = 1:max(nct, nrt)
    if k <= nct
      # Compute the transformation for the k-th column and
      # place the k-th diagonal in s[k].
      # Compute 2-norm of k-th column without under/overflow.
      s[k] = 0.0
      @inbounds for i = k:m
        s[k] = hypot(s[k],A[i,k]);
      end
      if s[k] != 0.0
        if A[k,k] < 0.0
          s[k] = -s[k];
        end
        for i = k:m
          A[i,k] /= s[k]
        end
        A[k,k] += 1.0;
      end
      s[k] = -s[k]
    end
    @inbounds for j = k+1:n
      if k <= nct && s[k] != 0.0
        # apply transformation
        t = 0.0
        @inbounds for i = k:m
          t += A[i,k]*A[i,j]
        end
        t = -t/A[k,k]
        @inbounds for i = k:m
          A[i,j] += t*A[i,k]
        end
      end
      # Place the k-th row of A into e for the
      # subsequent calculation of the row transformation.
      e[j] = A[k, j]
    end
    if k < nct
      # Place the transformation in U for subsequent back
      # multiplication.
      @inbounds for i = k:m
        U[i,k] = A[i,k]
      end
    end
    if k <= nrt
      # Compute the k-th row transformation and place the
      # k-th super-diagonal in e[k].
      # Compute 2-norm without under/overflow.
      e[k] = 0.0
      @inbounds for i = k+1:n
        e[k] = hypot(e[k],e[i])
      end
      if e[k] != 0.0
        if e[k+1] < 0.0
          e[k] = -e[k]
        end
        @inbounds for i = k+1:n
          e[i] /= e[k]
        end
        e[k+1] += 1.0
      end
      e[k] = -e[k]
      if k+1 <= m && e[k] != 0.0

        # Apply the transformation.
        for i = k+1:m
          work[i] = 0.0
        end
        @inbounds for j = k+1:n
          @inbounds for i = k+1:m
            work[i] += e[j]*A[i,j]
          end
        end
        @inbounds for j = k+1:n
          t = -e[j]/e[k+1]
          @inbounds for i = k+1:m
            A[i,j] += t*work[i]
          end
        end
      end

      # Place the transformation in V for subsequent
      # back multiplication.
      @inbounds for i = k+1:n
        V[i,k] = e[i]
      end
    end
  end

  # Set up the final bidiagonal matrix or order n.

  if nct < n
    s[nct+1] = A[nct+1, nct+1]
  end
  if nrt+1 < n
    e[nrt+1] = A[nrt+1, n]
  end
  e[n] = 0.0

  # Generate U
  @inbounds @simd for j = nct+1:n
    for i = 1:m
      U[i,j] = 0.0
    end
    U[j,j] = 1.0
  end

  @inbounds @simd for k = nct:-1:1
    if s[k] != 0.0
      @inbounds for j = k+1:n
        t = 0
        @inbounds for i = k:m
          t += U[i,k]*U[i,j]
        end
        t = -t/U[k,k]
        @inbounds for i = k:m
          U[i,j] += t*U[i,k]
        end
      end
      @inbounds for i = k:m
        U[i,k] = -U[i,k]
      end
      U[k,k] = 1.0 + U[k,k]
      @inbounds for i = 1:k-1
        U[i,k] = 0.0
      end
    else
      @inbounds for i = 1:m
        U[i,k] = 0.0
      end
      U[k,k] = 1.0
    end
  end

  # Generate V
  @inbounds for k = n:-1:1
    if k <= nrt && e[k] != 0.0
      @inbounds for j = k+1:n
        t = 0.0
        @inbounds for i = k+1:n
          t += V[i,k]*V[i,j]
        end
        t = -t/V[k+1,k]
        @inbounds for i = k+1:n
          V[i,j] += t*V[i,k]
        end
      end
    end
    @inbounds for i = 1:n
      V[i,k] = 0.0
    end
    V[k,k] = 1.0
  end

  # Main iteration loop for singular values
  p = n
  pp = p-1
  iter = 0

  @inbounds while p > 0
    # Here is where a test for too many iterations would go.

    # This section of the program inspects for
    # negligible elements in the s and e arrays.  On
    # completion the variables kase and k are set as follows.

    # kase = 1     if s(p) and e[k-1] are negligible and k<p
    # kase = 2     if s(k) is negligible and k<p
    # kase = 3     if e[k-1] is negligible, k<p, and
    #              s(k), ..., s(p) are not negligible (qr step).
    # kase = 4     if e(p-1) is negligible (convergence).

    for k = p -1:-1:0
      if k == 0
        break
      end

      if abs(e[k]) <= EPS_VAL*(abs(s[k]) + abs(s[k+1]))
        e[k] = 0.0
        break
      end
    end

    if k == p-1
        kase = 4
    else
      ks = 0
      @inbounds for ks = p:-1:k
        if ks == k
          break
        end

        t = (ks != p+1 ? abs(e[ks]) : 0.0) + (ks != k+1 ? abs(e[ks-1]) : 0.0)
        if abs(s[ks]) <= EPS_VAL * t
          s[ks] = 0.0;
          break
        end
      end
      if ks == k
        kase = 3
      elseif ks == p
        kase = 1
      else
        kase = 2
        k = ks
      end
    end
    k += 1

    # perform the task indicated by kase
    iter, p, k = svd_singular_val_creation!(Val{kase}, e, s, U, V, p, pp, iter, k, m, n)

  end
  return SVD(U, V, s, transp)
end


text:
using ODE

# Definition of parameters
const a = 1.0
const b = 0.1
const c = 1.5
const d = 0.75

function diffeq(t,u)
    du = similar(u)
    # u1: rabbits
    # u2: foxes
    du[1] = a*u[1] - b*u[1]*u[2]
    du[2] = -c*u[2] + d*b*u[1]*u[2]

    return du

end

# Simulation time
const t = collect(range(0.0,stop=18.0,length=1000))

# Initial condition
u0 = zeros(2)
u0[1] = 10
u0[2] = 5

# Solve ODE
(sol_t,sol_u) = ode45(diffeq,u0,tpoints=:specified)

#= Visualization of Results

using PyPlot
plot(sol_t,sol_u)
show()

=#


text:
function toindices(times, fs)
    round.(Int, times.*fs) .+ 1
end

function totimes(indices, fs)
    (indices .- 1)./fs
end

"""
    findduplicates(itr)

"""
function findduplicates(itr)
    dups = Dict{eltype(itr), Vector{Int}}()

    for (i, v) in enumerate(itr)
        push!(get!(dups, v, Int[]), i)
    end
    filter!(((k,v),) -> length(v) !== 1, dups)

    return dups
end

function matchevents(pred::AbstractVector, act::AbstractVector; Tthresh=0.5*median(diff(act)))
    mdiff_pred = median(diff(pred))
    mdiff_act = median(diff(act))
    err(a,b) = abs(a - b)/b
    if err(mdiff_pred, mdiff_act) > 10 || err(mdiff_act, mdiff_pred) > 10
        throw(error("""detected a large difference in event frequency. Are predicted events
            in the same units as actual events (eg frames vs sec)?"""))
    end
    missed = max(0, length(act) - length(pred))

    _pred = copy(pred)
    tree = KDTree(reshape(act, 1, length(act)), Chebyshev())
    idxs, dists = nn(tree, reshape(pred, 1, length(pred)))

    dups = findduplicates(idxs)
    delidxs = Int[]
    if !isempty(dups)
        foreach(dups) do (k, v)
            append!(delidxs, v[setdiff(eachindex(v), argmin(dists[v]))])
        end
        sort!(delidxs)
        deleteat!(idxs, delidxs)
        deleteat!(dists, delidxs)
        deleteat!(_pred, delidxs)
    end

    dists .= _pred - act[idxs]

    return setdiff(eachindex(pred), delidxs), idxs, dists, missed
end


text:
abstract type AbstractEquationOfState <: AbstractPhysicsComponent end

Base.@kwdef struct BarotropicFluid{FT} <: AbstractEquationOfState
    ρₒ::FT
    cₛ::FT
end

Base.@kwdef struct DryIdealGas{FT} <: AbstractEquationOfState
    R::FT
    pₒ::FT
    γ::FT
end

# Energy Prognostic Cases, perhaps abstract type

Base.@kwdef struct TotalEnergy{FT} <: AbstractEquationOfState
    γ::FT
end

Base.@kwdef struct DryEuler{FT} <: AbstractEquationOfState
    γ::FT
end

Base.@kwdef struct LinearizedTotalEnergy{FT} <: AbstractEquationOfState
    γ::FT
end

Base.@kwdef struct LinearizedDryEuler{FT} <: AbstractEquationOfState
    γ::FT
end

# Linearized Equation of States
linearize(eos::TotalEnergy) = LinearizedTotalEnergy(eos.γ)
linearize(eos::DryEuler) = LinearizedDryEuler(eos.γ)

"""
  Thermodynamic relationships
"""

@inline function calc_pressure(eos::BarotropicFluid, state)
    cₛ = eos.cₛ
    ρₒ = eos.ρₒ
    ρ = state.ρ

    return (cₛ * ρ)^2 / (2 * ρₒ)
end

@inline function calc_sound_speed(eos::BarotropicFluid, state)
    cₛ = eos.cₛ
    ρₒ = eos.ρₒ
    ρ = state.ρ

    return cₛ * sqrt(ρ / ρₒ)
end

@inline function calc_pressure(eos::DryIdealGas, state)
    R = eos.R
    pₒ = eos.pₒ
    γ = eos.γ
    ρθ = state.ρθ

    return pₒ * (R / pₒ * ρθ)^γ
end

@inline function calc_sound_speed(eos::DryIdealGas, state)
    R = eos.R
    pₒ = eos.pₒ
    γ = eos.γ
    ρ = state.ρ

    return sqrt(γ * calc_pressure(eos, state) / ρ)
end

"""
Modified Maciek's world
"""
@inline function calc_pressure(eos::TotalEnergy, state, aux)
    γ = eos.γ
    Φ = aux.Φ
    ρ = state.ρ
    ρe = state.ρe
    ρu = state.ρu
    return (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ - ρ * Φ)
end

@inline function calc_sound_speed(eos::TotalEnergy, state, aux)
    γ = eos.γ
    ρ = state.ρ
    return sqrt(γ * calc_pressure(eos, state, aux) / ρ)
end

@inline function calc_pressure(eos::LinearizedTotalEnergy, state, aux)
    γ = eos.γ
    Φ = aux.Φ
    ρ = state.ρ
    ρe = state.ρe
    return (γ - 1) * (ρe - ρ * Φ)
end

@inline function calc_sound_speed(eos::LinearizedTotalEnergy, state, aux)
    γ = eos.γ
    ρ = state.ρ
    return sqrt(γ * calc_pressure(eos, state, aux) / ρ)
end

@inline function calc_pressure(eos::DryEuler, state, aux)
    γ = eos.γ
    ρ = state.ρ
    ρe = state.ρe
    ρu = state.ρu
    return (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ)
end

@inline function calc_sound_speed(eos::DryEuler, state, aux)
    γ = eos.γ
    ρ = state.ρ
    return sqrt(γ * calc_pressure(eos, state, aux) / ρ)
end

@inline function calc_pressure(eos::LinearizedDryEuler, state, aux)
    γ = eos.γ
    ρe = state.ρe
    return (γ - 1) * ρe
end

@inline function calc_sound_speed(eos::LinearizedDryEuler, state, aux)
    γ = eos.γ
    ρ = state.ρ
    return sqrt(γ * calc_pressure(eos, state, aux) / ρ)
end

"""
  Maciek's world
"""
function pressure(ρ, ρu, ρe, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ - ρ * Φ)
    else
        (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ)
    end
end

function totalenergy(ρ, ρu, p, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        return p / (γ - 1) + dot(ρu, ρu) / 2ρ + ρ * Φ
    else
        return p / (γ - 1) + dot(ρu, ρu) / 2ρ
    end
end

function soundspeed(ρ, p)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    sqrt(γ * p / ρ)
end

@inline function linearized_pressure(ρ, ρe, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        (γ - 1) * (ρe - ρ * Φ)
    else
        (γ - 1) * ρe
    end
end

"""
  Base extensions
"""
function info(::AbstractEquationOfState)
    error("Not implemented!")
end

function info(::BarotropicFluid)
    println("The equation of state (eos) is:")
    printstyled("eos = (cₛ ρ)^2 / (2 ρₒ) \n", color = 82)
    println("The sound speed is:")
    printstyled("soundspeed = cₛ sqrt(ρ / ρₒ)  \n", color = 82)
    println("ρ : density ")
    println("cₛ: reference soundspeed ")
    println("ρₒ: reference density ")
end

function info(::DryIdealGas)
    error("Not implemented!")
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
export get_iterates_algoparams

function get_alg_color(algoname)
    res = nothing

    if algoname == "ISTA"
        res = "rgb,1:red,0.0;green,0.3843;blue,0.9922"
    elseif algoname == "FISTA"
        res = "rgb,1:red,0.7451;green,0.0;blue,0.0"
    elseif algoname == "T1"
        res = "rgb,1:red,0.0;green,0.502;blue,0.3804"
    elseif algoname == "T2"
        res = "black"
    elseif algoname == "MFISTA"
        res = "rgb,255:red,231;green,84;blue,121"
    elseif algoname == "ISTA - bt"
        res = distinguishable_colors(9, transform=deuteranopic)[1]
    elseif algoname == "FISTA - bt"
        res = distinguishable_colors(9, transform=deuteranopic)[5]
    elseif algoname == "T1 - bt"
        res = distinguishable_colors(9, transform=deuteranopic)[7]
    elseif algoname == "T2 - bt"
        res = distinguishable_colors(9, transform=deuteranopic)[8]
    end

    return res
end

function get_iterates_algoparams(algoname)
    res = Dict{Any, Any}()

    FigsIteratesAlgoCurveswidth = "very thick"

    res[FigsIteratesAlgoCurveswidth] = nothing

    if algoname == "ISTA"
        res = Dict(
            "smooth" => nothing,
            FigsIteratesAlgoCurveswidth => nothing,
            "mark size" => "1pt",
            "color" => get_alg_color(algoname),
            "mark" => "*",
            "mark options" => ""
        )
    elseif algoname == "FISTA"
        res = Dict(
            "smooth" => nothing,
            FigsIteratesAlgoCurveswidth => nothing,
            "mark size" => "1pt",
            "color" => get_alg_color(algoname),
            "mark" => "triangle*",
            "mark options" => ""
        )
    elseif algoname == "T1"
        res = Dict(
            "smooth" => nothing,
            FigsIteratesAlgoCurveswidth => nothing,
            # "fill opacity" => "0",
            "mark size" => "1.5pt",
            "color" => get_alg_color(algoname),
            "mark" => "*",
            "mark options" => ""
        )
    elseif algoname == "T2"
        res = Dict(
            "smooth" => nothing,
            FigsIteratesAlgoCurveswidth => nothing,
            # "fill opacity" => "0",
            "mark size" => "1.5pt",
            "color" => get_alg_color(algoname),
            "mark" => "star",
            # "mark options" => ""
        )
    elseif algoname == "MFISTA"
        res = Dict(
            "smooth" => nothing,
            FigsIteratesAlgoCurveswidth => nothing,
            "mark size" => "1pt",
            "color" => get_alg_color(algoname),
            "mark" => "pentagon*",
            "mark options" => ""
        )
    end

    if !haskey(res, "color")
        (res["color"] = get_alg_color(algoname))
        res["mark options"] = "{fill=white}"
    end

    return res
end

function get_suboptimality_algoparams(algoname)
    res = Dict{Any, Any}(
        "mark" => "none",
        "color" => get_alg_color(algoname),
    )

    if algoname == "ISTA"
        merge!(res, Dict(
            "dotted" => nothing
        ))
    elseif algoname == "FISTA"
        merge!(res, Dict(
            "densely dotted" => nothing
        ))
    elseif algoname == "T1"
        merge!(res, Dict(
            "dashed" => nothing
        ))
    elseif algoname == "T2"
        # merge!(res, Dict(
        # ))
    elseif algoname == "MFISTA"
        # merge!(res, Dict(
        #     "dashed" => nothing
        # ))
    end

    return res
end

function get_legend(algoname)
    res = algoname

    if algoname == "ISTA"
        res = "Proximal Gradient"
    elseif algoname == "FISTA"
        res = "Accel. Proximal Gradient"
    elseif algoname == "T1"
        res = "Prov. Alg -- \$\\mathsf{T}^1\$"
    elseif algoname == "T2"
        res = "Prov. Alg -- \$\\mathsf{T}^2\$"
    elseif algoname == "MFISTA"
        res = "Monotone APG"
    end

    return res
end


text:
"""
    distance_lat_lon_euclidean(coords1::Location2D, coords2::Location2D)

Compute the Euclidean distance in metres between two lat-long coordinates close to each other.
"""
function distance_lat_lon_euclidean(coords1::Location2D, coords2::Location2D)
    deglen = 110.25
    x = coords1[1]- coords2[1]
    y = (coords1[2] - coords2[2])*cos(coords2[1])
    return deglen*sqrt(x^2 + y^2)*1000.0
end

Distances.evaluate(::EuclideanLatLongMetric, coords1::Location2D, coords2::Location2D) = distance_lat_lon_euclidean(coords1, coords2)

"""
    get_location2D_list(node_attribs_file::String)

Returns the list of parsed 2D locations for the road graph nodes.
"""
function get_location2D_list(node_attribs_file::String)

    node_attribs_dict = JSON.parsefile(node_attribs_file)
    nnodes = length(keys(node_attribs_dict))

    location_list = Vector{Location2D}(undef, nnodes)
    for (k, v) in node_attribs_dict
        location_list[parse(Int64,k)] = Location2D(parse(Float64, v["y"]), parse(Float64, v["x"]))
    end

    return location_list
end

function compute_bidir_astar_euclidean(g::LightGraphs.AbstractGraph, s::Int64, t::Int64, distmx::AbstractMatrix{Float64}, node_locations::Vector{Location2D}, nbrs_to_exclude::Set{Int64}=Set{Int64}())

    if s == t
        return (path=[s], dist=0.0)
    end

    metric = EuclideanLatLongMetric()
    fwd_heuristic(u, t) = Distances.evaluate(metric, node_locations[u], node_locations[t])
    bwd_heuristic(u, s) = Distances.evaluate(metric, node_locations[s], node_locations[u])

    bidir_astar = ShortestPaths.BidirAStar(fwd_heuristic, bwd_heuristic, nbrs_to_exclude)

    return ShortestPaths.shortest_paths(g, s, t, distmx, bidir_astar)
end

"""
    get_non_dominated_transit_points(env::CoordinatedMAPFEnv, gstart::Int64, gend::Int64, avoid_vertex_set::Set{Int64}, v::AerialMAPFState, flight_heur::Bool)

Returns the ground transit points that are non-dominated on both time and flight distance.
"""
function get_non_dominated_transit_points(env::CoordinatedMAPFEnv, gstart::Int64, gend::Int64, avoid_vertex_set::Set{Int64}, v::AerialMAPFState, flight_heur::Bool)

    time_dist_set = NamedTuple{(:idx, :timeval, :distval), Tuple{Int64, Float64, Float64}}[]
    vgtg = env.ground_transit_graph.vertices[v.ground_transit_idx]

    for seq = gstart:gend-1
        # As soon as you hit a conflict vertex, stop
        if (seq in avoid_vertex_set)
            break
        end
        seq_gtg = env.ground_transit_graph.vertices[seq]

        ddiff = flight_heur ? distance_lat_lon_euclidean(env.location_list[vgtg.road_vtx_id], env.location_list[seq_gtg.road_vtx_id]) : env.ground_transit_weights[(v.ground_transit_idx, seq_gtg.idx)]
        tval = max(v.timeval + ddiff/AvgSpeed, seq_gtg.timeval)

        push!(time_dist_set, (idx=seq_gtg.idx, timeval=tval, distval=v.distval+ddiff))
    end # seq = gstart:gend-1

    if isempty(time_dist_set)
        return []
    end

    # Already sorted in increasing time
    non_dom_idxs = [1]

    for (i, tup) in enumerate(time_dist_set)
        dom = false
        for ndi in non_dom_idxs
            if tup.distval >= time_dist_set[ndi].distval
                dom=true
                break
            end
        end

        if dom==false
            push!(non_dom_idxs, i)
        end
    end

    non_dom_set = time_dist_set[non_dom_idxs]

    return non_dom_set
end

"""
    get_best_transit_connection(env::CoordinatedMAPFEnv, gstart::Int64, gend::Int64, v::AerialMAPFState, flight_heur::Bool)

Returns the lowest-cost ground transit connection point.
"""
function get_best_transit_connection(env::CoordinatedMAPFEnv, gstart::Int64, gend::Int64, v::AerialMAPFState, flight_heur::Bool)

    best_weighted_cost = Inf

    vgtg = env.ground_transit_graph.vertices[v.ground_transit_idx]

    result = (idx=0, timeval=0.0, distval=0.0)

    for seq = gstart:gend-1
        # As soon as you hit a conflict vertex, stop
        seq_gtg = env.ground_transit_graph.vertices[seq]
        ddiff = flight_heur ? distance_lat_lon_euclidean(env.location_list[vgtg.road_vtx_id], env.location_list[seq_gtg.road_vtx_id]) : env.ground_transit_weights[(v.ground_transit_idx, seq_gtg.idx)]
        tval = max(v.timeval + ddiff/AvgSpeed, seq_gtg.timeval)
        dval = v.distval + ddiff
        tdiff = tval - seq_gtg.timeval

        weighted_cost = env.alpha_weight_distance*ddiff + (1.0-env.alpha_weight_distance)*tdiff
        if weighted_cost < best_weighted_cost
            best_weighted_cost = weighted_cost
            result = (idx=seq_gtg.idx, timeval=tval, distval=dval)
        end
    end # seq = gstart:gend-1

    if best_weighted_cost < Inf
        return result
    else
        return nothing
    end
end

function get_tasks_with_valid_path(road_graph::LightGraphs.AbstractGraph, road_graph_wts::AbstractMatrix{Float64}, location_list::Vector{Location2D}, num_tasks::Int64, rng::AbstractRNG)

    nvg = nv(road_graph)
    task_list = Vector{AgentTask}(undef, num_tasks)

    for t = 1:num_tasks

        path_exists = false

        start = -1
        goal = -1

        while ~path_exists
            start = rand(rng, 1:nvg)
            goal = rand(rng, 1:nvg)

            bidir_sp = compute_bidir_astar_euclidean(road_graph, start, goal, road_graph_wts, location_list)

            path_exists = (isempty(bidir_sp.path) == false)

        end
        task_list[t] = AgentTask((start, goal))
    end

    return task_list
end

"""
    Computes sum-of-costs of the paths.
"""
function compute_total_cost(env::CoordinatedMAPFEnv, ground_paths::Vector{AgentPathInfo}, aerial_paths::Vector{AgentPathInfo})

    total_dist = sum(ap.total_dist for ap in aerial_paths) + sum(gp.total_dist for gp in ground_paths)
    total_time = sum(ap.total_time for ap in aerial_paths) + sum(gp.total_time for gp in ground_paths)

    return env.alpha_weight_distance*total_dist + (1.0-env.alpha_weight_distance)*total_time
end

function count_cars_not_coordinating(ground_paths::Vector{AgentPathInfo})

    not_coord_ids = Set{Int64}()

    for (gid, gpath) in enumerate(ground_paths)

        is_coord = false
        for gps in gpath.path_states
            if ~isempty(gps.coord_agent_id)
                is_coord = true
                break
            end
        end
        if ~is_coord
            push!(not_coord_ids, gid)
        end
    end
    return not_coord_ids
end

function get_increasing_dist_ordering(agent_paths::Vector{AgentPathInfo})

    id_dist_tuples = [(id, path.total_dist) for (id, path) in enumerate(agent_paths)]
    sort!(id_dist_tuples, by=x->x[2], rev=true)

    return [idd[1] for idd in id_dist_tuples]
end


text:
using Sundials, Test, ModelingToolkit, LinearAlgebra, IncompleteLU
import AlgebraicMultigrid

const N = 32
const xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a

function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end
p = (3.4, 1., 10., step(xyd_brusselator))

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0,(0.,11.5),p)

prob_ode_brusselator_2d_mtk = ODEProblem(modelingtoolkitize(prob_ode_brusselator_2d),[],(0.0,11.5),jac=true,sparse=true);

u0 = prob_ode_brusselator_2d_mtk.u0
p  = prob_ode_brusselator_2d_mtk.p
const jaccache = prob_ode_brusselator_2d_mtk.f.jac(u0,p,0.0)
const W = I - 1.0*jaccache

prectmp = ilu(W, τ = 50.0)
const preccache = Ref(prectmp)

function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
  if jok
    prob_ode_brusselator_2d_mtk.f.jac(jaccache,u,p,t)
    jcurPtr[] = true

    # W = I - gamma*J
    @. W = -gamma*jaccache
    idxs = diagind(W)
    @. @view(W[idxs]) = @view(W[idxs]) + 1

    # Build preconditioner on W
    preccache[] = ilu(W, τ = 5.0)
  end
end

function precilu(z,r,p,t,y,fy,gamma,delta,lr)
  ldiv!(z,preccache[],r)
end

prectmp2 = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W, presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,1))), postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,1)))))
const preccache2 = Ref(prectmp2)
function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
  if jok
    prob_ode_brusselator_2d_mtk.f.jac(jaccache,u,p,t)
    jcurPtr[] = true

    # W = I - gamma*J
    @. W = -gamma*jaccache
    idxs = diagind(W)
    @. @view(W[idxs]) = @view(W[idxs]) + 1

    # Build preconditioner on W
    preccache2[] = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W, presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,1))), postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,1)))))
  end
end

function precamg(z,r,p,t,y,fy,gamma,delta,lr)
  ldiv!(z,preccache2[],r)
end

sol1 = solve(prob_ode_brusselator_2d,CVODE_BDF(linear_solver=:GMRES),save_everystep=false);
sol2 = solve(prob_ode_brusselator_2d,CVODE_BDF(linear_solver=:GMRES,prec=precilu,psetup=psetupilu,prec_side=1),save_everystep=false);
sol3 = solve(prob_ode_brusselator_2d,CVODE_BDF(linear_solver=:GMRES,prec=precamg,psetup=psetupamg,prec_side=1),save_everystep=false);
@test sol1.destats.nf > sol2.destats.nf
@test sol1.destats.nf > sol3.destats.nf


text:
using LinearAlgebra
using ProximalOperators
using Test

@testset "HuberLoss" begin

f = HuberLoss(1.5, 0.7)

predicates_test(f)

@test ProximalOperators.is_smooth(f) == true
@test ProximalOperators.is_quadratic(f) == false
@test ProximalOperators.is_set(f) == false

x = randn(10)
x = 1.6*x/norm(x)

call_test(f, x)
prox_test(f, x, 1.3)
grad_fx, fx = gradient_test(f, x)

@test abs(fx - f(x)) <= 1e-12
@test norm(0.7*1.5*x/norm(x) - grad_fx, Inf) <= 1e-12

x = randn(10)
x = 1.4*x/norm(x)

call_test(f, x)
prox_test(f, x, 0.9)
grad_fx, fx = gradient_test(f, x)

@test abs(fx - f(x)) <= 1e-12
@test norm(0.7*x - grad_fx, Inf) <= 1e-12

end


text:
@testset "Scheduling" begin

    @testset "Scheduling: FCFS" begin
        n = 3
        process_id = Any[1, 2, 3]
        burst_times = Any[10, 5, 8]
        @test fcfs(n, process_id, burst_times) == (Any[1, 2, 3], Any[10, 5, 8], Any[0, 10, 15], Any[10, 15, 23], 8.333333333333334, 16.0)
    end
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
module CommunalHelperDreamBooster

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/DreamBooster" DreamBooster(
    x::Integer,
    y::Integer,
    hidePath::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Dream Booster (Communal Helper)" => Ahorn.EntityPlacement(
        DreamBooster,
        "line",
    ),
    "Dream Booster (Hidden Path) (Communal Helper)" => Ahorn.EntityPlacement(
        DreamBooster,
        "line",
        Dict{String,Any}(
            "hidePath" => true,
        ),
    ),
)

Ahorn.nodeLimits(entity::DreamBooster) = 1, 1

function Ahorn.selection(entity::DreamBooster)
    x, y = Ahorn.position(entity)
    endX, endY = Int.(entity.data["nodes"][1])

    return [
        Ahorn.Rectangle(x - 9, y - 9, 18, 18),
        Ahorn.Rectangle(endX - 9, endY - 9, 18, 18),
    ]
end

const dreamBoosterSprite = "objects/CommunalHelper/boosters/dreamBooster/idle00"

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamBooster, room::Maple.Room)
    x, y = Ahorn.position(entity)
    node = entity.data["nodes"][1]

    Ahorn.drawLines(ctx, [(x, y), node], (1.0, 1.0, 1.0, 1.0), thickness=1)
    Ahorn.drawSprite(ctx, dreamBoosterSprite, x, y)
end

end


text:
#=

prec=precision(midpt); # bits of significand (signficant radix-2 digits
prec, typeof(midpt)
# (57,ArbFloats.ArbFloat{57})

1 + round(Int,log2(ufp2(midpt))) -  round(Int,log2(ulp2(midpt)))
# 57

prec == (1 + round(Int,log2(ufp2(midpt))) -  round(Int,log2(ulp2(midpt))))
# true

=#

"""
   x.midpoint -> (significand, exponentOf2)
                  [0.5,1.0)     2^expo
   x.radius   -> (radial significand, radial exponentOf2)
"""
function frexp{P}(x::ArfFloat{P})
    exponentOf2 = x.exponentOf2
    significandOf = deepcopy(x)
    significandOf.exponentOf2 = 0
    return significandOf, exponentOf2
end

function ldexp{P}(s::ArfFloat{P}, e::Int)
    z = deepcopy(s)
    z.exponentOf2 = e
    return z
end

function frexp{P}(x::ArbFloat{P})
    significandOf, exponentOf2 = frexp(ArfFloat{P}(x))
    return ArbFloat{P}(significandOf), exponentOf2
end

function ldexp{P}(s::ArbFloat{P}, e::Int)
    z = deepcopy(s)
    z.exponentOf2 = e
    return z
end
ldexp{P}(x::Tuple{ArbFloats.ArbFloat{P}, Int}) = ldexp(x[1],x[2])

#=

      gte_bits2digs(b_bits)    ~>  d_digits
      d_digits | d digits suffice to encode and recover b bits without error

      gte_digs2bits(d_digits)  ~>  b_bits
      b_bits |  b bits suffice to encode and recover d digits without error

      lte_bits2digs(b_bits)    ~>  c_digits
      c_digits | b bits are necessary encode and recover c digits without error

      lte_digs2bits(d_digits)  ~>  a_bits
      a_bits |  d digits are necessary encode and recover a bits without error

      lte_digs2bits( gte_bits2digs(b_bits  ) ) == b_bits
      lte_bits2digs( gte_digs2bits(d_digits) ) == d_digits

const log2_log10 = log(10,2)  # 0.3010299956639812 ~= log(2)/log(10)
const log10_log2 = log(2,10)  # 3.321928094887362  ~= log(10)/log(2)

lte_bits2digs(nbits::Int) = floor(Int, nbits * log2_log10)
lte_digs2bits(ndigs::Int) = floor(Int, ndigs * log10_log2)
=#

lte_bits2digs(nbits::Int) = floor( Int, nbits * 0.30102999566398125 ) # log(10,2) RoundUp
gte_bits2digs(nbits::Int)  = ceil( Int, nbits * 0.3010299956639812  ) # log(10,2) RoundDown
lte_digs2bits(ndigs::Int) = floor( Int, ndigs * 3.3219280948873626  ) # log(2,10) RoundUp
gte_digs2bits(ndigs::Int)  = ceil( Int, ndigs * 3.3219280948873617  ) # log(2,10) RoundDown

"""
logarithm_base(x)
"""
function log_base(x::Real, base::Int)
   z = if base == 2
           log2(x)
        elseif base == 10
           log10(x)
        else
           log(x) / log(base)
        end
   return z
end
log_base{P}(x::ArbFloat{P}, base::Int) = ArbFloats.logbase(x,base)
log_base{T<:ArbFloat}(x::T, base::Int) = ArbFloats.logbase(x,base)

#=
    position_first_place: the radix position of the most significant nonzero bit|digit

    pfp{T<:AbstractFloat}(x::T, base::Int=2)
    pfp{T<:AbstractFloat}(x::T)              == pfp(x,  2)  ==  pfp2{T<:AbstractFloat}(x::T)
    pfp{T<:AbstractFloat}(x::T, base=10)     == pfp(x, 10)  ==  pfp10{T<:AbstractFloat}(x::T)

    position_last_place: the radix position of the least significant nonzero bit|digit

    plp{T<:AbstractFloat}(x::T, base::Int=2)
    plp{T<:AbstractFloat}(x::T)              == p;p(x,  2)  ==  plp2{T<:AbstractFloat}(x::T)
    plp{T<:AbstractFloat}(x::T, base=10)     == plp(x, 10)  ==  p;p10{T<:AbstractFloat}(x::T)

    unit_first_place: the radix *value) of the most significant nonzero bit|digit

    pfp{T<:AbstractFloat}(x::T, base::Int=2)
    pfp{T<:AbstractFloat}(x::T)              == pfp(x,  2)  ==  pfp2{T<:AbstractFloat}(x::T)
    pfp{T<:AbstractFloat}(x::T, base=10)     == pfp(x, 10)  ==  pfp10{T<:AbstractFloat}(x::T)

    unit_last_place: the radix *value* of the least significant nonzero bit|digit

    plp{T<:AbstractFloat}(x::T, base::Int=2)
    plp{T<:AbstractFloat}(x::T)              == plp(x,  2)  ==  plp2{T<:AbstractFloat}(x::T)
    plp{T<:AbstractFloat}(x::T, base=10)     == plp(x, 10)  ==  plp10{T<:AbstractFloat}(x::T)

=#

"""
position_first_place
determine the position of the most significant nonzero bit|digit
"""
function pfp{T<:Real}(x::T, base::Int=2)
   z = 0 # if x==0.0
   if notzero(x)
       z = floor( Int, log_base(abs(x), base) )
   end
   return z
end
pfp{P}(x::ArbFloat{P}, base::Int=2) =
    return iszero(x) ? 0 : floor( Int, log_base(abs(smartvalue(x)), base) )
"""
binary position_first_place
determine the position of the most significant nonzero bit
"""
pfp2{T<:Real}(x::T) = (x==zero(T) ? 0 : floor( Int, log2(abs(x)) ))
pfp2{P}(x::ArbFloat{P}) = (iszero(x) ? 0 : floor( Int, log2(abs(smartvalue(x))) ))
"""
decimal position_first_place
determine the position of the most significant nonzero digit
"""
pfp10{T<:Real}(x::T) = (x==zero(T) ? 0 : floor( Int, log10(abs(x)) ))
pfp10{P}(x::ArbFloat{P}) = (iszero(x) ? 0 : floor( Int, log10(abs(smartvalue(x))) ))

"""
ufp is unit_first_place
the float value given by a 1 at the position of
  the most significant nonzero bit|digit in _x_
"""
function ufp(x::AbstractFloat, base::Int=2)
   z = pfp(x, base)
   b = convert(Float64, base)
   return b^z
end
function ufp{P}(x::ArbFloat{P}, base::Int=2)
   z = pfp(x, base)
   return Float64(base)^z
end
ufp(x::Integer, base::Int=2) = ufp(Float64(x), base)
"ufp2 is unit_first_place in base 2"
ufp2{T<:Real}(x::T) = 2.0^pfp2(x)
ufp2{P}(x::ArbFloat{P}) = 2.0^pfp2(x)
ufp2(x::Integer) = ufp2(Float64(x))
"ufp10 is unit_first_place in base 10"
ufp10{T<:Real}(x::T) = 10.0^pfp10(x)
ufp10{P}(x::ArbFloat{P}) = 10.0^pfp10(x)
ufp10(x::Integer) = ufp10(Float64(x))
"""
ulp   is unit_last_place
the float value given by a 1 at the position of
  the least significant nonzero bit|digit in _x_
"""
function ulp(x::Real, precision::Int, base::Int)
   unitfp  = base==2 ? ufp2(x) : (base==10 ? ufp10(x) : throw(DomainError()))
   twice_u = 2.0^(1-precision)
   return twice_u * unitfp
end
ulp{T<:AbstractFloat}(x::T, base::Int=2)  =
    ulp(x, 1+Base.significand_bits(T), base)
ulp{P}(x::ArbFloat{P}, base::Int=2)  =
    ulp(x, P, base)
ulp(x::Integer, base::Int=2) = ulp(Float64(x), base)

"""ulp2  is unit_last_place base 2"""
function ulp2(x::Real, precision::Int)
   unitfp  = ufp2(x)
   twice_u = 2.0^(1-precision)
   return (*)(promote(twice_u, unitfp)...)
end
function ulp2{P}(x::ArbFloat{P})
   unitfp  = ufp2(x)
   twice_u = 2.0^(1-P)
   return (*)(promote(twice_u, unitfp)...)
end
ulp2{T<:AbstractFloat}(x::T)  = ulp2(x, 1+Base.significand_bits(T))
ulp2(x::Integer) = ulp2(Float64(x))

"""ulp10 is unit_last_place base 10"""
function ulp10(x::Real, bitprecision::Int)
    unit10fp = ufp10(x)
    digitprecision = lte_bits2digs(bitprecision)
    twice_u10 = 10.0^(1-digitprecision)
    return twice_u10 * unit10fp
end
function ulp10{P}(x::ArbFloat{P})
    unit10fp = ufp10(x)
    digitprecision = lte_bits2digs(P)
    twice_u10 = 10.0^(1-digitprecision)
    return twice_u10 * unit10fp
end
ulp10{T<:AbstractFloat}(x::T) = ulp10( x, (1+Base.significand_bits(T)) )
ulp10(x::Integer) = ulp10(Float64(x))

function eps{T<:ArbFloat}(x::T)
    ieps = internal_eps(x)
    return T(ieps)
end
eps{P}(x::ArbFloat{P}) = ArbFloat{P}( internal_eps(x) )
eps{T<:ArbFloat}(::Type{T}) = T(internal_eps(T))

internal_eps{T<:ArbFloat}(::Type{T}) = ldexp(0.5,1-precision(T)) # for intertype workings
internal_eps{P}(::Type{ArbFloat{P}}) = ldexp(0.5,1-P) # for intertype workings
function internal_eps{P}(x::ArbFloat{P})              # for intratype workings
    m,r = midpoint(x), radius(x)
    z =
        if iszero(m)
            if iszero(r)
               ldexp(0.5, 1-P)
            else
               ufp2(r)
            end
        elseif iszero(r)
            ulp2(m)
        else
            max( ulp2(m), ufp2(r) )
        end
   return z
end

function nextfloat{P}(x::ArbFloat{P})
    x + ulp2(x)
end

function prevfloat{P}(x::ArbFloat{P})
    x - ulp2(x)
end


text:
"""
    food_step!(model::ABM)

Simulates growth of food on the grid
"""
function food_step!(model::ABM)
    # randomise the update order for updating tiles
    shuffle!(model.rng, model.food_update_order)
    # update tiles
    for pos in model.food_update_order
        # if tile is empty
        if model.food[pos] == 0.0
            # regenerate food if the `random_spread_chance` occurs
            rand(model.rng) <= model.food_data.random_spread_chance &&
                (model.food[pos] += model.food_data.regen_rate; continue)

            # iterate through neighboring tiles
            for offset in model.food_data.neighborhood
                # bounds check
                all((0, 0) .< Tuple(pos + offset) .< size(model.food)) || continue
                # spread to this tile with probability proportional to amount of food
                # for neighboring dead tiles, this doesn't do anything
                if rand(model.rng) <=
                   model.food_data.spread_multiplier * model.food[pos+offset] /
                   model.food_data.food_cap
                    model.food[pos+offset] += model.food_data.regen_rate
                    break
                end
            end
        else
            # the tile is a food tile, so regenerate food till the cap
            model.food[pos] =
                min(model.food[pos] + model.food_data.regen_rate, model.food_data.food_cap)
        end
    end
end

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text: # TODO: do bit hacking to make Clongdouble into Cdouble and do math with that (for now)


text:
AbstractSimplyConnectedRegion = AbstractConnectedRegion{1}

"""
	(type) SimplyConnectedRegion
Representation of a simply connected region in the extended complex plane.
	SimplyConnectedRegion(p::Union{AbstractClosedCurve,AbstractClosedPath})
Construct an open simply connected region by specifying its boundary. The region is "to the left" of the orientation of the boundary.
"""
struct InteriorSimplyConnectedRegion{T<:AbstractJordan} <: AbstractConnectedRegion{1}
	boundary::T
end

struct ExteriorSimplyConnectedRegion{T<:AbstractJordan} <: AbstractConnectedRegion{1}
	boundary::T
end

SimplyConnectedRegion = Union{InteriorSimplyConnectedRegion{T},ExteriorSimplyConnectedRegion{T}} where T<:AbstractJordan

"""
	interior(C)
Construct the region interior to the closed curve or path `C`. If `C` is bounded, the bounded enclosure is chosen regardless of the orientation of `C`; otherwise, the region "to the left" is the interior.
"""
function interior(C::AbstractJordan)
	if isfinite(C) && winding(1/C,0) > 0
		C = reverse(C)
	end
	InteriorSimplyConnectedRegion(C)
end

"""
	exterior(C)
Construct the region exterior to  the closed curve or path `C`. If `C` is bounded, the bounded enclosure is chosen regardless of the orientation of `C`; otherwise, the region "to the right" is the exterior.
"""
function exterior(C::AbstractJordan)
	if isfinite(C)
		if winding(1/C,0) < 0
			C = reverse(C)
		end
	else
		if C isa AbstractClosedPath && (length(filter(isinf,vertices(C))) > 1)
			@error "Disconnected exterior"
		end
		C = reverse(C)
	end
	ExteriorSimplyConnectedRegion(C)
end

boundary(R::SimplyConnectedRegion) = R.boundary
innerboundary(R::InteriorSimplyConnectedRegion) = nothing
innerboundary(R::ExteriorSimplyConnectedRegion) = R.boundary
outerboundary(R::InteriorSimplyConnectedRegion) = R.boundary
outerboundary(R::ExteriorSimplyConnectedRegion) = nothing
in(z::Number,R::InteriorSimplyConnectedRegion) = isinside(z,outerboundary(R))
in(z::Number,R::ExteriorSimplyConnectedRegion) = isoutside(z,innerboundary(R))
isfinite(R::InteriorSimplyConnectedRegion) = isfinite(outerboundary(R))
isfinite(R::ExteriorSimplyConnectedRegion) = false

function show(io::IO,R::InteriorSimplyConnectedRegion)
	print(IOContext(io,:compact=>true),"Region interior to ",R.boundary)
end

function show(io::IO,R::ExteriorSimplyConnectedRegion)
	print(IOContext(io,:compact=>true),"Region exterior to ",R.boundary)
end

function show(io::IO,::MIME"text/plain",R::SimplyConnectedRegion)
	show(io,R)
end

"""
	!(R::SimplyConnectedRegion)
Compute the region complementary to `R`. This is not quite set complementation, as neither region includes its boundary. The complement is always simply connected in the extended plane.
"""
!(R::InteriorSimplyConnectedRegion) = exterior(reverse(boundary(R)))
!(R::ExteriorSimplyConnectedRegion) = interior(reverse(boundary(R)))

"""
	isapprox(R1::SimplyConnectedRegion,R2::SimplyConnectedRegion; tol=<default>)
Determine whether `R1` and `R2` represent the same region, up to tolerance `tol`. Equivalently, determine whether their boundaries are the same.
"""
isapprox(R1::S,R2::T;tol=DEFAULT[:tol]) where {S,T<:SimplyConnectedRegion} = isapprox(boundary(R1),boundary(R2),tol=tol)

# disks
AbstractDisk = SimplyConnectedRegion{T} where T<:Circle
"""
	disk(C::Circle)
Construct the disk interior to `C`.
"""
disk(C::Circle) = interior(C)
"""
	disk(center::Number,radius::Real)
Construct the disk with the given `center` and `radius`.
"""
disk(center::Number,radius::Real) = interior(Circle(center,radius))
unitdisk = disk(complex(0.0),1.0)
function show(io::IO,::MIME"text/plain",R::AbstractDisk)
	side = in(Inf,R) ? "exterior" : "interior"
	print(io,"Disk $side to:\n   ",R.boundary)
end

# half-planes
AbstractHalfplane = SimplyConnectedRegion{T} where T<:Line
"""
	halfplane(L::Line)
Construct the half-plane to the left of `L`.
"""
halfplane(L::Line) = interior(L)
"""
	halfplane(a,b)
Construct the half-plane to the left of the line from `a` to `b`.
"""
halfplane(a::Number,b::Number) = interior(Line(a,b))
upperhalfplane = halfplane(Line(0.0,direction=1.0))
lowerhalfplane = halfplane(Line(0.0,direction=-1.0))
lefthalfplane = halfplane(Line(0.0,direction=1.0im))
righthalfplane = halfplane(Line(0.0,direction=-1.0im))
function show(io::IO,::MIME"text/plain",R::AbstractHalfplane)
	print(io,"Half-plane to the left of:\n   ",R.boundary)
end

"""
	(type) PolygonalRegion
Representation of a simply connected region bounded by a plane.
"""
PolygonalRegion = SimplyConnectedRegion{Polygon}