text:
@testset "Events/Notation Declarations" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    DC = E.DataContent
    ME = E.MarkupError

    function NotationDeclaration(name, public_id, system_id, location)
        return E.NotationDeclaration(name, E.ExternalIdentifier(public_id, system_id, location), location)
    end

    function NotationDeclaration(name, location)
        return E.NotationDeclaration(name, nothing, location)
    end

    @testset "Events/Notation Declarations (Positive)" begin
        # Basic tests ... fixed public/system identifier.
        #
        @test (evaluate("<!NOTATION a PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("a", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION abc PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("abc", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION é PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("é", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION aé PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("aé", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION éb PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("éb", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION aéb PUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("aéb", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION\na\nPUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("a", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION\na\tPUBLIC \"hello.not\" \"salut.not\">")
               == [ NotationDeclaration("a", "hello.not", "salut.not", L.Location("a buffer", -1)) ])

        # Basic tests ... fixed notation name ...
        #
        @test (evaluate("<!NOTATION a SYSTEM \"\">") == [ NotationDeclaration("a", nothing, "", L.Location("a buffer", -1)) ])
        @test (evaluate("<!NOTATION a\nSYSTEM\t\"\"	>")
               == [ NotationDeclaration("a", nothing, "", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a SYSTEM \"   \">") == [ NotationDeclaration("a", nothing, "   ", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC \"\" \"\">") == [ NotationDeclaration("a", "", "", L.Location("a buffer", -1)) ])
        @test (evaluate("<!NOTATION a\nPUBLIC\r\"\"\t\"\"	>")
               == [ NotationDeclaration("a", "", "", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC \"   \" \"   \">")
               == [ NotationDeclaration("a", "   ", "   ", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a SYSTEM \"hello.not\">")
               == [ NotationDeclaration("a", nothing, "hello.not", L.Location("a buffer", -1)) ])

        # ... fixed notation name with system identifier ...
        #
        @test (evaluate("<!NOTATION a SYSTEM \"hello.not\">")
               == [ NotationDeclaration("a", nothing, "hello.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a SYSTEM 'hello.not'>")
               == [ NotationDeclaration("a", nothing, "hello.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a SYSTEM '\"hello.not\"'>")
               == [ NotationDeclaration("a", nothing, "\"hello.not\"", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a SYSTEM \"'hello.not'\">")
               == [ NotationDeclaration("a", nothing, "'hello.not'", L.Location("a buffer", -1)) ])

        # ... fixed notation name with public/system identifier ...
        #
        @test (evaluate("<!NOTATION a PUBLIC \"salut.not\" \"hello.not\">")
               == [ NotationDeclaration("a", "salut.not", "hello.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC 'salut.not' 'hello.not'>")
               == [ NotationDeclaration("a", "salut.not", "hello.not", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC '\"salut.not\"' '\"hello.not\"'>")
               == [ NotationDeclaration("a", "\"salut.not\"", "\"hello.not\"", L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC \"'salut.not'\" \"'hello.not'\">")
               == [ NotationDeclaration("a", "'salut.not'", "'hello.not'", L.Location("a buffer", -1)) ])

        # ... fixed notation name with public identifier ...
        #
        @test (evaluate("<!NOTATION a PUBLIC \"salut.not\">")
               == [ NotationDeclaration("a", "salut.not", nothing, L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC 'salut.not'>")
               == [ NotationDeclaration("a", "salut.not", nothing, L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC '\"salut.not\"'>")
               == [ NotationDeclaration("a", "\"salut.not\"", nothing, L.Location("a buffer", -1)) ])

        @test (evaluate("<!NOTATION a PUBLIC \"'salut.not'\">")
               == [ NotationDeclaration("a", "'salut.not'", nothing, L.Location("a buffer", -1)) ])
    end

    @testset "Events/Notation Declarations (Negative ... invalid or absent notation name.)" begin
        @test (evaluate("<!NOTATION")
               == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION>")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)),
                           DC(">", false, L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION >")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)),
                           DC(">", false, L.Location("a buffer", -1)) ])

        # This one is a little weird: we're picking up "PUBLIC" as the notation name. I'll have to figure out later if
        # that's allowed or not, but there doesn't appear to be anything in the specification that forbids it.
        #
        events = evaluate("<!NOTATION PUBLIC")
        @test length(events) == 3
        @test (events == [ ME("ERROR: A notation declaration must specify an external identifier.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("PUBLIC", L.Location("a buffer", -1)) ])

        # This one is a little weird: we're picking up "SYSTEM" as the notation name. I'll have to figure out later if
        # that's allowed or not, but there doesn't appear to be anything in the specification that forbids it.
        #
        events = evaluate("<!NOTATION SYSTEM")
        @test length(events) == 3
        @test (events == [ ME("ERROR: A notation declaration must specify an external identifier.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("SYSTEM", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION \"\"")
        @test length(events) == 3
        @test (events == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)),
                           DC("\"", false, L.Location("a buffer", -1)),
                           DC("\"", false, L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION ''")
        @test length(events) == 3
        @test (events == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)),
                           DC("'", false, L.Location("a buffer", -1)),
                           DC("'", false, L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION <")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting a notation name.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/Notation Declarations (Negative ... missing/invalid notation value.)" begin
        events = evaluate("<!NOTATION n>")
        @test length(events) == 2
        @test (events == [ ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n >")
        @test length(events) == 2
        @test (events == [ ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n\n>")
        @test length(events) == 2
        @test (events == [ ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n PUBLIC")
        @test length(events) == 4
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n PUBLIC ")
        @test length(events) == 4
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n PUBLIC <")
        @test length(events) == 5
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n SYSTEM")
        @test length(events) == 4
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n SYSTEM ")
        @test length(events) == 4
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n SYSTEM <")
        @test length(events) == 5
        @test (events == [ ME("ERROR: Expecting a quoted string.", L.Location("a buffer", -1)),
                           ME("ERROR: A notation declaration must specify an external identifier.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/Notation Declarations (Negative ... missing TAGC.)" begin
        events = evaluate("<!NOTATION n PUBLIC \"salut.not\" \"hello.not\"")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting '>' to end a notation declaration.",
                              L.Location("a buffer", -1)),
                           NotationDeclaration("n", "salut.not", "hello.not", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n PUBLIC \"salut.not\" \"hello.not\" ")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", "salut.not", "hello.not", L.Location("a buffer", -1)) ])

        events = evaluate("<!NOTATION n PUBLIC \"salut.not\" \"hello.not\" <")
        @test length(events) == 3
        @test (events == [ ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", "salut.not", "hello.not", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/Notation Declarations (Negative ... missing white space between public and system identifier.)" begin
        events = evaluate("<!NOTATION n PUBLIC \"salut.not\"\"hello.not\"")
        @test length(events) == 3
        @test (events == [ ME("ERROR: Expecting white space following a public identifier.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting '>' to end a notation declaration.", L.Location("a buffer", -1)),
                           NotationDeclaration("n", "salut.not", "hello.not", L.Location("a buffer", -1)) ])
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
function evaluate()
  periodic_boundary()
  reorder()

  Threads.@threads for i = 1:Npart
    kernel_summation(i)
  end

  if gradient == "IntegralApproach"
    Threads.@threads for i = 1:Npart
      sph[i].c11 = 0
      # neighbour loop
      s = index_list[i]

      for t = (s - Nngb) : (s + Nngb)
        t = order_periodic(t)
        k = index_order[t]

        sph[i].c11 += correction_matrix(i, k)
      end
    end
  end

  Threads.@threads for i = 1:Npart
    # refresh
    sph[i].F1 = 0
    sph[i].dU = 0
    sph[i].dalpha = 0
    divv = 0

    # neighbour loop
    s = index_list[i]

    for t = (s - Nngb) : (s + Nngb)
      t = order_periodic(t)
      k = index_order[t]
      sph[i].F1 += force(i, k)
      sph[i].dU += dU(i, k)
      if time_dependent_viscosity
        divv += divergence_v(i, k)
      end
    end
    if time_dependent_viscosity
      sph[i].dalpha = dalpha_visc(i, divv)
    end
  end
end

function periodic_boundary()
  for i = 1:Npart
    if sph[i].x < x_min
      sph[i].x += lbox
    elseif sph[i].x > x_max
      sph[i].x -= lbox
    end
  end
end

function reorder()
  x_array = []
  for i = 1:Npart
    push!(x_array, sph[i].x)
  end
  global index_order = sortperm(x_array)
  global index_list = zeros(Int, Npart)
  for i = 1:Npart
    j = index_order[i]
    index_list[j] = i
  end
end

function order_periodic(j)
  if j < 1
    j += Npart
  elseif j > Npart
    j -= Npart
  end
  return j
end

function kernel_summation(i)
  # search neighbour, adjust smoothing length and estimate density
  evaluate_hsml(i)
  density(i, volume_element)
  # EoS
  eos(i, volume_element)
end

function evaluate_hsml(i)
  h0 = eta_hsml * (sph[i].m / sph[i].rho)
  CHA = 1
  iter = 0
  while CHA > TOL
    sph[i].hsml = h0
    density(i, "mass")
    h1 = h0 - g(i) / dgdh(i)
    CHA = abs(h0 - h1) / 2 / (h0 + h1)
    h0 = h1
    iter += 1
    if iter > 1e4
      sph[i].hsml = h0
      println("iteration loop of hsml evaluation exceeds 1e4.")
      @printf("i = %d, h0 = %.7f, CHA = %.7f, g(i) = %.7f, rho = %.7f\n", i, h0, CHA, g(i), density(i, "mass"))
      break
    end
  end
end

function g(i)
  return sph[i].hsml - eta_hsml * sph[i].m / sph[i].rho
end

function dgdh(i)
  drhodh = 0
  sph[i].dydh = 0
  s = index_list[i]
  for t = (s - Nngb) : (s + Nngb)
    t = order_periodic(t)
    k = index_order[t]
    drhodh += sph[k].m * dWdh(i, k)
    sph[i].dydh += Z(k) * dWdh(i, k)
  end
  dgdh = 1 + eta_hsml * sph[i].m / sph[i].rho^2 * drhodh
  return dgdh
end

function density(i::Int, Z_input)
  y = 0

  # particle "i" is sth particle from left
  s = index_list[i]

  # take kernel summation for 2Nngb particles around particle "i"
  for t = (s - Nngb) : (s + Nngb)
    # consider tth particle from left
    t = order_periodic(t)
    # k is the index number of tth particle from left
    k = index_order[t]
    if Z_input == "mass"
      y += sph[k].m * W(i, k)
    elseif Z_input == "U"
      y += sph[k].U * (gamma - 1) * W(i, k)
    else
      @assert 0
    end
  end

  if Z_input == "mass"
    sph[i].rho = y
  elseif Z_input == "U"
    sph[i].P = y
  else
    @assert 0
  end
end

function eos(i::Int, Z_input)
  if sph[i].U < 0
    @show i, sph[i].x, sph[i].p / sph[i].m
  end
  @assert sph[i].U > 0
  if Z_input == "mass"
    sph[i].P = (gamma - 1) * sph[i].rho / sph[i].m * sph[i].U
  elseif Z_input == "U"
    sph[i].rho = sph[i].P / (gamma - 1) * sph[i].m / sph[i].U
  else
    @assert 0
  end
end

ext: jl


text:
struct IntrinsicGenerator{T}
    triangle ::Triangle{T}
    path     ::Vector{Int}
    lengths  ::Vector{T}
end

IntrinsicGenerator(t::Triangle{T}, path, lens) where T =
    IntrinsicGenerator{T}(t, path, lens)

function Base.show(io::IO, gen::IntrinsicGenerator{T}) where T
    compact = get(io, :compact, false)
    if !compact
        print(io, "IntrinsicGenerator{$T}(")
        showcompact(io, gen.triangle)
        print(io, ", path: $(gen.path), lengths: $(gen.lengths))")
    else
        print(io, "IntrinsicGenerator(")
        showcompact(io, gen.triangle)
        print(io, ", n = $(length(gen.path)), l = $(round(sum(gen.lengths), 5)))")
    end
end

struct IntrinsicPersistenceResult{T, D, I<:IntrinsicMetric{T, D}}
    generators ::Vector{IntrinsicGenerator{T}}
    points     ::Vector{SVector{D, T}}
    subset     ::Vector{Int}
    metric     ::I
end

function Base.show(io::IO, res::IntrinsicPersistenceResult{T}) where T
    println(io, "IntrinsicPersistenceResult{$T}:")
    println(io, "  Number of points: $(length(res.points)) ",
            "($(length(res.subset)) selected)")
    print(io, "  Generators ($(length(res.generators))):")
    for gen in res.generators
        print(io, "\n    ")
        showcompact(io, gen)
    end
end

function Base.convert(::Type{PersistenceBarcode},
                      res::IntrinsicPersistenceResult{T}) where T
    pairs = PersistencePair{T, Void}[]
    for gen in res.generators
        p1 = maximum(gen.lengths)
        p2 = sum(gen.lengths)
        push!(pairs, PersistencePair(p1, p2))
    end
    PersistenceBarcode(PersistencePair{T, Void}[], pairs)
end

# ============================================================================ #
# Plots
# ============================================================================ #
@recipe function f(res::IntrinsicPersistenceResult)
    allpts = res.points
    pts    = allpts[res.subset]
    dim    = length(eltype(pts))
    2 ≤ dim ≤ 3 || throw(ArgumentError("Can't plot $dim-d points!"))

    @series begin
        seriestype := :scatter
        label --> "points"

        xs = getindex.(pts, 1)
        ys = getindex.(pts, 2)
        if dim == 2
            xs, ys
        else
            markersize --> 0.5
            zs = getindex.(pts, 3)
            xs, ys, zs
        end
    end

    @series begin
        seriestype := :scatter
        label --> "all points"
        alpha --> 0.3

        idxs = setdiff(1:size(allpts, 2), res.subset)

        xs = getindex.(allpts[idxs], 1)
        ys = getindex.(allpts[idxs], 2)
        if dim == 2
            xs, ys
        else
            markersize --> 0.5
            zs = getindex.(allpts[idxs], 2)
            xs, ys, zs
        end
    end

    for (i, gen) in enumerate(res.generators)
        @series begin
            l = sum(gen.lengths)

            seriestype := :path
            linewidth --> 5
            label --> "gen$(i); l = $(round(l, 3))"

            path = copy(gen.path)
            push!(path, path[1])

            xs = getindex.(allpts[path], 1)
            ys = getindex.(allpts[path], 2)
            if dim == 2
                xs, ys
            else
                zs = getindex.(allpts[path], 3)
                xs, ys, zs
            end
        end
    end
end


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin
    f(x) = 5x^2 - 16x + 3
    f′(x) = 10x - 16
    g(x) = x^2 + 2
    g′(x) = 2x
    @test abs(newtonroot(f, f′, x₀=4) - 3) < 1e-7
    @test abs(newtonroot(f, f′, x₀=0) - 0.2) < 1e-7
    @test abs(newtonroot(f, x₀=4) - 3) < 1e-7
    @test abs(newtonroot(f, x₀=0) - 0.2) < 1e-7
    @test abs(newtonroot(f, f′, x₀=BigFloat(0.3)) - 0.2) < 1e-7
    @test abs(newtonroot(f, x₀=BigFloat(0.3)) - 0.2) < 1e-7
    @test isnan(newtonroot(g, g′, x₀=0))
    @test isnan(newtonroot(f, f′, x₀=4, maxiter=1))
    @test abs(newtonroot(f, f′, x₀=10, tolerance=1E+4) - 3) > 1e-1
end


text:
using Measures

# 1 function to be moved to Measures.jl
Base.convert(::Type{T} , x::Measure) where T<:Real  = T(x.value)

struct DiscreteSizeScale <: Gadfly.ScaleElement
    f::Function
    levels::Union{Nothing, AbstractVector}
    order::Union{Nothing, AbstractVector}
    preserve_order::Bool
end
DiscreteSizeScale(f; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteSizeScale(f, levels, order, preserve_order)

element_aesthetics(scale::DiscreteSizeScale) = [:size]

default_discrete_sizes(n::Int) = range(0.45mm, 1.8mm, length=n)

"""
    size_discrete2(f=default_discrete_sizes;
                    levels=nothing, order=nothing, preserve_order=true)

A discrete size scale that maps the categorical values in the `size`
aesthetic to x-axis units or `Measure` units (from Measures.jl).  `f` is a function that produces a vector of size units.
`levels` are the categorical levels, and level order will be respected.  `order` is
a vector of integers giving a permutation of the levels default order.  If
`preserve_order` is `true` levels are ordered as they appear in the data.
"""
size_discrete2(f::Function=Gadfly.current_theme().discrete_sizemap; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteSizeScale(f, levels, order, preserve_order)

function Scale.apply_scale(scale::DiscreteSizeScale, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)

    d = []
    for (aes, data) in zip(aess, datas)
        data.size === nothing && continue
        append!(d, skipmissing(data.size))
    end
    levelset = unique(d)

    if scale.levels == nothing
        scale_levels = [levelset...]
        scale.preserve_order || sort!(scale_levels)
    else
        scale_levels = scale.levels
    end
    scale.order == nothing || permute!(scale_levels, scale.order)

    sizes = scale.f(length(scale_levels))

    size_map = Dict(s=>string(label) for (s, label) in zip(sizes, scale_levels))
    labeler(xs) = [size_map[x] for x in xs]
    key_vals = OrderedDict(s=>i for (i,s) in enumerate(sizes))

    for (aes, data) in zip(aess, datas)
        data.size === nothing && continue
        ds = discretize([d for d in skipmissing(data.size)], scale_levels)
        vals = sizes[ds.index]
        aes.size = discretize_make_ia(vals, sizes)
        aes.size_label = labeler
        aes.size_key_vals = key_vals
    end
end

area_transform = ContinuousScaleTransform(a->sqrt(a/π), r->π*r*r, identity_formatter)

struct ContinuousSizeScale <: Gadfly.ScaleElement
    # A function of the form f(p) where 0 ≤ p ≤ 1, that returns a Measure.
    f::Function
    trans::ContinuousScaleTransform
    minvalue::Maybe(Compose.MeasureOrNumber)
    maxvalue::Maybe(Compose.MeasureOrNumber)
    format::Union{Nothing, Symbol}
end

Scale.element_aesthetics(scale::ContinuousSizeScale) = [:size]

default_continuous_sizes(p::Float64; min=0mm, max=2mm) = min + p*(max-min)

"""
    Scale.size_radius(f=default_continuous_sizes; minvalue=0.0, maxvalue=nothing)

A scale for continuous sizes. Values in the `size` aesthetic are mapped either to:
- x-axis units (if `maxvalue=nothing`), or
-  `Measure` units (from Measures.jl).
If `maxvalue` is specified, the continuous sizes are converted to a proportion
(between `minvalue` and `maxvalue`), and then mapped to absolute sizes using the function `f(p)` where `0≤p≤1`.
"""
size_radius(f::Function=Gadfly.current_theme().continuous_sizemap; minvalue=0.0, maxvalue=nothing) =
     ContinuousSizeScale(f, identity_transform, minvalue, maxvalue, nothing)

"""
    Scale.size_area(f=default_continuous_sizes; minvalue=0.0, maxvalue=nothing)

Similar to [`Scale.size_radius`](@ref), except that the values in the `size` aesthetic are
scaled to area rather than radius, before mapping to x-axis units or `Measure` units.
"""
function size_area(f::Function=Gadfly.current_theme().continuous_sizemap; minvalue=0.0, maxvalue=nothing)
    (isa(minvalue, Measure) || isa(maxvalue, Measure)) &&
        throw(ArgumentError("Scale.size_area maps the size variable to absolute size via the function `f`. See `?Scale.size_radius` for more info."))
     return ContinuousSizeScale(f, area_transform, minvalue, maxvalue, nothing)
end

function apply_scale(scale::ContinuousSizeScale, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)

    sdata = reduce(vcat, [d.size for d in datas if d.size≠nothing])

    showvals = length(sdata)>1
    dmax = maximum(skipmissing(sdata))

    strict_span = false
    (smin, smax) =
        if scale.maxvalue===nothing
            promote(scale.minvalue, dmax)
        else
            strict_span = true
            promote(scale.minvalue, scale.maxvalue)
        end
    ticks = Gadfly.optimize_ticks(smin, smax, strict_span=strict_span)[1]
    Δ = scale.trans.f(ticks[end])-scale.trans.f(ticks[1])
    labels = scale.trans.label(ticks)

# Transform ticks e.g. to areas/porportions/sizes
    keyvals = if scale.maxvalue===nothing
            showvals = false
            scale.trans.f.(ticks.-smin)
        else
            p = (scale.trans.f.(ticks) .- scale.trans.f(ticks[1]))./Δ
            scale.f.(p)
        end

    labeldict = Dict(k=>v for (k,v) in zip(keyvals, labels))
    key_vals= OrderedDict(s=>i for (i,s) in enumerate(keyvals))
    labeler(xs) = [labeldict[x] for x in xs]

    for (aes, data) in zip(aess, datas)
            data.size===nothing && continue

            ds = if scale.maxvalue === nothing
                    scale.trans.f.(data.size.-smin)
                else
                    p = (scale.trans.f.(data.size).-scale.trans.f(ticks[1]))./Δ
                    scale.f.(p)
                end
        aes.size = ds
        showvals && (aes.size_key_vals = key_vals)
        aes.size_label = labeler
    end
end


text:
using SymSemiseparableMatrices
using Test
using LinearAlgebra
import SymSemiseparableMatrices: spline_kernel, spline_kernel_matrix

W = ones(3,100)
U = ones(3,100)
V = ones(3,100)
X = ones(100,2)
for (u,v,x) in zip(eachcol(U),eachcol(V),eachrow(X))
end

for (u,v) in zip(reverse(eachcol(U)),eachcol(V))
    println(dot(u,v))
end

d = rand(10)
for ds in eachrow(d)
    println(typeof(ds[1]))
end

A = rand(5,5)

for j=1:size(A,2),i = 1:size(A,1)
    println((i,j))
end


text:
using Documenter
using EDF

makedocs(modules=[EDF],
         sitename="EDF.jl",
         authors="Beacon Biosignals, Inc.",
         pages=["API" => "index.md"])

deploydocs(repo="github.com/beacon-biosignals/EDF.jl.git")


text:
#####
##### Gradient AD implementation using ForwardDiff
#####

import .ForwardDiff

struct ForwardDiffLogDensity{L, C} <: ADGradientWrapper
    ℓ::L
    gradientconfig::C
end

function Base.show(io::IO, ℓ::ForwardDiffLogDensity)
    print(io, "ForwardDiff AD wrapper for ", ℓ.ℓ,
          ", w/ chunk size ", length(ℓ.gradientconfig.seeds))
end

_default_chunk(ℓ) = ForwardDiff.Chunk(dimension(ℓ))

# defined to make the tag match
_logdensity_closure(ℓ) = x -> logdensity(ℓ, x)

function _default_gradientconfig(ℓ, chunk::ForwardDiff.Chunk)
    ForwardDiff.GradientConfig(_logdensity_closure(ℓ), zeros(dimension(ℓ)), chunk)
end

function _default_gradientconfig(ℓ, chunk::Integer)
    _default_gradientconfig(ℓ, ForwardDiff.Chunk(chunk))
end

"""
    ADgradient(:ForwardDiff, ℓ; chunk, gradientconfig)
    ADgradient(Val(:ForwardDiff), ℓ; chunk, gradientconfig)

Wrap a log density that supports evaluation of `Value` to handle `ValueGradient`, using
`ForwardDiff`.

Keywords are passed on to `ForwardDiff.GradientConfig` to customize the setup. In
particular, chunk size can be set with a `chunk` keyword argument (accepting an integer or a
`ForwardDiff.Chunk`).
"""
function ADgradient(::Val{:ForwardDiff}, ℓ;
                    chunk = _default_chunk(ℓ),
                    gradientconfig = _default_gradientconfig(ℓ, chunk))
    ForwardDiffLogDensity(ℓ, gradientconfig)
end

function logdensity_and_gradient(fℓ::ForwardDiffLogDensity, x::AbstractVector)
    @unpack ℓ, gradientconfig = fℓ
    buffer = _diffresults_buffer(ℓ, x)
    result = ForwardDiff.gradient!(buffer, _logdensity_closure(ℓ), x, gradientconfig)
    _diffresults_extract(result)
end


text:
# Helmholtz operators

"""
    helmholtz!(v,w,α)

Evaluate the discrete Helmholtz operator (iα - L) on `w` and return it as `v`.
The data `w` can be of type dual/primary nodes or edges; `v` must be of the same type.
However, both have to be of complex data type.
"""
function helmholtz!(out::Nodes{Dual,NX, NY,T}, w::Nodes{Dual,NX, NY,T}, α::Number) where {NX, NY, T<:ComplexF64}
    @inbounds for y in 2:NY-1, x in 2:NX-1
        out[x,y] = im*α*w[x,y] - (w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1])
    end
    out
end

function helmholtz!(out::Nodes{Primal,NX, NY,T}, w::Nodes{Primal,NX, NY,T}, α::Number) where {NX, NY, T<:ComplexF64}
    @inbounds for y in 2:NY-2, x in 2:NX-2
        out[x,y] = im*α*w[x,y] - (w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1])
    end
    out
end

function helmholtz!(out::Edges{Dual,NX, NY,T}, w::Edges{Dual,NX, NY,T}, α::Number) where {NX, NY, T<:ComplexF64}
  @inbounds for y in 2:NY-1, x in 2:NX-2
      out.u[x,y] = im*α*w.u[x,y] - (w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1])
  end
  @inbounds for y in 2:NY-2, x in 2:NX-1
      out.v[x,y] = im*α*w.v[x,y] - (w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1])
  end
  out
end

function helmholtz!(out::Edges{Primal,NX, NY, T}, w::Edges{Primal,NX, NY, T}, α::Number) where {NX, NY, T<:ComplexF64}
  @inbounds for y in 2:NY-2, x in 2:NX-1
      out.u[x,y] = im*α*w.u[x,y] - (w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1])
  end
  @inbounds for y in 2:NY-1, x in 2:NX-2
      out.v[x,y] = im*α*w.v[x,y] - (w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1])
  end
  out
end

"""
    helmholtz(w,α)

Evaluate the discrete Helmholtz operator (iα - L) on `w`. The data `w` can be of complex type dual/primary
nodes or edges. The returned result is of the same type as `w`.

# Example
"""
function helmholtz(w::Nodes{C,NX,NY,T},α::Number) where {C<:CellType,NX,NY,T<:ComplexF64}
  helmholtz!(Nodes(C,w), w, α)
end

function helmholtz(w::Edges{C,NX,NY,T},α::Number) where {C<:CellType,NX,NY,T<:ComplexF64}
  helmholtz!(Edges(C,w), w, α)
end

"""
    plan_helmholtz(dims::Tuple,α::Number,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [factor=1.0],[dx=1.0])

Constructor to set up an operator for evaluating the discrete Helmholtz operator on
complex dual or primal nodal data of dimension `dims`. If the optional keyword
`with_inverse` is set to `true`, then it also sets up the inverse Helmholtz operator
(the lattice Green's function, LGF). These can then be applied, respectively, with
`*` and `\\` operations on data of the appropriate size. The optional
parameter `factor` is a scalar used to multiply the result of the operator and
divide the inverse. The optional parameter
`dx` is used in adjusting the uniform value of the LGF to match the behavior
of the continuous analog at large distances; this is set to 1.0 by default.

Instead of the first argument, one can also supply `w::Nodes` to specify the
size of the domain.

# Example

"""
function plan_helmholtz end

"""
    plan_helmholtz!(dims::Tuple,α::Number,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [factor=1.0],[dx=1.0])

Same as [`plan_helmholtz`](@ref), but forms an operator that works in-place on the
data it operates on.
"""
function plan_helmholtz! end

struct Helmholtz{NX, NY, R, DX, inplace}
    α::Number
    factor :: ComplexF64
    conv::Union{CircularConvolution{NX, NY},Nothing}
end

for (lf,inplace) in ((:plan_helmholtz,false),
                     (:plan_helmholtz!,true))
    @eval function $lf(dims::Tuple{Int,Int},α::Number;
                   with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor::Number = 1.0, dx = 1.0, nthreads = length(Sys.cpu_info()))
        NX, NY = dims
        if !with_inverse
            return Helmholtz{NX, NY, false, dx, $inplace}(α,convert(ComplexF64,factor),nothing)
        end
        lgfh_table = load_lgf_helmholtz(NX+1,α)
        G = view(lgfh_table, 1:NX, 1:NY)
        Helmholtz{NX, NY, true, dx, $inplace}(α,convert(ComplexF64,factor),CircularConvolution(G, fftw_flags,dtype=ComplexF64))
    end

    @eval function $lf(nx::Int, ny::Int,α::Number;
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor::Number = 1.0, dx = 1.0, nthreads = length(Sys.cpu_info()))
        $lf((nx, ny), α, with_inverse = with_inverse, fftw_flags = fftw_flags, factor = factor, dx = dx, nthreads = nthreads)
    end

    @eval function $lf(::GridData{NX,NY,T},α::Number;
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor::Number = 1.0, dx = 1.0, nthreads = length(Sys.cpu_info())) where {NX,NY,T<:ComplexF64}
        $lf((NX,NY), α, with_inverse = with_inverse, fftw_flags = fftw_flags, factor = factor, dx = dx, nthreads = nthreads)
    end
end

function Base.show(io::IO, L::Helmholtz{NX, NY, R, DX, inplace}) where {NX, NY, R, DX, inplace}
    nodedims = "(nx = $NX, ny = $NY)"
    inverse = R ? " (and inverse)" : ""
    isinplace = inplace ? " in-place" : ""
    print(io, "Discrete$isinplace Helmholtz operator$inverse on a $nodedims grid with factor $(L.factor) and spacing $DX")
end

mul!(out::Nodes{C,NX,NY,T}, L::Helmholtz, s::Nodes{C,NX,NY,T}) where {C<:CellType,NX,NY,T<:ComplexF64} = (helmholtz!(out, s, L.α); output .*= L.factor)
*(L::Helmholtz{MX,MY,R,DX,false}, s::Nodes{C,NX,NY,T}) where {MX,MY,R,DX,C <: CellType,NX,NY,T<:ComplexF64} =
      L.factor*helmholtz(s,L.α)
function (*)(L::Helmholtz{MX,MY,R,DX,true}, s::Nodes{C,NX,NY,T}) where {MX,MY,R,DX,C <: CellType,NX,NY,T<:ComplexF64}
    helmholtz!(s,deepcopy(s),L.α)
end

function ldiv!(out::Nodes{C,NX, NY,T},
                   L::Helmholtz{MX, MY, true, DX, inplace},
                   s::Nodes{C, NX, NY,T}) where {C <: CellType, NX, NY, T<:ComplexF64, MX, MY, DX, inplace}

    mul!(out.data, L.conv, s.data)
    inv_factor = 1.0/L.factor

    # Adjust the behavior at large distance to match continuous kernel
    #out.data .-= (sum(s.data)/2π)*(GAMMA+log(8)/2-log(DX))
    out.data .*= inv_factor
    out
end

\(L::Helmholtz{MX,MY,R,DX,false},s::Nodes{C,NX,NY,T}) where {MX,MY,R,DX,C <: CellType,NX,NY,T<:ComplexF64} =
  ldiv!(Nodes(C,s), L, s)

\(L::Helmholtz{MX,MY,R,DX,true},s::Nodes{C,NX,NY,T}) where {MX,MY,R,DX,C <: CellType,NX,NY,T<:ComplexF64} =
  ldiv!(s, L, deepcopy(s))

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using OMEinsum
using StatsBase, Random

function test()
    for i=1:50
        ranka = rand(1:8)
        rankb = rand(1:8)
        ta = [1:ranka...]
        rankab = rand(1:min(ranka, rankb))
        tb = sample(ta, rankab; replace=false)
        tout = setdiff(ta, tb)
        for k=1:rankb-rankab
            push!(tb, ranka+k)
            push!(tout, ranka+k)
        end
        shuffle!(tb)
        shuffle!(tout)
        A = randn(fill(2, ranka)...)
        B = randn(fill(2, rankb)...)
        OMEinsum.batched_contract(Val((ta...,)), A, Val((tb...,)), B, Val((tout...,)))
    end
end

@time test()

function batched_contract2(iAs, A::AbstractArray, iBs, B::AbstractArray, iOuts)
    pA, iAps, iAbs, iAss, pB, iBps, iBbs, iBss, pOut = OMEinsum.analyse_batched_dim(iAs, iBs, iOuts)
    sAb, sAs, sAp, sBs, sBb, sBp, sAB = OMEinsum.analyse_batched_size(iAs, iAps, iAbs, iAss, size(A), iBs, iBps, iBbs, iBss, size(B))

    A, B = OMEinsum.align_eltypes(A, B)
    Apr = reshape(OMEinsum.tensorpermute(A, pA), sAb, sAs, sAp)
    Bpr = reshape(OMEinsum.tensorpermute(B, pB), sBs, sBb, sBp)
    AB = OMEinsum._batched_gemm('N','N', Apr, Bpr)
    AB = OMEinsum.tensorpermute(reshape(AB, sAB...), (pOut...,))
end

function test2()
    for i=1:10
        ranka = rand(1:8)
        rankb = rand(1:8)
        ta = [1:ranka...]
        rankab = rand(1:min(ranka, rankb))
        tb = sample(ta, rankab; replace=false)
        tout = setdiff(ta, tb)
        for k=1:rankb-rankab
            push!(tb, ranka+k)
            push!(tout, ranka+k)
        end
        shuffle!(tb)
        shuffle!(tout)
        A = randn(fill(2, ranka)...)
        B = randn(fill(2, rankb)...)
        batched_contract2(ta, A, tb, B, tout)
    end
end

@time test2()


text:
using Pkg, Random, Distributions, StatsBase, Parameters
using Pkg
ENV["PYTHON"] = "python3.7"
Pkg.build("PyCall")

Random.seed!(123)

# structure representing a node in an isolation tree
@with_kw mutable struct It_node
	l::Union{It_node, Nothing}
	r::Union{It_node, Nothing}
	split_attr::Union{Int, Nothing}
	split_val::Union{Float64, Nothing}
	level::Int
	mass::Int = 0
end

# !!!TESTED!!!
function get_random_itree(data_sub)

    # random_itree: define auxiliary function to implement recursion.
    function random_itree(x_in; current_height, lim)
		if current_height >= lim || size(x_in, 1) <= 1 # Base case check
        	return It_node(l=nothing, r=nothing, split_attr=nothing, split_val=nothing, level=current_height)
        else
        # Randomly select an attribute q.
        q = rand(1:size(x_in, 2))
        # Randomly select a split point p between min and max values of attribute q in X.
		min_q_x_in = minimum(x_in[:, q])
		max_q_x_in = maximum(x_in[:, q])
		if min_q_x_in == max_q_x_in
			p = min_q_x_in
		else
        	p = rand(Uniform(minimum(x_in[:, q]), maximum(x_in[:, q])))
		end
        # Get left and right subtrees.
        xl = x_in[x_in[:, q] .< p, :]
        xr = x_in[x_in[:, q] .>= p, :]
        # Recursive case
        return It_node(l=random_itree(xl, current_height=current_height+1, lim=lim),
                   r=random_itree(xr, current_height=current_height+1, lim=lim),
                   split_attr=q,
                   split_val=p,
                   level=current_height)
        end
    end

    # Build itree
    return random_itree(data_sub, current_height=0, lim=10)
end

# !!!TESTED!!!
function get_n_random_itrees(n, subs_size, data)
    random_itrees = Array{It_node, 1}(undef, n)
    # TODO: parallelize!
    for k = 1:n
        # Get a random sample of training examples to build next random itree.
        data_sub = data[sample(1:size(data, 1), subs_size, replace=false), :]
        random_itrees[k] = get_random_itree(data_sub)  # Get next random itree
    end
    return random_itrees, subs_size
end

# !!!TESTED!!!
function get_lowest_common_node_mass(itree, x1, x2)

    # If node is a leaf, return its mass.
    if itree.split_val == nothing
    	return itree.mass
    end

    # If x1 and x2 are in different subtrees, return current node's mass.
	if (x1[itree.split_attr] < itree.split_val) != (x2[itree.split_attr] < itree.split_val)
    	return itree.mass
    end

    # If both examples if left subtree, make recursive call.
	if (x1[itree.split_attr] < itree.split_val) && (x2[itree.split_attr] < itree.split_val)
    	return get_lowest_common_node_mass(itree.l, x1, x2)
    end

    # If both examples in right subtree, make recursive call.
	if (x1[itree.split_attr] >= itree.split_val) && (x2[itree.split_attr] >= itree.split_val)
    	return get_lowest_common_node_mass(itree.r, x1, x2)
    end
end

# !!!TESTED!!!
function mass_based_dissimilarity(x1, x2, itrees, subs_size)
    # In each i-tree, find lowest nodes containing both examples and accumulate masses.
    sum_masses = 0
    for i = 1:length(itrees)
    	sum_masses += get_lowest_common_node_mass(itrees[i], x1, x2)/subs_size
    end

    return (1/length(itrees)) * sum_masses  # Divide by number of space partitioning models.

end

# !!!TESTED!!!
function get_node_masses(itrees, data)

    # traverse: traverse itree with example and increment masses of visited nodes
    function traverse(example, it_node)

        # base case - in leaf
        if it_node.l == nothing && it_node.r == nothing
        	it_node.mass += 1

        # if split attribute value lower than split value
        elseif example[it_node.split_attr] < it_node.split_val
        	it_node.mass += 1
        	traverse(example, it_node.l)  # Traverse left subtree.

        # if split attribute value greater or equal to split value
        else
        	it_node.mass += 1
        	traverse(example, it_node.r)  # Traverse right subtree.
        end
    end

    # compute_masses: compute masses of nodes in itree
    function compute_masses(itree, data)
		for example_idx in 1:size(data, 1)
			traverse(data[example_idx, :], itree)
        end
    end

    # TODO: parallelize!
    for itree in itrees  # Go over itrees and set masses of nodes.
        compute_masses(itree, data)
    end
end

function get_dist_func(num_itrees, data)
    itrees, subs_size = get_n_random_itrees(num_itrees, size(data, 1), data)
    get_node_masses(itrees, data)
	res_func = (i1, i2) -> mass_based_dissimilarity(data[i1.+1, :], data[i2.+1, :], itrees, subs_size)
	return (m, i1, i2, kwargs...) -> res_func.(i1, i2)
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
#===============================================================================

	Code for finding the terminal state of the HJB for a Hugget Model

===============================================================================#

using LinearAlgebra, SparseArrays, Plots

function terminal()
	# Define model parameters

	global s = 2
	global ρ = 0.05
	z1 = .1
	z2 = .2
	global z = [z1 z2]
	la1 = 0.6
	la2 = 0.8
	λ = [la1 la2]

	r0 = 0.03
	rmin = 0.01
	rmax = 0.04

	H= 1000
	global amin = -0.15
	global amax = 4
	a = LinRange(amin,amax,H)'
	global a = convert(Array, a)'
	global da = (amax-amin)/(H-1);

	global aa = [a a]
	global zz = ones(H,1)*z

	maxit= 100
	crit = 1e-6
	Delta = 100

	dVf, dVb, c = [zeros(H,2) for i in 1:3]

	global Aswitch = [-sparse(I, H, H)*λ[1] sparse(I, H, H)*λ[1];sparse(I, H, H)*λ[2] -sparse(I, H, H)*λ[2]]

	Ir = 40
	crit_S = 1e-8

	# INITIAL GUESS
	r = r0

	v0=zeros(H,2)
	v0[:,1] = (z[1] .+ r.*a).^(1-s)/(1-s)/ρ
	v0[:,2] = (z[2] .+ r.*a).^(1-s)/(1-s)/ρ
	dist =[]
	r_r=[]
	rmin_r=[]
	rmax_r=[]
	V_n =[]
	g_r=[]
	adot=[]
	V_r=[]
	S=[]

	for ir=1:Ir
		push!(r_r,r)
		push!(rmin_r, rmin)
		push!(rmax_r, rmax)

		global v = v0

		for n=1:maxit
		    V = v
		    push!(V_n,V)
		    #forward difference
		    dVf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
		    dVf[H,:] = (z .+ r.*amax).^(-s) #will never be used, but impose state constraint a<=amax just in case
		    #backward difference
		    dVb[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
		    dVb[1,:] = (z .+ r.*amin).^(-s) #state constraint boundary condition

		    I_concave = dVb .> dVf #indicator whether value function is concave (problems arise if this is not the case)

		    # consumption and savings with forward difference
		    cf = dVf.^(-1/s)
		    ssf = zz + r.*aa - cf
		    # consumption and savings with backward difference
		    cb = dVb.^(-1/s)
		    ssb = zz + r.*aa - cb
		    # consumption and derivative of value function at steady state
		    c0 = zz + r.*aa
		    dV0 = c0.^(-s)

		    #= dV_upwind makes a choice of forward or backward differences based on
		     the sign of the drift    =#
		    If = ssf .> 0 #positive drift → forward difference
		    Ib = ssb .< 0 #negative drift → backward difference
		    I0 = (1.0 .-If-Ib) #at steady state

		    dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0  #important to include third term
		    c = dV_Upwind.^(-1/s)
		    u = c.^(1-s)/(1-s)

		    # CONSTRUCT MATRIX
		    X = - min.(ssb,0)/da
		    Y = - max.(ssf,0)/da + min.(ssb,0)/da
		    Z = max.(ssf,0)/da

		    A1=spdiagm(0 => Y[:,1], -1 => X[2:H,1], 1 =>Z[1:H-1,1])
		    A2=spdiagm(0=> Y[:,2], -1 =>X[2:H,2], 1 => Z[1:H-1,2])

			global A = [A1 spzeros(H,H);spzeros(H,H) A2] + Aswitch

		    B = (1/Delta + ρ)*sparse(I, 2*H, 2*H) - A

		    u_stacked = [u[:,1];u[:,2]]
		    V_stacked = [V[:,1];V[:,2]]

		    b = u_stacked + V_stacked/Delta
		    V_stacked = B\b #SOLVE SYSTEM OF EQUATIONS

		    V = [V_stacked[1:H] V_stacked[H+1:2*H]]

		    V_change = V - v
		    v = V

		    push!(dist,findmax(abs.(V_change))[1])
		    if dist[n]<crit
		        println("Value Function Converged, Iteration = ")
		        println(n)
		        break
		    end
		end

	#===============================================================================

		Code for finding the terminal state of the KFE

	===============================================================================#
		AT = A'
		b = zeros(2*H,1)

		# need to fix one value, otherwise matrix is singular
		i_fix = 1
		b[i_fix]=.1
		row = [zeros(1,i_fix-1) 1 zeros(1,2*H-i_fix)]
		AT[i_fix,:] = row

		# Solve linear system
		gg = AT\b
		g_sum = gg'*ones(2*H,1)*da
		gg = gg./g_sum;

		global g = [gg[1:H] gg[H+1:2*H]]

		check1 = g[:,1]'*ones(H,1)*da
		check2 = g[:,2]'*ones(H,1)*da

		push!(g_r, g)
		push!(adot,zz + r.*aa - c)
		push!(V_r,V)

		push!(S,g[:,1]'*a*da + g[:,2]'*a*da)

			# UPDATE INTEREST RATE
			if S[ir][1] > crit_S
			    println("Excess Supply")
			    rmax = r
			    r = 0.5*(r+rmin)
			elseif S[ir][1]< -crit_S
			    println("Excess Demand")
			    rmin = r
			    r = 0.5*(r+rmax)
			elseif abs(S[ir][1])<crit_S
			    println("Equilibrium Found, Interest rate =")
			    println(r)
			    break
			end
	end

return g, v, r

end


text:
using LabelledArrays

export logistic3

"""
    eom_logistic3(u, p, t)

Equations of motion for a system consisting of three coupled logistic map
representing the response of two independent dynamical variables to the
forcing from a common driver. The dynamical influence goes in the directions
Z → X and Z → Y.

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x(t+1) = (x(t)(r - r_1 x(t) - z(t) + σ_x η_x)) \\mod 1 \\\\
y(t+1) = (y(t)(r - r_2 y(t) - z(t) + σ_y η_y)) \\mod 1 \\\\
z(t+1) = (z(t)(r - r_3 z(t) + σ_z η_z)) \\mod 1
\\end{aligned}
```

Dynamical noise may be added to each of the dynamical variables by tuning the
parameters `σz`, `σx` and `σz`. Default values for the parameters
`r₁`, `r₂` and `r₃` are set such that the system exhibits chaotic behaviour,
with `r₁ = r₂ = r₃ = 4`.

## References

1. Runge, Jakob. Causal network reconstruction from time series: From theoretical
    assumptions to practical estimation, Chaos 28, 075310 (2018);
    doi: 10.1063/1.5025050
"""
function eom_logistic3(u, p, t)
    r₁, r₂, r₃, σx, σy, σz = (p...,)
    x, y, z = (u...,)

    # Independent dynamical noise for each variable.
    ηx = rand()
    ηy = rand()
    ηz = rand()

    dx = (x*(r₁ - r₁*x - z + σx*ηx)) % 1
    dy = (y*(r₂ - r₂*y - z + σy*ηy)) % 1
    dz = (z*(r₃ - r₃*z + σz*ηz)) % 1
    return SVector{3}(dx, dy, dz)
end

function logistic3(u₀, r₁, r₂, r₃, σx, σy, σz)
    p = @LArray [r₁, r₂, r₃, σx, σy, σz] (:r₁, :r₂, :r₃, :σx, :σy, :σz)
    DiscreteDynamicalSystem(eom_logistic3, u₀, p)
end

"""
    logistic3(;u₀ = rand(3), r = 4,
        σx = 0.05, σy = 0.05, σz = 0.05) → DiscreteDynamicalSystem

Initialise a dynamical system consisting of three coupled logistic map
representing the response of two independent dynamical variables to the
forcing from a common driver. The dynamical influence goes in the directions
``Z \\to X`` and ``Z \\to Y``.

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x(t+1) = (x(t)(r - r_1 x(t) - z(t) + σ_x η_x)) \\mod 1 \\\\
y(t+1) = (y(t)(r - r_2 y(t) - z(t) + σ_y η_y)) \\mod 1 \\\\
z(t+1) = (z(t)(r - r_3 z(t) + σ_z η_z)) \\mod 1
\\end{aligned}
```

Dynamical noise may be added to each of the dynamical variables by tuning the
parameters `σz`, `σx` and `σz`. Default values for the parameters
`r₁`, `r₂` and `r₃` are set such that the system exhibits chaotic behaviour,
with `r₁ = r₂ = r₃ = 4`.

## References

1. Runge, Jakob. Causal network reconstruction from time series: From theoretical
    assumptions to practical estimation, Chaos 28, 075310 (2018);
    doi: 10.1063/1.5025050
"""
logistic3(;u₀ = rand(3), r₁ = 4, r₂ = 4, r₃ = 4,
    σx = 0.05, σy = 0.05, σz = 0.05) = logistic3(u₀, r₁, r₂, r₃, σx, σy, σz)


text:
println("\n\n\nStarting runtests.jl $(join(ARGS, " ")) ...")
# addprocs(3)
using FactCheck
importall FunctionalData
FactCheck.setstyle(:compact)

shouldtest(f, a) = length(ARGS) == 0 || a == ARGS[1] ? facts(f, a) : nothing
shouldtestcontext(f, a) = length(ARGS) < 2 || a == ARGS[2] ? facts(f, a) : nothing

type A
    a
    b
end

shouldtest("exports") do
    # Check whether all exported symbols really exist
    for a in names(FunctionalData)
        # @fact eval(a) --> eval("FunctionalData.$a")
    end
end

if VERSION < v"0.5-"
    readstring = readall
end

function checkcodeexamples(filename)
    absname = @p functionloc shouldtest | fst | dirname | joinpath "../doc/"*filename
    mdlines = @p readstring absname | lines

    startasserts = true
    insidecode = false

    for line in mdlines
        isempty(strip(line)) && continue
        if startswith(line, "```jl")
            insidecode = true && startasserts
        elseif startswith(line, "```")
            insidecode = false
        elseif insidecode
            # println(line)
            if contains(line, "=>")
                line = @p split line "=>" | concat "@fact (" fst(_) ")  -->  (" snd(_) ")"
            end
            eval(parse(line))
        end
    end
end

shouldtest("doc") do
    shouldtestcontext("computing") do
        checkcodeexamples("computing.md")
    end
    shouldtestcontext("dataflow") do
        checkcodeexamples("dataflow.md")
    end
    shouldtestcontext("helpers") do
        checkcodeexamples("helpers.md")
    end
    shouldtestcontext("io") do
        checkcodeexamples("io.md")
    end
    shouldtestcontext("lensize") do
        checkcodeexamples("lensize.md")
    end
    shouldtestcontext("output") do
        checkcodeexamples("output.md")
    end
    shouldtestcontext("pipeline") do
        checkcodeexamples("pipeline.md")
    end
end

shouldtest("views") do
    a = [1,2,3]
    @fact FD.view(a,1)  -->  [1]
    @fact FD.view(a,3)  -->  [3]
    a = [1 2 3]
    @fact FD.view(a,1)  -->  row([1])
    @fact FD.view(a,3)  -->  row([3])
    a = [1 2 3; 4 5 6]
    @fact FD.view(a,1)  -->  col([1,4])
    @fact FD.view(a,3)  -->  col([3,6])
    v = FD.view(a,2)
    v[2] = 10
    @fact a --> [1 2 3; 4 10 6]
    a = UInt8[1 2 3; 4 5 6]
    @fact FD.view(a,1)  -->  col(UInt8[1,4])
    @fact FD.view(a,3)  -->  col(UInt8[3,6])
    @fact FD.view(a,2:3)  --> part(a,2:3)
end

shouldtest("lensize") do
    @fact siz(1)      -->  transpose([1 1])
    @fact siz([1])    -->  transpose([1 1])
    @fact siz([1,2])  -->  transpose([2 1])
    @fact siz([1;2])  -->  transpose([2 1])
    @fact siz([1 2])  -->  transpose([1 2])
    @fact siz(transpose([1 2]))  -->  transpose([2 1])

    @fact siz3(rand(1))  --> [1 1 1]'
    @fact siz3(rand(1,2))  --> [1 2 1]'
    @fact siz3(rand(1,2,3))  --> [1 2 3]'

    @fact len(1)        -->  1
    @fact len([1])      -->  1
    @fact len([1,2])    -->  2
    @fact len([1;2])    -->  2
    @fact len([1 2])    -->  2
    @fact len([1,2,3])  -->  3
    @fact len([1 1 1 ;2 2 3])  -->  3
    @fact len("adsf")   -->   4
    @fact len('a')      -->  1
    @fact len(['a',1])  -->  2

end

shouldtest("basics") do
    shouldtestcontext("arraylike") do
        @fact size(FunctionalData.arraylike([1],2)) --> (1,2)
        @fact size(FunctionalData.arraylike([1 2],2)) --> (1,2,2)
        @fact size(FunctionalData.arraylike(1,2)) --> (2,)
    end
    shouldtestcontext("ones") do
        @fact onessiz([2 3 4]') --> ones(2,3,4)
        @fact zerossiz([2 3 4]') --> zeros(2,3,4)
        @fact size(randsiz([2 3 4]')) --> (2,3,4)
        @fact size(randnsiz([2 3 4]')) --> (2,3,4)
    end
    shouldtestcontext("shones") do
        @fact shonessiz([2 3 4]') --> ones(2,3,4)
        @fact shzerossiz([2 3 4]') --> zeros(2,3,4)
        @fact size(randsiz([2 3 4]')) --> (2,3,4)
        @fact size(randnsiz([2 3 4]')) --> (2,3,4)
    end
    shouldtestcontext("repeat") do
        @fact repeat('a',0) --> ""
        @fact repeat('a',1) --> "a"
        @fact repeat('a',3) --> "aaa"
        @fact repeat("a",0) --> ""
        @fact repeat("a",1) --> "a"
        @fact repeat("a",3) --> "aaa"
        @fact repeat(1,3)  --> [1,1,1]
        @fact repeat([1],3) --> [1,1,1]
        @fact repeat([1 2]',3) --> [1 1 1; 2 2 2]
    end
if VERSION >= v"0.4-"
    shouldtestcontext("minimum") do
        @fact minimum(Float32)+maximum(Float32) --> 0f0
        @fact minimum(Float64)+maximum(Float64) --> 0.0
        @fact minimum(Int)+maximum(Int) --> -1
        @fact maximum() --> maximum(Float64)
        @fact minimum() --> minimum(Float64)
    end
end

end

type _somedummytype
  a
  b
end

shouldtest("accessors") do
    shouldtestcontext("at") do
        @fact at([1,2,3],1) --> 1
        @fact at([1 2 3],1) --> col([1])
        @fact at([1;2;3],1) --> 1
        @fact at((1,2,3),1) --> 1
        @fact at([1 2 3; 4 5 6],1) --> col([1 4])
        @fact at(cat(3,[1 1],[2 2],[3 3]),1) --> [1 1]

        @fact at("asdf",1) --> 'a'
        @fact at(Any["aa",1],1) --> "aa"
        @fact at(['a','b'],2) --> 'b'
        @fact at(['a','b'],1) --> 'a'

        @fact size(at(rand(2,3,4),(1,))) --> ()
        @fact size(at(rand(2,3,4),(1:2,))) --> (2,)
        @fact size(at(rand(2,3,4),([1,2],1:2))) --> (2,2)

        d = Dict(:a => 1, :b => Dict(:c => 2, :d => Dict(:e => 3)))
        @fact at(d,:a) --> 1
        @fact at(d,:b,:c) --> 2
        @fact at(d,:b,:c) --> 2
        @fact at(d,:b,:d,:e) --> 3
    end
    shouldtestcontext("atend") do
        @fact atend(1:10,1) --> 10
        @fact atend(1:10,2) --> 9
    end

    shouldtestcontext("setat") do
        a = [1,2,3]
        setat!(a, 1, 10)
        @fact at(a,1) --> 10
        a = [1 2 3]
        setat!(a, 1, 10)
        @fact at(a,1) --> col([10])
        setat!(a, 1, col([10]))
        @fact at(a,1) --> col([10])
        a = [1 2 3; 4 5 6]
        setat!(a, 1, col([10,11]))
        @fact at(a,1) --> col([10,11])
        a = cat(3,[1 1],[2 2],[3 3])
        setat!(a,1,[10 11])
        @fact at(a,1) --> [10 11]

        a = Any["aa",1]
        setat!(a, 1, "bb")
        @fact at(a,1) --> "bb"
    end

    # shouldtestcontext("at!") do
    #     @fact at!([1,2,3],1)   1
    #     @fact at!([1 2 3],1)    slicedim([1 2 3],2,1)
    #     @fact at!([1;2;3],1)    1
    #     @fact at!([1 2 3; 4 5 6],1)    slicedim([1 2 3; 4 5 6],2,1)
    #     @fact at!("asdf",1)    'a'
    #     @fact at!(['a','b'],1)    'a'
    #     @fact at!([1,[2 3],'a'],2)    [2 3]

    #     D = [1 2 3; 4 5 6]
    #     b = at!(D,2)
    #     b[2] = 1000
    #     @fact D    [1 2 3; 4 100 6]
    # end

    shouldtestcontext("part") do
        @fact part([1,2,3],[1])  --> [1]
        @fact part([1,2,3],1:2) --> [1,2]
        @fact part([1,2,3],[1,2]) --> [1,2]
        @fact part([1 2 3],[1,3]) --> [1 3]
        @fact part([1 2 3],[1,3]) --> slicedim([1 2 3],2,[1,3])
        @fact part([1;2;3],[1,3]) --> [1;3]
        @fact part([1 2 3; 4 5 6],[1,3]) --> [1 3;4 6]
        @fact part([1 2 3; 4 5 6],[1 2; 3 2]) --> [3,5]
    end

    shouldtestcontext("dict") do
        d = Dict(:a => 1, :b => 2)
        @fact part(d, :a) --> Dict(:a => 1)
        @fact part(d, :a, :b) --> d
        @fact values(d, :a) --> [1]
        @fact values(d, :a, :b) --> [1,2]
        @fact values(A(1,2), :a, :b)  -->  [1,2]
    end

    shouldtestcontext("trimmedpart") do
        @fact trimmedpart(collect(1:10), -1:3) --> [1,2,3]
        @fact trimmedpart(collect(1:10), 1:3) --> [1,2,3]
        @fact trimmedpart(collect(1:10), 8:13) --> [8,9,10]
        @fact trimmedpart(collect(1:10), 13:15) --> []
        @fact trimmedpart(1:10, [1,3,30,2,-10]) --> [1,3,2]
    end

    shouldtestcontext("fst") do
        @fact fst([1 2 3]) --> col([1])
        @fact fst(1:3) --> 1
        @fact fst([1 2 3; 4 5 6]) --> col([1,4])
        @fact fst('a') --> 'a'
        @fact fst("asdf") --> 'a'
    end

    shouldtestcontext("last") do
        @fact last([1 2 3]) --> col([3])
        @fact last(1:3) --> 3
        @fact last([1 2 3; 4 5 6]) --> col([3,6])
        @fact last('a') --> 'a'
        @fact last("asdf") --> 'f'
    end

    shouldtestcontext("drop") do
        @fact drop([1,2,3],1) --> [2,3]
        @fact drop(Any["test",2,"asdf"],1) --> Any[2,"asdf"]
        @fact drop(Any["test",2,"asdf"],2) --> ["asdf"]

        @fact drop(1:3,1) --> 2:3
        @fact drop([1 2 3],1) --> [2 3]
        @fact drop([1 2 3; 4 5 6],1) --> [2 3; 5 6]
        @fact drop([1 2 3; 4 5 6],2) --> col([3; 6])
    end

    shouldtestcontext("dropat") do
        @fact dropat(1:10,3:9) --> [1,2,10]
    end

    shouldtestcontext("take") do
        @fact take(1:3,1) --> 1:1
        @fact take([1 2 3],1) --> col([1])
        @fact take([1 2 3; 4 5 6],1) --> col([1; 4])
        @fact take([1 2 3; 4 5 6],2) --> [1 2; 4 5]
        @fact take("asdf",1) --> "a"
        @fact take("asdf",2) --> "as"

        @fact last(1:3,1) --> 3:3
        @fact last([1 2 3],1) --> col([3])
        @fact last([1 2 3; 4 5 6],1) --> col([3; 6])
        @fact last([1 2 3; 4 5 6],2) --> [2 3; 5 6]
        @fact last("asdf",1) --> "f"
        @fact last("asdf",2) --> "df"
    end

    shouldtestcontext("takelast") do
        @fact takelast("asdf",1) --> "f"
        @fact takelast("asdf",2) --> "df"
        @fact takelast("asdf",10) --> "asdf"
    end

    shouldtestcontext("takewhile") do
        @fact takewhile(1:10,x->x<5) --> 1:4
        @fact (@p takewhile (1:10) isless 5) --> 1:4
        @fact (@p takewhile (1:10) isless 50) --> 1:10
        @fact (@p takewhile (1:10) isless 0) --> []
    end

    shouldtestcontext("droplast") do
        @fact droplast(1:3,1) --> 1:2
        @fact droplast([1 2 3],1) --> [1 2]
        @fact droplast([1 2 3; 4 5 6],1) --> [1 2 ; 4 5 ]
        @fact droplast([1 2 3; 4 5 6],2) --> col([1; 4])
        @fact droplast([]) --> []
        @fact droplast([1]) --> []
        @fact droplast([],2) --> []
    end

    shouldtestcontext("dropwhile") do
        @fact dropwhile(1:10,x->x<5) --> 5:10
        @fact (@p dropwhile (1:10) isless 5) --> 5:10
        @fact (@p dropwhile (1:10) isless 50) --> []
        @fact (@p dropwhile (1:10) isless 0) --> 1:10
    end

    shouldtestcontext("cut") do
        a = [1,2,3,4,5]
        b = [1 2 3; 4 5 6]
        @fact cut(a,3)  -->  ([1,2,4,5],[3])
        @fact cut(b,1)  -->  ([2 3; 5 6], col([1,4]))
        @fact cut(b,2:3)  -->  (col([1,4]), [2 3; 5 6])
    end

    shouldtestcontext("partition") do
        @fact partition(1:3,1)  --> Any[1:3]
        a = partition(1:3,2)
        @fact (a == Any[1:2, 3:3] || a == Any[1:1, 2:3])  -->  true  # Julia v0.3 and v0.4 work differently
        @fact partition(1:3,3)  --> Any[1:1, 2:2, 3:3]
        @fact partition(1:3,4)  --> Any[1:1, 2:2, 3:3]
    end
    shouldtestcontext("partsoflen") do
        @fact partsoflen(1:4,2)  -->  Any[1:2, 3:4]
        @fact partsoflen(1:4,3)  -->  Any[1:3, 4:4]
    end
    shouldtestcontext("extract") do
        @fact extract(_somedummytype(1,2), :a)  -->  1
        @fact extract(_somedummytype(1,2), :b)  -->  2
        @fact extract([_somedummytype(1,2), _somedummytype(3,4)], :b)  -->  [2,4]
        d1 = Dict(:a => 1)
        d2 = Dict(:b => 2)
        @fact extract(d1, :a)  -->  1
        @fact extract(d1, :b)  -->  nothing
        @fact extract([d1,d2], :a, 10)  -->  [1, 10]
    end
    shouldtestcontext("extractnested") do
        a = [Dict(:a => 1, :b => Dict(:c => 1)), Dict(:a => 2, :b => Dict(:c => 3))]
        @fact extractnested(a,:b,:c) --> Any[1,3]
    end
    shouldtestcontext("fieldvalues") do
        @fact fieldvalues(A(1,2)) --> [1,2]
        @fact dict(A(1,2)) --> Dict(:a => 1, :b => 2)
    end
    shouldtestcontext("isnil") do
        @fact isnil(Void) --> true
        @fact isnil(nothing) --> true
        @fact isnil(Nullable()) --> true
        @fact isnil(Nullable(1)) --> false
        @fact isnil(1) --> false
        @fact isnil("asdf") --> false
    end
end

shouldtest("computing") do
    shouldtestcontext("fold") do
        @fact fold([1,2,3], max)  -->  3
        @fact fold(["1","2","3"], concat)  -->  "123"
    end
    shouldtestcontext("sort") do
        @fact FunctionalData.sort([1,2,3], id) --> [1,2,3]
        @fact FunctionalData.sort([1 2 3], id) --> [1 2 3]
        @fact FunctionalData.sort([1,2,3], x->-x) --> [3,2,1]
        @fact FunctionalData.sort([1 2 3], x->-x) --> [3 2 1]
        local D = [Dict(:id => x, :a => string(x)) for x in 1:3]
        @fact FunctionalData.sort(D, :id) --> D
        @fact FunctionalData.sort("dcba", x->convert(Int,x)) --> "abcd"
        @fact FunctionalData.sort("dcba", x->convert(Int,x); rev = true) --> "dcba"
        @fact FunctionalData.sortrev("dcba", x->convert(Int,x)) --> "dcba"
    end
    shouldtestcontext("groupdict") do
        a = [1,2,3,2,3,3]
        @fact (@p groupdict a id) --> Dict(1 => Any[1], 2 => Any[2,2], 3 => Any[3,3,3])
    end
     shouldtestcontext("groupby") do
        a = [1,2,3,2,3,3]
        @fact (@p groupby a id) --> Any[Any[1],Any[2,2],Any[3,3,3]]
        a = [2 1 3 2 3 3; 20 10 30 20 30 30]
        @fact (@p groupby a getindex 2) --> Any[[1 10]', [2 2; 20 20], [3 3 3; 30 30 30]]
    end
    shouldtestcontext("filter") do
        @fact filter([1,2,3],x->isodd(x)) --> [1,3]
        @fact filter([1,2,3],x->iseven(x)) --> [2]
        @fact (@p filter Any[1,2,3] unequal 3) --> [1,2]
        @fact (@p filter [1,2,3] unequal 3) --> [1,2]
        @fact (@p select [1,2,3] unequal 3) --> [1,2]
        @fact (@p reject [1,2,3] unequal 3) --> [3]
    end
    shouldtestcontext("uniq") do
        @fact uniq([10 20 30])  -->  [10 20 30]
        @fact uniq([20 20 10], id)  -->  [20 10]
        @fact uniq([20 -10 10], abs)  -->  [20 -10]
        @fact uniq([20 10 -10], abs)  -->  [20 10]
        @fact (@p uniq [20 10 -10] abs)  -->  [20 10]
        @fact (@p uniq [20 10 -10] getindex 1)  -->  [20 10 -10]
    end
    shouldtestcontext("map") do
        @fact map([1 2 3; 4 5 6], x->[size(x,1)]) -->   [2 2 2]
        @fact map([1 2 3; 4 5 6], x->Any[size(x,1)]) -->   Any[2 2 2]
        @fact map([1 2 3; 4 5 6], x->[size(x,1);size(x,1)]) --> [2 2 2; 2 2 2]
        @fact map([1 2 3; 4 5 6], x->[size(x,1);size(x,1)]) --> [2 2 2; 2 2 2]
        @fact map([1 2 3; 4 5 6], x->[size(x,1) size(x,1)]) --> cat(3,[2 2],[2 2],[2 2])
        @fact map((Dict(1 => 2)), x->(fst(x), 10*snd(x))) --> Dict(1 => 20)
        @fact map((Dict(1 => 2)), x->nothing) --> Dict()
        @fact (@p map Dict(1 => 2) x->(fst(x), 10*snd(x))) --> Dict(1 => 20)
        @fact (@p map Dict(1 => 2) x->nothing) --> Dict()

        @fact map(Dict(1 => 2, 3 => 4), x->(fst(x),10*fst(x)+snd(x))) --> Dict(1=>12,3=>34)

        @fact mapkeys((Dict(1 => 2)), x -> 2x) --> Dict(2 => 2)
        @fact mapvalues((Dict(1 => 2)), x -> 2x) --> Dict(1 => 4)
        d = Dict(:a => 1, :b => 2)
        @fact (@p mapmap d id | sort) --> [1,2]
    end
    shouldtestcontext("mapi") do
        @fact reduce(&, mapi(0:9, (x,i)->x+1==i)) --> true
    end
    shouldtestcontext("worki") do
        a = zeros(Int,3)
        worki(a,(x,i)->a[i]=2*i)
        @fact a  -->  [2,4,6]
    end
    shouldtestcontext("map2") do
        @fact map2(1:3, 10:12, (+))  -->  [11,13,15]
    end
    shouldtestcontext("mapmap") do
        @fact mapmap([[1 2]; [3 4]], x -> x+1)  -->  [[2 3]; [4 5]]
    end
    shouldtestcontext("map!") do
        f(x) = x*2
        f!(x) = x[:] = f(x)
        f!r(x) = f(x)
        f2!(r,x) = r[:] = vcat(f(x),f(x))
        a = [1 2 3; 4 5 6]
        @fact map(a,f) --> [2 4 6; 8 10 12]
        a = [1 2 3; 4 5 6]
        map!(a,f!)
        @fact a --> [2 4 6; 8 10 12]
        a = [1 2 3; 4 5 6]
        map!r(a,f)
        @fact a --> [2 4 6; 8 10 12]
        a = [1 2 3; 4 5 6]
        r = zeros(Int, 4, 3)
        map2!(a, r, f2!)
        @fact r --> 2*vcat(a,a)
        r = zeros(Int, 4, 3)
        r = map2!(a, x->vcat(x,x)*2, f2!)
        @fact r --> 2*vcat(a,a)
    end
    shouldtestcontext("shmap") do
        a = row(collect(1:3))
        r = shmap(a, x->x+1)
        @fact r --> a + 1
        a = rand(10,round(Int,1e3))
        r = shmap(a, x->x+1)
        @fact r --> a+1
    end
    shouldtestcontext("shmap!") do
        a = share(row(collect(1:3)))
        orig = copy(a)
        shmap!(a, x->x[:] = x+1)
        @fact a --> orig + 1
        a = share(rand(10,round(Int,1e3)))
        orig = copy(a)
        shmap!(a, x->x[:] = x+1)
        @fact a --> orig + 1
    end
    shouldtestcontext("shmap!r") do
        a = share(row(collect(1:3)))
        orig = copy(a)
        shmap!r(a, x-> x+1)
        @fact a --> orig + 1
        # a = share(rand(10,round(Int,1e3)))
        # orig = copy(a)
        # shmap!r(a, x-> x+1)
        # @fact a --> orig + 1
    end
    shouldtestcontext("shmap2!") do
        a = row(collect(1:3))
        r = shmap2!(a, x->x+1, (r,x)->r[:] = x+1)
        @fact r --> a + 1
        r = shzerossiz(siz(a))
        shmap2!(a, r, (r,x)->r[:] = x+1)
        @fact r --> a + 1
        a = rand(10,round(Int,1e3))
        r = shmap2!(a, x->x+1, (r,x)->r[:] = x+1)
        @fact r --> a + 1
    end
    shouldtestcontext("pmap") do
        a = row(collect(1:3))
        r = pmap(a, x->x+1)
        @fact r --> a + 1
        @fact eltype(r) --> Int
        r = pmapvec(a, x->x+1)
        @fact eltype(r) --> Any
        a = rand(2,10)
        r = pmap(a, x->x+1)
        @fact r --> a+1
    end
    shouldtestcontext("lmap") do
        a = row(collect(1:3))
        r = lmap(a, x->x+1)
        @fact r --> a + 1
        @fact eltype(r) --> Int
        r = lmapvec(a, x->x+1)
        @fact eltype(r) --> Any
        a = rand(2,10)
        r = lmap(a, x->x+1)
        @fact r --> a+1
    end
    shouldtestcontext("amap") do
        a = row(1:30)
        r = amap(a, x->x+1)
        @fact r --> a + 1
        a = rand(2,10)
        r = amap(a, x->x+1)
        @fact r --> a+1
        @fact amap2(1:10, 1:10, +) --> collect(2*(1:10))
        f = (VERSION < v"0.5-") ? utf8 : String
        @fact amap2("abc", 1:3, (x,y)->"$x$y") --> map(f, ["a1","b2","c3"])
        # @fact amapvec2(1:10, 1:10, +) --> unstack(2*(1:10))
    end

    shouldtestcontext("table") do
        adder(x,y) = x+y
        pass(x,y) = [x,y]
        passarray(x,y) = col([x,y])
        @fact table(id,[1,2,3])  -->  [1,2,3]
        @fact table([1,2,3],id)  -->  [1,2,3]
        @fact table(pass,[1,2],1:3)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact table([1,2],1:3,pass)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact ltable(id,[1,2,3])  -->  [1,2,3]
        @fact ltable([1,2,3],id)  -->  [1,2,3]
        @fact ltable(pass,[1,2],1:3)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact ltable([1,2],1:3,pass)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact ptable(id,[1,2,3])  -->  [1,2,3]
        @fact ptable([1,2,3],id)  -->  [1,2,3]
        @fact ptable(pass,[1,2],1:3)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact ptable([1,2],1:3,pass)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact tableany(id,[1,2,3])  -->  Any[1,2,3]
        @fact tableany([1,2,3],id)  -->  Any[1,2,3]
        @fact tableany(pass,[1,2],1:3)  -->  reshape(Any[[1,1], [2,1], [1,2], [2,2], [1,3], [2,3]], 2, 3)
        @fact tableany([1,2],1:3,pass)  -->  reshape(Any[[1,1], [2,1], [1,2], [2,2], [1,3], [2,3]], 2, 3)
        @fact table(passarray,[1,2],1:3)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact table([1,2],1:3,passarray)  -->  cat(3, [1 2; 1 1], [1 2; 2 2], [1 2; 3 3])
        @fact tableany(passarray,[1,2],1:3)  -->  reshape(map(Any[[1,1], [2,1], [1,2], [2,2], [1,3], [2,3]],col), 2, 3)
        @fact tableany([1,2],1:3,passarray)  -->  reshape(map(Any[[1,1], [2,1], [1,2], [2,2], [1,3], [2,3]],col), 2, 3)
        @fact table(adder,[1 2; 3 4],1:3)  -->  cat(3, [2 3; 4 5], [3 4; 5 6], [4 5; 6 7])
        @fact table([1 2; 3 4],1:3,adder)  -->  cat(3, [2 3; 4 5], [3 4; 5 6], [4 5; 6 7])
        # @fact size(ptableany((x,y)->myid(), 1:3, 1:4, nworkers = 2)) --> (3,4) # FIXME
        # @fact size(ptableany(1:3, 1:4, (x,y)->myid(), nworkers = 2)) --> (3,4) # FIXME
    end

    shouldtestcontext("tee") do
        a = Any[]
        @fact tee(1:3, x->push!(a,x+1))  -->  1:3
        @fact a  -->  Any[2:4]
        b = Any[]
        pushadd(x,param) = push!(b,x+param)
        c = @p tee 1 pushadd 10
        @fact b  -->  [11]
        @fact c  -->  1
    end
    shouldtestcontext("*") do
        @fact (sum*abs)([-1,0,1])  -->  2
        @fact (join*split)("a b c")  -->  "abc"
        @fact (last*join*split)("a b c")  -->  'c'
    end
    shouldtestcontext("apply") do
        f(a,b) = a+b
        g() = 1
        @fact apply(f,1,2) --> 3
        @fact apply(g) --> 1
    end
    shouldtestcontext("minelem") do
        d = [Dict(:a => 1), Dict(:a => 2), Dict(:a => 3)]
        @fact minelem(d,x->at(x,:a)) --> Dict(:a => 1)
        @fact (@p minelem d at :a) --> Dict(:a => 1)
        @fact maxelem(d,x->at(x,:a)) --> Dict(:a => 3)
        @fact (@p maxelem d at :a) --> Dict(:a => 3)
        @fact extremaelem(d,x->at(x,:a)) --> [Dict(:a => 1), Dict(:a => 3)]
        @fact (@p extremaelem d at :a) --> [Dict(:a => 1), Dict(:a => 3)]
        @fact (@p map [10,11,12] x->@p minelem [10,0,12,13] y->(abs(y-x))) --> [10,10,12]
    end
    shouldtestcontext("isany") do
        @fact isany(zeros(3), x->x==0) --> true
        @fact areall(zeros(3), x->x==0) --> true
        @fact isany(zeros(3), x->x!=0) --> false
        @fact areall(zeros(3), x->x!=0) --> false
        @fact (@p zeros 3 | areall unequal 0) --> false
        @fact (@p zeros 3 | isany unequal 0) --> false
    end
end

shouldtest("dataflow") do
    shouldtestcontext("unflatten") do
        a = Any[[1],[2,3,4],[5,6]]
        @fact unflatten(flatten(a,),a)  -->  a
    end
    shouldtestcontext("reshape") do
        @fact size(reshape(rand(9)))  -->  (3,3)
    end
    shouldtestcontext("rowcol") do
        @fact row(1)  -->  ones(Int, 1, 1)
        @fact row([1,2,3])  -->  [1 2 3]
        @fact row(1,2,3)    -->  [1 2 3]
        @fact col(1)  -->  ones(Int, 1, 1)
        @fact col([1,2,3])  -->  [1 2 3]'
        @fact col(1,2,3)    -->  [1 2 3]'
    end
    shouldtestcontext("stack") do
        @fact stack(Any[1,2]) --> [1,2]
        @fact stack(Any[[1 2],[3 4]]) --> cat(3,[1 2], [3 4])
        @fact stack(Any[zeros(2,3,4),ones(2,3,4)]) --> cat(4,zeros(2,3,4),ones(2,3,4))
    end
    shouldtestcontext("flatten") do
        @fact flatten(Any[[1],[2]]) --> [1,2]
        @fact flatten(Any[row([1]),col([2])]) --> [1 2]
        @fact flatten(Any[[1 2],[2 3]]) --> [1 2 2 3]
        @fact flatten(Any[[1 2]',[2 3]']) --> [1 2; 2 3]
        @fact flatten(Any[[1, 2],[2, 3]]) --> [1,2,2,3]
        @fact flatten(Any[[1; 2],[2; 3]]) --> [1,2,2,3]
        @fact flatten(Any[Any[[1; 2],[2; 3]]]) --> Any[[1; 2],[2; 3]]
        @fact flatten(Char['a','b','c']) --> "abc"
        @fact flatten(Char['a' 'b' 'c']) --> "abc"
        @fact flatten(["a","b","c"]) --> "abc"
        @fact flatten(["a" "b" "c"]) --> "abc"
        @fact flatten(Any["a","b","c"]) --> "abc"
        @fact flatten(Any["a" "b" "c"]) --> "abc"
    end
    shouldtestcontext("concat") do
        @fact concat([1,2]) --> [1,2]
        @fact concat([[1,2];]) --> [1,2]
        @fact concat([[1,2];3]) --> [1,2,3]
        @fact concat(ones(2,3),zeros(2,4)) --> hcat(ones(2,3),zeros(2,4))
        VERSION.minor == 3 && begin  # to avoid false 0.4 deprecation warnings
            @fact concat([[1,2];3],4) --> [1,2,3,4]
            @fact concat([[1,2],3],[4,5]) --> [1,2,3,4,5]
            @fact concat([[[1,2],3];[4,5]]) --> [1,2,3,4,5]
        end
    end
    shouldtestcontext("unstack") do
        @fact unstack(cat(3,[1 2],[2 3])) --> Any[[1 2],[2 3]]
        @fact unstack(stack(Any[[1 2],[2 3]])) --> Any[[1 2],[2 3]]
        @fact unstack(stack(Any[zeros(2,3,4),ones(2,3,4)])) --> Any[zeros(2,3,4),ones(2,3,4)]
        @fact unstack([1 2 3; 4 5 6]) --> Any[col([1,4]), col([2,5]), col([3,6])]
        @fact unstack((1,2,3)) --> Any[1,2,3]
    end
    shouldtestcontext("riffle") do
        @fact riffle([1,2,3],0) --> [1,0,2,0,3]
        @fact riffle(1,0) --> 1
        @fact riffle([1 2 3; 4 5 6],zeros(2,1)) --> [1 0 2 0 3; 4 0 5 0 6]
        @fact riffle([1 2 3; 4 5 6],[8;9]) --> [1 8 2 8 3; 4 9 5 9 6]
        @fact riffle("abc",'_') --> "a_b_c"
        @fact riffle("abc","_") --> "a_b_c"
        @fact riffle("abc",", ") --> "a, b, c"
        @fact riffle([1,2,3],0) --> [1,0,2,0,3]
    end
    shouldtestcontext("matrix") do
        a = Any[ones(2,3), zeros(2,3)]
        @fact size(matrix(a))  -->  (6,2)
        @fact size(matrix(zeros(2,3,4))) --> (6,4)
        @fact unmatrix(matrix(a),a)  -->  a
    end
    shouldtestcontext("lines") do
        @fact lines("line1\nline2\r\nline3") --> ["line1","line2","line3"]
        @fact unlines(lines("line1\nline2\r\nline3")) --> "line1\nline2\nline3"
    end
    shouldtestcontext("findsub") do
        a = [0 1 -1; 1 0 0]
        @fact findsub(a)  -->  [2 1 1; 1 2 3]
    end
    shouldtestcontext("subtoind") do
        @fact subtoind([1 1]', rand(2,3)) --> 1
        @fact subtoind([2 3]', rand(2,3)) --> 6
        @fact subtoind([1 2]', rand(2,3)) --> 3
        @fact subtoind([1 1 1]', rand(2,3,4)) --> 1
        @fact subtoind([2 3 4]', rand(2,3,4)) --> 24
        @fact subtoind([1 1 2]', rand(2,3,4)) --> 7
    end
    shouldtestcontext("randsample") do
        @fact size(randsample(1:10,5))  -->  (5,)
        @fact size(randsample(rand(2,10),20))  -->  (2,20)
        @fact randsample("aaa",5)  -->  "aaaaa"
    end
    shouldtestcontext("flip") do
        @fact flip([])  -->  []
        @fact flip("abc")  -->  "cba"
        @fact flip(1:10)  -->  10:-1:1
        b = [1,2,3]
        flip!(b)
        @fact b --> [3,2,1]
        @fact flip(Pair(1,2))  -->  Pair(2,1)
        @fact flip(Dict(1=>2,3=>4))  -->  Dict(2=>1, 4=>3)
    end
    shouldtestcontext("flipdims") do
        @fact size(flipdims(rand(2,3,4),1,1))  -->  (2,3,4)
        @fact size(flipdims(rand(2,3,4),1,3))  -->  (4,3,2)
        @fact size(flipdims(rand(2,3,4),2,1))  -->  (3,2,4)
    end
end

shouldtest("unzip") do
    @fact unzip([(1,1),(2,2)])  -->  ([1,2],[1,2])
    @fact unzip([(1,2), "ab"])  -->  (Any[1,'a'], Any[2,'b'])
end

shouldtest("pipeline") do
    shouldtestcontext("general") do
        add(x,y) = x.+y
        minus(x,y) = x.-y

        x = @p add 1 2
        @fact x --> 3

        @fact (@p add 1 2) --> 3
        @fact (@p add 1 2 | minus 2) --> 1
        @fact (@p add 1 2 | minus _ 2) --> 1
        @fact (@p add 1 2 | minus 3 _) --> 0

        @fact (@p map [1 2 3] add 1)  -->  [2 3 4]
        @fact (@p map [1 2 3] minus 1)  -->  [0 1 2]

        x = @p linspace 1 5 5 | map add 1
        @fact x  -->  [2,3,4,5,6]

        x = @p add (1:5) 1 | map add 1 | map minus 1
        @fact x  -->  [2,3,4,5,6]

        x = @p map (1:5) add 1
        @fact x  -->  [2,3,4,5,6]

        x = @p id [1] | map _ (x->x+_+_+_)
        @fact x  -->  row([4])

        @fact square(2)  --> 4
        @fact power(2,3)  --> 8
    end

    add2(a,b) = a+b

    shouldtestcontext("map2") do
        x = @p map2 1:3 10:12 add2
        @fact x  -->  [11,13,15]

        o = ones(2,3)
        z = zeros(2,3)
        @fact (@p map2 o z plus) --> ones(2,3)
        @fact (@p map3 o z z (a,b,c)->a+b+c) --> ones(2,3)
        @fact (@p map4 o z z z (a,b,c,d)->a+b+c+d) --> ones(2,3)
        @fact (@p map5 o z z z z (a,b,c,d,e)->a+b+c+d+e) --> ones(2,3)

        @fact (@p mapvec2 (1:2) [3,4] (a,b)->"$a$b") --> ["13","24"]
        @fact (@p map2 (1:2) [3,4] (a,b)->"$a$b") --> "1324"

        Z = ones(Int,3)
        add3(a,b,c) = a+b+c
        x = @p map3 1:3 10:12 100Z add3
        @fact x  -->  [111,113,115]

        add4(a,b,c,d) = a+b+c+d
        x = @p map4 1:3 10:12 100Z 2000Z add4
        @fact x  -->  [2111,2113,2115]

        add5(a,b,c,d,e) = a+b+c+d+e
        x = @p map5 1:3 10:12 100Z 2000Z 30000Z add5
        @fact x  -->  [32111,32113,32115]
    end
end

shouldtest("io") do
    shouldtestcontext("readwrite") do
        filename = tempname()
        write("line1\nline2\nline3",filename)
        @fact readstring(filename)  -->  "line1\nline2\nline3"
        @fact lines(readstring(filename))  -->  lines("line1\nline2\nline3")

        @p readstring filename | lines
        @fact (@p readstring filename | lines | unlines | lines | unlines)  -->  "line1\nline2\nline3"
    end

    # matname = string(tempname(),".mat")
    # data = ["a"-->1:3,"b"-->1,"c"-->[1:3],"d"-->"test"]
    # #FIXME writemat(matname, data)

    shouldtestcontext("filenames") do
        d = mktempdir()
        mk(a) = mkpath(joinpath(d,a))
        t(a...) = touch(joinpath(d,a...))
        mk("adir")
        mk("bdir")
        t("a")
        t("b")
        t("adir","a")
        t("adir","b")
        @fact filenames(d)  -->  ["a","b"]
        @fact  dirnames(d)  -->  ["adir","bdir"]
        rm(d, recursive = true)
    end
end

s = getstats()
if s["nNonSuccessful"] == 0
    print("   PASSED!")
else
    print("   FAILED: $(s["nFailures"]) failures and $(s["nErrors"]) errors")
end
println("   ( $(s["nNonSuccessful"]+s["nSuccesses"]) tests for runtests.jl $(join(ARGS, " ")))")
exit(s["nNonSuccessful"])


text:
## JL-ML SPGL Comparison
export ex_cs

"""
An example of using implicit JOLI operators to solve a system using GenSPGL

Use: x, r, g, info, SNR = ex_cs()
"""
function ex_cs()
    #n = 1024
    #k = 50

    # Create a sparse vector
    #inds = randperm(n)
    #x0 = zeros(n)
    #x0[inds[1:10]] = randn(10)

    # Load sparse vector x and restriction mask
    fid = matopen(Pkg.dir("GenSPGL")*"/src/Examples/pres/x0.mat")
    x0 = vec(read(fid, "x0"))
    r_inds = Int.(vec(read(fid, "r_inds")))
    n = length(vec(x0))
    close(fid)

    # Sparsity Promoting Transform
    S = joDCT(n)

    # Restriction Operator
    R = joRestriction(n, r_inds)

    # Modelling op
    A = R*S

    # Create data
    b = A*x0

    # Set BLAS Threads
    BLAS.set_num_threads(8)

    # Solve
    opts = spgOptions(verbosity = 1)
    x, r, g, info = spgl1(A, b, options = opts)
    #t = spgl1(A, b, options = opts)

    # Solve again silently for timing
    opts = spgOptions(verbosity = 0)
    t_jl = @elapsed spgl1(A, b, options = opts)

    # Calc SNR of recovery
    SNR = -20*log10(norm(x0-x)/norm(x0));

    # Write results
    fw = matopen("x_jl.mat", "w")
    write(fw, "x_jl", x)
    write(fw, "t_jl", t_jl)
    close(fw)

    return x, r, g, info, SNR
    #return t
end


text:
"""
Compute coeffs that are required to solve the memory variable eq.37 from Robertsson et. al.
"""
function get_memcoeff(mod::Medium)
	tau_sigma = mod[:tau_sigma]
	tau_epsilon = mod[:tau_epsilon]

	@assert size(tau_epsilon) == size(tau_sigma)

	nsls, nz, nx = size(tau_sigma)
	# compute factor to convert from relaxed (ω -> 0) to unrelaxed (ω -> ∞)

	K_relaxation_factor=relaxed_conversion_factor(tau_sigma, tau_epsilon)

	memcoeff2=inv.(tau_sigma) .* (1.0 .- tau_epsilon .* inv.(tau_sigma))
	for ix in 1:nx
		for iz in 1:nz
			for i in 1:nsls
				memcoeff2[i,iz,ix] *= inv(K_relaxation_factor[iz,ix])
			end
		end
	end

	memcoeff1=inv.(tau_sigma)

	return memcoeff1, memcoeff2
end


text:
module Transforms

using AutoHashEquals
using LinearAlgebra

export applytrans, applytrans_grad
export AbstractTransform
export Shift, Updim

abstract type AbstractTransform end

TransformChain = Tuple{Vararg{AbstractTransform}}

@auto_hash_equals struct Shift <: AbstractTransform
    offset :: Array{Float64,1}
end

fromdims(self::Shift) = length(self.offset)
todims(self::Shift) = length(self.offset)
codegen(::Type{Shift}, transform::Expr, out::Symbol, ::Symbol) = :($out .+= $transform.offset)
codegen_grad(::Type{Shift}, ::Expr, ::Symbol, ::Symbol, ::Symbol) = :()

@auto_hash_equals struct Updim{N,D} <: AbstractTransform
    value :: Float64
end

fromdims(self::Updim{N}) where N = N
todims(self::Updim{N}) where N = N + 1

function codegen(::Type{Updim{N,D}}, transform::Expr, out::Symbol, _) where {N,D}
    quote
        $out[$(D+1):$(N+1)] .= $out[$D:$N]
        $out[$D] = $transform.value
    end
end

function codegen_grad(::Type{Updim{1,D}}, transform::Expr, pt::Symbol, out::Symbol, ::Symbol) where D
    quote
        $out[$(D+1):2, :] .= $out[$D:1, :]
        $out[$D, :] .= 0.0
        $out[1,2] = $out[2,1]
        $out[2,2] = -$out[1,1]
    end
end

function flatten_transforms!(ret, transform, rootexpr)
    if transform <: AbstractTransform
        push!(ret, (transform, rootexpr))
    elseif transform <: Tuple
        for (i, trf) in enumerate(transform.parameters)
            flatten_transforms!(ret, trf, :($rootexpr[$i]))
        end
    elseif transform <: NamedTuple
        for (i, trf) in enumerate(transform.parameters[2].parameters)
            flatten_transforms!(ret, trf, :($rootexpr[$i]))
        end
    end
end

function flatten_transforms(transforms)
    ret = Tuple{DataType,Expr}[]
    flatten_transforms!(ret, transforms, :transforms)
    ret
end

@generated function applytrans(transforms, points::Vector{Float64},
                               workspace::Vector{Float64}, output::Vector{Float64})
    flat = flatten_transforms(transforms)
    trfcode = [codegen(dtype, code, :output, :workspace) for (dtype, code) in flat]

    quote
        output .= 0.0
        output[1:length(points)] = points
        $(trfcode...)
        output
    end
end

@generated function applytrans_grad(transforms, points::Vector{Float64},
                                    ptworkspace::Vector{Float64}, ptoutput::Vector{Float64},
                                    mxworkspace::Matrix{Float64}, mxoutput::Matrix{Float64})
    flat = flatten_transforms(transforms)
    pt_trfcode = [codegen(dtype, code, :ptoutput, :ptworkspace) for (dtype, code) in flat]
    mx_trfcode = [codegen_grad(dtype, code, :ptoutput, :mxoutput, :mxworkspace) for (dtype, code) in flat]
    trfcode = [:($mx; $pt) for (pt, mx) in zip(pt_trfcode, mx_trfcode)]

    quote
        copyto!(ptoutput, points)
        copyto!(mxoutput, I)
        $(trfcode...)
        mxoutput
    end
end

end # module

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
function f(x,θ)
  α,β,γ,δ = θ
  x1,x2=x
   dx1 = α*x1 - β*x1*x2
   dx2 = δ*x1*x2 - γ*x2
  [dx1,dx2]
end

function solve(f,x0::AbstractVector,θ,dt,N)
  X = hcat([zero(x0) for i=1:N]...)
  X[:,1]=x0
  for t=1:N-1
     X[:,t+1]=X[:,t]+dt*f(X[:,t],θ)
  end
  X
end

θ0 = [0.1,0.2,0.3,0.2]
x0 = [1.0,1.0]
dt = 0.1
N  = 1000
X=solve(f,x0,θ0,dt,N)

using Plots

p=plot(X[1,:],xlabel="t",label="x",color=:blue)
p=plot!(X[2,:],xlabel="t",label="y",color=:red)

K=100
X0 = [x0 .+ 0.1*randn(2) for k=1:K]
Xens=[X=solve(f,X0[i],θ0,dt,N) for i=1:K]

for i=1:K
  p=plot!(Xens[i][1,:],label="",color=:blue)
  p=plot!(Xens[i][2,:],label="",color=:red)
end

savefig("LV_ensemble.svg")

using StatsBase
Xm=mean(Xens)
Xstd = std(Xens)

struct GaussNum{T<:Real}
  μ::T
  σ::T
end
import Base: +, -, *, zero
+(x::GaussNum{T},a::T) where T =GaussNum(x.μ+a,x.σ)
+(a::T,x::GaussNum{T}) where T =GaussNum(x.μ+a,x.σ)
-(x::GaussNum{T},a::T) where T =GaussNum(x.μ-a,x.σ)
-(a::T,x::GaussNum{T}) where T =GaussNum(x.μ-a,x.σ)
*(x::GaussNum{T},a::T) where T =GaussNum(x.μ*a,a*x.σ)
*(a::T,x::GaussNum{T}) where T =GaussNum(x.μ*a,a*x.σ)

# TODO
# sin(x::GaussNum)= @uncertain sin(x)

+(x1::GaussNum{T},x2::GaussNum{T}) where T =GaussNum(x1.μ+x2.μ, sqrt(x1.σ.^2 + x2.σ.^2))
-(x1::GaussNum{T},x2::GaussNum{T}) where T =GaussNum(x1.μ-x2.μ, sqrt(x1.σ.^2 + x2.σ.^2))
*(x1::GaussNum{T},x2::GaussNum{T}) where T =GaussNum(x1.μ*x2.μ, sqrt((x2.μ*x1.σ).^2 + (x1.μ * x2.σ).^2))
zero(::Type{GaussNum{T}}) where T =GaussNum(zero(T),zero(T))
zero(x::AbstractVector{T}) where T =[zero(T) for i=1:length(x)]

function MV(x::AbstractArray{GaussNum{T}}) where T
  M=similar(x,T)
  V=similar(x,T)
  for i=1:length(x)
    M[i]=x[i].μ
    V[i]=x[i].σ
  end
  (M,V)
end

GX=solve(f,[GaussNum(1.0,0.1),GaussNum(1.0,0.1)],[0.1,0.2,0.3,0.2],0.1,1000)
M,V=MV(GX)
plot(M')
plot(M[1,1:30:end],errorbar=V[1,1:30:end],label="x",color=:blue)
plot!(M[2,1:30:end],errorbar=V[2,1:30:end],label="y",color=:red)

savefig("LV_GaussNum.svg")

GX=solve(f,[GaussNum(1.0,0.1),GaussNum(1.0,0.1)],[GaussNum(0.1,0.1),0.2,0.3,0.2],0.1,1000)
M,V=MV(GX)
plot(M')
plot(M[1,1:30:end],errorbar=V[1,1:30:end],label="x",color=:blue)
plot!(M[2,1:30:end],errorbar=V[2,1:30:end],label="y",color=:red)

savefig("LV_GaussNum2.svg")

using Measurements
MX=solve(f,[1.0±0.1,1.0±0.1],[0.1,0.2,0.3,0.2],0.1,1000)
plot(MX[1,1:30:end],label="x",color=:blue)
plot!(MX[2,1:30:end],label="y",color=:red)

savefig("LV_Measurements.svg")

MX=solve(f,[1.0±0.1,1.0±0.1],[0.1±0.01,0.2±0.01,0.3±0.01,0.2±0.01],0.1,1000)
plot(MX[1,1:30:end],label="x",color=:blue)
plot!(MX[2,1:30:end],label="y",color=:red)

savefig("LV_Measurements2.svg")

# Plot receipe
# plot(Vector{GaussNum})

using LinearAlgebra
function solve(f,x0::AbstractVector,sqΣ0, θ,dt,N,Nr)
   n = length(x0)
   n2 = 2*length(x0)
   Qp = sqrt(n)*[I(n) -I(n)]

  X = hcat([zero(x0) for i=1:N]...)
  S = hcat([zero(x0) for i=1:N]...)
  X[:,1]=x0
  Xp = x0 .+ sqΣ0*Qp
  sqΣ = sqΣ0
  Σ = sqΣ* sqΣ'
  S[:,1]= diag(Σ)
  for t=1:N-1
    if rem(t,Nr)==0
      Xp .= X[:,t] .+ sqΣ * Qp
    end
    for i=1:n2 # all quadrature points
      Xp[:,i].=Xp[:,i] + dt*f(Xp[:,i],θ)
    end
    mXp=mean(Xp,dims=2)
    X[:,t+1]=mXp
    Σ=Matrix((Xp.-mXp)*(Xp.-mXp)'/n2)
    S[:,t+1]=sqrt.(diag(Σ))
    # @show Σ

    sqΣ = cholesky(Σ).L

  end
  X,S
end

## Extension to arbitrary

QX,QS=solve(f,[1.0,1.0],(0.1)*I(2),θ0,0.1,1000,1)
plot(QX[1,1:30:end],label="x",color=:blue,errorbar=QS[1,1:30:end])
plot!(QX[2,1:30:end],label="y",color=:red,errorbar=QS[2,1:30:end])

savefig("LV_Quadrics.svg")


text:
using Mmap

########################
# fractions of formula #
########################

factorial_div(i::Integer, j::Integer) = factorial(big(i)) / factorial(big(j))

z(x::Real, p::Real) = sqrt(2.)*(x + p*im)
z(x::Vector{<:Real}, p::Vector{<:Real}) = z.(x, p')

gaussian_function(x::Real, p::Real) = exp(-0.5 * abs2(z(x,p))) / π
gaussian_function(x::Vector{<:Real}, p::Vector{<:Real}) = gaussian_function.(x, p')

neg_one_to_power_of(i::Integer) = (i % 2 == 0) ? 1 : -1

function coefficient_of_wave_function(m::Integer, n::Integer)
    if n ≥ m
        # adjust index bases for number state on `m`
        return neg_one_to_power_of(m-1) * sqrt(factorial_div(m, n))
    else
        # adjust index bases for number state on `n`
        return neg_one_to_power_of(n-1) * sqrt(factorial_div(n, m))
    end
end

function z_to_power(m::Integer, n::Integer, x::Real, p::Real)
    if n ≥ m
        return conj(z(x, p'))^(n - m)
    else
        return z(x, p')^(m - n)
    end
end

function z_to_power(m::Integer, n::Integer, x::Vector{<:Real}, p::Vector{<:Real})
    z_to_power.(m, n, x, p')
end

function laguerre(m::Integer, n::Integer, x::Real, p::Real)
    if n ≥ m
        # adjust index bases for number state on `m`
        return laguerre(m-1, n - m, abs2(z(x, p)))
    else
        # adjust index bases for number state on `n`
        return laguerre(n-1, m - n, abs2(z(x, p)))
    end
end

function laguerre(m::Integer, n::Integer, x::Vector{<:Real}, p::Vector{<:Real})
    return laguerre.(m, n, x, p')
end

#########
# utils #
#########

function save_𝐰(bin_name::String, 𝐰::Array{ComplexF64,4})
    global my_artifacts
    hash = create_my_artifact() do artifact_dir
        path = joinpath(artifact_dir, bin_name)
        mem = open(path, "w+")
        write(mem, 𝐰)
        close(mem)

        return path
    end
    bind_my_artifact!(my_artifacts[], bin_name, hash)

    @info "Save Wₘₙₓₚ to $bin_name"

    return my_artifact_path(hash)
end

function load_𝐰(
    bin_path::String,
    x_range::AbstractRange,
    p_range::AbstractRange,
    dim::Integer
)
    mem = open(bin_path)
    𝐰 = Mmap.mmap(
        mem,
        Array{ComplexF64,4},
        (dim, dim, length(x_range), length(p_range))
    )
    close(mem)

    @info "Load Wₘₙₓₚ from $bin_path"

    return 𝐰
end

function gen_wigner_bin_name(
    x_range::AbstractRange,
    p_range::AbstractRange,
    dim::Integer,
)
    return string(hash((x_range, p_range, dim)))
end


text:
# Map from SE(3) embedded in R12
if !(@isdefined LinkSpherePosition)
    @computed struct LinkSpherePosition{R12,R3,S,n} <: TaskMap{R12,R3,S}
        point_cache::fulltype(FramePointCache{n,S})
    end
end
struct LinkSpherePositionT{R12,R3,S,n} <: TaskMapT{R12,R3,S}
    base_map::LinkSpherePosition{R12,R3,S,n}
end
LinkSpherePosition{R12,R3}(point_cache::FramePointCache{n,S}) where {n,S} =
    LinkSpherePosition{R12,R3,S,n}(point_cache)

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::LinkSpherePosition{R12,R3}) =
    point_world_position(xme, task_map.point_cache.point.v)

function task_jacobian_emb(xme, task_map::LinkSpherePosition{R12,R3}, arg...)
    Jfc = task_map.point_cache.world_position_jacobian
    point_frame_position = task_map.point_cache.point.v
    update_point_world_position_jacobian!(Jfc, xme, point_frame_position)
    Jf = Jfc.data[1]
end

task_jacobian_chart(xm, task_map::LinkSpherePosition{R12,R3}, CM::Chart{1,R12}, CN::Chart{1,R3}) =
    task_jacobian_emb(xm, task_map)

function task_jacobian_emb_dot(xme, vme, task_map::LinkSpherePosition{R12,R3}, arg...)
    Jfc_dot = task_map.point_cache.world_position_jacobian_dot
    point_frame_position = task_map.point_cache.point.v
    update_point_world_position_jacobian_dot!(Jfc_dot, xme, vme, point_frame_position)
    Jf_dot = Jfc_dot.data[1]
end
task_jacobian_chart_dot(xm, vm, task_map::LinkSpherePosition{R12,R3},
    CM::Chart{1,R12}, CN::Chart{1,R3}) = task_jacobian_emb_dot(xm, vm, task_map)


text:
module RocketCollectLatestObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "CollectLatestObservable" begin

    println("Testing: collectLatest")

    run_testset([
        (
            source = collectLatest([ of(1), from(1:5) ]),
            values = @ts([ [1, 1], [1, 2], [1, 3], [1, 4], [1, 5], c ]),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest([ of(1), from(1:5) ], (v) -> v .^ 2),
            values = @ts([ [1, 1], [1, 4], [1, 9], [1, 16], [1, 25], c ]),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest(Int, Int, [ of(1), from(1:5) ], sum),
            values = @ts([ 2, 3, 4, 5, 6, c ]),
            source_type = Int
        ),
        (
            source = collectLatest([ of(1) |> async(0), from(1:5) ]),
            values = @ts([ [1, 5] ] ~ c),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest([ of(1) |> async(0), from(1:5) |> async(0) ]),
            values = @ts([ [1, 1] ] ~ [ [1, 2] ] ~ [ [1, 3] ] ~ [ [1, 4] ] ~ [ [1, 5] ] ~ c),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest([ from(1:5), of(2.0) ]),
            values = @ts([ [5, 2], c ]),
            source_type = Vector{Union{Float64, Int}}
        ),
        (
            source = collectLatest([ from(1:5) |> async(0), of(2.0) ]),
            values = @ts([ [1, 2.0] ] ~ [ [2, 2.0] ] ~ [ [3, 2.0] ] ~ [ [4, 2.0] ] ~ [ [5, 2.0] ] ~ c),
            source_type = Vector{Union{Float64, Int}}
        ),
        (
            source = collectLatest([ completed(Int), of(1) ]),
            values = @ts(c),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest([ completed(Int), completed(Int) ]),
            values = @ts(c),
            source_type = Vector{Int}
        ),
        (
            source = collectLatest([ faulted(Float64, "err"), completed(Int) ]),
            values = @ts(e("err")),
            source_type = Vector{Union{Float64, Int}}
        ),
        (
            source = collectLatest([ completed(Int), faulted(Float64, "err") ]),
            values = @ts(c),
            source_type = Vector{Union{Float64, Int}}
        ),
        (
            source = collectLatest([ faulted(Float64, "err1"), faulted(Float64, "err2") ]),
            values = @ts(e("err1")),
            source_type = Vector{Float64}
        ),
    ])

    somenumbers    = Union{Float64, Int}[ 1, 1, 1, 1, 1.0, 1.0, 1.0, 1.0 ]
    expectedoutput = somenumbers

    collected100 = collectLatest(map(n -> of(n), somenumbers))

    run_testset([
        (
            source = collected100,
            values = @ts([ [ 1, 1, 1, 1, 1.0, 1.0, 1.0, 1.0 ], c ]),
            source_type = Vector{Union{Float64, Int}}
        )
    ])

    @testset begin
        s1 = Subject(Int)
        s2 = Subject(Float64)
        s3 = Subject(String)
        s4 = of(5)
        s5 = from(1:10)

        latest = collectLatest([ s1, s2, s3, s4, s5 ])
        values = Vector{Any}()
        subscription = subscribe!(latest, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        ))

        @test values == [ ]

        next!(s1, 1)

        @test values == [ ]

        next!(s2, 2.0)

        @test values == [ ]

        next!(s3, "Hello")

        @test values == [ [1, 2.0, "Hello", 5, 10] ]

        next!(s1, 2)

        @test values == [ [1, 2.0, "Hello", 5, 10] ]

        next!(s3, "Hello, world!")

        @test values == [ [1, 2.0, "Hello", 5, 10]  ]

        complete!(s1);

        next!(s2, 3.0)

        @test values == [ [1, 2.0, "Hello", 5, 10], [2, 3.0, "Hello, world!", 5, 10]  ]

        complete!(s2)

        @test values == [ [1, 2.0, "Hello", 5, 10], [2, 3.0, "Hello, world!", 5, 10] ]

        next!(s3, "Something")

        @test values == [ [1, 2.0, "Hello", 5, 10], [2, 3.0, "Hello, world!", 5, 10], [2, 3.0, "Something", 5, 10] ]

        complete!(s3)

        @test values == [ [1, 2.0, "Hello", 5, 10], [2, 3.0, "Hello, world!", 5, 10], [2, 3.0, "Something", 5, 10], "completed" ]

        unsubscribe!(subscription)
    end

end

end


text:
@testset "fsm_stream_reassembler_holes.jl" begin

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""

        @test !eof(rs.output)
    end
    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        @test rs.cur_idx == 0
        push_substring!(rs, "a", 0)
        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1, true)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "a", 0)
        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"
        @test eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        push_substring!(rs, "ab", 0)

        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "d", 3)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "c", 2)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "ab", 0)

        @test rs.cur_idx == 4
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "d", 3)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "abc", 0)

        @test rs.cur_idx == 4
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "b", 1)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "d", 3)
        @test rs.cur_idx == 0
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "a", 0)
        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"
        @test !eof(rs.output)

        push_substring!(rs, "c", 2)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "cd"
        @test !eof(rs.output)

        push_substring!(rs, "", 4, true)

        @test rs.cur_idx == 4
        @test read!(rs.output) == ""
        @test eof(rs.output)
    end

end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge