text:
module RocketLimitSubscribersOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: limit_subscribers()" begin

    println("Testing: operator limit_subscribers()")

    run_proxyshowcheck("LimitSubscribers", limit_subscribers(), args = (check_subscription = true, ))

    run_testset([
        (
            source = from(1:5) |> limit_subscribers(),
            values = @ts([ 1, 2, 3, 4, 5, c ])
        ),
        (
            source = completed() |> limit_subscribers(),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> limit_subscribers(),
            values = @ts(e(1))
        ),
        (
            source = never() |> limit_subscribers(),
            values = @ts()
        )
    ])

    @testset begin
        values1 = Int[]
        values2 = Int[]
        values3 = Int[]

        guard   = LimitSubscribersGuard(1, true)
        subject = Subject(Int)
        source  = subject |> limit_subscribers(guard)

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        @test values1 == [ ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        next!(subject, 1)

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription2 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values2, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values2, -2)
        ))

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 2)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        subscription3 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values3, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values3, -3)
        ))

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ ]

        next!(subject, 3)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3 ]

        source  = subject |> limit_subscribers(guard)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3, -3 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
        unsubscribe!(subscription3)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3, -3 ]
    end

    @testset begin
        values1 = Int[]
        values2 = Int[]
        values3 = Int[]

        guard   = LimitSubscribersGuard(2, true)
        subject = Subject(Int)
        source  = subject |> limit_subscribers(guard)

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        @test values1 == [ ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        next!(subject, 1)

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription2 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values2, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values2, -2)
        ))

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 2)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        subscription3 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values3, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values3, -3)
        ))

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        next!(subject, 3)

        @test values1 == [ 1, 1,2  ]
        @test values2 == [ 2, 3 ]
        @test values3 == [ 3 ]

        source  = subject |> limit_subscribers(guard)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2, 3, -2 ]
        @test values3 == [ 3, -3 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
        unsubscribe!(subscription3)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2, 3, -2 ]
        @test values3 == [ 3, -3 ]
    end

    @testset begin
        values1 = Int[]
        values2 = Int[]
        values3 = Int[]

        guard   = LimitSubscribersGuard(1, false)
        subject = Subject(Int)
        source  = subject |> limit_subscribers(guard)

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        @test values1 == [ ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        next!(subject, 1)

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription2 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values2, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values2, -2)
        ))

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 2)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        subscription3 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values3, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values3, -3)
        ))

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ ]

        next!(subject, 3)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3 ]

        source  = subject |> limit_subscribers(guard)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3 ]

        subscribe!(source, void())

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3, -3 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
        unsubscribe!(subscription3)

        @test values1 == [ 1, 1, -1 ]
        @test values2 == [ 2, -2 ]
        @test values3 == [ 3, -3 ]
    end

    @testset begin
        values1 = Int[]
        values2 = Int[]
        values3 = Int[]

        guard   = LimitSubscribersGuard(2, false)
        subject = Subject(Int)
        source  = subject |> limit_subscribers(guard)

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        @test values1 == [ ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        next!(subject, 1)

        @test values1 == [ 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription1 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values1, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values1, -1)
        ))

        next!(subject, 1)

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        subscription2 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values2, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values2, -2)
        ))

        @test values1 == [ 1, 1 ]
        @test values2 == [ ]
        @test values3 == [ ]

        next!(subject, 2)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        unsubscribe!(subscription1)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        subscription3 = subscribe!(source, lambda(
            on_next     = (d) -> push!(values3, d),
            on_error    = (e) -> error(e),
            on_complete = () -> push!(values3, -3)
        ))

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2 ]
        @test values3 == [ ]

        next!(subject, 3)

        @test values1 == [ 1, 1,2  ]
        @test values2 == [ 2, 3 ]
        @test values3 == [ 3 ]

        source  = subject |> limit_subscribers(guard)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2, 3 ]
        @test values3 == [ 3 ]

        subscribe!(source, void())

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2, 3, -2 ]
        @test values3 == [ 3 ]

        unsubscribe!(subscription1)
        unsubscribe!(subscription2)
        unsubscribe!(subscription3)

        @test values1 == [ 1, 1, 2 ]
        @test values2 == [ 2, 3, -2 ]
        @test values3 == [ 3 ]
    end

end

end


text:
import BioSimulator: execute_jump!, get_site, initialize_datastructs!, get_neighbor_class, get_ptype, coordinates

@testset "Execute" begin
  @testset "Predator-Prey Model" begin
    ##### reaction definition #####
    intermediate = @def_reactions begin
      A + B --> A + A, k1 # predation / predator reproduction
      B + 0 --> 0 + B, k2 # prey movement
      A + 0 --> 0 + A, k3 # predator movement
      B + 0 --> B + B, k4 # prey reproduction
      A --> 0, k5         # predator death
      B --> 0, k6         # prey death
    end k1 k2 k2 k3 k4 k5 k6

    ##### assign parameter values #####
    params = [2.0, 1.0, 1.0, 2.0, 1.0, 1.0]

    ##### generate the full reaction list #####
    model = @enumerate_with_sclass intermediate VonNeumann() 2 params

    ##### initial configuration #####
    initial_sites = [(0,1), (1,0), (1,1), (1,2), (1,3), (2,0), (2,1), (2,2), (3,0), (3,1)]
    initial_sites = [initial_sites[j][i] for i in 1:2, j in 1:length(initial_sites)]
    ptypes        = [1, 1, 1, 2, 2, 1, 2, 2, 2, 2]

    lattice = Lattice(initial_sites, ptypes)

    initialize_datastructs!(lattice, model)
    enum = model.enumeration

    # Test Reaction:
    execute_jump!(lattice, model, 9)

    # aggregate coordinates
    test_coord = [tuple(coordinates(x)...) for x in lattice.site]

    # check new open sites
    @test (2, -1) in test_coord

    # check type changes
    @test get_ptype(get_site(lattice, (2,  0))) == 1
    @test get_ptype(get_site(lattice, (2, -1))) == 2

    # check neighbor class changes
    # (2, 1): 13 -> 9
    # (1, 0): 4  -> 2
    # (3, 0): 5  -> 3
    # (2, 0): 9  -> 13
    # (2, -1): 2 -> 1
    points_to_check = [(2,1), (1,0), (3,0), (2,0), (2,-1)]
    list = [get_site(lattice, x) for x in points_to_check]

    @test get_neighbor_class(list[1]) == 9
    @test get_neighbor_class(list[2]) == 2
    @test get_neighbor_class(list[3]) == 3
    @test get_neighbor_class(list[4]) == 13
    @test get_neighbor_class(list[5]) == 1

    # check sample class changes
    # (2, 1): [14] -> [14, 5]
    # (1, 0): [10, 13] -> [13, 11]
    # (3, 0): [6, 14] -> [14, 7]
    # (2, 0): [2,9,13] -> []
    # (2, -1): [] -> [12, 13]
    particle_sample_classes = [
      (2, 1)  => [5, 14],
      (1, 0)  => [11, 13],
      (3, 0)  => [7, 14],
      (2, 0)  => Int[],
      (2, -1) => [12, 13]
      ]

    for (coord, expected_sample_classes) in particle_sample_classes
      x = get_site(lattice, coord)

      l = get_ptype(x)
      k = get_neighbor_class(x)

      test_sample_classes = sort(get(enum.pair_to_classes, (l,k), Int[]))

      @test test_sample_classes == expected_sample_classes
    end

    # check invariants
    @test issorted(lattice.site)
    @test issorted(lattice.coord_order, by = coordinates)
  end
end

ext: jl


text:
module YaoBlocks

using LinearAlgebra

include("utils.jl")
# include("traits.jl")

include("abstract_block.jl")

# concrete blocks
include("routines.jl")
include("primitive/primitive.jl")
include("composite/composite.jl")

include("algebra.jl")
include("blocktools.jl")
include("layout.jl")
include("treeutils/treeutils.jl")

include("autodiff/autodiff.jl")
export AD

include("deprecations.jl")

end # YaoBlocks


text:
for idx = (12):-1:6
  data = readdlm(string("posit2-", idx, ".csv"))

  print(sum(data)/1000,",,")
  print(sum(data .== 0.0)/1000,",")
  println(sum(data .<= 5.0)/1000)
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
export pswap, singlet_block, basis_rotor
using YaoBlocks
using YaoArrayRegister: mulrow!, u1rows!

function Yao.apply!(reg::ArrayReg, pb::PutBlock{N,2,RotationGate{2,T,G}}) where {N,T,G<:SWAPGate}
    mask1 = bmask(pb.locs[1])
    mask2 = bmask(pb.locs[2])
    mask12 = mask1|mask2
    a, c, b_, d = mat(Rx(pb.content.theta))
    e = exp(-im/2*pb.content.theta)
    state = statevec(reg)
    for b in basis(reg)
        if b&mask1==0
            i = b+1
            i_ = b ⊻ mask12 + 1
            if b&mask2==mask2
                u1rows!(state, i, i_, a, b_, c, d)
            else
                mulrow!(state, i, e)
                mulrow!(state, i_, e)
            end
        end
    end
    return reg
end

mutable struct Bag{N}<:TagBlock{AbstractBlock, N}
    content::AbstractBlock{N}
end

Yao.content(bag) = bag.content
Yao.chcontent(bag::Bag, content) = Bag(content)
Yao.mat(bag::Bag) = mat(bag.content)
Yao.apply!(reg::AbstractRegister, bag::Bag) = apply!(reg, bag.content)
YaoBlocks.PreserveStyle(::Bag) = YaoBlocks.PreserveAll()
setcontent!(bag::Bag, content) = (bag.content = content; bag)

function YaoBlocks.print_annotation(io::IO, bag::Bag)
    printstyled(io, "[⊞] "; bold=true, color=:blue)
end

"""parametrized swap gate."""
function pswap(nbit::Int, i::Int, j::Int)
    put(nbit, (i,j)=>rot(SWAP, 0.0))
end

"""block for generating singlets."""
function singlet_block(nbit::Int, i::Int, j::Int)
    unit = chain(nbit)
    push!(unit, put(nbit, i=>chain(X, H)))
    push!(unit, control(nbit, -i, j=>X))
end

basis_rotor(::ZGate) = I2Gate()
basis_rotor(::XGate) = Ry(-0.5π)
basis_rotor(::YGate) = Rx(0.5π)

basis_rotor(basis::PauliGate, nbit, locs) = repeat(nbit, basis_rotor(basis), locs)


text:
#
# Create multiscene plots for Makie, PyPlot and VTKView
#

import Pkg
Pkg.activate(mktempdir())
Pkg.add(["GLMakie","PyPlot","ExtendableGrids"])
Pkg.Registry.add(Pkg.RegistrySpec(url="https://github.com/j-fu/PackageNursery"))
Pkg.add(Pkg.PackageSpec(name="VTKView", uuid="63b67fce-66f7-11ea-1808-1361d8c55bf4"))

import PyPlot
import GLMakie
import VTKView

using ExtendableGrids

include("../plotting.jl")

function test_multiscene()
    plotting_multiscene(Plotter=PyPlot)
    plotting__multiscene(Plotter=VTKView)
    plotting__multiscene(Plotter=GLMakie)
end


text:
nsamps_to_thin(nsamps::Int, nmcmc::Int) = max(1, div(nmcmc, nsamps))

monitor1 = [:theta__Z, :theta__v, :theta__alpha,
            :omega, :r, :theta__lam, :W_star, :theta__eta,
            :theta__W, :theta__delta, :theta__sig2]

monitor2 = [:theta__y_imputed, :theta__gam]

function _convertStateToDicParam(s::State, c::Constants, d::Data)::DICparam
  I = d.I
  N = d.N
  J = d.J
  beta = c.beta

  p = [[prob_miss(s.y_imputed[i][n, j], beta[:, i])
        for n in 1:N[i], j in 1:J] for i in 1:I]

  mu = [[mus(i, n, j, s, c, d)
         for n in 1:N[i], j in 1:J] for i in 1:I]

  sig = [fill(sqrt(s.sig2[i]), N[i]) for i in 1:I]

  y = deepcopy(s.y_imputed)

  return DICparam(p, mu, sig, y)
end

function initDicStream(data::Data)
  param = DICparam(p=deepcopy(data.y),
                   mu=deepcopy(data.y),
                   sig=[zeros(Float64, data.N[i]) for i in 1:data.I],
                   y=deepcopy(data.y))
  return MCMC.DICstream{DICparam}(param)
end

"""
printFreq: defaults to 0 => prints every 10%. turn off printing by
           setting to -1.
thins: defaults to `[2, nsamps_to_thin(10, nmcmc)]`.
       i.e. Keep every other sample for monitor 1 (so the total number of
       samples is `nmcmc/2`). And keep 10 samples total for monitor 2.
"""
function fit_fs!(init::StateFS, c::ConstantsFS, d::DataFS;
                 nmcmc::Int, nburn::Int,
                 tuners::Union{Nothing, TunersFS}=nothing,
                 monitors=[monitor1, monitor2],
                 fix::Vector{Symbol}=Symbol[],
                 thins::Vector{Int}=[2, nsamps_to_thin(10, nmcmc)],
                 thin_dden::Int=1,
                 printFreq::Int=0,
                 computeDIC::Bool=false, computeLPML::Bool=false,
                 computedden::Bool=false,
                 sb_ibp::Bool=false,
                 use_repulsive::Bool=true, Z_marg_lamgam::Bool=true,
                 verbose::Int=1, time_updates::Bool=false, Z_thin::Int=1,
                 seed::Int=-1)

  # Set random seed if needed
  if seed >= 0
    Random.seed!(seed)
  end

  # We don't want to use noisy distribution.
  # Assert that all eps == 0
  if any(init.theta.eps .!= 0)
    msg = "WARNING: `init.theta.eps`: $(init.theta.eps) contains non-zeros values! "
    msg *= "Setting all eps to 0!"
    println(msg)

    # Set eps == 0
    init.theta.eps .= 0
  end
  flush(stdout)

  @assert Z_thin >= 0
  if Z_thin == 0
    Z_thin = d.data.J
  end

  if verbose >= 1
    fixed_vars_str = join(fix, ", ")
    if fixed_vars_str == ""
      fixed_vars_str = "nothing"
    end
    println("fixing: $fixed_vars_str")
    println("Use stick-breaking IBP: $(sb_ibp)")
    println("Z_marg_lamgam: $(Z_marg_lamgam)")
    println("use_repulsive: $(use_repulsive)")
    println("Z_thin: $(Z_thin)")
    flush(stdout)
  end

  @assert printFreq >= -1
  if printFreq == 0
    numPrints = 10
    printFreq = Int(ceil((nburn + nmcmc) / numPrints))
  end

  # Initialize tuners if needed
  if tuners == nothing
    tuners = TunersFS(Tuners(d.data.y, c.constants.K),
                      init.theta,
                      d.data.X)
  end

  function printMsg(iter::Int, msg::String)
    if printFreq > 0 && iter % printFreq == 0
      print(msg)
    end
  end

  # Loglike
  loglike = Float64[]

  # Instantiate (but not initialize) CPO stream
  if computeLPML
    cpoStream = MCMC.CPOstream{Float64}()
  end

  # DIC
  if computeDIC
    dicStream = initDicStream(d.data)
    loglikeDIC(param::DICparam) = computeLoglikeDIC(d.data, param)

    convertStateToDicParam(s::State)::DICparam = let
      _convertStateToDicParam(s, c.constants, d.data)
    end
  end

  # Cache for data density
  dden = Matrix{Vector{Float64}}[]

  function update!(s::StateFS, iter::Int, out)
    update_state_feature_select!(s, c, d, tuners,
                                 ll=loglike, fix=fix,
                                 use_repulsive=use_repulsive,
                                 Z_marg_lamgam=Z_marg_lamgam,
                                 sb_ibp=sb_ibp, time_updates=time_updates,
                                 Z_thin=Z_thin)

    if computedden && iter > nburn && (iter - nburn) % thin_dden == 0
      # NOTE: `datadensity(s, c, d)` returns an (I x J) matrix of vectors of
      # length g.
      append!(dden, [datadensity(s.theta, c.constants, d.data)])
    end

    if computeLPML && iter > nburn
      # Inverse likelihood for each data point
      like = [[compute_like(i, n, s.theta, c.constants, d.data)
               for n in 1:d.data.N[i]] for i in 1:d.data.I]

      # Update (or initialize) CPO
      MCMC.updateCPO(cpoStream, vcat(like...))

      # Add to printMsg
      printMsg(iter, " -- LPML: $(MCMC.computeLPML(cpoStream))")
    end

    if computeDIC && iter > nburn
      # Update DIC
      MCMC.updateDIC(dicStream, s.theta, updateParams,
                     loglikeDIC, convertStateToDicParam)

      # Add to printMsg
      printMsg(iter, " -- DIC: $(MCMC.computeDIC(dicStream, loglikeDIC,
                                                 paramMeanCompute))")

      DICg = MCMC.DIC_gelman(loglike[(nburn+1):end])
      printMsg(iter, " -- DIC_gelman: $(DICg)")
    end

    printMsg(iter, "\n")
    flush(stdout)
  end

  # Loglike of initial state
  init_ll = compute_marg_loglike(init.theta, c.constants,
                                 d.data, c.constants.temper)

  if isinf(init_ll)
    println("Warning: Initial state yields likelihood of zero.")
    msg = """
    It is likely the case that the initialization of missing values
    is not consistent with the provided missing mechanism. The MCMC
    will almost certainly reject the initial values and sample new
    ones in its place.
    """
    println(join(split(msg), " "))
  else
    println("")
  end

  out, lastState = MCMC.gibbs(init, update!, monitors=monitors,
                              thins=thins, nmcmc=nmcmc, nburn=nburn,
                              printFreq=printFreq,
                              loglike=loglike,
                              printlnAfterMsg=false)

  mega_out = Dict{Symbol, Any}(:samples => out,
                               :lastState => lastState,
                               :loglike => loglike,
                               :c => c,
                               :init => init)

  if computeDIC || computeLPML
    LPML = computeLPML ? MCMC.computeLPML(cpoStream) : NaN
    Dmean, pD = computeDIC ? MCMC.computeDIC(dicStream,
                                             loglikeDIC,
                                             paramMeanCompute,
                                             return_Dmean_pD=true) : (NaN, NaN)
    DICg = MCMC.DIC_gelman(loglike[(nburn+1):end])

    metrics = Dict(:LPML => LPML,
                   :DIC => Dmean + pD,
                   :Dmean => Dmean,
                   :pD => pD,
                   :DICg => DICg)
    println()
    println("metrics:")
    for (k, v) in metrics
      mega_out[k] = v
      println("$k => $v")
    end
  end
  flush(stdout)

  if computedden
    mega_out[:dden] = dden
  end

  mega_out[:nburn] = nburn
  mega_out[:Z_thin] = Z_thin
  mega_out[:m] = Matrix{Bool}.(d.data.m)

  return mega_out
end


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

# Data analysis, export SIR populations and parameters

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

"""
    get_populations(agents, params)
Get number of agents in each state:
    - S: Suceptible
    - I: Infected
    - R: Recovered
    - Q: Quarantained (not implemented)
    - D: Decesead (not implemented)
"""
function get_populations_old(agents, params)
    # species = num_S, num_I, num_R, num_Q
    species = zeros(4)

    for ag in agents
        if ag.state == "S"
            species[1] += 1.0
        elseif ag.state == "I"
            species[2] += 1.0
        elseif ag.state == "R"
            species[3] += 1.0
        # else
        #     species[4] += 1.0
        end
    end
    #  return species ./ params.num_agents
    return species
end
##=================####==============##

"""
    get_populations(agents)
It returns a dictionary with the number
of agents in each state. Might be:
    - S: Suceptible
    - I: Infected
    - R: Recovered
    - Q: Quarantained (not implemented)
    - D: Decesead (not implemented)
"""
function get_populations(agents)
    map(x -> x.state, agents) |> countmap
end

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

"""
    export_parameters(params, out_path)
Export parameters to file
"""
function export_parameters(params, out_path)

    filename = "parameters_N_$(params.num_agents)_rep_$(params.repetition).txt"

    open(joinpath(out_path, filename), "w") do io
        println(io, "num_agents   :\t", params.num_agents   )
        println(io, "p_link       :\t", params.p_link       )
        println(io, "p_infected_t0:\t", params.p_infected_t0)
        println(io, "attack_rate  :\t", params.attack_rate  )
        println(io, "repetition   :\t", params.repetition   )
    end

end
##=================####==============##


text:
using EquationsSolver
using Test
using Symbolics

@testset "check_eqs" begin
    @variables x, y
    res = EquationsSolver.check_eqs(x + y ~ 1)
    @test isequal(res, [x + y - 1])
    res = EquationsSolver.check_eqs([x + y ~ 1])
    @test isequal(res, [x + y - 1])
end

@testset "get_vars" begin
    @variables x, y
    res = EquationsSolver.get_vars(x + y)
    @test EquationsSolver.isin(x, res) == true
    @test EquationsSolver.isin(x, res) == true
end

@testset "get_all_vars" begin
    @variables x, y, z
    expr1 = EquationsSolver.toexpr([x - y ~ 0, x + y ~ 1])
    expr2 = EquationsSolver.toexpr([x - y ~ 0, x + y - z ~ 1])
    res1 = EquationsSolver.get_all_vars(expr1)
    res2 = EquationsSolver.get_all_vars(expr2)
    @test isequal(res1, [x, y])
    @test isequal(res2, [x, y, z])
end

@testset "check_vars" begin
    @variables x, y
    eqs = [
        x - y - 0,
        x + y - 1
    ]
    vars = Dict(x => 1.0, y => 1.0)
    @test EquationsSolver.check_vars(eqs, vars)
end

@testset "LinearProblem" begin
    @variables x, y
    eqs = [
        x - y ~ 0,
        x + y ~ 1
    ]
    vars = Dict(x => 1.0, y => 1.0)
    res = LinearProblem(eqs, vars)
    @test typeof(res) == LinearProblem
end

@testset "solve" begin
    @variables x, y
    eqs = [
        x + y ~ 3.0,
        x + 2y ~ 5.0
    ]
    vars = Dict(x => 1.0, y => 1.0)
    res = EquationsSolver.solve(LinearProblem(eqs, vars))
    @test res == Dict(x => 1.0, y => 2.0)
end

@testset "one_cal" begin
    @variables x
    eqs = [
        x + 5 - exp(x)
    ]
    vars = Dict(x => 2.0)
    jac = Symbolics.jacobian(eqs, [x])
    res = EquationsSolver.one_cal([1.0], vars, eqs, jac)
    @test length(res) == 1
end

@testset "get_dict" begin
    @variables x, y, z
    value = [1.0, 2.0, 3.0]
    vars = Dict(x => 2.0, y => 1.0, z => 3.0)
    res = EquationsSolver.get_dict(value, vars, [x, y, z])
    @test res == Dict(x => 1.0, y => 2.0, z => 3.0)
end

@testset "_solve" begin
    @variables x, y, z
    eqs = [
        x + 5 - exp(x)
    ]
    vars = Dict(x => 2.0)
    jac = Symbolics.jacobian(eqs, [x])
    res = EquationsSolver._solve(vars, eqs, jac, 1000, 1.0E-6, [x])
    @test length(res) == 1

    eqs = [
        x + y + z - 3,
        x + 4y + 9z - 14,
        x + 2y + 3z - 6
    ]
    vars = Dict(x => 1.0, y => 0.0, z => 0.0)
    jac = Symbolics.jacobian(eqs, [x, y, z])
    res = EquationsSolver._solve(vars, eqs, jac, 10000, 1.0E-6, [x, y, z])
    @test round.(collect(values(res)), digits=1) == [1.0,1.0,1.0]
end

@testset "solve-NLProblem" begin
    @variables x
    eqs = [
        x + 5 ~ exp(x)
    ]
    vars = Dict(x => 2.0)
    pro = NLProblem(eqs,vars)
    res = solve(pro)
    println(res)
    @test length(res) == 1
end

@testset "solve-LinearProblem" begin
    @variables x,y,z
    eqs = [
        x + y + z ~ 3,
        x + 4y + 9z ~ 14,
        x + 2y + 3z ~ 6
    ]
    vars = Dict(x => 1.0, y => 0.0, z => 0.0)
    pro = LinearProblem(eqs,vars)
    res = solve(pro)
    @show res
    @test round.(collect(values(res)), digits=1) == [1.0,1.0,1.0]
end

@testset "NLProblem_1: x , y" begin
    @variables x, y
    eqs = [
        cos(x^2 + 0.4 * y) + x^2 + y^2 ~ 1.6
        1.5 * x^2 - 1 / 0.36 * y^2 ~ 1.0
    ]
    vars = Dict(x => 2.0, y => 1.0)
    pro = NLProblem(eqs, vars)
    res = solve(pro)
    @show res
    @test round(res[x], digits=2) == 1.04 &&
          round(res[y], digits=2) == 0.47
end

@testset "NLProblem_2: x , y" begin
    @variables x, y
    eqs = [
        4 * x + 0.1 * exp(x) ~ y + 1
        4 * y + x^2 / 8 ~ x
    ]
    vars = Dict(x => 2.0, y => 1.0)
    pro = NLProblem(eqs, vars)
    res = solve(pro)
    @test round(res[x], digits=2) == 0.23 &&
          round(res[y], digits=2) == 0.06
end

@testset "NLProblem: x , y , z" begin
    @variables x, y, z
    eqs = [
        x^2 + y^2 + z^2 ~ 1.0
        2 * x^2 + y^2 - 4 * z ~ 0
        3 * x^2 - 4 * y^2 + z^2 ~ 0
    ]
    vars = Dict(x => 2.0, y => 1.0, z => 1.0)
    pro = NLProblem(eqs, vars)
    res = solve(pro)
    @test round(res[x], digits=2) == 0.70 &&
          round(res[y], digits=2) == 0.63 &&
          round(res[z], digits=2) == 0.34
end

@testset "NLProblem: x , y , z" begin
    @variables x, y, z
    eqs = [
        x^2 + y^2 + z^2 ~ 1.0
        2 * x^2 + y^2 - 4 * z ~ 0
        3 * x^2 - 4 * y^2 + z^2 ~ 0
    ]
    vars = Dict(x => 2.0, y => 1.0, z => 1.0)
    pro = NLProblem(eqs, vars)
    res = solve(pro)
    @test round(res[x], digits=2) == 0.70 &&
          round(res[y], digits=2) == 0.63 &&
          round(res[z], digits=2) == 0.34
end

@testset "solve-LinearProblem" begin
    N = 20
    @variables x[N]
    eqs = Vector{Equation}([])
    for i in 1:N
        rand_num = rand(1:20, Vector{Int}, 3)
        s = sum(rand_num)
        if i == 1
            push!(eqs, sum([rand_num[j-i+1] * x[j] for j in i:i+2]) ~ s)
        elseif i == N
            push!(eqs, sum([rand_num[j-i+3] * x[j] for j in i-2:i]) ~ s)
        else
            push!(eqs, sum([rand_num[j-i+2] * x[j] for j in i-1:i+1]) ~ s)
        end
    end
    vars=Dict([x[i]=>0.0 for i in 1:N])
    pro = LinearProblem(eqs,vars)
    res = solve(pro)
    @test round.(collect(values(res)), digits=1) == [1.0 for i in 1:N]
end

@testset "solve-LinearProblem by NLProblem" begin
    N = 100
    @variables x[N]
    eqs = Vector{Equation}([])
    for i in 1:N
        rand_num = rand(1:20, Vector{Int}, 3)
        s = sum(rand_num)
        if i == 1
            push!(eqs, sum([rand_num[j-i+1] * x[j] for j in i:i+2]) ~ s)
        elseif i == N
            push!(eqs, sum([rand_num[j-i+3] * x[j] for j in i-2:i]) ~ s)
        else
            push!(eqs, sum([rand_num[j-i+2] * x[j] for j in i-1:i+1]) ~ s)
        end
    end
    vars=Dict([x[i]=>0.0 for i in 1:N])
    pro = NLProblem(eqs,vars)
    res = solve(pro)
    @test round.(collect(values(res)), digits=1) == [1.0 for i in 1:N]
end


text:
"""
    Plate(normal::AbstractVector{T}, width::T, [, origin::AbstractVector{T}=zeros()])

A plate defined by all the points ``\\mathbf x`` that satify ``|(\\mathbf x - \\mathbf o) \\cdot \\mathbf n| < w /2`` where ``\\mathbf n`` is the unit normal, ``\\mathbf o`` is the origin, and ``w`` is the width.
"""
struct Plate{T,Dim} <: Shape{T,Dim}
    normal::SVector{Dim,T} # a unit vector which is orthogonal to the plate
    width::T # the width
    origin::SVector{Dim,T}
end

# Constructors
Plate(normal::AbstractVector{T}, width::T, origin::AbstractVector{T} = zeros(T,length(normal))) where {T} = Plate{T,length(normal)}(normal ./ norm(normal), width, origin)

Plate(normal::NTuple{Dim,T}, width::T, origin::AbstractVector{T} = zeros(T,Dim)) where {T,Dim} = Plate{T,Dim}(normal ./ norm(normal), width, origin)

name(shape::Plate) = "Plate"

volume(shape::Plate) = Inf
outer_radius(hs::Plate{T}) where T = T(Inf)

# import Base.issubset
# function issubset(inner_rect::Plate{T}, outer_rect::Plate{T}) where T
#     all(topright(inner_rect) .<= topright(outer_rect)) &&
#     all(bottomleft(inner_rect) .>= bottomleft(outer_rect))
# end

import Base.in
function in(x::AbstractVector{T}, p::Plate)::Bool where T
    abs(dot(x - p.origin, p.normal)) < p.width / T(2)
end

import Base.(==)
function ==(h1::Plate{T}, h2::Plate{T}) where T
    h1.origin == h2.origin &&
    h1.width == h2.width &&
    h1.normal == h2.normal
end

import Base.isequal
function isequal(h1::Plate{T}, h2::Plate{T}) where T
    isequal(h1.origin, h2.origin) &&
    isequal(h1.width, h2.width) &&
    isequal(h1.normal, h2.normal)
end

function iscongruent(h1::Plate{T}, h2::Plate{T}) where T
    (h1.normal  == h2.normal) && (h1.width  == h2.width)
end

function congruent(h::Plate{T}, x) where T
    Plate(h.normal, h.width, x)
end

# function boundary_functions(h::Plate{T,2}, scale::T = 10.0) where T
#     surface_vec = [-h.normal[2], h.normal[1]]
#
#     function x(t)
#         check_boundary_coord_range(t)
#         h.origin[1] + surface_vec[1] * width/2 + t * scale * surface_vec[1]
#     end
#
#     function y(t)
#         check_boundary_coord_range(t)
#         h.origin[2] + t * scale * surface_vec[2]
#     end
#
#     return x, y
# end


text:
include("init.jl")

@testset "MessageFederate test_message_federate_initialize" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    state = h.helicsFederateGetState(mFed)
    @test state == 0
    h.helicsFederateEnterExecutingMode(mFed)

    state = h.helicsFederateGetState(mFed)
    @test state == 2

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate test_message_federate_endpoint_registration" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    epid1 = h.helicsFederateRegisterEndpoint(mFed, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "ep2", "random")

    h.helicsFederateEnterExecutingMode(mFed)

    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed)

    endpoint_name = h.helicsEndpointGetName(epid1)
    @test endpoint_name == "TestA Federate/ep1"

    endpoint_name = h.helicsEndpointGetName(epid2)
    @test endpoint_name == "ep2"

    endpoint_name = h.helicsEndpointGetType(epid1)
    @test endpoint_name == ""

    endpoint_name = h.helicsEndpointGetType(epid2)
    @test endpoint_name == "random"

    epid_b = h.helicsFederateGetEndpoint(mFed, "ep2")
    type = h.helicsEndpointGetType(epid_b)
    @test type == "random"

    epid_c = h.helicsFederateGetEndpointByIndex(mFed, 0)
    name = h.helicsEndpointGetName(epid_c)
    @test name == "TestA Federate/ep1"

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate test_message_federate_send" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    epid1 = h.helicsFederateRegisterEndpoint(mFed, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "ep2", "random")

    h.helicsFederateSetTimeProperty(mFed, h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    h.helicsFederateEnterExecutingMode(mFed)

    data = "random-data"

    h.helicsEndpointSendEventRaw(epid1, "ep2", data, 1.0)

    granted_time = h.helicsFederateRequestTime(mFed, 2.0)
    @test granted_time == 1.0

    res = h.helicsFederateHasMessage(mFed)
    @test res == true

    res = h.helicsEndpointHasMessage(epid1)
    @test res == false

    res = h.helicsEndpointHasMessage(epid2)
    @test res == true

    message = h.helicsEndpointGetMessageObject(epid2)

    @test h.helicsMessageGetString(message) == "random-data"
    @test h.helicsMessageGetRawDataSize(message) == 11
    @test h.helicsMessageGetOriginalDestination(message) == ""
    @test h.helicsMessageGetOriginalSource(message) == "TestA Federate/ep1"
    @test h.helicsMessageGetSource(message) == "TestA Federate/ep1"
    @test h.helicsMessageGetTime(message) == 1.0

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate send_receive_2fed_multisend" begin

    broker = createBroker(2)
    mFed1, fedinfo1 = createMessageFederate(1, "A Federate")
    mFed2, fedinfo2 = createMessageFederate(1, "B Federate")

    epid1 = h.helicsFederateRegisterEndpoint(mFed1, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed2, "ep2", "random")

    h.helicsFederateSetTimeProperty(mFed1, h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    h.helicsFederateSetTimeProperty(mFed2, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    h.helicsFederateEnterExecutingModeAsync(mFed1)
    h.helicsFederateEnterExecutingMode(mFed2)
    h.helicsFederateEnterExecutingModeComplete(mFed1)

    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed1)
    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed2)

    h.helicsEndpointSetDefaultDestination(epid1, "ep2");

    h.helicsEndpointSendMessageRaw(epid1, "", "a")
    h.helicsEndpointSendMessageRaw(epid1, "", "a")
    h.helicsEndpointSendMessageRaw(epid1, "", "a")

    h.helicsFederateRequestTimeAsync(mFed1, 1.0)
    granted_time = h.helicsFederateRequestTime(mFed2, 1.0)
    complete_time = h.helicsFederateRequestTimeComplete(mFed1)

    @test granted_time == 1.0
    @test complete_time == 1.0

    res = h.helicsEndpointPendingMessages(epid2)
    @test res == 3

    res = h.helicsFederatePendingMessages(mFed2)
    @test res == 3

    @test h.helicsEndpointGetDefaultDestination(epid1) == "ep2"

    destroyFederate(mFed1, fedinfo1)
    destroyFederate(mFed2, fedinfo2)
    destroyBroker(broker)

end

@testset "MessageFederate message_object_tests" begin

    broker = createBroker(1)
    mFed1, fedinfo1 = createMessageFederate(1, "test")

    epid1 = h.helicsFederateRegisterEndpoint(mFed1, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed1, "ep2", "random");

    h.helicsFederateSetTimeProperty(mFed1, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    h.helicsFederateEnterExecutingMode(mFed1)

    @test h.helicsFederateGetState(mFed1) == h.HELICS_STATE_EXECUTION

    msg = h.helicsFederateCreateMessageObject(mFed1)
    h.helicsMessageSetDestination(msg, "ep2")
    h.helicsMessageGetDestination(msg) == "ep2"
    h.helicsMessageSetData(msg, repeat("a", 500))
    h.helicsMessageSetTime(msg, 0.0)

    h.helicsEndpointSendMessageObject(epid1, msg)
    time = h.helicsFederateRequestTime(mFed1, 1.0)
    @test time == 1.0

    @test h.helicsFederateHasMessage(mFed1) == true
    @test h.helicsEndpointHasMessage(epid1) == false
    @test h.helicsEndpointHasMessage(epid2) == true

    msg = h.helicsEndpointGetMessageObject(epid2)
    @test h.helicsMessageGetRawDataSize(msg) == 500

    rdata = h.helicsMessageGetRawDataPointer(msg)
    @test Char(unsafe_load(Ptr{Cchar}(rdata), 245)) == 'a'

    h.helicsFederateFinalize(mFed1)

    @test h.helicsFederateGetState(mFed1) == h.HELICS_STATE_FINALIZE

    h.helicsMessageSetFlagOption(msg, 7, true);
    @test h.helicsMessageCheckFlag(msg, 7) == true
    h.helicsMessageClearFlags(msg);
    @test h.helicsMessageCheckFlag(msg, 7) == false

    h.helicsEndpointSetDefaultDestination(epid1, "ep2")

    destroyFederate(mFed1, fedinfo1)
    destroyBroker(broker)

end


text:
module IPR

import ..ED
using ..SimLib
using ..SimLib: FArray
using SpinSymmetry: basissize
using SharedArrays: sdata

import Statistics

export ipr, inverse_participation_ratio, IPRData, IPRDataDescriptor, load_ipr, InverseParticipationRatio

## Data structure

"""
    struct IPRDataDescriptor <: ED.EDDerivedDataDescriptor

See [`EDDataDescriptor`](!ref) and [`EDDerivedDataDescriptor`](!ref).
"""
struct IPRDataDescriptor <: ED.EDDerivedDataDescriptor
    derivedfrom::ED.EDDataDescriptor
end

IPRDataDescriptor(args...; kwargs...) = IPRDataDescriptor(EDDataDescriptor(args...; kwargs...))

"""
    struct IPRData

Stores the actual inverse participation ratios (IPR) specified by the descriptor [`LSRDataDescriptor`](@ref).
Inverse participation ratio is defined as IPR(|ψ⟩) = 1/(∑ᵢ |⟨i|ψ⟩|⁴ ) relative to a basis |i⟩.
Usually |i⟩ is just the z-basis and then the IPR tells a story of localization.
IPR = 2^N for maximally entangled states
IPR = 1 for maximally localized states

The first dimension hold the participation ratio of the i-th eigenstate IPR(|ϕᵢ⟩).

The default save directory is "ipr".

`Statistics.mean` and `Statistics.std` are overloaded to act on the first dimension to conveniently compute
mean LSR and its variance.
"""
struct IPRData{N} <: SimLib.AbstractSimpleData
    descriptor::IPRDataDescriptor
    data::FArray{N}
end

IPRData(iprdd::IPRDataDescriptor) = IPRData(iprdd, FArray{4}(undef, ED.ed_size(iprdd), iprdd.shots, length(iprdd.fields), length(iprdd.ρs)))

## Saving/Loading

ED._default_folder(::IPRDataDescriptor) = "ipr"

"""
    load_ipr(edd)
    load_ipr(model[, diagtype][, location])
"""
load_ipr(args...; kwargs...) = load(IPRDataDescriptor(args...; kwargs...))

## Functions

ipr(mat) = ipr!(Vector{Float64}(undef, size(mat,1)), mat)

const inverse_participation_ratio = ipr

function ipr!(v, mat)
    I, J = size(mat)
    J == length(v) || throw(ArgumentError("Sizes not compatible: $(size(v)) and $(size(mat))"))
    @inbounds for j in 1:J
        tmp = zero(Float64)
        for i in 1:I
            tmp += abs2(mat[i,j])^2
        end
        v[j] = 1/tmp
    end
    v
end

### Task

mutable struct IPRTask <: ED.EDTask
    data
end

InverseParticipationRatio() = IPRTask(nothing)

function ED.initialize!(task::IPRTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, spectral_size)
end

function ED.compute_task!(task::IPRTask, evals, evecs, inds...)
    n = min(size(task.data,1), size(evecs, 2))
    @views ipr!(task.data[1:n, inds...], evecs[:, 1:n])
end

function ED.failed_task!(task::IPRTask, inds...)
    task.data[:, inds...] .= NaN64
end

function ED.assemble(task::IPRTask, edd)
    IPRData(IPRDataDescriptor(edd), sdata(task.data))
end

end #module


text:
eltype_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where N) where T = T
eltype_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where {N, T}) = or

ndims_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where T) where N = N
ndims_or(::Type{G}, or) where G <: (AbstractGeometry{N, T} where {N, T}) = or

Base.eltype(T::Type{<:AbstractGeometry}) = eltype_or(T, Any)

function Base.ndims(T::Type{<:AbstractGeometry})
    ndims_or(T, Any)
end

Base.eltype(x::AbstractGeometry{N, T}) where {N, T} = T
Base.ndims(x::AbstractGeometry{N, T}) where {N, T} = N

Base.eltype(::Type{AFG{T}}) where {T} = T
Base.eltype(::Type{FG}) where {FG <: AFG} = eltype(supertype(FG))


text:
# this object represents a slice of the hamiltonian at a single position.
# you can then do things like transfer_left(vector,ham[5],st.AL[5],st.AR[5])

struct MPOHamSlice{O<:MPOHamiltonian}
    ham::O
    i::Int
end

function Base.getproperty(x::MPOHamSlice,s::Symbol)
    if s == :imspaces
        return x.ham.imspaces[x.i,:]
    elseif s == :domspaces
        return x.ham.domspaces[x.i,:]
    elseif s == :odim
        return x.ham.odim
    elseif s == :Os
        return x.ham.Os[x.i,:,:]
    else
        return getfield(x,s)
    end
end

Base.getindex(x::MPOHamiltonian,a::AbstractVector{Int}) = [x[i] for i in a];
Base.getindex(x::MPOHamiltonian,a::Int) = MPOHamSlice(x,a);

Base.getindex(x::MPOHamSlice,args...) = Base.getindex(x.ham,x.i,args...);
isscal(x::MPOHamSlice,j,k) = isscal(x.ham,x.i,j,k);
Base.keys(x::MPOHamSlice) = keys(x.ham,x.i);
scalkeys(x::MPOHamSlice) = scalkeys(x.ham,x.i);
opkeys(x::MPOHamSlice) = opkeys(x.ham,x.i);
Base.contains(x::MPOHamSlice,j,k) = contains(x.ham,x.i,j,k);

#the usual mpoham transfer
function transfer_left(vec::Vector{V},ham::MPOHamSlice,A::V,Ab::V=A) where V<:MPSTensor
    toreturn = [TensorMap(zeros,eltype(A),_lastspace(Ab)'*ham.imspaces[i],_lastspace(A)') for i in 1:ham.odim]::Vector{V}

    for (j,k) in keys(ham)
        if isscal(ham,j,k)
            toreturn[k]+=ham.Os[j,k]*transfer_left(vec[j],A,Ab)
        else
            v = transfer_left(vec[j],ham[j,k],A,Ab)
            toreturn[k]+=transfer_left(vec[j],ham[j,k],A,Ab)
        end
    end

    return toreturn
end
function transfer_right(vec::Vector{V},ham::MPOHamSlice,A::V,Ab::V=A) where V<:MPSTensor
    toreturn = [TensorMap(zeros,eltype(A),_firstspace(A)*ham.domspaces[i],_firstspace(Ab)) for i in 1:ham.odim]::Vector{V}

    for (j,k) in keys(ham)
        if isscal(ham,j,k)
            toreturn[j]+=ham.Os[j,k]*transfer_right(vec[k],A,Ab)
        else
            toreturn[j]+=transfer_right(vec[k],ham[j,k],A,Ab)
        end
    end

    return toreturn
end


text:
""
function constraint_power_balance_shunt_dispatch(pm::_PM.AbstractWRModels, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    w   = var(pm, n, :w, i)
    p    = get(var(pm, n),    :p, Dict()); _PM._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); _PM._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); _PM._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); _PM._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); _PM._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); _PM._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); _PM._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); _PM._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); _PM._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); _PM._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")
    wbs = get(var(pm, n), :wbs, Dict()); _PM._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*w)
    cstr_q = JuMP.@constraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*w + sum(wbs[s] for s in bus_shunts_var))

    for s in bus_shunts_var
        _IM.relaxation_product(pm.model, w, bs[s], wbs[s])
    end

    if _IM.report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = cstr_q
    end
end

""
function constraint_power_balance_shunt_dispatch_soft(pm::_PM.AbstractWRModels, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    w = var(pm, n, :w, i)
    p_delta_abs = var(pm, n, :p_delta_abs, i)
    q_delta_abs = var(pm, n, :q_delta_abs, i)

    p = var(pm, n, :p)
    q = var(pm, n, :q)
    pg = var(pm, n, :pg)
    qg = var(pm, n, :qg)
    bs = var(pm, n, :bs)
    wbs = var(pm, n, :wbs)
    p_dc = var(pm, n, :p_dc)
    q_dc = var(pm, n, :q_dc)

    #p_delta = - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*w
    #q_delta = - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*w + sum(wbs[s] for s in bus_shunts_var)

    @constraint(pm.model,  p_delta_abs >= - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*w)
    @constraint(pm.model, -p_delta_abs <= - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*w)

    @constraint(pm.model,  q_delta_abs >= - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*w + sum(wbs[s] for s in bus_shunts_var))
    @constraint(pm.model, -q_delta_abs <= - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*w + sum(wbs[s] for s in bus_shunts_var))

    for s in bus_shunts_var
        _IM.relaxation_product(pm.model, w, bs[s], wbs[s])
    end
end

function constraint_ohms_yt_from_goc(pm::_PM.AbstractWRModels, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    w_fr = var(pm, n, :w, f_bus)
    wr   = var(pm, n, :wr, (f_bus, t_bus))
    wi   = var(pm, n, :wi, (f_bus, t_bus))

    JuMP.@constraint(pm.model, p_fr ==  (g/tm^2+g_fr)*w_fr + (-g*tr+b*ti)/tm^2*wr + (-b*tr-g*ti)/tm^2*wi )
    JuMP.@constraint(pm.model, q_fr == -(b/tm^2+b_fr)*w_fr - (-b*tr-g*ti)/tm^2*wr + (-g*tr+b*ti)/tm^2*wi )
end


text:
"""
    tfm = linmap(mat, img, initial_tfm=IdentityTransformation())

Construct a linear transformation from the values in `mat`.
`mat` can be any list of `N^2` values, where `N` is the dimensionality.
`img` is the array (or array of the same dimensionality) that `tfm` will be applied to.

If `initial_tfm` is supplied, it will be composed with the one from `mat`.

See also [`RegisterQD.aff`](@ref), which includes a translational component.
"""
function linmap(mat, img::AbstractArray{T,N}, initial_tfm=IdentityTransformation()) where {T,N}
    # img isn't used *except* to get the dimensionality in an inferrable way, so that
    # the (static) size of the array in `tfm` is known.
    mat = [mat...]
    mat = reshape(mat, N,N)
    lm = LinearMap(SMatrix{N,N}(mat))
    # The order below is subtle: you might think that `initial_tfm` should be applied
    # to `x` first, so that the order of composition should be `lm ∘ initial_tfm`.
    # But what we're going to be doing is calculating an `lm` that aligns
    #    imgw(x) = moving(initial_tfm(x))
    # to `fixed`. The final aligned image is
    #    imgw2(x) = imgw(lm(x))
    # Putting these two together, we see that
    #    imgw2(x) = moving(initial_tfm(lm(x)))
    # and thus the resulting transformation is
    #    initial_tfm ∘ lm
    return initial_tfm ∘ lm
end

"""
    tfm = aff(params, img, initial_tfm=IdentityTransformation())

Construct an affine transformation from the values in `params`.
`params` can be any list of `N + N^2` values, where `N` is the dimensionality;
the first `N` are for the translation,
and the last `N^2` are for the linear portion (see [`RegisterQD.linmap`](@ref)).

`img` is the array (or array of the same dimensionality) that `tfm` will be applied to.

If `initial_tfm` is supplied, it will be composed with the one from `params`.
"""
function aff(params, img::AbstractArray{T,N}, initial_tfm=IdentityTransformation()) where {T,N}
    # img isn't used *except* to get the dimensionality in an inferrable way, so that
    # the (static) sizes of the arrays in `tfm` are known.
    params = [params...]
    length(params) == (N+N^2) || throw(DimensionMismatch("expected $(N+N^2) parameters, got $(length(params))"))
    offs = Float64.(params[1:N])
    mat = Float64.(params[(N+1):end])
    mat = reshape(mat,N,N)
    # The order below is explained in `linmap`
    return initial_tfm ∘ AffineMap(SMatrix{N,N}(mat), SVector{N}(offs))
end

# This acts as the objective function during optimization, see qd_affine_coarse
# which defines `f(x)` in terms of this function.
# here params contains parameters of a linear map
# The "fast" part means that it handles translation using `best_shift`, so this is
# concerned only with the linear-map portion of the affine transformation.
function affine_mm_fast(params, mxshift, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation())
    tfm = arrayscale(linmap(params, moving, initial_tfm), SD)
    moving, fixed = warp_and_intersect(moving, fixed, tfm)
    bshft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity)
    return mm
end

# Similar to the above but includes the translation in `params`
# This is used for optimization of the complete affine transformation.
function affine_mm_slow(params, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation())
    tfm = arrayscale(aff(params, moving, initial_tfm), SD)
    moving, fixed = warp_and_intersect(moving, fixed, tfm)
    mm = mismatch0(fixed, moving; normalization=:intensity)
    return ratio(mm, thresh, Inf)
end

"""
    tfm, mm = qd_affine_coarse(fixed, moving, mxshift, linmins, linmaxs; kwargs...)

Compute an affine transformation `tfm` that coarsely minimizes the mismatch between
`fixed(x)` and `moving(tfm(x))`. `mm` is the total mismatch.
The translational component is handled at the level of integer-pixel shifts,
which is why this is "coarse" optimization.
"""
function qd_affine_coarse(fixed, moving, mxshift, linmins, linmaxs;
                          SD=I,
                          initial_tfm=IdentityTransformation(),
                          thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])),
                          minwidth=default_lin_minwidths(moving),
                          maxevals=5e4,
                          kwargs...)
    # Define the objective function for the coarse stage
    f(x) = affine_mm_fast(x, mxshift, fixed, moving, thresh, SD; initial_tfm=initial_tfm)
    # QuadDIRECT benefits from bounds on the search region, build those bounds
    upper = linmaxs
    lower = linmins
    # Try to find the global minimum (within the needed precision)
    root, x0 = _analyze(f, lower, upper;
                        minwidth=minwidth, print_interval=100, maxevals=maxevals, kwargs..., atol=0, rtol=1e-3)
    box = minimum(root)
    # Convert the coarse-minimum into a complete affine transformation
    params = position(box, x0)
    tfmcoarse0 = linmap(params, moving, initial_tfm)
    best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=arrayscale(tfmcoarse0, SD))
    tfmcoarse = tfmcoarse0 ∘ pscale(Translation(best_shft), SD)
    return tfmcoarse, mm
end

"""
    tfm, mm = qd_affine_fine(fixed, moving, linmins, linmaxs; initial_tfm=IdentityTransformation(), kwargs...)

Compute an affine transformation `tfm` that minimizes the mismatch between
`fixed(x)` and `moving(tfm(x))`. `mm` is the total mismatch.
This only allows small translations (just a couple of pixels).
As a consequence, any large translations must be supplied with reasonable accuracy in
`initial_tfm` (see [`RegisterQD.qd_affine_coarse`](@ref)).
"""
function qd_affine_fine(fixed, moving, linmins, linmaxs;
                        SD=I,
                        initial_tfm=IdentityTransformation(),
                        thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])),
                        minwidth_mat=default_lin_minwidths(fixed)./10,
                        maxevals=5e4,
                        kwargs...)
    # Define the objective function
    f(x) = affine_mm_slow(x, fixed, moving, thresh, SD; initial_tfm=initial_tfm)
    # Limit translations to only a couple of pixels
    upper_shft = fill(2.0, ndims(fixed))
    # Combine translation limits with the affine-transform limits
    upper = vcat(upper_shft, linmaxs)
    lower = vcat(-upper_shft, linmins)
    # Don't worry about translations that less than 1% of a pixel
    minwidth_shfts = fill(0.01, ndims(fixed))
    minwidth = vcat(minwidth_shfts, minwidth_mat)
    # Try to find the global minimum
    root, x0 = _analyze(f, lower, upper;
                        minwidth=minwidth, print_interval=100, maxevals=maxevals, kwargs...)
    box = minimum(root)
    # Convert the result to an affine transformation
    params = position(box, x0)
    tfmfine = aff(params, moving, initial_tfm)
    return tfmfine, value(box)
end

#TODO oops maybe dmax and ndax is too closely related to linmax and linmin?
# Response: dmax and ndmax are designed to allow you to control linmax and linmin.
# You supply one or the other, so I don't see a problem.
#TODO I think that this is a tad loquatious, and not enough examples. Permission to rework it?
"""
`tform, mm = qd_affine(fixed, moving, mxshift, linmins, linmaxs, SD=I; presmoothed=false, thresh, initial_tfm, kwargs...)`
`tform, mm = qd_affine(fixed, moving, mxshift, SD=I; presmoothed=false, thresh, initial_tfm, kwargs...)`
optimizes an affine transformation (linear map + translation) to minimize the mismatch between `fixed` and
`moving` using the QuadDIRECT algorithm.  The algorithm is run twice: the first step samples the search space
at a coarser resolution than the second.  `kwargs...` may contain any keyword argument that can be passed to
`QuadDIRECT.analyze`. It's recommended that you pass your own stopping criteria when possible
(i.e. `rtol`, `atol`, and/or `fvalue`).  If you provide `rtol` and/or `atol` they will apply only to the
second (fine) step of the registration; the user may not adjust these criteria for the coarse step.

`tform` will be centered on the origin-of-coordinates, i.e. (0,0) for a 2D image.  Usually it is more natural to consider rotations
around the center of the image.  If you would like `mxrot` and the returned rotation to act relative to the center of the image, then you must
move the origin to the center of the image by calling `centered(img)` from the `ImageTransformations` package.  Call `centered` on both the
fixed and moving image to generate the `fixed` and `moving` that you provide as arguments.  If you later want to apply the returned transform
to an image you must remember to call `centered` on that image as well.  Alternatively you can re-encode the transformation in terms of a
different origin by calling `recenter(tform, newctr)` where `newctr` is the displacement of the new center from the old center.

The `linmins` and `linmaxs` arguments set the minimum and maximum allowable values in the linear map matrix.
They can be supplied as NxN matrices or flattened vectors.  If omitted then a modest default search space is chosen.
`mxshift` sets the magnitude of the largest allowable translation in each dimension (It's a vector of length N).
This default search-space allows for very little rotation.
Alternatively, you can submit `dmax` or `ndmax` values as keyword functions, which will use diagonal or non-diagonal variation from the identity matrix
to generate less modest `linmins` and `linmaxs` arguments for you.

Use `presmoothed=true` if you have called [`qsmooth`](@ref) on `fixed` before calling `qd_affine`.
Do not smooth `moving`.

`kwargs...` can also include any other keyword argument that can be passed to `QuadDIRECT.analyze`.
It's recommended that you pass your own stopping criteria when possible (i.e. `rtol`, `atol`, and/or `fvalue`).

If you have a good initial guess at the solution, pass it with the `initial_tfm` kwarg to jump-start the search.

Use `SD` if your axes are not uniformly sampled, for example `SD = diagm(voxelspacing)` where `voxelspacing`
is a vector encoding the spacing along all axes of the image. `thresh` enforces a certain amount of sum-of-squared-intensity
overlap between the two images; with non-zero `thresh`, it is not permissible to "align" the images by shifting one entirely out of the way of the other.
"""
function qd_affine(fixed, moving, mxshift, linmins, linmaxs;
                   presmoothed=false,
                   SD=I,
                   thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])),
                   initial_tfm=IdentityTransformation(),
                   print_interval=100,
                   kwargs...)
    fixed, moving = float(fixed), float(moving)
    if presmoothed
        moving = qinterp(eltype(fixed), moving)
    end
    linmins = [linmins...]
    linmaxs = [linmaxs...]
    print_interval < typemax(Int) && print("Running coarse step\n")
    mw = default_lin_minwidths(moving)
    tfm_coarse, mm_coarse = qd_affine_coarse(fixed, moving, mxshift, linmins, linmaxs;
                                             SD = SD, minwidth=mw, initial_tfm=initial_tfm, thresh=thresh, print_interval=print_interval, kwargs...)
    print_interval < typemax(Int) && print("Running fine step\n")
    mw = mw./100
    linmins, linmaxs = scalebounds(linmins, linmaxs, 0.5) # should these be narrowed further?
    final_tfm, final_mm = qd_affine_fine(fixed, moving, linmins, linmaxs;
                                         SD = SD, minwidth_mat=mw, initial_tfm=tfm_coarse, thresh=thresh, print_interval=print_interval, kwargs...)
    return final_tfm, final_mm
end

function qd_affine(fixed, moving, mxshift;
                   dmax = 0.05, ndmax = 0.05,
                   kwargs...)
    minb, maxb = default_linmap_bounds(fixed; dmax = dmax, ndmax = ndmax)
    return qd_affine(fixed, moving, mxshift, minb, maxb; kwargs...)
end


text:
## Analyzing lowered code

function add_docexpr!(docexprs::AbstractDict{Module,V}, mod::Module, ex) where V
    docexs = get(docexprs, mod, nothing)
    if docexs === nothing
        docexs = docexprs[mod] = V()
    end
    push!(docexs, ex)
    return docexprs
end

function lookup_callexpr(frame, stmt)
    fargs = JuliaInterpreter.collect_args(frame, stmt)
    return Expr(:call, fargs...)
end

function assign_this!(frame, value)
    frame.framedata.ssavalues[frame.pc] = value
end

# This defines the API needed to store signatures using methods_by_execution!
# This default version is simple and only used for testing purposes.
# The "real" one is CodeTrackingMethodInfo in Revise.jl.
const MethodInfo = IdDict{Type,LineNumberNode}
add_signature!(methodinfo::MethodInfo, @nospecialize(sig), ln) = push!(methodinfo, sig=>ln)
push_expr!(methodinfo::MethodInfo, mod::Module, ex::Expr) = methodinfo
pop_expr!(methodinfo::MethodInfo) = methodinfo
add_dependencies!(methodinfo::MethodInfo, be::CodeEdges, src, isrequired) = methodinfo
add_includes!(methodinfo::MethodInfo, mod::Module, filename) = methodinfo

# This is not generally used, see `is_method_or_eval` instead
function hastrackedexpr(stmt; heads=LoweredCodeUtils.trackedheads)
    haseval = false
    if isa(stmt, Expr)
        haseval = matches_eval(stmt)
        if stmt.head === :call
            f = stmt.args[1]
            callee_matches(f, Core, :_typebody!) && return true, haseval
            callee_matches(f, Core, :_setsuper!) && return true, haseval
            f === :include && return true, haseval
        elseif stmt.head === :thunk
            any(s->any(hastrackedexpr(s; heads=heads)), (stmt.args[1]::Core.CodeInfo).code) && return true, haseval
        elseif stmt.head ∈ heads
            return true, haseval
        end
    end
    return false, haseval
end

function matches_eval(stmt::Expr)
    stmt.head === :call || return false
    f = stmt.args[1]
    return f === :eval ||
           (callee_matches(f, Base, :getproperty) && is_quotenode_egal(stmt.args[end], :eval)) ||
           (isa(f, GlobalRef) && f.name === :eval) || is_quotenode_egal(f, Core.eval)
end

function categorize_stmt(@nospecialize(stmt))
    ismeth, haseval, isinclude, isnamespace, istoplevel = false, false, false, false, false
    if isa(stmt, Expr)
        haseval = matches_eval(stmt)
        ismeth = stmt.head === :method
        istoplevel = stmt.head === :toplevel
        isnamespace = stmt.head === :export || stmt.head === :import || stmt.head === :using
        isinclude = stmt.head === :call && stmt.args[1] === :include
    end
    return ismeth, haseval, isinclude, isnamespace, istoplevel
end

"""
    isrequired, evalassign = minimal_evaluation!([predicate,] methodinfo, src::Core.CodeInfo, mode::Symbol)

Mark required statements in `src`: `isrequired[i]` is `true` if `src.code[i]` should be evaluated.
Statements are analyzed by `isreq, haseval = predicate(stmt)`, and `predicate` defaults
to `Revise.is_method_or_eval`.
`haseval` is true if the statement came from `@eval` or `eval(...)` call.
Since the contents of such expression are difficult to analyze, it is generally
safest to execute all such evals.
"""
function minimal_evaluation!(@nospecialize(predicate), methodinfo, src::Core.CodeInfo, mode::Symbol)
    edges = CodeEdges(src)
    # LoweredCodeUtils.print_with_code(stdout, src, edges)
    isrequired = fill(false, length(src.code))
    evalassign = false
    for (i, stmt) in enumerate(src.code)
        if !isrequired[i]
            isrequired[i], haseval = predicate(stmt)
            if haseval                              # line `i` may be the equivalent of `f = Core.eval`, so...
                isrequired[edges.succs[i]] .= true  # ...require each stmt that calls `eval` via `f(expr)`
                isrequired[i] = true
            end
        end
        if mode === :evalassign && isexpr(stmt, :(=))
            evalassign = isrequired[i] = true
            lhs = (stmt::Expr).args[1]
            if isa(lhs, Symbol)
                isrequired[edges.byname[lhs].succs] .= true  # mark any `const` statements or other "uses" in this block
            end
        end
    end
    # Check for docstrings
    if length(src.code) > 1 && mode !== :sigs
        stmt = src.code[end-1]
        if isexpr(stmt, :call) && (stmt::Expr).args[1] === Base.Docs.doc!
            isrequired[end-1] = true
        end
    end
    # All tracked expressions are marked. Now add their dependencies.
    # LoweredCodeUtils.print_with_code(stdout, src, isrequired)
    lines_required!(isrequired, src, edges;
                    norequire=mode===:sigs ? LoweredCodeUtils.exclude_named_typedefs(src, edges) : ())
    # LoweredCodeUtils.print_with_code(stdout, src, isrequired)
    add_dependencies!(methodinfo, edges, src, isrequired)
    return isrequired, evalassign
end
minimal_evaluation!(@nospecialize(predicate), methodinfo, frame::JuliaInterpreter.Frame, mode::Symbol) =
    minimal_evaluation!(predicate, methodinfo, frame.framecode.src, mode)

function minimal_evaluation!(methodinfo, frame, mode::Symbol)
    minimal_evaluation!(methodinfo, frame, mode) do @nospecialize(stmt)
        ismeth, haseval, isinclude, isnamespace, istoplevel = categorize_stmt(stmt)
        isreq = ismeth | isinclude | istoplevel
        return mode === :sigs ? (isreq, haseval) : (isreq | isnamespace, haseval)
    end
end

function methods_by_execution(mod::Module, ex::Expr; kwargs...)
    methodinfo = MethodInfo()
    docexprs = DocExprs()
    value, frame = methods_by_execution!(JuliaInterpreter.Compiled(), methodinfo, docexprs, mod, ex; kwargs...)
    return methodinfo, docexprs, frame
end

"""
    methods_by_execution!(recurse=JuliaInterpreter.Compiled(), methodinfo, docexprs, mod::Module, ex::Expr;
                          mode=:eval, disablebp=true, skip_include=mode!==:eval, always_rethrow=false)

Evaluate or analyze `ex` in the context of `mod`.
Depending on the setting of `mode` (see the Extended help), it supports full evaluation or just the minimal
evaluation needed to extract method signatures.
`recurse` controls JuliaInterpreter's evaluation of any non-intercepted statement;
likely choices are `JuliaInterpreter.Compiled()` or `JuliaInterpreter.finish_and_return!`.
`methodinfo` is a cache for storing information about any method definitions (see [`CodeTrackingMethodInfo`](@ref)).
`docexprs` is a cache for storing documentation expressions; obtain an empty one with `Revise.DocExprs()`.

# Extended help

The action depends on `mode`:

- `:eval` evaluates the expression in `mod`, similar to `Core.eval(mod, ex)` except that `methodinfo` and `docexprs`
  will be populated with information about any signatures or docstrings. This mode is used to implement `includet`.
- `:sigs` analyzes `ex` and extracts signatures of methods and docstrings (specifically, statements flagged by
  [`Revise.minimal_evaluation!`](@ref)), but does not evaluate `ex` in the traditional sense.
  It will selectively execute statements needed to form the signatures of defined methods.
  It will also expand any `@eval`ed expressions, since these might contain method definitions.
- `:evalmeth` analyzes `ex` and extracts signatures and docstrings like `:sigs`, but takes the additional step of
  evaluating any `:method` statements.
- `:evalassign` acts similarly to `:evalmeth`, and also evaluates assignment statements.

When selectively evaluating an expression, Revise will incorporate required dependencies, even for
minimal-evaluation modes like `:sigs`. For example, the method definition

    max_values(T::Union{map(X -> Type{X}, Base.BitIntegerSmall_types)...}) = 1 << (8*sizeof(T))

found in `base/abstractset.jl` requires that it create the anonymous function in order to compute the
signature.

The other keyword arguments are more straightforward:

- `disablebp` controls whether JuliaInterpreter's breakpoints are disabled before stepping through the code.
  They are restored on exit.
- `skip_include` prevents execution of `include` statements, instead inserting them into `methodinfo`'s
  cache. This defaults to `true` unless `mode` is `:eval`.
- `always_rethrow`, if true, causes an error to be thrown if evaluating `ex` triggered an error.
  If false, the error is logged with `@error`. `InterruptException`s are always rethrown.
  This is primarily useful for debugging.
"""
function methods_by_execution!(@nospecialize(recurse), methodinfo, docexprs, mod::Module, ex::Expr;
                               mode::Symbol=:eval, disablebp::Bool=true, always_rethrow::Bool=false, kwargs...)
    mode ∈ (:sigs, :eval, :evalmeth, :evalassign) || error("unsupported mode ", mode)
    lwr = Meta.lower(mod, ex)
    isa(lwr, Expr) || return nothing, nothing
    if lwr.head === :error || lwr.head === :incomplete
        error("lowering returned an error, ", lwr)
    end
    if lwr.head !== :thunk
        mode === :sigs && return nothing, nothing
        return Core.eval(mod, lwr), nothing
    end
    frame = JuliaInterpreter.Frame(mod, lwr.args[1])
    mode === :eval || LoweredCodeUtils.rename_framemethods!(recurse, frame)
    # Determine whether we need interpreted mode
    isrequired, evalassign = minimal_evaluation!(methodinfo, frame, mode)
    # LoweredCodeUtils.print_with_code(stdout, frame.framecode.src, isrequired)
    if !any(isrequired) && (mode===:eval || !evalassign)
        # We can evaluate the entire expression in compiled mode
        if mode===:eval
            ret = try
                Core.eval(mod, ex)
            catch err
                (always_rethrow || isa(err, InterruptException)) && rethrow(err)
                loc = location_string(whereis(frame)...)
                bt = trim_toplevel!(catch_backtrace())
                throw(ReviseEvalException(loc, err, Any[(sf, 1) for sf in stacktrace(bt)]))
            end
        else
            ret = nothing
        end
    else
        # Use the interpreter
        local active_bp_refs
        if disablebp
            # We have to turn off all active breakpoints, https://github.com/timholy/CodeTracking.jl/issues/27
            bp_refs = JuliaInterpreter.BreakpointRef[]
            for bp in JuliaInterpreter.breakpoints()
                append!(bp_refs, bp.instances)
            end
            active_bp_refs = filter(bp->bp[].isactive, bp_refs)
            foreach(disable, active_bp_refs)
        end
        ret = try
            methods_by_execution!(recurse, methodinfo, docexprs, frame, isrequired; mode=mode, kwargs...)
        catch err
            (always_rethrow || isa(err, InterruptException)) && (disablebp && foreach(enable, active_bp_refs); rethrow(err))
            loc = location_string(whereis(frame)...)
            sfs = []  # crafted for interaction with Base.show_backtrace
            frame = JuliaInterpreter.leaf(frame)
            while frame !== nothing
                push!(sfs, (Base.StackTraces.StackFrame(frame), 1))
                frame = frame.caller
            end
            throw(ReviseEvalException(loc, err, sfs))
        end
        if disablebp
            foreach(enable, active_bp_refs)
        end
    end
    return ret, lwr
end
methods_by_execution!(methodinfo, docexprs, mod::Module, ex::Expr; kwargs...) =
    methods_by_execution!(JuliaInterpreter.Compiled(), methodinfo, docexprs, mod, ex; kwargs...)

function methods_by_execution!(@nospecialize(recurse), methodinfo, docexprs, frame::Frame, isrequired::AbstractVector{Bool}; mode::Symbol=:eval, skip_include::Bool=true)
    isok(lnn::LineTypes) = !iszero(lnn.line) || lnn.file !== :none   # might fail either one, but accept anything

    mod = moduleof(frame)
    # Hoist this lookup for performance. Don't throw even when `mod` is a baremodule:
    modinclude = isdefined(mod, :include) ? getfield(mod, :include) : nothing
    signatures = []  # temporary for method signature storage
    pc = frame.pc
    while true
        JuliaInterpreter.is_leaf(frame) || (@warn("not a leaf"); break)
        stmt = pc_expr(frame, pc)
        if !isrequired[pc] && mode !== :eval && !(mode === :evalassign && isexpr(stmt, :(=)))
            pc = next_or_nothing!(frame)
            pc === nothing && break
            continue
        end
        if isa(stmt, Expr)
            head = stmt.head
            if head === :toplevel
                local value
                for ex in stmt.args
                    ex isa Expr || continue
                    value = methods_by_execution!(recurse, methodinfo, docexprs, mod, ex; mode=mode, disablebp=false, skip_include=skip_include)
                end
                isassign(frame, pc) && assign_this!(frame, value)
                pc = next_or_nothing!(frame)
            elseif head ∈ structheads
                if mode !== :sigs
                    pc = step_expr!(recurse, frame, stmt, true)  # This checks that they are unchanged
                else
                    pc = next_or_nothing!(frame)
                end
            # elseif head === :thunk && isanonymous_typedef(stmt.args[1])
            #     # Anonymous functions should just be defined anew, since there does not seem to be a practical
            #     # way to find them within the already-defined module.
            #     # They may be needed to define later signatures.
            #     # Note that named inner methods don't require special treatment.
            #     pc = step_expr!(recurse, frame, stmt, true)
            elseif head === :method
                empty!(signatures)
                ret = methoddef!(recurse, signatures, frame, stmt, pc; define=mode!==:sigs)
                if ret === nothing
                    # This was just `function foo end` or similar.
                    # However, it might have been followed by a thunk that defined a
                    # method (issue #435), so we still need to check for additions.
                    if !isempty(signatures)
                        file, line = whereis(frame.framecode, pc)
                        lnn = LineNumberNode(Int(line), Symbol(file))
                        for sig in signatures
                            add_signature!(methodinfo, sig, lnn)
                        end
                    end
                    pc = ret
                else
                    pc, pc3 = ret
                    # Get the line number from the body
                    stmt3 = pc_expr(frame, pc3)::Expr
                    lnn = nothing
                    if line_is_decl
                        sigcode = @lookup(frame, stmt3.args[2])::Core.SimpleVector
                        lnn = sigcode[end]
                        if !isa(lnn, LineNumberNode)
                            lnn = nothing
                        end
                    end
                    if lnn === nothing
                        bodycode = stmt3.args[end]
                        if !isa(bodycode, CodeInfo)
                            bodycode = @lookup(frame, bodycode)
                        end
                        if isa(bodycode, CodeInfo)
                            lnn = linetable(bodycode, 1)
                            if !isok(lnn)
                                lnn = nothing
                                if length(bodycode.code) > 1
                                    # This may be a kwarg method. Mimic LoweredCodeUtils.bodymethod,
                                    # except without having a method
                                    stmt = bodycode.code[end-1]
                                    if isa(stmt, Expr) && length(stmt.args) > 1
                                        stmt = stmt::Expr
                                        a = stmt.args[1]
                                        nargs = length(stmt.args)
                                        hasself = let stmt = stmt, slotnames::Vector{Symbol} = bodycode.slotnames
                                            any(i->LoweredCodeUtils.is_self_call(stmt, slotnames, i), 2:nargs)
                                        end
                                        if isa(a, Core.SlotNumber)
                                            a = bodycode.slotnames[a.id]
                                        end
                                        if hasself && (isa(a, Symbol) || isa(a, GlobalRef))
                                            thismod, thisname = isa(a, Symbol) ? (mod, a) : (a.mod, a.name)
                                            if isdefined(thismod, thisname)
                                                f = getfield(thismod, thisname)
                                                mths = methods(f)
                                                if length(mths) == 1
                                                    mth = first(mths)
                                                    lnn = LineNumberNode(Int(mth.line), mth.file)
                                                end
                                            end
                                        end
                                    end
                                end
                                if lnn === nothing
                                    # Just try to find *any* line number
                                    for lnntmp in linetable(bodycode)
                                        lnntmp = lnntmp::LineTypes
                                        if isok(lnntmp)
                                            lnn = lnntmp
                                            break
                                        end
                                    end
                                end
                            end
                        elseif isexpr(bodycode, :lambda)
                            bodycode = bodycode::Expr
                            lnntmp = bodycode.args[end][1]::LineTypes
                            if isok(lnntmp)
                                lnn = lnntmp
                            end
                        end
                    end
                    if lnn === nothing
                        i = codelocs(frame, pc3)
                        while i > 0
                            lnntmp = linetable(frame, i)
                            if isok(lnntmp)
                                lnn = lnntmp
                                break
                            end
                            i -= 1
                        end
                    end
                    if lnn !== nothing && isok(lnn)
                        for sig in signatures
                            add_signature!(methodinfo, sig, lnn)
                        end
                    end
                end
            elseif head === :(=)
                # If we're here, either isrequired[pc] is true, or the mode forces us to eval assignments
                pc = step_expr!(recurse, frame, stmt, true)
            elseif head === :call
                f = @lookup(frame, stmt.args[1])
                if f === Core.eval
                    # an @eval or eval block: this may contain method definitions, so intercept it.
                    evalmod = @lookup(frame, stmt.args[2])::Module
                    evalex = @lookup(frame, stmt.args[3])
                    value = nothing
                    for (newmod, newex) in ExprSplitter(evalmod, evalex)
                        if is_doc_expr(newex)
                            add_docexpr!(docexprs, newmod, newex)
                            newex = newex.args[4]
                        end
                        newex = unwrap(newex)
                        push_expr!(methodinfo, newmod, newex)
                        value = methods_by_execution!(recurse, methodinfo, docexprs, newmod, newex; mode=mode, skip_include=skip_include, disablebp=false)
                        pop_expr!(methodinfo)
                    end
                    assign_this!(frame, value)
                    pc = next_or_nothing!(frame)
                elseif skip_include && (f === modinclude || f === Base.include || f === Core.include)
                    # include calls need to be managed carefully from several standpoints, including
                    # path management and parsing new expressions
                    add_includes!(methodinfo, mod, @lookup(frame, stmt.args[2]))
                    assign_this!(frame, nothing)  # FIXME: the file might return something different from `nothing`
                    pc = next_or_nothing!(frame)
                elseif f === Base.Docs.doc! # && mode !== :eval
                    fargs = JuliaInterpreter.collect_args(frame, stmt)
                    popfirst!(fargs)
                    length(fargs) == 3 && push!(fargs, Union{})  # add the default sig
                    dmod::Module, b::Base.Docs.Binding, str::Base.Docs.DocStr, sig = fargs
                    if isdefined(b.mod, b.var)
                        tmpvar = getfield(b.mod, b.var)
                        if isa(tmpvar, Module)
                            dmod = tmpvar
                        end
                    end
                    # Workaround for julia#38819 on older Julia versions
                    if !isdefined(dmod, Base.Docs.META)
                        Base.Docs.initmeta(dmod)
                    end
                    m = get!(Base.Docs.meta(dmod), b, Base.Docs.MultiDoc())::Base.Docs.MultiDoc
                    if haskey(m.docs, sig)
                        currentstr = m.docs[sig]::Base.Docs.DocStr
                        redefine = currentstr.text != str.text
                    else
                        push!(m.order, sig)
                        redefine = true
                    end
                    # (Re)assign without the warning
                    if redefine
                        m.docs[sig] = str
                        str.data[:binding] = b
                        str.data[:typesig] = sig
                    end
                    assign_this!(frame, Base.Docs.doc(b, sig))
                    pc = next_or_nothing!(frame)
                else
                    # A :call Expr we don't want to intercept
                    pc = step_expr!(recurse, frame, stmt, true)
                end
            else
                # An Expr we don't want to intercept
                pc = step_expr!(recurse, frame, stmt, true)
            end
        else
            # A statement we don't want to intercept
            pc = step_expr!(recurse, frame, stmt, true)
        end
        pc === nothing && break
    end
    return isrequired[frame.pc] ? get_return(frame) : nothing
end


text:
module FieldArraysOperationsTests

using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using FillArrays

using Test
# using BenchmarkTools

np = 4
ncells = 10
d = 2
p = Point(1.0,2.0)
x = fill(p,np)
cell_to_x = Fill(x,ncells)

# ConstantField

cell_to_a = rand(ncells)
cell_to_f = lazy_map(ConstantField,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,map(i->fill(i,np),cell_to_a))

#print_op_tree(cell_to_fx)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# linear_combination -> Field

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = VectorValue{2,Float64}[(1,1),(4,2),(3,5),(1,2)]
f = linear_combination(b,a)

cell_to_a = Fill(a,ncells)
cell_to_b = fill(b,ncells)
cell_to_r = fill(f(x),ncells)

cell_to_f = lazy_map(linear_combination,cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
@test isa(cell_to_fx.maps.value,Fields.LinearCombinationMap)
test_array(cell_to_fx,cell_to_r)

cell_to_∇r = fill(∇(f)(x),ncells)
cell_to_∇f = lazy_map(∇,cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
@test isa(cell_to_∇fx.maps.value,Fields.LinearCombinationMap)
test_array(cell_to_∇fx,cell_to_∇r)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# linear_combination -> Field array

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = rand(length(avals),5)
f = linear_combination(b,a)

cell_to_a = Fill(a,ncells)
cell_to_b = fill(b,ncells)
cell_to_r = fill( evaluate(linear_combination(b,a),x),ncells)

cell_to_f = lazy_map(linear_combination,cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
@test isa(cell_to_fx.maps.value,Fields.LinearCombinationMap)
test_array(cell_to_fx,cell_to_r)

cell_to_∇r = fill(evaluate(Broadcasting(∇)(f),x),ncells)
cell_to_∇f = lazy_map(Broadcasting(∇),cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
@test isa(cell_to_∇fx.maps.value,Fields.LinearCombinationMap)
test_array(cell_to_∇fx,cell_to_∇r)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# transpose of an array of fields

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = rand(length(avals),5)
c = linear_combination(b,a)

cell_to_a = Fill(a,ncells)
cell_to_b = fill(b,ncells)
cell_to_c = lazy_map(linear_combination,cell_to_b,cell_to_a)

cell_to_f = lazy_map(transpose,cell_to_c)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
@test isa(cell_to_fx.maps.value,Fields.TransposeMap)
cell_to_r = fill(evaluate(transpose(c),x),ncells)
test_array(cell_to_fx,cell_to_r)

cell_to_∇f = lazy_map(Broadcasting(∇),cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
@test isa(cell_to_fx.maps.value,Fields.TransposeMap)
cell_to_∇r = fill(evaluate(transpose(∇.(c)),x),ncells)
test_array(cell_to_∇fx,cell_to_∇r)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# Composition between fields

mfun(g) = 3*g[1]
gfun(x) = 2*x
ffun(x) = mfun(gfun(x))

m = GenericField(mfun)
g = GenericField(gfun)
f = m∘g

cell_to_m = Fill(m,ncells)
cell_to_g = fill(g,ncells)
cell_to_f = lazy_map(∘,cell_to_m,cell_to_g)
cell_to_r = fill(f(x),ncells)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)

cell_to_∇f = lazy_map(∇,cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
cell_to_∇r = fill(evaluate(∇(f),x),ncells)
test_array(cell_to_∇fx,cell_to_∇r)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# Composition between field array and field

afun(x) = 2*x
bfun(x) = x[1] + 3
a = GenericField(afun)
b = fill(GenericField(bfun),4)
f = Broadcasting(∘)(b,a)

cell_to_a = fill(a,ncells)
cell_to_b = Fill(b,ncells)

cell_to_f = lazy_map(Broadcasting(∘),cell_to_b,cell_to_a)
cell_to_r = fill(evaluate(f,x),ncells)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)

cell_to_∇f = lazy_map(Broadcasting(∇),cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
cell_to_∇r = fill(evaluate(∇.(f),x),ncells)
test_array(cell_to_∇fx,cell_to_∇r)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# Operations

afun(x) = x[2]
bfun(x) = x[1] + 3
a = GenericField(afun)
b = GenericField(bfun)
f = Operation(+)(a,b)

cell_to_a = fill(a,ncells)
cell_to_b = Fill(b,ncells)
cell_to_r = fill(evaluate(f,x),ncells)
cell_to_∇r = fill(evaluate(∇(f),x),ncells)

cell_to_f = lazy_map(Operation(+),cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)
@test cell_to_fx.maps.value == Broadcasting(+)

cell_to_∇f = lazy_map(∇,cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
test_array(cell_to_∇fx,cell_to_∇r)
@test cell_to_fx.maps.value == Broadcasting(+)

T = GenericField{Nothing}
cell_to_f = lazy_map(Operation(+),T,cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)
@test cell_to_fx.maps.value == Broadcasting(+)

cell_to_∇f = lazy_map(∇,cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
test_array(cell_to_∇fx,cell_to_∇r)
@test cell_to_fx.maps.value == Broadcasting(+)

#using BenchmarkTools
#c = array_cache(cell_to_fx)
#@btime getindex!($c,$cell_to_fx,1)

# Operations with Broadcasting

afun(x) = x[2]
bfun(x) = x[1] + 3
a = GenericField(afun)
b = fill(GenericField(bfun),4)
f = Broadcasting(Operation(+))(a,b)

cell_to_a = fill(a,ncells)
cell_to_b = Fill(b,ncells)
cell_to_r = fill(evaluate(f,x),ncells)
cell_to_∇r = fill(evaluate(Broadcasting(∇)(f),x),ncells)

cell_to_f = lazy_map(Broadcasting(Operation(+)),cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)
#@test cell_to_fx.maps.value == Broadcasting(+)
@test cell_to_fx.maps.value == Fields.BroadcastingFieldOpMap(+)

cell_to_∇f = lazy_map(Broadcasting(∇),cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
test_array(cell_to_∇fx,cell_to_∇r)
#@test cell_to_∇fx.maps.value == Broadcasting(+)
@test cell_to_∇fx.maps.value == Fields.BroadcastingFieldOpMap(+)

T = GenericField{Nothing}
cell_to_f = lazy_map(Broadcasting(Operation(+)),T,cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)
#@test cell_to_fx.maps.value == Broadcasting(+)
@test cell_to_fx.maps.value == Fields.BroadcastingFieldOpMap(+)

# Operations with Broadcasting (product)

afun(x) = x[2]
bfun(x) = x[1] + 3
a = GenericField(afun)
b = fill(GenericField(bfun),4)
f = Broadcasting(Operation(*))(a,b)

cell_to_a = fill(a,ncells)
cell_to_b = Fill(b,ncells)
cell_to_r = fill(evaluate(f,x),ncells)
cell_to_∇r = fill(evaluate(Broadcasting(∇)(f),x),ncells)

cell_to_f = lazy_map(Broadcasting(Operation(*)),cell_to_b,cell_to_a)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
test_array(cell_to_fx,cell_to_r)

cell_to_∇f = lazy_map(Broadcasting(∇),cell_to_f)
cell_to_∇fx = lazy_map(evaluate,cell_to_∇f,cell_to_x)
test_array(cell_to_∇fx,cell_to_∇r)

# Integration

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = rand(length(avals),5)
d = linear_combination(b,a)
c = Broadcasting(Operation(*))(d,transpose(d))
w = rand(np)
fx = integrate(c,x,w)

cell_to_a = Fill(a,ncells)
cell_to_b = fill(b,ncells)
cell_to_w = Fill(w,ncells)
cell_to_d = lazy_map(linear_combination,cell_to_b,cell_to_a)
cell_to_dt = lazy_map(transpose,cell_to_d)
cell_to_c = lazy_map(Broadcasting(Operation(*)),cell_to_d,cell_to_dt)
cell_to_fx = lazy_map(integrate,cell_to_c,cell_to_x,cell_to_w)
cell_to_r = fill(fx,ncells)
test_array(cell_to_fx,cell_to_r)
@test isa(cell_to_fx.maps.value,IntegrationMap)

ϕ = GenericField(x->2*x)
fx = integrate(c,x,w,∇(ϕ))

cell_to_ϕ = fill(ϕ,ncells)
cell_to_J = lazy_map(∇,cell_to_ϕ)

cell_to_fx = lazy_map(integrate,cell_to_c,cell_to_x,cell_to_w,cell_to_J)
cell_to_r = fill(fx,ncells)
test_array(cell_to_fx,cell_to_r)
@test isa(cell_to_fx.maps.value,IntegrationMap)

# Reindex

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = VectorValue{2,Float64}[(1,1),(4,2),(3,5),(1,2)]
f = linear_combination(b,a)

cell_to_a = Fill(a,ncells)
cell_to_b = fill(b,ncells)
cell_to_c = lazy_map(linear_combination,cell_to_b,cell_to_a)

j_to_cell = [2,3,2,2,1]
j_to_c = lazy_map(Reindex(cell_to_c),j_to_cell)
j_to_x = Fill(x,length(j_to_cell))

j_to_cx = lazy_map(evaluate,j_to_c,j_to_x)
j_to_r = fill( f(x), length(j_to_cell) )
@test isa(j_to_cx.maps.value,Fields.LinearCombinationMap)
test_array(j_to_cx,j_to_r)

j_to_∇c = lazy_map(Broadcasting(∇),j_to_c)
j_to_∇cx = lazy_map(evaluate,j_to_∇c,j_to_x)
@test isa(j_to_∇cx.maps.value,Fields.LinearCombinationMap)
j_to_∇r = fill( ∇(f)(x), length(j_to_cell) )
test_array(j_to_∇cx,j_to_∇r)

# PosNegReindex

aval = 1.0
a = MockField(aval)
ipos_to_cell = [2,5,1]
cell_to_i = PosNegPartition(ipos_to_cell,ncells)
npos = length(cell_to_i.ipos_to_i)
nneg = length(cell_to_i.ineg_to_i)

ipos_to_a = fill(a,npos)
ineg_to_a = fill(0*a,nneg)

T = Union{eltype(ipos_to_a),eltype(ineg_to_a)}
cell_to_f = lazy_map(PosNegReindex(ipos_to_a,ineg_to_a),T,cell_to_i)
cell_to_fx = lazy_map(evaluate,cell_to_f,cell_to_x)
cell_to_r = [ zeros(np) for cell in 1:ncells]
cell_to_r[ipos_to_cell] = [ fill(aval,np) for cell in 1:npos]
#print_op_tree(cell_to_fx)

@test isa(cell_to_fx.maps.value,PosNegReindex)
test_array(cell_to_fx,cell_to_r)

#npp = np #(np,np)
#
#v = VectorValue{d}(1.0,1.0)
#f = MockField{d}(v)
#fx = fill(v,np)
#
#nf = 3
#nff = nf #(nf,nf)
#
#fa = fill(f,nf)
#
## Gradients
#
#∇fa = Broadcasting(∇)(fa)
#c = return_cache(∇fa,x)
## @santiagobadia : Not sure we can make it allocation-free
#evaluate!(c,∇fa,x)
## @btime evaluate!($c,$∇fa,$x)
#
#_∇fa = gradient.(fa)
#c = return_cache(∇fa,x)
## @santiagobadia : Whereas this one is allocation-free...
## @btime evaluate!($c,$_∇fa,$x)
#
#@test evaluate!(c,∇fa,x) == evaluate!(c,_∇fa,x)
#
## Transpose
#
#c = return_cache(fa,p)
## @btime evaluate!($c,$fa,$p);
#
#c = return_cache(fa,x)
## @btime evaluate!($c,$fa,$x);
#
#tfa = transpose(fa)
#c = return_cache(tfa,p)
## @btime evaluate!($c,$tfa,$p);
#tr = evaluate!(c,tfa,p)
#r = evaluate!(c,fa,p)
#@test transpose(tr) == r
#
#c = return_cache(tfa,x)
## @btime evaluate!($c,$tfa,$x);
#tr = evaluate!(c,tfa,x)
#r = evaluate!(c,fa,x)
#@test size(tr) == (size(r,1),1,size(r,2))
#@test reshape(permutedims(tr,[1,3,2]),size(r,1),size(r,2)) == r
#
## Broadcasted operations
#
#op = +
#b = Broadcasting(Operation(op))
#bpfa = b(fa,fa)
#@test bpfa == op(fa,fa)
#
#fa isa AbstractArray{<:Field}
#fa+fa
#+(fa,fa)
#
#c = return_cache(bpfa,x)
## @btime evaluate!($c,$bpfa,$x);
#
#op = ⋅
#b = Broadcasting(Operation(op))
#bpfa = b(fa,fa)
#@test bpfa == op.(fa,fa)
#
#c = return_cache(bpfa,x)
## @btime evaluate!($c,$bpfa,$x);
#
## transpose
#
#tfa = transpose(fa)
#tr = evaluate(tfa,x)
#r = evaluate(fa,x)
#
#c = return_cache(tfa,x)
## @btime evaluate!($c,$tfa,$x);
#
## field + field array
#
#op = ⋅
#b = Broadcasting(Operation(op))
#bpfa = b(f,fa)
#@test evaluate(bpfa,x) == broadcast(⋅,evaluate(f,x),evaluate(fa,x))
#
#c = return_cache(bpfa,x)
## @btime evaluate!(c,bpfa,x)
#
## column vector * row vector -> field array
#
#op = ⋅
#b = Broadcasting(Operation(op))
#bpfa = b(fa,tfa)
#evaluate(bpfa,x) == broadcast(⋅,evaluate(fa,x),evaluate(tfa,x))
#@test all(bpfa .== broadcast(⋅,fa,tfa))
#
#c = return_cache(bpfa,x)
## @btime evaluate!(c,bpfa,x)
#
## row vector * column vector -> field
#
#dopa = tfa*fa
#c = return_cache(dopa,x)
#evaluate!(c,dopa,x)
#@test evaluate(dopa,x)[1] == evaluate(tfa,p)⋅evaluate(fa,p)
## @btime evaluate!(c,dopa,x);
#
#dopa = fa⋅fa
#c = return_cache(dopa,x)
#evaluate!(c,dopa,x)
#@test evaluate(dopa,x)[1] == evaluate(tfa,p)⋅evaluate(fa,p)
## @btime evaluate!(c,dopa,x);
#
## linear combination
#
#va = [1.0,2.0,3.0]
#dopa = fa⋅va
#@test evaluate(va⋅fa,x) == evaluate(fa⋅va,x)
#c = return_cache(dopa,x)
## @btime evaluate!(c,dopa,x);
#
#@test linear_combination(fa,va) == linear_combination(va,fa)
#lc = linear_combination(fa,va)
#fax = evaluate(fa,x)
#lcx = evaluate(lc,x)
#for i in 1:length(x)
#  @test lcx[i] == fax[i,:]⋅va
#end
#
#c  =return_cache(lc,x)
## @btime evaluate!($c,$lc,$x);
#
#vb = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
#@test linear_combination(fa,vb) == linear_combination(vb,fa)
#lc = linear_combination(fa,vb)
#fax = evaluate(fa,x)
#lcx = evaluate(lc,x)
#for i in 1:length(x)
#  @test lcx[i,:] == transpose(transpose(fax[i,:])*vb)
#end
#
#c  =return_cache(lc,x)
## @btime evaluate!($c,$lc,$x);
#
#lc = linear_combination(fa,va)
#fap = evaluate(fa,p)
#lcx = evaluate(lc,p)
#@test lcx == fap⋅va
#
#c  =return_cache(lc,p)
## @btime evaluate!($c,$lc,$p);
#
#lc = linear_combination(fa,vb)
#fap = evaluate(fa,p)
#lcp = evaluate(lc,p)
#@test lcp == transpose(transpose(fap) * vb)
#
#c  =return_cache(lc,p)
## @btime evaluate!($c,$lc,$p);
#
## composition
#
#bm = Broadcasting(∘)
#fc = bm(fa,f);
#
#c = return_cache(fc,x)
## @btime evaluate!(c,fc,x)

end # module