text:
# Benchmark Catlab.Graphs against LightGraphs and MetaGraphs.
using BenchmarkTools
const SUITE = BenchmarkGroup()

using Random
import LightGraphs, MetaGraphs
const LG, MG = LightGraphs, MetaGraphs

using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphs
using Catlab.Graphs.BasicGraphs: TheoryGraph

# Helpers
#########

# `bench_iter_edges` and `bench_has_edge` adapted from LightGraphs:
# https://github.com/JuliaGraphs/LightGraphs.jl/blob/master/benchmark/core.jl

function bench_iter_edges(g)
  count = 0
  for e in edges(g)
    s, t = src(g,e), tgt(g,e)
    count += 1
  end
  count
end
function bench_iter_edges_vectorized(g)
  count = 0
  for (s,t) in zip(src(g), tgt(g))
    count += 1
  end
  count
end

function bench_has_edge(g)
  Random.seed!(1)
  n = nv(g)
  srcs, tgts = rand(1:n, n÷4), rand(1:n, n÷4)
  count = 0
  for (s, t) in zip(srcs, tgts)
    if has_edge(g, s, t)
      count += 1
    end
  end
  count
end

function bench_iter_neighbors(g)
  count = 0
  for v in vertices(g)
    count += length(neighbors(g, v))
  end
  count
end

@inline Graphs.nv(g::LG.AbstractGraph) = LG.nv(g)
@inline Graphs.vertices(g::LG.AbstractGraph) = LG.vertices(g)
@inline Graphs.edges(g::LG.AbstractGraph) = LG.edges(g)
@inline Graphs.src(g::LG.AbstractGraph, e::LG.AbstractEdge) = LG.src(e)
@inline Graphs.tgt(g::LG.AbstractGraph, e::LG.AbstractEdge) = LG.dst(e)
@inline Graphs.has_edge(g::LG.AbstractGraph, args...) = LG.has_edge(g, args...)
@inline Graphs.neighbors(g::LG.AbstractGraph, args...) = LG.neighbors(g, args...)

# Graphs
########

bench = SUITE["Graph"] = BenchmarkGroup()

n = 10000
bench["make-path"] = @benchmarkable begin
  g = Graph()
  add_vertices!(g, n)
  add_edges!(g, 1:(n-1), 2:n)
end

bench["make-path-lightgraphs"] = @benchmarkable begin
  g = LG.DiGraph()
  LG.add_vertices!(g, n)
  for v in 1:(n-1)
    LG.add_edge!(g, v, v+1)
  end
end

g = Graph(n)
add_edges!(g, 1:(n-1), 2:n)
lg = LG.DiGraph(g)

bench["iter-edges"] = @benchmarkable bench_iter_edges($g)
bench["iter-edges-vectorized"] = @benchmarkable bench_iter_edges_vectorized($g)
bench["iter-edges-lightgraphs"] = @benchmarkable bench_iter_edges($lg)
bench["has-edge"] = @benchmarkable bench_has_edge($g)
bench["has-edge-lightgraphs"] = @benchmarkable bench_has_edge($lg)
bench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($g)
bench["iter-neighbors-lightgraphs"] = @benchmarkable bench_iter_neighbors($lg)

n₀ = 2000
g₀ = Graph(n₀)
add_edges!(g₀, 1:(n₀-1), 2:n₀)
g = ob(coproduct(fill(g₀, 5)))
lg = LG.DiGraph(g)
bench["path-graph-components"] = @benchmarkable connected_components($g)
bench["path-graph-components-proj"] =
  @benchmarkable connected_component_projection($g)
bench["path-graph-components-lightgraphs"] =
  @benchmarkable LG.weakly_connected_components($lg)

g₀ = Graph(n₀)
add_edges!(g₀, fill(1,n₀-1), 2:n₀)
g = ob(coproduct(fill(g₀, 5)))
lg = LG.DiGraph(g)
bench["star-graph-components"] = @benchmarkable connected_components($g)
bench["star-graph-components-proj"] =
  @benchmarkable connected_component_projection($g)
bench["star-graph-components-lightgraphs"] =
  @benchmarkable LG.weakly_connected_components($lg)

# Symmetric graphs
##################

bench = SUITE["SymmetricGraph"] = BenchmarkGroup()

n = 10000
bench["make-path"] = @benchmarkable begin
  g = SymmetricGraph()
  add_vertices!(g, n)
  add_edges!(g, 1:(n-1), 2:n)
end

bench["make-path-lightgraphs"] = @benchmarkable begin
  g = LG.Graph()
  LG.add_vertices!(g, n)
  for v in 1:(n-1)
    LG.add_edge!(g, v, v+1)
  end
end

g = SymmetricGraph(n)
add_edges!(g, 1:(n-1), 2:n)
lg = LG.Graph(g)

bench["iter-edges"] = @benchmarkable bench_iter_edges($g)
bench["iter-edges-vectorized"] = @benchmarkable bench_iter_edges_vectorized($g)
bench["iter-edges-lightgraphs"] = @benchmarkable bench_iter_edges($lg)
bench["has-edge"] = @benchmarkable bench_has_edge($g)
bench["has-edge-lightgraphs"] = @benchmarkable bench_has_edge($lg)
bench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($g)
bench["iter-neighbors-lightgraphs"] = @benchmarkable bench_iter_neighbors($lg)

function lg_connected_components_projection(g)
  label = Vector{Int}(undef, LG.nv(g))
  LG.connected_components!(label, g)
end

n₀ = 2000
g₀ = SymmetricGraph(n₀)
add_edges!(g₀, 1:(n₀-1), 2:n₀)
g = ob(coproduct(fill(g₀, 5)))
lg = LG.Graph(g)
bench["path-graph-components"] = @benchmarkable connected_components($g)
bench["path-graph-components-proj"] =
  @benchmarkable connected_component_projection($g)
bench["path-graph-components-lightgraphs"] =
  @benchmarkable LG.connected_components($lg)
bench["path-graph-components-proj-lightgraphs"] =
  @benchmarkable lg_connected_components_projection($lg)

g₀ = SymmetricGraph(n₀)
add_edges!(g₀, fill(1,n₀-1), 2:n₀)
g = ob(coproduct(fill(g₀, 5)))
lg = LG.Graph(g)
bench["star-graph-components"] = @benchmarkable connected_components($g)
bench["star-graph-components-proj"] =
  @benchmarkable connected_component_projection($g)
bench["star-graph-components-lightgraphs"] =
  @benchmarkable LG.connected_components($lg)
bench["star-graph-components-proj-lightgraphs"] =
  @benchmarkable lg_connected_components_projection($lg)

# Weighted graphs
#################

bench = SUITE["WeightedGraph"] = BenchmarkGroup()

n = 10000
g = WeightedGraph{Float64}(n)
add_edges!(g, 1:(n-1), 2:n, weight=range(0, 1, length=n-1))
mg = MG.MetaDiGraph(g)

bench["sum-weights-vectorized"] = @benchmarkable sum(weight($g))
bench["sum-weights"] = @benchmarkable begin
  total = 0.0
  for e in edges($g)
    total += weight($g, e)
  end
  total
end
bench["sum-weights-metagraphs"] = @benchmarkable begin
  total = 0.0
  for e in MG.edges($mg)
    total += MG.get_prop($mg, e, :weight)
  end
  total
end

bench["increment-weights-vectorized"] = @benchmarkable begin
  $g[:weight] = $g[:weight] .+ 1.0
end
bench["increment-weights"] = @benchmarkable begin
  for e in edges($g)
    $g[e,:weight] += 1.0
  end
end
bench["increment-weights-metagraphs"] = @benchmarkable begin
  for e in MG.edges($mg)
    MG.set_prop!($mg, e, :weight, MG.get_prop($mg, e, :weight) + 1.0)
  end
end

# Labeled graphs
################

bench = SUITE["LabeledGraph"] = BenchmarkGroup()

@present TheoryLabeledGraph <: TheoryGraph begin
  Label::Data
  label::Attr(V,Label)
end
const LabeledGraph = ACSetType(TheoryLabeledGraph, index=[:src,:tgt])
const IndexedLabeledGraph = ACSetType(TheoryLabeledGraph, index=[:src,:tgt],
                                      unique_index=[:label])

function discrete_labeled_graph(n::Int; indexed::Bool=false)
  g = (indexed ? IndexedLabeledGraph{String} : LabeledGraph{String})()
  add_vertices!(g, n, label=("v$i" for i in 1:n))
  g
end

function discrete_labeled_metagraph(n::Int; indexed::Bool=false)
  mg = MG.MetaDiGraph()
  for i in 1:n
    MG.add_vertex!(mg, :label, "v$i")
  end
  if indexed; MG.set_indexing_prop!(mg, :label) end
  mg
end

n = 5000
bench["make-discrete"] = @benchmarkable discrete_labeled_graph($n)
bench["make-discrete-metagraphs"] = @benchmarkable discrete_labeled_metagraph($n)
bench["make-discrete-indexed"] =
  @benchmarkable discrete_labeled_graph($n, indexed=true)
bench["make-discrete-indexed-metagraphs"] =
  @benchmarkable discrete_labeled_metagraph($n, indexed=true)

n = 10000
g = discrete_labeled_graph(n)
mg = discrete_labeled_metagraph(n)
bench["iter-labels"] = @benchmarkable begin
  for v in vertices($g)
    label = $g[v,:label]
  end
end
bench["iter-labels-metagraphs"] = @benchmarkable begin
  for v in MG.vertices($mg)
    label = MG.get_prop($mg, v, :label)
  end
end

g = discrete_labeled_graph(n, indexed=true)
mg = discrete_labeled_metagraph(n, indexed=true)
Random.seed!(1)
σ = randperm(n)
bench["indexed-lookup"] = @benchmarkable begin
  for i in $σ
    @assert incident($g, "v$i", :label) == i
  end
end
bench["indexed-lookup-metagraphs"] = @benchmarkable begin
  for i in $σ
    @assert $mg["v$i", :label] == i
  end
end


text:
using ITensorNetworkAD
using Documenter

DocMeta.setdocmeta!(
  ITensorNetworkAD, :DocTestSetup, :(using ITensorNetworkAD); recursive=true
)

makedocs(;
  modules=[ITensorNetworkAD],
  authors="Matthew Fishman <mfishman@flatironinstitute.org>, Linjian Ma <lma16@illinois.edu>",
  repo="https://github.com/ITensor/ITensorNetworkAD.jl/blob/{commit}{path}#{line}",
  sitename="ITensorNetworkAD.jl",
  format=Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    canonical="https://itensor.github.io/ITensorNetworkAD.jl",
    assets=String[],
  ),
  pages=["Home" => "index.md"],
)

deploydocs(;
  repo="github.com/ITensor/ITensorNetworkAD.jl.git", devbranch="main", push_preview=true
)

ext: jl


text:
using Plots, LaTeXStrings
gr(size=(600,300),colorbar=false,legend=false,grid=false)

using FourierGPE

## initialize
L = (60.0,)
N = (512,)
sim = Sim(L,N)
@unpack_Sim sim

## simulation parameters
μ = 1.0
g = 0.01
γ = 0.0
tf = 20.0pi
Nt = 150
ti = 0.0
t = LinRange(ti,tf,Nt)

## Bogoliubov state
x = X[1]
k = K[1]

## Define Bog state
f(k) = 1 + 4/k^2
u(k) = 0.5*(sqrt(f(k))+1)/f(k)^(1/4)
v(k) = 0.5*(sqrt(f(k))-1)/f(k)^(1/4)
lam = 0.01
bog(x,k) = u(k)*exp(im*k*x) - conj(v(k))*exp(-im*k*x)
ψb(x,k) = sqrt(μ/g)*(complex(1) + lam*bog(x,k))

## choose an allowed k for periodic system
kb = k[20]
ψi = ψb.(x,kb)
ϕi = kspace(ψi,sim)

# Set time evolution and pack
# reltol = 1e-7
# alg = Vern7()

@pack_Sim! sim

## Evolve in k space
@time sol = runsim(sim)
y = abs2.(xspace(sol[end-1],sim)); plot(x,y)

## density movie
anim = @animate for i=1:Nt
    ψ = xspace(sol[i],sim)
    y = g*abs2.(ψ)
    plot(x,y,c=c3,fill=(0,0.5,c3))
    plot!(x,one.(x)*μ,c=:black,alpha=0.3)
    ylims!(0,1.1)
    xlims!(first(x),last(x))
    xlabel!(L"x/\xi"); ylabel!(L"gn(x)/\mu")
end

gif(anim, "./examples/bogoliubov_density.gif", fps = 20)


text:
module BayesLaplace

include("data/Data.jl")
using .Data

include("curvature/Curvature.jl")
using .Curvature

include("laplace.jl")

export Laplace, fit!, predict, plugin

end


text:
using Test, StaticArrays
using ClimateMachine.VariableTemplates: varsindex, @vars, varsindices
using StaticArrays

@testset "varsindex" begin
    struct TestMoistureModel{FT} end
    struct TestAtmosModel{FT}
        moisture::TestMoistureModel{FT}
    end

    function vars_state(::TestMoistureModel, FT)
        @vars begin
            ρq_tot::FT
            ρq_x::SVector{5, FT}
            ρq_liq::FT
            ρq_vap::FT
        end
    end
    function vars_state(m::TestAtmosModel, FT)
        @vars begin
            ρ::FT
            ρu::SVector{3, FT}
            ρe::FT
            moisture::vars_state(m.moisture, FT)
            S::SHermitianCompact{3, FT, 6}
        end
    end
    FT = Float64
    m = TestAtmosModel(TestMoistureModel{FT}())

    @test 1:1 === varsindex(vars_state(m, FT), :ρ)
    @test 2:4 === varsindex(vars_state(m, FT), :ρu)
    @test 5:5 === varsindex(vars_state(m, FT), :ρe)

    # Since moisture is defined recusively this will get all the fields
    moist = varsindex(vars_state(m, FT), :moisture)
    @test 6:13 === moist

    # To get the specific ones we can do
    @test 6:6 === varsindex(vars_state(m, FT), :moisture, :ρq_tot)
    @test 7:11 === varsindex(vars_state(m, FT), :moisture, :ρq_x)
    @test 12:12 === varsindex(vars_state(m, FT), :moisture, :ρq_liq)
    @test 13:13 === varsindex(vars_state(m, FT), :moisture, :ρq_vap)
    # or
    @test 6:6 === moist[varsindex(vars_state(m.moisture, FT), :ρq_tot)]
    @test 7:11 === moist[varsindex(vars_state(m.moisture, FT), :ρq_x)]
    @test 12:12 === moist[varsindex(vars_state(m.moisture, FT), :ρq_liq)]
    @test 13:13 === moist[varsindex(vars_state(m.moisture, FT), :ρq_vap)]

    @test 14:19 == varsindex(vars_state(m, FT), :S)

    @test (1,) === varsindices(vars_state(m, FT), :ρ)
    @test (2, 3, 4) === varsindices(vars_state(m, FT), :ρu)
    @test (1, 5) === varsindices(vars_state(m, FT), :ρ, :ρe)
    @test (12,) === varsindices(vars_state(m, FT), :(moisture.ρq_liq))
    let
        vars = ("ρe", "moisture.ρq_x", "moisture.ρq_vap")
        @test (5, 7, 8, 9, 10, 11, 13) === varsindices(vars_state(m, FT), vars)
    end
end


text:
using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays,
        FillArrays, SpecialFunctions, Test, SemiseparableMatrices, LazyBandedMatrices, BlockArrays
import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!, arguments
import LazyBandedMatrices: BroadcastBandedLayout, InvDiagTrav, BroadcastBandedBlockBandedLayout
import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout, adaptiveqr
import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout

@testset "Adaptive QR" begin
    @testset "banded" begin
        @testset "test partialqr!" begin
            A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
            l,u = bandwidths(A)
            F = AdaptiveQRData(A);
            @test bandwidths(F.data) == (1,2)
            n = 3
            partialqr!(F,n);
            F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
            @test F̃.factors ≈ F.data[1:n+l,1:n]
            @test F̃.τ  ≈ F.τ
            @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
            n  = 6
            partialqr!(F,n);
            F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
            @test F̃.factors ≈ F.data[1:n+l,1:n]
            @test F̃.τ  ≈ F.τ
            @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
        end

        @testset "AdaptiveQRFactors" begin
            A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
            F = qr(A);
            @test F.factors[1,1] ≈ -sqrt(2)
            @test F.factors[100,100] ≈ qrunblocked(A[1:101,1:100]).factors[100,100]
            @test F.τ[1] ≈ 1+sqrt(2)/2
            @test F.τ[100] ≈ qrunblocked(A[1:101,1:100]).τ[100]
            Q,R = F;
            @test bandwidths(R) == (0,2)
            @test size(R) == (∞,∞)
            @test R[1,1] == -sqrt(2)
        end

        @testset "col/rowsupport" begin
            A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
            F = qr(A);
            @test MemoryLayout(typeof(F.factors)) isa AdaptiveLayout{BandedColumns{DenseColumnMajor}}
            @test bandwidths(F.factors) == (1,2)
            @test colsupport(F.factors,1) ==  1:2
            @test colsupport(F.factors,5) ==  3:6
            @test rowsupport(F.factors,1) ==  1:3
            @test rowsupport(F.factors,5) ==  4:7
            Q,R = F;
            @test MemoryLayout(typeof(R)) isa TriangularLayout
            @test colsupport(R,1) == 1:1
            @test colsupport(R,5) == 3:5
            @test rowsupport(R,1) == 1:3
            @test rowsupport(R,5) == 5:7
            @test colsupport(Q,1) == 1:2
            @test colsupport(Q,5) == 3:6
            @test rowsupport(Q,1) == 1:3
            @test rowsupport(Q,5) == 4:7
        end

        @testset "Qmul" begin
            A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
            Q,R = qr(A);
            b = Vcat([1.,2,3],Zeros(∞))
            @test colsupport(Q,1:3) == colsupport(Q.factors,1:3) == 1:4
            @test lmul!(Q, Base.copymutable(b)).datasize[1] == 4
            @test lmul!(Q, Base.copymutable(b)).data[1:4] ≈ qr(A[1:4,1:3]).Q*[1,2,3]

            @test Q[1,1] ≈ -1/sqrt(2)
            @test Q[200_000,200_000] ≈ -1.0
            @test Q[1:101,1:100] ≈ qr(A[1:101,1:100]).Q[:,1:100]

            @test lmul!(Q, Base.copymutable(b)) == Q*b

            r = lmul!(Q', Base.copymutable(b))
            nr = length(r.data)
            @test qr(A[1:nr+1,1:nr]).Q'b[1:nr+1] ≈ r[1:nr+1]

            @test Q'*b == r

            @test factorize(A) isa typeof(qr(A))
            @test qr(A)\b == A\b
            @test (A*(A\b))[1:100] ≈ [1:3; Zeros(97)]

            @test Q * view(b,:) == Q*b
            @test Q' * view(b,:) == Q'*b
        end

        @testset "Bessel J" begin
            z = 1000; # the bigger z the longer before we see convergence
            A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
            b = Vcat([besselj(1,z)], Zeros(∞))
            F = qr(A)
            @test qr(A[1:3000,1:3000]).Q'b[1:3000] ≈ (F.Q'b)[1:3000]
            @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
            @test J[1:2000] ≈ [besselj(k,z) for k=0:1999]

            z = 10_000; # the bigger z the longer before we see convergence
            A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
            @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
            @test J[1:20_000] ≈ [besselj(k,z) for k=0:20_000-1]

            # Tridiagonal works too
            A =  LazyBandedMatrices.Tridiagonal(Ones(∞), -2*(0:∞)/z, Ones(∞))
            @test factorize(A) isa MatrixFactorizations.QR
            @time J = A \ [besselj(1,z); Zeros(∞)]
            @test J[1:20_000] ≈ [besselj(k,z) for k=0:20_000-1]
        end

        @testset "5-band" begin
            A = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))
            b = Vcat([3,4,5],Zeros(∞))
            x = qr(A) \ b
            @test x[1:2000] ≈ (A[1:2000,1:2000]\b[1:2000])
        end

        @testset "broadcast" begin
            A = BandedMatrix(0 => -2*(0:∞), 1 => Ones(∞), -1 => Ones(∞))
            B = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))

            AB = BroadcastArray(+,A,B)
            @test MemoryLayout(typeof(AB)) isa BroadcastBandedLayout{typeof(+)}

            C = cache(AB);
            resizedata!(C,103,100); resizedata!(C,203,200);
            @test C[103,104] ≈ 1.0
            C = qr(AB).factors.data.data;
            resizedata!(C,102,104); resizedata!(C,202,204);
            @test C[202,204] == AB[202,204]
            F = qr(AB);
            partialqr!(F.factors.data, 100); partialqr!(F.factors.data, 200);
            @test norm(F.factors.data.data.data[Base.OneTo.(F.factors.data.data.datasize)...]) ≤ 4000
            b = Vcat([3,4,5],Zeros(∞))
            @time x = qr(AB) \ b;
            @test x[1:300] ≈ AB[1:300,1:300] \ b[1:300]
        end

        @testset "triangular infqr" begin
            A = BandedMatrix(0 => 1:∞, 2 => Ones(∞))
            F = qr(A)
            @test F.Q[1:10,1:10] == Eye(10)
            @test F.R[1:10,1:10] == A[1:10,1:10]
        end

        @testset "diag special case" begin
            A = _BandedMatrix((1:∞)', ℵ₀, 0, 0)
            b = [[1,2,3]; zeros(∞)]
            @test A \ b == [ones(3); zeros(∞)]
        end

        @testset "Symmetric" begin
            A = Symmetric(BandedMatrix(0 => 1:∞, 1=> Ones(∞)))
            Ã = BandedMatrix(0 => 1:∞, 1=> Ones(∞), -1=> Ones(∞))
            @test qr(A).R[1:10,1:10] ≈ qr(Ã).R[1:10,1:10]
        end
    end

    @testset "almost-banded" begin
        @testset "one-band" begin
            A = Vcat(Ones(1,∞), BandedMatrix(0 => -Ones(∞), 1 => 1:∞))
            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
            V = view(A,1:10,1:10)
            @test MemoryLayout(typeof(V)) isa VcatAlmostBandedLayout
            @test A[1:10,1:10] isa AlmostBandedMatrix
            @test AlmostBandedMatrix(V) == Matrix(V) == A[1:10,1:10]

            C = cache(A);
            @test C[1000,1000] ≡ 999.0
            F = adaptiveqr(A);
            partialqr!(F.factors.data,2);
            @test F.factors.data.data[1:3,1:5] ≈ qr(A[1:3,1:5]).factors
            partialqr!(F.factors.data,3);
            @test F.factors.data.data[1:4,1:6] ≈ qr(A[1:4,1:6]).factors

            F = adaptiveqr(A);
            partialqr!(F.factors.data,10);
            @test F.factors[1:11,1:10] ≈ qr(A[1:11,1:10]).factors
            @test F.τ[1:10] ≈ qr(A[1:11,1:10]).τ
            partialqr!(F.factors.data,20);
            @test F.factors[1:21,1:20] ≈ qr(A[1:21,1:20]).factors

            @test adaptiveqr(A).R[1:10,1:10] ≈ qr(A[1:11,1:10]).R

            @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
            @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}

            @test (adaptiveqr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (qr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (A \ [ℯ; zeros(∞)])[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
        end

        @testset "two-bands" begin
            B = BandedMatrix(0 => -Ones(∞), 2 => (1:∞).* (2:∞))
            A = Vcat(Vcat(Ones(1,∞), ((-1).^(0:∞))'), B)
            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout

            @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
            @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
            u = qr(A) \ [1; zeros(∞)]
            x = 0.1
            @test (exp(1 - x)*(-1 + exp(2 + 2x)))/(-1 + exp(4)) ≈ dot(u[1:1000], x.^(0:999))
            u = qr(A) \ Vcat([ℯ,1/ℯ], zeros(∞))
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
            u = qr(A) \ Vcat(ℯ,1/ℯ, zeros(∞))
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
            u = qr(A) \ [ℯ; 1/ℯ; zeros(∞)]
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]

            A = Vcat(Ones(1,∞), ((-1.0).^(0:∞))', B)
            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
            u = A \ Vcat(ℯ,1/ℯ, zeros(∞))
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
            u = qr(A) \ [ℯ; 1/ℯ; zeros(∞)]
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]

            A = Vcat(Ones(1,∞), ((-1).^(0:∞))', B)
            u = A \ [ℯ; 1/ℯ; zeros(∞)]
            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
        end

        @testset "more bands" begin
            L = Vcat(Ones(1,∞), ((-1).^(0:∞))',
                     BandedMatrix(-1 => Ones(∞), 1 => Ones(∞), 2 => 4:2:∞, 3 => Ones(∞), 5 => Ones(∞)))
            F = qr(L).factors.data;
            resizedata!(F.data,13,19)
            @test F.data.data[2,8] == -1
            F = qr(L);
            partialqr!(F,10);
            @test F.factors[1:10,1:10] ≈ qr(L[1:13,1:10]).factors[1:10,1:10]
            @test qr(L).factors[1:10,1:10] ≈ qr(L[1:13,1:10]).factors[1:10,1:10]

            u = L \ [1; 2; zeros(∞)]
            @test L[1:1000,1:1000]*u[1:1000] ≈ [1; 2; zeros(998)]
        end
    end

    @testset "block-banded" begin
        Δ = BandedMatrix(1 => Ones(∞), -1 => Ones(∞))/2
        A = KronTrav(Δ - 2I, Eye(∞))
        F = qr(A);
        @test abs.(F.factors[1:15,1:10]) ≈ abs.(qr(A[1:15,1:10]).factors)

        @test (F.Q' * [1; zeros(∞)])[1:6] ≈ [-0.9701425001453321,0,-0.23386170701251197,0,0,-0.06193705069863463]
        @test (F.Q*[1;zeros(∞)])[1:6] ≈ [-0.9701425001453321,0,0.24253562503633297,0,0,0]

        u = F \ [1; zeros(∞)]
        @test blockisequal(axes(A,2),axes(u,1))
        @test (A*u)[1:10] ≈ [1; zeros(9)]

        x = 0.1
        θ = acos(x)
        @test dot(u[getindex.(Block.(1:50),1:50)], sin.((1:50) .* θ)/sin(θ)) ≈ 1/(x-2)

        B = KronTrav(Eye(∞), Δ - 2I)
        u = B \ [1; zeros(∞)]

        @test dot(u[getindex.(Block.(1:50),1)], sin.((1:50) .* θ)/sin(θ)) ≈ 1/(x-2)

        L = A+B;
        @test MemoryLayout(L) isa BroadcastBandedBlockBandedLayout{typeof(+)}
        V = view(L,Block.(1:400),Block.(1:400))
        @test blockbandwidths(V) == blockbandwidths(L)
        @test subblockbandwidths(V) == blockbandwidths(L)
        @time u = L \ [1;zeros(∞)]

        x,y = 0.1,0.2
        θ,φ = acos(x),acos(y)
        @test u[Block.(1:50)] isa PseudoBlockArray
        @test (sin.((1:50) .* φ)/sin(φ))' * InvDiagTrav(u[Block.(1:50)]) * sin.((1:50) .* θ)/sin(θ) ≈ 1/(x+y-4)
        @test (L*u)[1:10] ≈ [1; zeros(9)]

        X = KronTrav(Δ,Eye(∞))
        Y = KronTrav(Eye(∞),Δ)
        II = KronTrav(Eye(∞),Eye(∞))
        @test MemoryLayout(2X + Y - 4II) isa BroadcastBandedBlockBandedLayout
        u =  (2X + Y - 8II) \ [1; zeros(∞)]
        x,y = 0.1,0.2
        θ,φ = acos(x),acos(y)
        @test (sin.((1:100) .* φ)/sin(φ))' * InvDiagTrav(u[Block.(1:100)]) * sin.((1:100) .* θ)/sin(θ) ≈ 1/(2x+y-8)
    end

    @testset "SymTridiagonal Toeplitz" begin
        Δ = SymTridiagonal(Fill(-2,∞), Fill(1,∞))
        h = 0.01
        A = I - h*Δ
        b = [1; 2; 3; zeros(∞)]
        @test (qr(A) \ b) ≈ (ul(A) \ b)
    end
end


text:
using DiffEqBase
using Test, LinearAlgebra

@time begin
@time @testset "Number of Parameters Calculation" begin include("numargs_test.jl") end
@time @testset "Data Arrays" begin include("data_array_tests.jl") end
@time @testset "Extended Functions" begin include("extended_function_tests.jl") end
@time @testset "Callbacks" begin include("callbacks.jl") end
@time @testset "Plot Variables" begin include("plot_vars.jl") end
@time @testset "Problem Creation Tests" begin include("problem_creation_tests.jl") end
@time @testset "Remake tests" begin include("remake_tests.jl") end
@time @testset "Affine differential equation operators" begin include("affine_operators_tests.jl") end
@time @testset "TableTraits" begin include("tabletraits_tests.jl") end
@time @testset "Integrator interface" begin include("integrator_tests.jl") end
@time @testset "Export tests" begin include("export_tests.jl") end
@time @testset "High Level solve Interface" begin include("high_level_solve.jl") end
end


text:
# add Dash   # in package mode

using Dash

app = dash()                             # 1

app.layout = html_div() do               # 2
    html_h1("Hello from Dash!"),
    html_div("Dash: A web application framework for your data."),
    dcc_graph(
        id = "example-graph-1",
        figure = (
            data = [
                (x = ["Chocolate", "Strawberry", "Vanilla"], y = [21, 19, 33], type = "bar", name = "Male"),
                (x = ["Chocolate", "Strawberry", "Vanilla"], y = [38, 18, 12], type = "bar", name = "Female"),
            ],
            layout = (title = "Flavor Preferences by Gender", barmode="group")
        )
    )
end

run_server(app, "0.0.0.0", debug=true)   # 3


text:
#=
minimize the condition number of positive definite matrix
M(x) = M_0 + sum_i x_i*M_i
subject to F(x) = F_0 + sum_i x_i*F_i in S_+

original formulation:
minimize    gamma
subject to  mu >= 0
            F(x) in S_+
            M(x) - mu*I in S_+
            mu*gamma*I - M(x) in S_+
introduce nu = inv(mu), y = x/mu:
minimize    gamma
subject to  nu >= 0
            nu*F_0 + sum_i y_i*F_i in S_+
            nu*M_0 + sum_i y_i*M_i - I in S_+
            gamma*I - nu*M_0 - sum_i y_i*M_i in S_+
we make F_0 and M_0 positive definite to ensure existence of a feasible solution

see section 3.2 "Linear Matrix Inequalities in System and Control Theory" by
S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan
=#

struct ConditionNumJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    side::Int
    len_y::Int
    use_linmatrixineq::Bool # use linmatrixineq cone, else PSD formulation
end

function build(inst::ConditionNumJuMP{T}) where {T <: Float64}
    (side, len_y) = (inst.side, inst.len_y)

    rand_pd() = (Mh = randn(side, side); Symmetric(Mh * Mh'))
    Mi = [rand_pd() for i in 1:len_y]
    M0 = rand_pd()
    # make some F_i matrices pos def
    Fi = [(rand() > 0.5 || i <= 2) ? rand_pd() :
        Symmetric(randn(side, side)) for i in 1:len_y]
    F0 = rand_pd() + I

    model = JuMP.Model()
    JuMP.@variables(model, begin
        gamma
        nu >= 0
        y[1:len_y]
    end)
    JuMP.@objective(model, Min, gamma)

    if inst.use_linmatrixineq
        lmiT = Hypatia.LinMatrixIneqCone{T}
        JuMP.@constraints(model, begin
            vcat(nu, y) in lmiT([F0, Fi...])
            vcat(-1, nu, y) in lmiT([I, M0, Mi...])
            vcat(gamma, -nu, -y) in lmiT([I, M0, Mi...])
        end)
    else
        S1 = Symmetric(nu * F0 + sum(y[i] * Fi[i] for i in eachindex(y)))
        S2 = Symmetric(nu * M0 + sum(y[i] * Mi[i] for i in eachindex(y)) - I)
        S3 = Symmetric(gamma * Matrix(I, side, side) - nu * M0 -
            sum(y[i] * Mi[i] for i in eachindex(y)))

        JuMP.@constraints(model, begin
            S1 in JuMP.PSDCone()
            S2 in JuMP.PSDCone()
            S3 in JuMP.PSDCone()
        end)
    end

    return model
end


text:
function control_comparison(
      parent_folder::String,
      output_folder::String,
      objective::String;
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      print_results::Bool=false,
      results_folders::Array=[],
      output_to_files::Bool=true,
      override_param::Dict=Dict(),
      start_vals=Dict{String, Dict}("sn"=>Dict()),
      regularize_f::Float64=0.0
   )

   # Run the opf for the base case then all dc configurations
   results_dict_allplots = []
   n_subnets = 0
   subnet_array = []
   idx_sorted = []
   series_labels = ["control X(ω)","control P, Q","control P, Q, X(ω)"]
   if length(results_folders) == 0
      # Base case opf
      (results_dict,n_subnets,subnet_array,idx_sorted,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas=gen_areas,
         area_transfer=area_transfer,
         gen_zones=gen_zones,
         zone_transfer=zone_transfer,
         enum_branches=enum_branches,
         plot_best_x=plot_best_x,
         print_results=print_results,
         suffix="",
         output_to_files=output_to_files,
         start_vals=start_vals,
         override_param=override_param,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)

      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(series_output_folder)/all_results_dict.json")
      open("$series_output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end
      push!(results_folders, series_output_folder)

      # fixed frequency case
      fix_f_override=true
      direct_pq=true
      master_subnet=1
      suffix="_fixf"
      (results_dict,n_subnets_f,subnet_array_f,idx_sorted_f,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas=gen_areas,
         area_transfer=area_transfer,
         gen_zones=gen_zones,
         zone_transfer=zone_transfer,
         enum_branches=enum_branches,
         plot_best_x=plot_best_x,
         print_results=print_results,
         suffix=suffix,
         fix_f_override=fix_f_override,
         direct_pq=direct_pq,
         master_subnet=master_subnet,
         output_to_files=output_to_files,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)

      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(series_output_folder)/all_results_dict.json")
      open("$series_output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end
      push!(results_folders, series_output_folder)

      # indirect PQ case
      fix_f_override=false
      direct_pq=false
      master_subnet=1
      suffix="_indirPQ"
      (results_dict,n_subnets_pq,subnet_array_pq,idx_sorted_pq,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas=gen_areas,
         area_transfer=area_transfer,
         gen_zones=gen_zones,
         zone_transfer=zone_transfer,
         enum_branches=enum_branches,
         plot_best_x=plot_best_x,
         print_results=print_results,
         suffix=suffix,
         fix_f_override=fix_f_override,
         direct_pq=direct_pq,
         master_subnet=master_subnet,
         output_to_files=output_to_files,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)
      push!(results_folders, series_output_folder)

      if output_to_files
         results_dict_string = JSON.json(results_dict)
         println("Saving results at $(series_output_folder)/all_results_dict.json")
         open("$series_output_folder/all_results_dict.json", "w") do f
            write(f, results_dict_string)
         end
      end
   else
      # Use results data that has already been saved
      results_dict_allplots = []
      for (i,output_folder) in enumerate(results_folders)

         results_dict_parsed = JSON.parsefile("$(output_folder)/all_results_dict.json")
         if i == 1
            (n_subnets, subnet_array, idx_sorted) = get_sort_indices(results_dict_parsed, objective, gen_areas, area_transfer, gen_zones, zone_transfer)
         end

         push!(results_dict_allplots, results_dict_parsed)
      end
   end
   if output_to_files
      println("Saving results at $(output_folder)/results_dict_allplots.json")
      results_dict_string = JSON.json(results_dict_allplots)
      if !isdir(output_folder)
         mkpath(output_folder)
      end
      open("$(output_folder)/results_dict_allplots.json", "w") do f
         write(f, results_dict_string)
      end
   end
   plot_results_dicts_bar(results_dict_allplots,n_subnets,subnet_array,idx_sorted,output_folder,plot_best_x,series_labels,color_palette=:Dark2_8)

end


text:
using Documenter, XCB

makedocs(;
    modules=[XCB],
    format=Documenter.HTML(prettyurls = true),
    pages=[
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "API" => "api.md",
        "Troubleshooting" => "troubleshooting.md",
        "Developer documentation" => [
            "WindowAbstractions Interface" => "interface.md",
            "Utility tools" => "utility.md",
        ]
    ],
    repo="https://github.com/JuliaGL/XCB.jl/blob/{commit}{path}#L{line}",
    sitename="XCB.jl",
    authors="serenity4 <cedric.bel@hotmail.fr>, SimonDanisch"
)

deploydocs(
    repo = "github.com/JuliaGL/XCB.jl.git",
)


text:
module NEC

abstract type StructuralElement end
const Structure = Vector{StructuralElement}
abstract type AbstractAntenna end

export
    Wire,

    Excitation,
    VoltageSource,
    LinearWave,

    AbstractAntenna,
    Yagi,
    Moxon,
    Dipole,

    RadiationPattern,
    f2b

export
    get_lambda,
    divide

export
    AWG1,
    AWG2,
    AWG3,
    AWG4,
    AWG5,
    AWG6,
    AWG7,
    AWG8,
    AWG9,
    AWG10,
    AWG11,
    AWG12,
    AWG13,
    AWG14,
    AWG15,
    AWG16,
    AWG17,
    AWG18,
    AWG19,
    AWG20

# functions I might need
include("math.jl")

# structural elements
include("wire.jl")

# basic radiation things
include("excitation.jl")
include("radiation_pattern.jl")

# pre-made antenna types
include("yagi.jl")
include("moxon.jl")
include("dipole.jl")

include("backend.jl")

end # module


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

"""
top_level_mk_body: lower domain nodes in the AST to parfor nodes.
"""
function top_level_mk_body(ast::Array{Any,1}, depth, state)
    len  = length(ast)
    body = Any[]
    fuse_number = 1

    # Process the top-level expressions of a function and do fusion and useless assignment elimination.
    for i = 1:len
        # Record the top-level statement number in the processing state.
        state.top_level_number = i
        @dprintln(2,"Processing top-level ast #",i," depth=",depth)

        # Convert the current expression.
        new_exprs = from_expr(ast[i], depth, state, true)
        assert(isa(new_exprs,Array))
        # If conversion of current statement resulted in anything.
        if length(new_exprs) != 0
            append!(body, new_exprs)
        else
            @dprintln(3,"This statement was filtered: ", ast[i])
        end
    end
    return body
end

"""
fusion_pass: perform fusion of adjacent parfor nodes in the AST.
"""
function fusion_pass(ast::Array{Any,1}, depth, state)
    len  = length(ast)
    body = Any[]
    fuse_number = 1

    # Process the top-level expressions of a function and do fusion and useless assignment elimination.
    for i = 1:len
        # Record the top-level statement number in the processing state.
        state.top_level_number = i
        @dprintln(2,"Processing top-level ast #",i," depth=",depth)

        # If this isn't the first statement...
        if length(body) != 0
            last_node = body[end]  # Get the last statement inserted.
            @dprintln(3, "Should fuse?")
            @dprintln(3, "new = ", ast[i])
            @dprintln(3, "last = ", last_node)

            # See if the previous expression is a parfor.
            is_last_parfor = isParforAssignmentNode(last_node) || isBareParfor(last_node)
            # See if the new expression is a parfor.
            is_new_parfor  = isParforAssignmentNode(ast[i])    || isBareParfor(ast[i])
            @dprintln(3,"is_new_parfor = ", is_new_parfor, " is_last_parfor = ", is_last_parfor)

            # If both are parfors then try to fuse them.
            if is_new_parfor && is_last_parfor
                @dprintln(3,"Starting fusion ", fuse_number)
                fuse_number = fuse_number + 1
                fuse_ret = fuse(body, length(body), ast[i], state)
                if fuse_ret>0
                    # 2 means combination of old and new parfors has no output and both are dead
                    if fuse_ret==2
                        # remove last parfor and don't add anything new
                        pop!(body)
                    end
                    # If fused then the new node is consumed and no new node is added to the body.
                    continue
                end
            end

            push!(body, ast[i])
        else
            push!(body, ast[i])
        end
    end
    return body
end

"""
fusion_pass: perform fusion on the AST but callable as an external pass.  Just takes the
function name and function lambda.
"""
function fusion_pass(function_name::String, LambdaVarInfo, body::Expr)
    #LambdaVarInfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast)
    @assert body.head==:body "invalid body expr head"
    max_label = getMaxLabel(0, body.args)
    lives = computeLiveness(body, LambdaVarInfo)
    input_arrays = getArrayParams(LambdaVarInfo)
    state = expr_state(function_name, lives, max_label, input_arrays)
    genEquivalenceClasses(LambdaVarInfo, body, state)
    @dprintln(3,"new fusion lambda = ")
    printLambda(3, LambdaVarInfo, body)
    @dprintln(3,"new fusion equivalence classes")
    print_correlations(3, state)
    body = fusion_pass(body.args, 1, state)
    body = Expr(:body, body...)
    # return CompilerTools.LambdaHandling.LambdaVarInfoToLambda(LambdaVarInfo, body, ParallelAccelerator.ParallelIR.AstWalk)
    return LambdaVarInfo, body
end

# Remove the pre-statements from parfor nodes and expand them into the top-level expression array.
function top_level_expand_pre(body, state)

    expanded_body = Any[]

    for i = 1:length(body)
        if isParforAssignmentNode(body[i])
            parfor_assignment = body[i]
            @dprintln(3,"Expanding a parfor assignment node = ", parfor_assignment)

            the_parfor = getParforNode(parfor_assignment)
            @dprintln(3,"the_parfor = ", the_parfor)
            lhs = getLhsFromAssignment(parfor_assignment)
            rhs = getRhsFromAssignment(parfor_assignment)

            # Add all the pre-parfor statements to the expanded body.
            append!(expanded_body, the_parfor.preParFor)
            append!(expanded_body, the_parfor.hoisted)
            the_parfor.preParFor = Any[]
            the_parfor.hoisted = Any[]

            # Add just the parfor to the expanded body.  The post-parfor part removed below.
            assert(typeof(rhs) == Expr)
            rhs.typ = typeof(0) # a few lines down you can see that the last post-statement of 0 is added.
            push!(expanded_body, rhs)

            # All the post parfor part to the expanded body.
            # The regular part of the post parfor is just added.
            # The part that indicates the return values creates individual assignment statements for each thing returned.
            if isFusionAssignment(parfor_assignment)
                append!(expanded_body, the_parfor.postParFor[1:end-1])
                for j = 4:length(parfor_assignment.args)
                    push!(expanded_body, mk_assignment_expr(parfor_assignment.args[j], the_parfor.postParFor[end][j-3], state))
                end
            else
                append!(expanded_body, the_parfor.postParFor[1:end-1])
                push!(expanded_body, mk_assignment_expr(lhs, the_parfor.postParFor[end], state))
            end
            the_parfor.postParFor = Any[]
            push!(the_parfor.postParFor, 0)
            #createInstructionCountEstimate(the_parfor, state)
        elseif isBareParfor(body[i])
            rhs = body[i]
            the_parfor = rhs.args[1]

            # Add all the pre-parfor statements to the expanded body.
            append!(expanded_body, the_parfor.preParFor)
            append!(expanded_body, the_parfor.hoisted)
            the_parfor.preParFor = Any[]
            the_parfor.hoisted= Any[]

            # Add just the parfor to the expanded body.  The post-parfor part removed below.
            assert(typeof(rhs) == Expr)
            rhs.typ = typeof(0) # a few lines down you can see that the last post-statement of 0 is added.
            push!(expanded_body, rhs)

            the_parfor.postParFor = Any[]
            push!(the_parfor.postParFor, 0)
            #createInstructionCountEstimate(the_parfor, state)
        else
            push!(expanded_body, body[i])
        end
    end
    return expanded_body
end

function top_level_mk_task_graph(body, state, new_lives, loop_info)
    task_start = time_ns()

    rr = ReplacedRegion[]
    # TODO: another pass of alias analysis to re-use dead but uniquely allocated arrays
    # AliasAnalysis.analyze_lambda_body(fake_body, state.param, state.meta2_typed, new_lives)

    # Create a mapping between the minimized basic block numbers and indices in body that are in those correponding basic blocks.
    map_reduced_bb_num_to_body = Dict{Int,Array{Int,1}}()
    for i = 1:length(body)
        # Get the basic block number for the first top level number associated with this entry in body.
        bb_num = CompilerTools.LivenessAnalysis.find_bb_for_statement(i, new_lives)
        if bb_num == nothing
            if typeof(body[i]) != LabelNode
                @dprintln(0,"statement that couldn't be found in liveness analysis ", body[i])
                throw(string("find_bb_for_statement should not fail for non-LabelNodes"))
            end
            continue
        end
        # If not previously in the map then initialize it with the current body index.  Otherwise, add the current body index
        # as also mapping to its containing basic block.
        if !haskey(map_reduced_bb_num_to_body, bb_num)
            map_reduced_bb_num_to_body[bb_num] = [i]
        else
            map_reduced_bb_num_to_body[bb_num] = [map_reduced_bb_num_to_body[bb_num]; i]
        end
    end

    @dprintln(3,"map_reduced_bb_num_to_body = ", map_reduced_bb_num_to_body)

    bbs_in_task_graph_loops = Set()
    bbs = new_lives.basic_blocks
    #tgsections = TaskGraphSection[]

    if put_loops_in_task_graph
        # For each loop that loop analysis identified.
        for one_loop in loop_info.loops
            # If the loop is "simple", i.e., has no just a body and a back-edge and no conditionals or nested loops.
            if length(one_loop.members) == 2
                # This is sort of sanity checking because the way Julia creates loops, these conditions should probably always hold.
                head_bb = bbs[one_loop.head]
                back_bb = bbs[one_loop.back_edge]

                if length(head_bb.preds) == 2 &&
                    length(head_bb.succs) == 1 &&
                    length(back_bb.preds) == 1 &&
                    length(back_bb.succs) == 2
                    before_head = getNonBlock(head_bb.preds, one_loop.back_edge)
                    assert(typeof(before_head) == CompilerTools.LivenessAnalysis.BasicBlock)
                    @dprintln(3,"head_bb.preds = ", head_bb.preds, " one_loop.back_edge = ", one_loop.back_edge, " before_head = ", before_head)
                    # assert(length(before_head) == 1)
                    after_back  = getNonBlock(back_bb.succs, one_loop.head)
                    assert(typeof(after_back) == CompilerTools.LivenessAnalysis.BasicBlock)
                    #assert(length(after_back) == 1)

                    head_indices = map_reduced_bb_num_to_body[one_loop.head]
                    head_first_parfor = nothing
                    for j = 1:length(head_indices)
                        if isParforAssignmentNode(body[head_indices[j]]) || isBareParfor(body[head_indices[j]])
                            head_first_parfor = j
                            break
                        end
                    end

                    back_indices = map_reduced_bb_num_to_body[one_loop.back_edge]
                    back_first_parfor = nothing
                    for j = 1:length(back_indices)
                        if isParforAssignmentNode(body[back_indices[j]]) || isBareParfor(body[back_indices[j]])
                            back_first_parfor = j
                            break
                        end
                    end

                    if head_first_parfor != nothing || back_first_parfor != nothing
                        new_bbs_for_set = Set(one_loop.head, one_loop.back_edge, before_head.label, after_back.label)
                        assert(length(intersect(bbs_in_task_graph_loops, new_bbs_for_set)) == 0)
                        bbs_in_task_graph_loops = union(bbs_in_task_graph_loops, new_bbs_for_set)

                        before_indices = map_reduced_bb_num_to_body[before_head.label]
                        before_first_parfor = nothing
                        for j = 1:length(before_indices)
                            if isParforAssignmentNode(body[before_indices[j]]) || isBareParfor(body[body_indices[j]])
                                before_first_parfor = j
                                break
                            end
                        end

                        after_indices = map_reduced_bb_num_to_body[after_back.label]
                        after_first_parfor = nothing
                        for j = 1:length(after_indices)
                            if isParforAssignmentNode(body[after_indices[j]]) || isBareParfor(body[after_indices[j]])
                                after_first_parfor = j
                                break
                            end
                        end

                        #          bb_live_info = new_lives.basic_blocks[bb_num]
                        #          push!(replaced_regions, (first_parfor, last_parfor, makeTasks(first_parfor, last_parfor, body, bb_live_info)))

                    end
                else
                    @dprintln(1,"Found a loop with 2 members but unexpected head or back_edge structure.")
                    @dprintln(1,"head = ", head_bb)
                    @dprintln(1,"back_edge = ", back_bb)
                end
            end
        end
    end

    for i in map_reduced_bb_num_to_body
        bb_num = i[1]
        body_indices = i[2]
        bb_live_info = new_lives.basic_blocks[new_lives.cfg.basic_blocks[bb_num]]

        if !in(bb_num, bbs_in_task_graph_loops)
            if task_graph_mode == SEQUENTIAL_TASKS
                # Find the first parfor in the block.
                first_parfor = nothing
                for j = 1:length(body_indices)
                    if isParforAssignmentNode(body[body_indices[j]]) || isBareParfor(body[body_indices[j]])
                        first_parfor = body_indices[j]
                        break
                    end
                end

                # If we found a parfor in the block.
                if first_parfor != nothing
                    # Find the last parfor in the block...it might be the same as the first.
                    last_parfor = nothing
                    for j = length(body_indices):-1:1
                        if isParforAssignmentNode(body[body_indices[j]]) || isBareParfor(body[body_indices[j]])
                            last_parfor = body_indices[j]
                            break
                        end
                    end
                    assert(last_parfor != nothing)

                    # Remember this section of code as something to transform into task graph format.
                    #push!(tgsections, TaskGraphSection(first_parfor, last_parfor, body[first_parfor:last_parfor]))
                    #@dprintln(3,"Adding TaskGraphSection ", tgsections[end])

                    push!(rr, ReplacedRegion(first_parfor, last_parfor, bb_num, makeTasks(first_parfor, last_parfor, body, bb_live_info, state, task_graph_mode)))
                end
            elseif task_graph_mode == ONE_AT_A_TIME
                for j = 1:length(body_indices)
                    if taskableParfor(body[body_indices[j]])
                        # Remember this section of code as something to transform into task graph format.
                        cur_start = cur_end = body_indices[j]
                        #push!(tgsections, TaskGraphSection(cur_start, cur_end, body[cur_start:cur_end]))
                        #@dprintln(3,"Adding TaskGraphSection ", tgsections[end])

                        push!(rr, ReplacedRegion(body_indices[j], body_indices[j], bb_num, makeTasks(cur_start, cur_end, body, bb_live_info, state, task_graph_mode)))
                    end
                end
            elseif task_graph_mode == MULTI_PARFOR_SEQ_NO
                cur_start = nothing
                cur_end   = nothing
                stmts_in_batch = Int64[]

                for j = 1:length(body_indices)
                    if taskableParfor(body[body_indices[j]])
                        if cur_start == nothing
                            cur_start = cur_end = body_indices[j]
                        else
                            cur_end = body_indices[j]
                        end
                        push!(stmts_in_batch, body_indices[j])
                    else
                        if cur_start != nothing
                            @dprintln(3,"Non-taskable parfor ", stmts_in_batch, " ", body[body_indices[j]])
                            in_vars, out, locals = io_of_stmts_in_batch = getIO(stmts_in_batch, bb_live_info.statements)
                            @dprintln(3,"in_vars = ", in_vars)
                            @dprintln(3,"out_vars = ", out)
                            @dprintln(3,"local_vars = ", locals)

                            cur_in_vars, cur_out, cur_locals = io_of_stmts_in_batch = getIO([body_indices[j]], bb_live_info.statements)
                            @dprintln(3,"cur_in_vars = ", cur_in_vars)
                            @dprintln(3,"cur_out_vars = ", cur_out)
                            @dprintln(3,"cur_local_vars = ", cur_locals)

                            if isempty(intersect(out, cur_in_vars))
                                @dprintln(3,"Sequential statement doesn't conflict with batch.")
                                push!(stmts_in_batch, body_indices[j])
                                cur_end = body_indices[j]
                            else
                                # Remember this section of code (excluding current statement) as something to transform into task graph format.
                                #push!(tgsections, TaskGraphSection(cur_start, cur_end, body[cur_start:cur_end]))
                                #@dprintln(3,"Adding TaskGraphSection ", tgsections[end])

                                push!(rr, ReplacedRegion(cur_start, cur_end, bb_num, makeTasks(cur_start, cur_end, body, bb_live_info, state, task_graph_mode)))

                                cur_start = cur_end = nothing
                                stmts_in_batch = Int64[]
                            end
                        end
                    end
                end

                if cur_start != nothing
                    # Remember this section of code as something to transform into task graph format.
                    #push!(tgsections, TaskGraphSection(cur_start, cur_end, body[cur_start:cur_end]))
                    #@dprintln(3,"Adding TaskGraphSection ", tgsections[end])

                    push!(rr, ReplacedRegion(cur_start, cur_end, bb_num, makeTasks(cur_start, cur_end, body, bb_live_info, state, task_graph_mode)))
                end
            else
                throw(string("Unknown Parallel IR task graph formation mode."))
            end
        end
    end

    @dprintln(3,"Regions prior to sorting.")
    @dprintln(3,rr)
    # We replace regions in reverse order of index so that we don't mess up indices that we need to replace later.
    sort!(rr, by=x -> x.end_index, rev=true)
    @dprintln(3,"Regions after sorting.")
    @dprintln(3,rr)

    printBody(3,body)

    @dprintln(2, "replaced_regions")
    for i = 1:length(rr)
        @dprintln(2, rr[i])

        if ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE
            # New body starts with the pre-task graph portion
            new_body = body[1:rr[i].start_index-1]

            # Then adds calls for each task
            for j = 1:length(rr[i].tasks)
                cur_task = rr[i].tasks[j]
                @dprintln(3,"cur_task = ", cur_task, " type = ", typeof(cur_task))

                process_cur_task(cur_task, new_body, state)
            end

            # Then appends the post-task graph portion
            append!(new_body, body[rr[i].end_index+1:end])

            body = new_body

            @dprintln(3,"new_body after region ", i, " replaced")
            printBody(3,body)
        elseif ParallelAccelerator.client_intel_task_graph
            # new body starts with the pre-task graph portion
            new_body = body[1:rr[i].start_index-1]

            # then adds calls for each task
            for j = 1:length(rr[i].tasks)
                cur_task = rr[i].tasks[j]
                @dprintln(3,"cur_task = ", cur_task, " type = ", typeof(cur_task))
                if typeof(cur_task) == TaskInfo
                    @dprintln(3,"Inserting call to insert_divisible_task")
                    @dprintln(3,cur_task.function_sym, " type = ", typeof(cur_task.function_sym))

                    in_len  = length(cur_task.input_symbols)
                    mod_len = length(cur_task.modified_inputs)
                    io_len  = length(cur_task.io_symbols)
                    red_len = length(cur_task.reduction_vars)
                    @dprintln(3, "inputs, modifieds, io_sym, reductions = ", cur_task.input_symbols, " ", cur_task.modified_inputs, " ", cur_task.io_symbols, " ", cur_task.reduction_vars)

                    dims = length(cur_task.loopNests)
                    if dims > 0
                        itn = InsertTaskNode()

                        if ParallelAccelerator.client_intel_task_graph_mode == 0
                            itn.task_options = TASK_STATIC_SCHEDULER | TASK_AFFINITY_XEON
                        elseif ParallelAccelerator.client_intel_task_graph_mode == 1
                            itn.task_options = TASK_STATIC_SCHEDULER | TASK_AFFINITY_PHI
                        elseif ParallelAccelerator.client_intel_task_graph_mode == 2
                            itn.task_options = 0
                        else
                            throw(string("Unknown task graph mode option ", ParallelAccelerator.client_intel_task_graph_mode))
                        end

                        if task_graph_mode == SEQUENTIAL_TASKS
                            # intentionally do nothing
                        elseif task_graph_mode == ONE_AT_A_TIME
                            itn.task_options |= TASK_FINISH
                        elseif task_graph_mode == MULTI_PARFOR_SEQ_NO
                            # intentionally do nothing
                        else
                            throw(string("Unknown task_graph_mode."))
                        end

                        # Fill in pir_range
                        # Fill in the min and max grain sizes.
                        itn.ranges.dim = dims
                        itn.host_grain_size.dim = dims
                        itn.phi_grain_size.dim = dims
                        for l = 1:dims
                            # Note that loopNest is outer-dimension first
                            push!(itn.ranges.lower_bounds, Domain_IR.sub_expr(cur_task.loopNests[dims - l + 1].lower, 1))
                            push!(itn.ranges.upper_bounds, Domain_IR.sub_expr(cur_task.loopNests[dims - l + 1].upper, 1))
                            push!(itn.host_grain_size.sizes, 2)
                            push!(itn.phi_grain_size.sizes, 2)
                        end

                        # Fill in the arg metadata.
                        for l = 1:in_len
                            if isArrayType(cur_task.input_symbols[l].typ)
                                @dprintln(3, "is array")
                                push!(itn.args, pir_arg_metadata(cur_task.input_symbols[l], ARG_OPT_IN, create_array_access_desc(cur_task.input_symbols[l])))
                            else
                                @dprintln(3, "is not array")
                                push!(itn.args, pir_arg_metadata(cur_task.input_symbols[l], ARG_OPT_IN))
                            end
                        end
                        for l = 1:mod_len
                            @dprintln(3, "mod_len loop: ", l, " ", cur_task.modified_inputs[l])
                            if isArrayType(cur_task.modified_inputs[l].typ)
                                @dprintln(3, "is array")
                                push!(itn.args, pir_arg_metadata(cur_task.modified_inputs[l], ARG_OPT_OUT, create_array_access_desc(cur_task.modified_inputs[l])))
                            else
                                @dprintln(3, "is not array")
                                push!(itn.args, pir_arg_metadata(cur_task.modified_inputs[l], ARG_OPT_OUT))
                            end
                        end
                        for l = 1:io_len
                            @dprintln(3, "io_len loop: ", l, " ", cur_task.io_symbols[l])
                            if isArrayType(cur_task.io_symbols[l].typ)
                                @dprintln(3, "is array")
                                push!(itn.args, pir_arg_metadata(cur_task.io_symbols[l], ARG_OPT_INOUT, create_array_access_desc(cur_task.io_symbols[l])))
                            else
                                @dprintln(3, "is not array")
                                push!(itn.args, pir_arg_metadata(cur_task.io_symbols[l], ARG_OPT_INOUT))
                            end
                        end
                        for l = 1:red_len
                            @dprintln(3, "red_len loop: ", l, " ", cur_task.reduction_vars[l])
                            push!(itn.args, pir_arg_metadata(cur_task.reduction_vars[l], ARG_OPT_ACCUMULATOR))
                        end

                        # Fill in the task function.
                        itn.task_func = cur_task.task_func
                        itn.join_func = cur_task.join_func

                        @dprintln(3,"InsertTaskNode = ", itn)

                        insert_task_expr = TypedExpr(Int, :insert_divisible_task, itn)
                        push!(new_body, insert_task_expr)
                    else
                        throw(string("insert sequential task not implemented yet"))
                    end
                else
                    push!(new_body, cur_task)
                end
            end

            # Insert call to wait on the scheduler to complete all tasks.
            #push!(new_body, TypedExpr(Cint, :call, TopNode(:ccall), QuoteNode(:pert_wait_all_task), Type{Cint}, ()))

            # Then appends the post-task graph portion
            append!(new_body, body[rr[i].end_index+1:end])

            body = new_body

            @dprintln(3,"new_body after region ", i, " replaced")
            printBody(3,body)
        elseif run_as_task_decrement()
            # new body starts with the pre-task graph portion
            new_body = body[1:rr[i].start_index-1]

            # then adds calls for each task
            for j = 1:length(rr[i].tasks)
                cur_task = rr[i].tasks[j]
                assert(typeof(cur_task) == TaskInfo)

                @dprintln(3,"Inserting call to function")
                @dprintln(3,cur_task, " type = ", typeof(cur_task))
                @dprintln(3,cur_task.function_sym, " type = ", typeof(cur_task.function_sym))

                in_len = length(cur_task.input_symbols)
                out_len = length(cur_task.output_symbols)

                real_out_params = Any[]

                for k = 1:out_len
                    this_param = cur_task.output_symbols[k]
                    atype = Array{getType(this_param, state.LambdaVarInfo), 1}
                    temp_param_array = createStateVar(state, string(this_param.name, "_out_array"), atype, ISASSIGNED)
                    push!(real_out_params, temp_param_array)
                    new_temp_array = mk_alloc_array_expr(this_param.typ, atype, 1)
                    push!(new_body, mk_assignment_expr(temp_param_array, new_temp_array, state))
                end

                #push!(new_body, TypedExpr(Void, :call, cur_task.function_sym, cur_task.input_symbols..., real_out_params...))
                #push!(new_body, TypedExpr(Void, :call, TopNode(cur_task.function_sym), cur_task.input_symbols..., real_out_params...))
                push!(new_body, TypedExpr(Void, :call, mk_parallelir_ref(cur_task.function_sym), TypedExpr(pir_range_actual, :call, :pir_range_actual), cur_task.input_symbols..., real_out_params...))

                for k = 1:out_len
                    push!(new_body, mk_assignment_expr(cur_task.output_symbols[k], mk_arrayref1(real_out_params[k], 1, false, state), state))
                end
            end

            # then appends the post-task graph portion
            append!(new_body, body[rr[i].end_index+1:end])

            body = new_body

            @dprintln(3,"new_body after region ", i, " replaced")
            printBody(3,body)
        end
    end

    #    throw(string("debugging task graph"))
    if print_times
    @dprintln(1,"Task formation time = ", ns_to_sec(time_ns() - task_start))
    end

    return body
end

function recreateFromLoophead(new_body, stmt :: Expr, LoopEndDict :: Dict{LHSVar,Array{Any,1}}, state, newLambdaVarInfo)
    # Only handle 1D loophead right now.

    if stmt.head == :loophead
        assert(length(stmt.args) == 3)
        loop_id    = toLHSVar(stmt.args[1])
        loop_start = stmt.args[2]
        loop_end   = stmt.args[3]
        @dprintln(3, "recreateFromLoophead ", loop_id, " ", loop_start, " ", loop_end)
    else
        assert(length(stmt.args) == 1)
        loop_id = toLHSVar(stmt.args[1])
        if haskey(LoopEndDict, loop_id)
            append!(new_body, LoopEndDict[loop_id])
        else
            throw(string("loophead information missing for loop ", loop_id))
        end
        return nothing
    end

    assert(isa(loop_id, LHSVar))
    uniq = get_unique_num()

    if VERSION >= v"0.5.0-dev+4449"
        label_head = next_label(state)
        label_end  = next_label(state)

        num_vars = 6
        parfor_nest_level = 1
        loop_nest_level = 1

        steprange_last_var  = Symbol(string("#recreate_steprange_last_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 0))
        steprange_first_var = Symbol(string("#recreate_steprange_first_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 1))
        steprange_step_var  = Symbol(string("#recreate_steprange_step_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 2))
        recreate_temp_var   = Symbol(string("#recreate_temp_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 3))
        recreate_ssa5_var   = Symbol(string("#recreate_ssa5_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 4))
        recreate_ssa6_var   = Symbol(string("#recreate_ssa6_", uniq, "_", parfor_nest_level, "_", (loop_nest_level-1) * num_vars + 5))

        steprange_last_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(steprange_last_var  , Int64, ISASSIGNED, newLambdaVarInfo)
        steprange_first_rhsvar = CompilerTools.LambdaHandling.addLocalVariable(steprange_first_var , Int64, ISASSIGNED, newLambdaVarInfo)
        steprange_step_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(steprange_step_var  , Int64, ISASSIGNED, newLambdaVarInfo)
        recreate_temp_rhsvar   = CompilerTools.LambdaHandling.addLocalVariable(recreate_temp_var   , Int64, ISASSIGNED, newLambdaVarInfo)
        recreate_ssa5_rhsvar   = CompilerTools.LambdaHandling.addLocalVariable(recreate_ssa5_var   , Int64, ISASSIGNED, newLambdaVarInfo)
        recreate_ssa6_rhsvar   = CompilerTools.LambdaHandling.addLocalVariable(recreate_ssa6_var   , Int64, ISASSIGNED, newLambdaVarInfo)

        steprange_last_lhsvar  = toLHSVar(steprange_last_rhsvar)
        steprange_first_lhsvar = toLHSVar(steprange_first_rhsvar)
        steprange_step_lhsvar  = toLHSVar(steprange_step_rhsvar)
        recreate_temp_lhsvar   = toLHSVar(recreate_temp_rhsvar)
        recreate_ssa5_lhsvar   = toLHSVar(recreate_ssa5_rhsvar)
        recreate_ssa6_lhsvar   = toLHSVar(recreate_ssa6_rhsvar)

        push!(new_body, mk_assignment_expr(deepcopy(steprange_last_lhsvar), Expr(:call, GlobalRef(Base,:steprange_last), deepcopy(loop_start), 1, deepcopy(loop_end)), newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(deepcopy(steprange_first_lhsvar), deepcopy(loop_start), newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(deepcopy(steprange_step_lhsvar), 1, newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(deepcopy(recreate_temp_lhsvar), deepcopy(steprange_first_rhsvar), newLambdaVarInfo))

#        push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "steprange_last_var = ", deepcopy(steprange_last_rhsvar)))
#        push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "steprange_first_var = ", deepcopy(steprange_first_rhsvar)))
#        push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "steprange_step_var = ", deepcopy(steprange_step_rhsvar)))

        push!(new_body, LabelNode(label_head))

#        push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "after label_head"))

        push!(new_body, mk_gotoifnot_expr(
               boxOrNot(GlobalRef(Base, :Bool),
                   Expr(:call, GlobalRef(Base, :not_int),
                           Expr(:call, GlobalRef(Base, :or_int),
                                   Expr(:call, GlobalRef(Base, :and_int),
                                           Expr(:call, GlobalRef(Base, :not_int),
                                               Expr(:call, GlobalRef(Base, :(===)), deepcopy(steprange_first_rhsvar), deepcopy(steprange_last_rhsvar))
                                           ),
                                           Expr(:call, GlobalRef(Base, :not_int),
                                               Expr(:call, GlobalRef(Base, :(===)),
                                                   Expr(:call, GlobalRef(Base, :slt_int),
                                                       0,
                                                       deepcopy(steprange_step_rhsvar)
                                                   ),
                                                   Expr(:call, GlobalRef(Base, :slt_int),
                                                       deepcopy(steprange_first_rhsvar),
                                                       deepcopy(steprange_last_rhsvar)
                                                   )
                                               )
                                           )
                                       ),
                               Expr(:call, GlobalRef(Base, :(===)), deepcopy(recreate_temp_rhsvar),
                                   boxOrNot(Int64, Expr(:call, GlobalRef(Base, :add_int), deepcopy(steprange_last_rhsvar), deepcopy(steprange_step_rhsvar)))
                               )
                           )
                   )
               )
               , label_end))

        push!(new_body, mk_assignment_expr(deepcopy(recreate_ssa5_lhsvar), deepcopy(recreate_temp_rhsvar), newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(deepcopy(recreate_ssa6_lhsvar), boxOrNot(Int64, Expr(:call, GlobalRef(Base, :add_int), deepcopy(recreate_temp_rhsvar), deepcopy(steprange_step_rhsvar))), newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(CompilerTools.LambdaHandling.toLHSVar(deepcopy(loop_id), newLambdaVarInfo), deepcopy(recreate_ssa5_rhsvar), newLambdaVarInfo))
        push!(new_body, mk_assignment_expr(deepcopy(recreate_temp_lhsvar), deepcopy(recreate_ssa6_rhsvar), newLambdaVarInfo))

        for_loop_end = Any[]

        push!(for_loop_end, GotoNode(label_head))
        push!(for_loop_end, LabelNode(label_end))

        LoopEndDict[loop_id] = for_loop_end
    else

        label_after_first_unless   = next_label(state)
    #    label_before_second_unless = next_label(state)
        label_after_second_unless  = next_label(state)
    #    label_last                 = next_label(state)

        gensym2_var = string("#recreate_gensym2_", uniq, "_", 0)
        gensym2_sym = Symbol(gensym2_var)
        gensym0_var = string("#recreate_gensym0_", uniq, "_", 1)
        gensym0_sym = Symbol(gensym0_var)
        pound_s1_var = string("#recreate_pound_s1_", uniq, "_", 2)
        pound_s1_sym = Symbol(pound_s1_var)
        gensym3_var = string("#recreate_gensym3_", uniq, "_", 3)
        gensym3_sym = Symbol(gensym3_var)
        gensym4_var = string("#recreate_gensym4_", uniq, "_", 4)
        gensym4_sym = Symbol(gensym4_var)
        gensym2_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(gensym2_sym, Int64, ISASSIGNED, newLambdaVarInfo)
        gensym0_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(gensym0_sym, StepRange{Int64,Int64}, ISASSIGNED, newLambdaVarInfo)
        pound_s1_rhsvar = CompilerTools.LambdaHandling.addLocalVariable(pound_s1_sym, Int64, ISASSIGNED, newLambdaVarInfo)
        gensym3_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(gensym3_sym, Int64, ISASSIGNED, newLambdaVarInfo)
        gensym4_rhsvar  = CompilerTools.LambdaHandling.addLocalVariable(gensym4_sym, Int64, ISASSIGNED, newLambdaVarInfo)

        gensym2_lhsvar  = toLHSVar(gensym2_rhsvar)
        gensym0_lhsvar  = toLHSVar(gensym0_rhsvar)
        pound_s1_lhsvar = toLHSVar(pound_s1_rhsvar)
        gensym3_lhsvar  = toLHSVar(gensym3_rhsvar)
        gensym4_lhsvar  = toLHSVar(gensym4_rhsvar)

        #push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "ranges = ", toRHSVar(:ranges, pir_range_actual, state.LambdaVarInfo)))
        #push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "this_nest.lower = ", this_nest.lower))
        #push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "this_nest.step  = ", this_nest.step))
        #push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "this_nest.upper = ", this_nest.upper))

        push!(new_body, mk_assignment_expr(deepcopy(gensym2_lhsvar), TypedExpr(Int64, :call, GlobalRef(Base,:steprange_last), loop_start, 1, loop_end), state))
        push!(new_body, mk_assignment_expr(deepcopy(gensym0_lhsvar), TypedExpr(StepRange{Int64,Int64}, :new, StepRange{Int64,Int64}, loop_start, 1, deepcopy(gensym2_rhsvar)), state))
        push!(new_body, mk_assignment_expr(deepcopy(pound_s1_lhsvar), TypedExpr(Int64, :call, GlobalRef(Base, :getfield), deepcopy(gensym0_rhsvar), QuoteNode(:start)), state))
        push!(new_body, mk_gotoifnot_expr(TypedExpr(Bool, :call, GlobalRef(ParallelAccelerator.ParallelIR, :first_unless), deepcopy(gensym0_rhsvar), deepcopy(pound_s1_rhsvar)), label_after_second_unless))
        push!(new_body, LabelNode(label_after_first_unless))

    #       push!(new_body, Expr(:call, GlobalRef(Base,:println), GlobalRef(Base,:STDOUT), " in label_after_first_unless section"))

        push!(new_body, mk_assignment_expr(deepcopy(gensym3_lhsvar), deepcopy(pound_s1_rhsvar), state))
        push!(new_body, mk_assignment_expr(deepcopy(gensym4_lhsvar), TypedExpr(Int64, :call, GlobalRef(ParallelAccelerator.ParallelIR, :assign_gs4), deepcopy(gensym0_rhsvar), deepcopy(pound_s1_rhsvar)), state))
        push!(new_body, mk_assignment_expr(deepcopy(loop_id), deepcopy(gensym3_rhsvar), state))
        push!(new_body, mk_assignment_expr(deepcopy(pound_s1_lhsvar), deepcopy(gensym4_rhsvar), state))

        for_loop_end = Any[]

    #    push!(for_loop_end, LabelNode(label_before_second_unless))
        push!(for_loop_end, mk_gotoifnot_expr(TypedExpr(Bool, :call, GlobalRef(ParallelAccelerator.ParallelIR, :second_unless), deepcopy(gensym0_rhsvar), deepcopy(pound_s1_rhsvar)), label_after_first_unless))
        push!(for_loop_end, LabelNode(label_after_second_unless))
    #    push!(for_loop_end, LabelNode(label_last))

        LoopEndDict[loop_id] = for_loop_end
    end

    return nothing
end

#function got_here_1(x)
#println("got here! ", x)
#end

run_nested_task_directly = true

function red_alloc(elem_type, length, value)
    elem_type[deepcopy(value) for i = 1:length]
end

function process_cur_task(cur_task::TaskInfo, new_body, state)
    @dprintln(3,"process_cur_task = ", cur_task.function_sym, " type = ", typeof(cur_task.function_sym), " function = ", state.function_name)
    range_var = string(Base.function_name(cur_task.task_func),"_range_var")
    range_sym = Symbol(range_var)
    range_rhsvar = CompilerTools.LambdaHandling.addLocalVariable(range_sym, pir_range_actual, CompilerTools.LambdaHandling.ISASSIGNED, state.LambdaVarInfo)
    range_lhsvar = toLHSVar(range_rhsvar)

    run_task_directly = haskey(ENV,"PROSPECT_RUN_TASK_DIRECTLY") || (state.in_nested && run_nested_task_directly)
    @dprintln(3,"run_task_directly = ", run_task_directly, " state.in_nested = ", state.in_nested, " run_nested_task_directly = ", run_nested_task_directly)

    if run_task_directly
        @dprintln(3,"Inserting call to task function directly ", range_sym, " ", range_rhsvar)
    elseif haskey(ENV,"PROSPECT_CALL_ISF")
        @dprintln(3,"Inserting call to isf ", range_sym, " ", range_rhsvar)
#    elseif haskey(ENV,"PROSPECT_CALL_DYNAMIC_CHECK")
#        @dprintln(3,"Inserting call dynamic check ", range_sym, " ", range_rhsvar)
    else
        @dprintln(3,"Inserting call dynamic check ", range_sym, " ", range_rhsvar)
#        @dprintln(3,"Inserting call to jl_threading_run in function ", state.function_name, " ", range_sym, " ", range_rhsvar)
    end
    @dprintln(3,"Cur LambdaVarInfo = ", state.LambdaVarInfo)

    in_len  = length(cur_task.input_symbols)
    mod_len = length(cur_task.modified_inputs)
    io_len  = length(cur_task.io_symbols)
    red_len = length(cur_task.reduction_vars)
    @dprintln(3, "inputs, modifieds, io_sym, reductions = ", cur_task.input_symbols, " ", cur_task.modified_inputs, " ", cur_task.io_symbols, " ", cur_task.reduction_vars)

    dims = length(cur_task.loopNests)
    if dims > 0
        @dprintln(3,"dims > 0")
        assert(dims <= 3)
        cstr_params = Any[]
        for l = 1:dims
            # Note that loopNest is outer-dimension first
            # Should this still be outer-dimension first?  FIX FIX FIX
            push!(cstr_params, cur_task.loopNests[dims - l + 1].lower)
            push!(cstr_params, cur_task.loopNests[dims - l + 1].upper)
        end
        @dprintln(3, "cstr_params = ", cstr_params)
        cstr_expr = ParallelAccelerator.ParallelIR.pir_range_actual
        #cstr_expr = mk_parallelir_ref(:pir_range_actual, Any)
        whole_range_expr = mk_assignment_expr(deepcopy(range_lhsvar), TypedExpr(pir_range_actual, :call, cstr_expr, cstr_params...), state)
        @dprintln(3,"whole_range_expr = ", whole_range_expr)
        push!(new_body, whole_range_expr)

        #    push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "whole_range = ", range_rhsvar))

        real_args_build = Any[]
        args_type = Expr(:tuple)
        red_var_rhsvars = Any[]

        # Fill in the arg metadata.
        for l = 1:in_len
            @dprintln(3,cur_task.input_symbols[l].name, " ", typeof(cur_task.input_symbols[l].name))
            push!(real_args_build, cur_task.input_symbols[l].name)
            push!(args_type.args,  cur_task.input_symbols[l].typ)
        end
        for l = 1:mod_len
            @dprintln(3,cur_task.modified_inputs[l].name, " ", typeof(cur_task.modified_inputs[l].name))
            push!(real_args_build, cur_task.modified_inputs[l].name)
            push!(args_type.args,  cur_task.modified_inputs[l].typ)
        end
        for l = 1:io_len
            @dprintln(3,cur_task.io_symbols[l].name, " ", typeof(cur_task.io_symbols[l].name))
            push!(real_args_build, cur_task.io_symbols[l].name)
            push!(args_type.args,  cur_task.io_symbols[l].typ)
        end
        for l = 1:red_len
            push!(red_var_rhsvars, cur_task.reduction_vars[l].name)
            push!(real_args_build, cur_task.reduction_vars[l].name)
            push!(args_type.args,  cur_task.reduction_vars[l].typ)
        end

        @dprintln(3,"task_func = ", cur_task.task_func)
        #@dprintln(3,"whole_iteration_range = ", whole_iteration_range)
        @dprintln(3,"real_args_build = ", real_args_build)

        if cur_task.ret_types != Void
            red_array_lhsvar = createVarForTask(cur_task, "_red_array", Array{Any,1}, state)
#            if haskey(ENV,"PROSPECT_CALL_DYNAMIC_CHECK")
                push!(new_body, mk_assignment_expr(deepcopy(red_array_lhsvar), Expr(:call, GlobalRef(ParallelAccelerator.ParallelIR, :red_alloc), Any, TypedExpr(Int, :call, GlobalRef(Base.Threads, :nthreads)), Expr(:call, GlobalRef(Core, :tuple), red_var_rhsvars...)), state))
#            else
#                push!(new_body, mk_assignment_expr(deepcopy(red_array_lhsvar), mk_alloc_array_expr(Any, Array{Any,1}, TypedExpr(Int, :call, GlobalRef(Base.Threads, :nthreads))), state))
#            end

            red_output_tuple_typ = Tuple{cur_task.ret_types...}
            red_output_tuple_lhsvar = createVarForTask(cur_task, "_red_output_tuple", red_output_tuple_typ, state)
        end

        #push!(new_body, Expr(:call, GlobalRef(ParallelAccelerator.ParallelIR, :got_here_1), TypedExpr(pir_range_actual, :call, cstr_expr, cstr_params...)))
        task_globalref = GlobalRef(ParallelAccelerator.ParallelIR, Base.function_name(cur_task.task_func))
        #push!(new_body, Expr(:call, GlobalRef(ParallelAccelerator.ParallelIR, :got_here_1), range_lhsvar))
        if run_task_directly
          if cur_task.ret_types != Void
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...,
                                         deepcopy(red_array_lhsvar))
          else
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...)
          end
        elseif haskey(ENV,"PROSPECT_CALL_ISF")
          if cur_task.ret_types != Void
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         ParallelAccelerator.ParallelIR.isf,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...,
                                         deepcopy(red_array_lhsvar))
          else
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         ParallelAccelerator.ParallelIR.isf,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...)
          end
        #elseif haskey(ENV,"PROSPECT_CALL_DYNAMIC_CHECK")
        else
          if cur_task.ret_types != Void
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         ParallelAccelerator.ParallelIR.dynamic_check,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...,
                                         deepcopy(red_array_lhsvar))
          else
            insert_task_expr = TypedExpr(Void,
                                         :call,
                                         ParallelAccelerator.ParallelIR.dynamic_check,
                                         task_globalref,
                                         deepcopy(range_lhsvar),
                                         real_args_build...)
          end
#############   This code won't prevent nested parallelism.  When nested parallelism works then this code can be re-enabled.
#############   In the meantime, in the above section, we use a function that uses a global to make sure that only function is
#############   run in threaded mode at a time.
#        else
#            tup_var = string(Base.function_name(cur_task.task_func),"_tup_var")
#            tup_sym = Symbol(tup_var)
#            tup_rhsvar = CompilerTools.LambdaHandling.addLocalVariable(tup_sym, SimpleVector, CompilerTools.LambdaHandling.ISASSIGNED, state.LambdaVarInfo)
#            tup_lhsvar = toLHSVar(tup_rhsvar)
#
#            if cur_task.ret_types != Void
#                push!(new_body, mk_assignment_expr(deepcopy(tup_lhsvar), mk_svec_expr(GlobalRef(ParallelAccelerator.ParallelIR, :isf), task_globalref, deepcopy(range_lhsvar), real_args_build..., deepcopy(red_array_lhsvar)), state))
#            else
#                push!(new_body, mk_assignment_expr(deepcopy(tup_lhsvar), mk_svec_expr(GlobalRef(ParallelAccelerator.ParallelIR, :isf), task_globalref, deepcopy(range_lhsvar), real_args_build...), state))
#            end
#
#            svec_args = mk_svec_expr(Any)
#            insert_task_expr = TypedExpr(Void,
#                                         :call,
#                                         GlobalRef(Core, :ccall),
#                                         QuoteNode(:jl_threading_run),
#                                         GlobalRef(Base.Threads,:Void),
#                                         svec_args,
#                                         tup_lhsvar, 0)
        end
        push!(new_body, insert_task_expr)

        if cur_task.ret_types != Void
            if run_task_directly || haskey(ENV,"PROSPECT_CALL_ISF")
                for l = 1:red_len
                    push!(new_body,
                      mk_assignment_expr(
                        deepcopy(cur_task.reduction_vars[l].name),
                        TypedExpr(cur_task.ret_types[l], :call, GlobalRef(Core, :typeassert),
                            Expr(:call, GlobalRef(Base, :getfield), TypedExpr(red_output_tuple_typ, :call, GlobalRef(Base, :arrayref), deepcopy(red_array_lhsvar), TypedExpr(Int, :call, GlobalRef(Base.Threads,:threadid))), l), cur_task.ret_types[l]), state))
#                    push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "index threadid of red_array = ", deepcopy(cur_task.reduction_vars[l].name)))
                end
            else
                red_loop_index_lhsvar = createVarForTask(cur_task, "_reduction_loop_index", Int, state)
                unique_node_id = get_unique_num()
                #push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "whole reduction array = ", deepcopy(red_array_lhsvar)))

                # After the jl_threading_run call, we store the first element of the reduction array into their destinations.
                for l = 1:red_len
                    push!(new_body, mk_assignment_expr(
                                       deepcopy(cur_task.reduction_vars[l].name),
                                       TypedExpr(cur_task.ret_types[l], :call, GlobalRef(Core, :typeassert),
                                         Expr(:call, GlobalRef(Base, :getfield), TypedExpr(red_output_tuple_typ, :call, GlobalRef(Base, :arrayref), deepcopy(red_array_lhsvar), 1), l), cur_task.ret_types[l]), state))
#                    push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "index 1 of red_array = ", deepcopy(cur_task.reduction_vars[l].name)))
                end

                push!(new_body, Expr(:loophead, deepcopy(red_loop_index_lhsvar), 2, TypedExpr(Int, :call, GlobalRef(Base, :arraylen), deepcopy(red_array_lhsvar))))
                for l = 1:red_len
                    @dprintln(3, "reduction_vars[l] = ", cur_task.reduction_vars[l])
                    reduction_code = callDelayedFuncWith(cur_task.join_func[l].reductionFunc,
                                                    deepcopy(cur_task.reduction_vars[l].name),
                                                    TypedExpr(cur_task.ret_types[l], :call, GlobalRef(Core, :typeassert),
                                                      Expr(:call, GlobalRef(Base, :getfield), TypedExpr(red_output_tuple_typ, :call, GlobalRef(Base, :arrayref), deepcopy(red_array_lhsvar), deepcopy(red_loop_index_lhsvar)), l), cur_task.ret_types[l]))
                    @dprintln(3, "Adding reduction code for reduction variable ", l, " with ", length(reduction_code), " statements.")

                    for stmt in reduction_code
                        @dprintln(3, "reduction stmt = ", stmt)
                        if isBareParfor(stmt)
                            @dprintln(3, "reduction stmt is parfor so recreating loops")
                            recreateLoops(new_body, stmt.args[1], state, state.LambdaVarInfo)
                        else
                            @dprintln(3, "reduction stmt is normal")
                            push!(new_body, stmt)
                        end
                    end
#                    push!(new_body, TypedExpr(Any, :call, :println, GlobalRef(Base,:STDOUT), "next reduction from red_array = ", deepcopy(cur_task.reduction_vars[l].name)))
                end

                push!(new_body, Expr(:loopend, deepcopy(red_loop_index_lhsvar)))

                @dprintln(3,"LambdaVarInfo after reduction code = ", state.LambdaVarInfo)
            end
        end
    else
        throw(string("insert sequential task not implemented yet"))
    end
end

function createVarForTask(cur_task, suffix, typ, state)
    new_var_name = string(Base.function_name(cur_task.task_func),suffix)
    new_var_sym = Symbol(new_var_name)
    new_var_rhsvar = CompilerTools.LambdaHandling.addLocalVariable(new_var_sym, typ, CompilerTools.LambdaHandling.ISASSIGNED, state.LambdaVarInfo)
    return toLHSVar(new_var_rhsvar)
end

function run_reduction_func(red_array :: Array{Any,1}, join_func)
    @dprintln(3,"run_reduction_fun ", red_array, " ", length(red_array), " ", join_func, " type = ", typeof(join_func))
    if length(red_array) == 1
        return red_array[1]
    else
        return foldl((x,y) -> join_func(x,y), red_array)
    end
end

function process_cur_task(cur_task::Any, new_body, state)
    push!(new_body, cur_task)
end

# TOP_LEVEL
# sequence of expressions
# ast = [ expr, ... ]
function top_level_from_exprs(ast::Array{Any,1}, depth, state)
    # Perform lowering from domain nodes to parfors.
    main_proc_start = time_ns()
    body = top_level_mk_body(ast, depth, state)
    @dprintln(3,"Body after lowering but before fusion.")
    printBody(3, body)
    if print_times
    @dprintln(1,"Lowering time = ", ns_to_sec(time_ns() - main_proc_start))
    end

    # Perform fusion.
    fusion_start = time_ns()
    body = fusion_pass(body, depth, state)
    @dprintln(3,"Body after first pass before task graph creation.")
    printBody(3, body)
    if print_times
    @dprintln(1,"Fusion time = ", ns_to_sec(time_ns() - fusion_start))
    end

    # TASK GRAPH

    expand_start = time_ns()

    body = top_level_expand_pre(body, state)

    if print_times
    @dprintln(1,"Expanding parfors time = ", ns_to_sec(time_ns() - expand_start))
    end

    @dprintln(3,"expanded_body = ")
    @dprintln(3,"lvi = ", state.LambdaVarInfo)
    printBody(3, body)

    @dprintln(3,"LambdaVarInfo = ", state.LambdaVarInfo)
    new_lives = CompilerTools.LivenessAnalysis.from_lambda(state.LambdaVarInfo, body, pir_live_cb, state.LambdaVarInfo)
    @dprintln(1,"Starting loop analysis.")
    loop_info = CompilerTools.Loops.compute_dom_loops(new_lives.cfg)
    @dprintln(1,"Finished loop analysis.")

    if hoist_allocation == 1
        body = hoistAllocation(body, new_lives, loop_info, state)
        new_lives = CompilerTools.LivenessAnalysis.from_lambda(state.LambdaVarInfo, body, pir_live_cb, state.LambdaVarInfo)
        @dprintln(1,"Starting loop analysis again.")
        loop_info = CompilerTools.Loops.compute_dom_loops(new_lives.cfg)
        @dprintln(1,"Finished loop analysis.")
    end

    @dprintln(3,"new_lives = ", new_lives)
    @dprintln(3,"loop_info = ", loop_info)

    if ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE || ParallelAccelerator.getTaskMode() > 0 || run_as_task()
        @dprintln(3, "Entering top_level_mk_task_graph.")
        body = top_level_mk_task_graph(body, state, new_lives, loop_info)
    end  # end of task graph formation section

    @dprintln(3,"body after mk_task_graph = ")
    printBody(3, body)

    expanded_body = Any[]

    max_label   = getMaxLabel(0, body)
    LoopEndDict = Dict{LHSVar,Array{Any,1}}()

    for i = 1:length(body)
        @dprintln(3,"Flatten index ", i, " ", body[i], " type = ", typeof(body[i]))
            # Convert loophead and loopend into Julia loops.
        if ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE &&
            typeof(body[i]) == Expr &&
            (body[i].head == :loophead || body[i].head == :loopend)
            recreateFromLoophead(expanded_body, body[i], LoopEndDict, state, state.LambdaVarInfo)
        else
            push!(expanded_body, body[i])
        end
    end

    body = expanded_body
    # CompilerTools.OptFramework.cleanupBodyLabels(expanded_body)
    return body
end


text:
using CUDA
using Random

function training_runs(individuals, number_individuals, brains, environments; env_seed::Int, number_rounds::Int, threads::Int, blocks::Int, shared_memory::Int)

    rewards = CUDA.fill(0.0f0, number_individuals)
    individuals_gpu = CUDA.CuArray(individuals)
    environment_seeds = CUDA.fill(env_seed, number_individuals)
    CUDA.@cuda threads = threads blocks = blocks shmem = shared_memory kernel_eval_fitness(
        individuals_gpu ,
        rewards,
        environment_seeds,
        number_rounds,
        brains,
        environments
    )

    CUDA.synchronize()

    return Array(rewards)

end

function validation_runs(individual, individual_size, number_validation_runs, brains, environments; threads::Int, blocks::Int, shared_memory::Int)

    rewards = CUDA.fill(0.0f0, number_validation_runs)

    individuals = fill(0.0f0, number_validation_runs, individual_size)

    # Use same individual for each validation run
    for i = 1:number_validation_runs
        for j = 1:individual_size
            individuals[i, j] = individual[j]
        end
    end

    individuals_gpu = CuArray(individuals)
    environment_seeds = CuArray(1:number_validation_runs)
    number_rounds = 1
    CUDA.@cuda threads = threads blocks = blocks shmem = shared_memory kernel_eval_fitness(
        individuals_gpu,
        rewards,
        environment_seeds,
        number_rounds,
        brains,
        environments
    )

    CUDA.synchronize()

    return Array(rewards)

end