text:
export editcost
export editcosts
export editerrs
export alignops
export alignerrs

"""
    editcost(tar, src, subcost=1)

Minimum edit distance discribes how many operations to transform `src` into `tar`.
Editing operations are `Insertion(i)`, `Deletion(d)` and `Substitution(s)`, e.g.\n
    src: I N T E * N T I O N
    tar: * E X E C U T I O N
         d s s   i s          (Note: the transformation direction is src -> tar)
`subcost`=2 is so called Levenshtein distance. `subcost`=1 is often used in speech recognition tasks.
"""
function editcost(tar::V, src::V, subcost=1) where V
    T = length(tar)+1
    S = length(src)+1
    d = zeros(Int, T, S)
    d[1,2:S] = 1:S-1
    d[2:T,1] = 1:T-1
    for t = 2:T
        for s = 2:S
            c = ifelse(src[s-1]!=tar[t-1], subcost, 0)
            d[t,s] = min(d[t-1,s]+1, d[t,s-1]+1, d[t-1,s-1]+c)
            # -----------ins---------del---------sub----------
        end
    end
    return d[T,S]
end

"""
    editcosts(tar, src, subcost=1) -> d, ins,del,sub

`d` = `ins` + `del` + `sub`, where `d` is the total edit distance from `src` to `tar`,
`ins` is insertion operations, `del` is deletion operations and `sub` is substitution
operations. `subcost`=2 is so called Levenshtein distance. e.g.\n
    src: I N T E * N T I O N
    tar: * E X E C U T I O N
         d s s   i s          (Note: the transformation direction is src -> tar)

# Example
    julia> d,ins,del,sub = editcosts("ABCDEfghi","abcABCdE")
    (8, 4, 3, 1)
"""
function editcosts(tar::V, src::V, subcost=1) where V
    T = length(tar)+1
    S = length(src)+1
    d = zeros(Int, T, S)
    row = zeros(Int, T, S)
    col = zeros(Int, T, S)

    # if tar=[], Deletions on src
    for t=1, s=2:S
        d[t,s] = s - 1
        row[t,s] = t
        col[t,s] = s - 1
    end
    # if src=[], Insertions on src
    for t=2:T, s=1
        d[t,s] = t - 1
        row[t,s] = t - 1
        col[t,s] = s
    end

    # forward
    for t = 2:T
        for s = 2:S
            c = ifelse(src[s-1]!=tar[t-1], subcost, 0)
            way3 = [d[t-1,s]+1, d[t,s-1]+1, d[t-1,s-1]+c]
            #-------ins---------del---------sub----------
            indx = argmin(way3)
            d[t,s] = way3[indx]
            if indx==1
                row[t,s], col[t,s] = t-1, s   # ins
            elseif indx==2
                row[t,s], col[t,s] = t, s-1   # del
            elseif indx==3
                row[t,s], col[t,s] = t-1, s-1 # sub
            end
        end
    end

    # trace-back
    trace = [];
    push!(trace, (T, S))
    push!(trace, (row[T,S], col[T,S]))
    while trace[end] != (0,0)
        t,s = trace[end]
        push!(trace, (row[t,s], col[t,s]))
    end

    ins, del, sub = 0, 0, 0
    for (t,s) in trace[end-1:-1:1]
        r, c = row[t,s], col[t,s]
        if r==0 || c==0 continue end
        𝕥 = t!=1 ? t-1 : t  # edge case
        𝕤 = s!=1 ? s-1 : s  # edge case
        if d[t,s]==d[r,c]
            continue
        elseif d[t,s]==d[𝕥,s]+1
            ins += 1
        elseif d[t,s]==d[t,𝕤]+1
            del += 1
        elseif d[t,s]==d[r,c]+subcost
            sub += 1
        end
    end
    return d[T,S], ins,del,sub
end

"""
    editerrs(correct, estimated, subcost=1) -> d, ins,del,sub

`d` = `ins` + `del` + `sub`, where `d` is the total edit errors made from editing
`correct` to `estimated`, `ins` is insertion(i) errors, `del` is deletion(d) errors
and `sub` is substitution(s) errors, e.g.\n
    correct   : A * C D
    estimated : a b C
                s i   d  (Note: the edit direction is correct -> estimated)
`subcost`=1 is mostly used in speech transcription task to calculate word error rate.
`subcost`=2 is called Levenshtein distance.
# Example
    julia> d,ins,del,sub = editerrs("ABCDEfghi","abcABCdE")
    (8, 3, 4, 1)
    julia> d,ins,del,sub = editcosts("ABCDEfghi","abcABCdE")
    (8, 4, 3, 1)
"""
function editerrs(correct::V, estimated::V, subcost=1) where V
    d, ins,del,sub = editcosts(estimated, correct, subcost)
    return d, ins,del,sub
end

"""
    alignops(tar, src; subcost=1, show=true, cn=false) -> marks::Array{String,2}

align `src` and `tar`, `subcost` is the substitution cost, `show` decides whether
to display, `cn` stands for chinese. `marks` is the aligned results recorded in a
3×N Array{String,2}, of which the 1st row is `tar`, 2nd is `src` and 3rd is the
operations to transform `src` into `tar`. `subcost`=2 is Levenshtein distance.

# Example
    julia> r = alignops(
           ["are","are","you","ok","bu"], # tar
           ["ah","r","are","u","OKAY"]);  # src
    ───────────────────
    are     ah      sub
            r       del
    are     are
    you     u       sub
    ok      OKAY    sub
    bu              ins
"""
function alignops(tar::V, src::V; subcost=1, show=true, cn=false) where V
    T = length(tar) + 1
    S = length(src) + 1
    d =   zeros(Int, T, S)
    row = zeros(Int, T, S)
    col = zeros(Int, T, S)

    Tmax = 0
    Smax = 0
    Tlen = zeros(Int, T-1)
    Slen = zeros(Int, S-1)
    for t=1, s=2:S
        𝕤 = s - 1
        d[t,s] = 𝕤
        row[t,s] = t
        col[t,s] = 𝕤
        Slen[𝕤] = length(src[𝕤])
        if Smax < Slen[𝕤] Smax=Slen[𝕤] end
    end
    for t=2:T, s=1
        𝕥 = t - 1
        d[t,s] = 𝕥
        row[t,s] = 𝕥
        col[t,s] = s
        Tlen[𝕥] = length(tar[𝕥])
        if Tmax < Tlen[𝕥] Tmax=Tlen[𝕥] end
    end

    # forward
    for t = 2:T
        for s = 2:S
            c = ifelse(src[s-1]!=tar[t-1], subcost, 0)
            way3 = [d[t-1,s]+1, d[t,s-1]+1, d[t-1,s-1]+c]
            #-------ins---------del---------sub----------
            indx = argmin(way3)
            d[t,s] = way3[indx]
            if indx==1
                row[t,s], col[t,s] = t-1, s   # ins
            elseif indx==2
                row[t,s], col[t,s] = t, s-1   # del
            elseif indx==3
                row[t,s], col[t,s] = t-1, s-1 # sub
            end
        end
    end
    # trace-back
    trace = [];
    push!(trace, (T, S))
    push!(trace, (row[T,S], col[T,S]))
    while trace[end] != (0,0)
        t,s = trace[end]
        push!(trace, (row[t,s], col[t,s]))
    end

    steps = length(trace)-2
    marks = Array{String,2}(undef, 3, steps)
    INS = cn ? "插" : "ins"
    DEL = cn ? "删" : "del"
    SUB = cn ? "替" : "sub"
    NIL = cn ? "无" : "nil"
    _V_ = "\t"

    cn ? println(repeat("─",floor(Int,(Tmax+Smax+12)*3/2))) : println(repeat("─",Tmax+Smax+12))
    for (i,(t,s)) in enumerate(trace[end-1:-1:1])
        r, c = row[t,s], col[t,s]
        if r==0 || c==0 continue end
        𝕥 = t!=1 ? t-1 : t  # edge case
        𝕤 = s!=1 ? s-1 : s  # edge case
        𝕚 = i - 1
        if d[t,s]==d[r,c]
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = NIL
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = src[𝕤] * repeat(" ", Smax - Slen[𝕤])
            show && println(TAR, _V_, SRC) # Correct
        elseif d[t,s]==d[𝕥,s]+1
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = ""
            marks[3,𝕚] = INS
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = repeat(" ", Smax)
            show && println(TAR, _V_, SRC, _V_, "ins"|>blue!) # Insertion
        elseif d[t,s]==d[t,𝕤]+1
            marks[1,𝕚] = ""
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = DEL
            TAR = repeat(" ", Tmax)
            SRC = src[𝕤] * repeat(" ", Smax - length(src[𝕤]))
            show && println(TAR, _V_, SRC, _V_, "del"|>red!) # Deletion
        elseif d[t,s]==d[r,c]+subcost
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = SUB
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = src[𝕤] * repeat(" ", Smax - Slen[𝕤])
            show && println(TAR, _V_, SRC, _V_, "sub"|>yellow!) # Substitution
        end
    end
    return marks
end

"""
    alignerrs(tar, src; subcost=1, show=true, cn=false) -> marks::Array{String,2}

align `src`(estimated) and `tar`(desired), `subcost` is the substitution cost, `show`
decides whether to display, `cn` stands for chinese. `marks` is the aligned results
recorded in a 3×N Array{String,2}, of which the 1st row is `tar`, the 2nd is `src` and
the 3rd is the errors to align `src` to `tar`. `subcost`=2 is Levenshtein distance.

# Example
    julia> r = alignerrs(
           ["are","are","you","ok","bu"],
           ["ah" ,"r"  ,"are","u" ,"OKAY"]);
    ───────────────────
    are     ah      sub
            r       ins
    are     are
    you     u       sub
    ok      OKAY    sub
    bu              del
"""
function alignerrs(tar::V, src::V; subcost=1, show=true, cn=false) where V
    T = length(tar) + 1
    S = length(src) + 1
    d =   zeros(Int, T, S)
    row = zeros(Int, T, S)
    col = zeros(Int, T, S)

    Tmax = 0
    Smax = 0
    Tlen = zeros(Int, T-1)
    Slen = zeros(Int, S-1)
    for t=1, s=2:S
        𝕤 = s - 1
        d[t,s] = 𝕤
        row[t,s] = t
        col[t,s] = 𝕤
        Slen[𝕤] = length(src[𝕤])
        if Smax < Slen[𝕤] Smax=Slen[𝕤] end
    end
    for t=2:T, s=1
        𝕥 = t - 1
        d[t,s] = 𝕥
        row[t,s] = 𝕥
        col[t,s] = s
        Tlen[𝕥] = length(tar[𝕥])
        if Tmax < Tlen[𝕥] Tmax=Tlen[𝕥] end
    end

    # forward
    for t = 2:T
        for s = 2:S
            c = ifelse(src[s-1]!=tar[t-1], subcost, 0)
            way3 = [d[t-1,s]+1, d[t,s-1]+1, d[t-1,s-1]+c]
            #-------ins---------del---------sub----------
            indx = argmin(way3)
            d[t,s] = way3[indx]
            if indx==1
                row[t,s], col[t,s] = t-1, s   # ins op on src
            elseif indx==2
                row[t,s], col[t,s] = t, s-1   # del op on src
            elseif indx==3
                row[t,s], col[t,s] = t-1, s-1 # sub op on src
            end
        end
    end
    # trace-back
    trace = [];
    push!(trace, (T, S))
    push!(trace, (row[T,S], col[T,S]))
    while trace[end] != (0,0)
        t,s = trace[end]
        push!(trace, (row[t,s], col[t,s]))
    end

    # ops on src to transfer into tar, can be viewed as the
    # error edits that leads tar turning into src, so ins/del op
    # from src to tar are equivalent to del/ins op from tar to src
    steps = length(trace)-2
    marks = Array{String,2}(undef, 3, steps)
    INS = cn ? "删" : "del"
    DEL = cn ? "插" : "ins"
    SUB = cn ? "替" : "sub"
    NIL = cn ? "无" : "nil"
    _V_ = "\t"

    cn ? println(repeat("─",floor(Int,(Tmax+Smax+12)*3/2))) : println(repeat("─",Tmax+Smax+12))
    for (i,(t,s)) in enumerate(trace[end-1:-1:1])
        r, c = row[t,s], col[t,s]
        if r==0 || c==0 continue end
        𝕥 = t!=1 ? t-1 : t  # edge case
        𝕤 = s!=1 ? s-1 : s  # edge case
        𝕚 = i - 1
        if d[t,s]==d[r,c]
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = NIL
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = src[𝕤] * repeat(" ", Smax - Slen[𝕤])
            show && println(TAR, _V_, SRC) # Correct
        elseif d[t,s]==d[𝕥,s]+1
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = ""
            marks[3,𝕚] = INS
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = repeat(" ", Smax)
            show && println(TAR, _V_, SRC, _V_, "del"|>red!) # Deletion errors
        elseif d[t,s]==d[t,𝕤]+1
            marks[1,𝕚] = ""
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = DEL
            TAR = repeat(" ", Tmax)
            SRC = src[𝕤] * repeat(" ", Smax - length(src[𝕤]))
            show && println(TAR, _V_, SRC, _V_, "ins"|>blue!) # Insertion errors
        elseif d[t,s]==d[r,c]+subcost
            marks[1,𝕚] = tar[𝕥]
            marks[2,𝕚] = src[𝕤]
            marks[3,𝕚] = SUB
            TAR = tar[𝕥] * repeat(" ", Tmax - Tlen[𝕥])
            SRC = src[𝕤] * repeat(" ", Smax - Slen[𝕤])
            show && println(TAR, _V_, SRC, _V_, "sub"|>yellow!) # Substitution errors
        end
    end
    return marks
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Autogenerated wrapper script for Glib_jll for x86_64-apple-darwin
export libgio, libglib, libgmodule, libgobject, libgthread

using Libiconv_jll
using Libffi_jll
using Gettext_jll
using PCRE_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("Glib")
JLLWrappers.@declare_library_product(libgio, "@rpath/libgio-2.0.0.dylib")
JLLWrappers.@declare_library_product(libglib, "@rpath/libglib-2.0.0.dylib")
JLLWrappers.@declare_library_product(libgmodule, "@rpath/libgmodule-2.0.0.dylib")
JLLWrappers.@declare_library_product(libgobject, "@rpath/libgobject-2.0.0.dylib")
JLLWrappers.@declare_library_product(libgthread, "@rpath/libgthread-2.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Libffi_jll, Gettext_jll, PCRE_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libgio,
        "lib/libgio-2.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libglib,
        "lib/libglib-2.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgmodule,
        "lib/libgmodule-2.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgobject,
        "lib/libgobject-2.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgthread,
        "lib/libgthread-2.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
module TestDay13

import AOC2021.Day13.part1, AOC2021.Day13.part2
using Test

test_input = """
6,10
0,14
9,10
0,3
10,4
4,11
6,0
6,12
4,1
0,13
10,12
3,4
3,0
8,4
1,10
2,14
8,10
9,0

fold along y=7
fold along x=5"""

@testset "Day 13" begin
    @testset "part 1" begin
        @test part1(test_input) == 17
    end
    # @testset "part 2" begin
    # @test part2(test_input) == true
    # end
end

end # module


text:
#
# Boundary MPS contraction methods
# Approximately contract a 2D tensor network with boundary MPS.
# The boundary MPS will be used as projectors to be inserted
# into the network.
#

struct ContractionAlgorithm{algorithm} end

ContractionAlgorithm(s::AbstractString) = ContractionAlgorithm{Symbol(s)}()

macro ContractionAlgorithm_str(s)
  return :(ContractionAlgorithm{$(Expr(:quote, Symbol(s)))})
end

Base.@kwdef struct BoundaryMPS
  top::Vector{MPS} = MPS[]
  bottom::Vector{MPS} = MPS[]
  left::Vector{MPS} = MPS[]
  right::Vector{MPS} = MPS[]
end

struct BoundaryMPSDir{dir} end

BoundaryMPSDir(s::AbstractString) = BoundaryMPSDir{Symbol(s)}()

macro BoundaryMPSDir_str(s)
  return :(BoundaryMPSDir{$(Expr(:quote, Symbol(s)))})
end

function contract_approx(
  tn::Matrix{ITensor},
  alg::ContractionAlgorithm"boundary_mps",
  dirs::BoundaryMPSDir"top_to_bottom";
  cutoff,
  maxdim,
)
  nrows, ncols = size(tn)
  boundary_mps = Vector{MPS}(undef, nrows - 1)
  x = MPS(tn[1, :])
  boundary_mps[1] = orthogonalize(x, ncols)
  for nrow in 2:(nrows - 1)
    A = MPO(tn[nrow, :])
    x = contract(A, x; cutoff=cutoff, maxdim=maxdim)
    boundary_mps[nrow] = orthogonalize(x, ncols)
  end
  return BoundaryMPS(; top=boundary_mps)
end

function contract_approx(
  tn::Matrix{ITensor},
  alg::ContractionAlgorithm"boundary_mps",
  dirs::BoundaryMPSDir"bottom_to_top";
  cutoff,
  maxdim,
)
  tn = rot180(tn)
  tn = reverse(tn; dims=2)
  boundary_mps_top = contract_approx(
    tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim
  )
  return BoundaryMPS(; bottom=reverse(boundary_mps_top.top))
end

function contract_approx(
  tn::Matrix{ITensor},
  alg::ContractionAlgorithm"boundary_mps",
  dirs::BoundaryMPSDir"left_to_right";
  cutoff,
  maxdim,
)
  tn = rotr90(tn)
  boundary_mps = contract_approx(
    tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim
  )
  return BoundaryMPS(; left=boundary_mps.top)
end

function contract_approx(
  tn::Matrix{ITensor},
  alg::ContractionAlgorithm"boundary_mps",
  dirs::BoundaryMPSDir"right_to_left";
  cutoff,
  maxdim,
)
  tn = rotr90(tn)
  boundary_mps = contract_approx(
    tn, alg, BoundaryMPSDir"bottom_to_top"(); cutoff=cutoff, maxdim=maxdim
  )
  return BoundaryMPS(; right=boundary_mps.bottom)
end

function contract_approx(
  tn::Matrix{ITensor},
  alg::ContractionAlgorithm"boundary_mps",
  dirs::BoundaryMPSDir"all";
  cutoff,
  maxdim,
)
  mps_top = contract_approx(
    tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim
  )
  mps_bottom = contract_approx(
    tn, alg, BoundaryMPSDir"bottom_to_top"(); cutoff=cutoff, maxdim=maxdim
  )
  mps_left = contract_approx(
    tn, alg, BoundaryMPSDir"left_to_right"(); cutoff=cutoff, maxdim=maxdim
  )
  mps_right = contract_approx(
    tn, alg, BoundaryMPSDir"right_to_left"(); cutoff=cutoff, maxdim=maxdim
  )
  return BoundaryMPS(;
    top=mps_top.top, bottom=mps_bottom.bottom, left=mps_left.left, right=mps_right.right
  )
end

function contract_approx(
  tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps"; dirs, cutoff, maxdim
)
  return contract_approx(tn, alg, BoundaryMPSDir(dirs); cutoff=cutoff, maxdim=maxdim)
end

# Compute the truncation projectors for the network,
# contracting from top to bottom
function contract_approx(
  tn::Matrix{ITensor}; alg="boundary_mps", dirs="all", cutoff=1e-8, maxdim=maxdim_arg(tn)
)
  return contract_approx(
    tn, ContractionAlgorithm(alg); dirs=dirs, cutoff=cutoff, maxdim=maxdim
  )
end

function contraction_cache_top(tn, boundary_mps::Vector{MPS})
  tn_cache_top = Vector{typeof(tn)}(undef, length(boundary_mps))
  for n in 1:length(boundary_mps)
    tn_cache_top[n] = contraction_cache_top(tn, boundary_mps, n)
  end
  return tn_cache_top
end

function contraction_cache_bottom(tn, boundary_mps::Vector{MPS})
  tn_top = rot180(tn)
  tn_top = reverse(tn_top; dims=2)
  boundary_mps_top = reverse(boundary_mps)
  tn_cache_top = contraction_cache_top(tn_top, boundary_mps_top)
  tn_cache = reverse.(tn_cache_top; dims=2)
  tn_cache = rot180.(tn_cache)
  return tn_cache
end

function contraction_cache_left(tn, boundary_mps::Vector{MPS})
  tn_top = rotr90(tn)
  tn_cache_top = contraction_cache_top(tn_top, boundary_mps)
  tn_cache = rotl90.(tn_cache_top)
  return tn_cache
end

function contraction_cache_right(tn, boundary_mps::Vector{MPS})
  tn_bottom = rotr90(tn)
  tn_cache_bottom = contraction_cache_bottom(tn_bottom, boundary_mps)
  tn_cache = rotl90.(tn_cache_bottom)
  return tn_cache
end

function contraction_cache(tn, boundary_mps)
  cache_top = contraction_cache_top(tn, boundary_mps.top)
  cache_bottom = contraction_cache_bottom(tn, boundary_mps.bottom)
  cache_left = contraction_cache_left(tn, boundary_mps.left)
  cache_right = contraction_cache_right(tn, boundary_mps.right)
  return (top=cache_top, bottom=cache_bottom, left=cache_left, right=cache_right)
end

function boundary_mps_top(tn)
  return [MPS(tn[n][n, :]) for n in 1:length(tn)]
end

function boundary_mps_bottom(tn)
  tn_top = rot180.(tn)
  tn_top = reverse.(tn_top; dims=2)
  boundary_mps = boundary_mps_top(tn_top)
  return reverse(boundary_mps)
end

function boundary_mps_left(tn)
  tn_top = rotr90.(tn)
  return boundary_mps_top(tn_top)
end

function boundary_mps_right(tn)
  tn_bottom = rotr90.(tn)
  return boundary_mps_bottom(tn_bottom)
end

function boundary_mps(tn::NamedTuple)
  _boundary_mps_top = boundary_mps_top(tn.top)
  _boundary_mps_bottom = boundary_mps_bottom(tn.bottom)
  _boundary_mps_left = boundary_mps_left(tn.left)
  _boundary_mps_right = boundary_mps_right(tn.right)
  return BoundaryMPS(;
    top=_boundary_mps_top,
    bottom=_boundary_mps_bottom,
    left=_boundary_mps_left,
    right=_boundary_mps_right,
  )
end

function boundary_mps(tn::Matrix{ITensor}; cutoff, maxdim)
  #TODO
  # Contract in every direction
  combiner_gauge = combiners(linkinds, tn)
  tnᶜ = insert_gauge(tn, combiner_gauge)
  boundary_mpsᶜ = contract_approx(tnᶜ; maxdim=maxdim, cutoff=cutoff)
  tn_cacheᶜ = contraction_cache(tnᶜ, boundary_mpsᶜ)
  tn_cache = insert_gauge(tn_cacheᶜ, combiner_gauge)
  return boundary_mps(tn_cache)
end


text:
## Bayesian Linear Regression
## ==========================

using Sigma
using Distributions

Σ_1 = 2 * eye(5,5)
# Zero mean
μ = MvNormal(zeros(5), Σ_1)
Σ_2 = eye(5,5)
Σ_prior = Wishart(5,Σ_2)
w = MvNormal(rand(μ), rand(Σ_prior))
x = [Uniform(-1,1) for j = 1:5, i = 1:500]
τ = Gamma(0.5,2)
ϵ = [Normal(0,1/rand(τ)) for i = 1:500]
y = [rand(w)' * rand(x[i,j]) + ϵ[i] for j=1:5, i = 1:500]

dist(w, y = x[])

# Nubmer of random variables is
num_rvs = 5*5 + 5*500 + 1 * 500


text:
# Required packages

using Vann
using PyPlot
using DataFrames

# Model choices

snow_choice = TinBasic
hydro_choice = Gr4j

# Load data

path_inputs = Pkg.dir("Vann", "data", "roykenes")

date, tair, prec, q_obs, frac = load_data(path_inputs)

q_obs[q_obs .< 0] = NaN

# Compute potential evapotranspiration

epot = epot_zero(date)

# Initilize model

tstep = 3.0

time = date[1]

st_snow = eval(Expr(:call, snow_choice, tstep, time, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time))

# Run calibration

param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs)

println(param_snow)
println(param_hydro)

# Reinitilize model

st_snow = eval(Expr(:call, snow_choice, tstep, time, param_snow, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time, param_hydro))

# Run model with best parameter set

q_sim = run_model(st_snow, st_hydro, tair, prec, epot)

# Store results in data frame

q_obs = round(q_obs, 2)
q_sim = round(q_sim, 2)

df_res = DataFrame(date = date, q_sim = q_sim, q_obs = q_obs)

# Compute performance

perf_nse = nse(df_res[:q_sim], df_res[:q_obs])
perf_nse = round(perf_nse, 2)

perf_kge = kge(df_res[:q_sim], df_res[:q_obs])
perf_kge = round(perf_kge, 2)

# Plot results using rcode

fig = plt[:figure](figsize = (12,7))

plt[:style][:use]("ggplot")

plt[:plot](df_res[:date], df_res[:q_obs], linewidth = 1.2, color = "k", label = "Observed")
plt[:plot](df_res[:date], df_res[:q_sim], linewidth = 1, color = "r", label = "Simulated")
plt[:legend]()
plt[:title]("KGE = $perf_kge , NSE = $perf_nse")

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
using FFTW
using LinearAlgebra
using Statistics

export DistributionFunction

"""
    DistributionFunction( grid1, grid2 )

"""
struct DistributionFunction

    xgrid :: OneDGrid
    vgrid :: OneDGrid
    values :: Array{AbstractFloat, 2}

    function DistributionFunction( xgrid, vgrid )

        nx = xgrid.len
        nv = vgrid.len
        values = zeros(nx, nv)

        new( xgrid, vgrid, values )

    end

end

export landau!

"""
    landau!( f, α, kx )

Initialize the distribution function for the Landau damping test case

```math
f(x,v,t=0) = \\frac{1}{\\sqrt{2\\pi}} ( 1 + \\cos (kx x) ) \\exp (-\\frac{v^2}{2})
```

"""
function landau!( f :: DistributionFunction, α, kx)

    nx = f.xgrid.len
    nv = f.vgrid.len
    x  = f.xgrid.points
    v  = f.vgrid.points
    f.values .= (1 .+ α*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.^2))

end

export two_stream_instability!

function two_stream_instability!(f; eps = 0.01, xi = 0.90, v0 = 2.4)

    nx = f.xgrid.len
    nv = f.vgrid.len
    xg = f.xgrid.points
    vg = f.vgrid.points

    for (i,x) in enumerate(xg), (j,v) in enumerate(vg)
        f.values[i,j] = ((1 + eps*((cos(4π*x) + cos(6π*x)) / 1.2
                .+ cos(2π*x))) * (1/sqrt(2π)) * ((2-2xi)/(3-2xi))
                    * (1 + 5 * v^2 / (1-xi)) * exp(-.5 * v^2))
    end
end

max_stars_repo_name: JuliaVlasov/VlasovSolvers

max_issues_repo_name: JuliaVlasov/VlasovSolvers

max_forks_repo_name: JuliaVlasov/VlasovSolvers


text:
module test_voxelplots_magicavoxel_parser

using Test
using VoxelPlots.MagicaVoxel
using .MagicaVoxel: Voxel, Model, Size, Material, VoxData, ChunkError, ChunkStream, ChunkTree
using .MagicaVoxel: DEFAULT_PALETTE, DEFAULT_MATERIALS
using .MagicaVoxel: parse_chunk, parse_material, parse_vox_file, chunk_to_data, placeholder
using Colors: RGBA

function resource(block, filename)
    path = normpath(@__DIR__, "resources", filename)
    f = open(path)
    tree = ChunkTree([], UInt8[])
    stream = ChunkStream(f, tree)
    block(stream)
    close(stream)
end

# https://github.com/davidedmonds/dot_vox/blob/master/src/parser.rs#L209

resource("valid_size.bytes") do f
    chunk = Size(24, 24, 24)
    @test parse_chunk(f) == (:SIZE, chunk)

    seekstart(f)
    @test chunk_to_data(chunk) == read(f)
end

resource("valid_voxels.bytes") do f
    chunk = [Voxel(0, 0, 0, 225), Voxel(0, 1, 1, 215), Voxel(1, 0, 1, 235), Voxel(1, 1, 0, 5)]
    @test parse_chunk(f) == (:XYZI, chunk)

    seekstart(f)
    @test chunk_to_data(chunk) == read(f)
end

resource("valid_palette.bytes") do f
    (chunk_id, chunk) = parse_chunk(f)
    @test chunk isa Vector{RGBA}
    @test length(chunk) == 256
    @test chunk[1] == RGBA(1, 1, 1, 1)
    @test chunk[end-1] == RGBA(17/255, 17/255, 17/255, 1)
    @test chunk[end] == RGBA(0, 0, 0, 0)

    seekstart(f)
    @test chunk_to_data(chunk) == read(f)
end

resource("valid_material.bytes") do f
    chunk = Material(0, (_type = "_diffuse", _weight = "1", _rough = "0.1", _spec = "0.5", _ior = "0.3"))
    @test parse_material(f) == chunk

    seekstart(f)
    @test chunk_to_data(chunk)[12+1:end] == read(f)
end

resource("default_palette.bytes") do f
    chunk = MagicaVoxel.build_chunk(Val{:RGBA}(), f, 0, 0)
    @test first(chunk) == first(DEFAULT_PALETTE)
    @test MagicaVoxel.toRGBA(0xff99ccff) == RGBA(1, 0.8, 0.6, 1)
end

# https://github.com/davidedmonds/dot_vox/blob/master/src/lib.rs#L162

@test_throws ChunkError("Not a valid MagicaVoxel .vox file") MagicaVoxel.load(normpath(@__DIR__, "resources", "not_a.vox"))
@test_throws ChunkError("Unable to read file")               MagicaVoxel.load(normpath(@__DIR__, "resources", "not_here.vox"))

resource("placeholder.vox") do f
    chunk = parse_vox_file(f)
    @test chunk isa VoxData
    @test chunk.version == 150
    @test chunk.models[1].voxels == [Voxel(0, 0, 0, 225), Voxel(0, 1, 1, 215), Voxel(1, 0, 1, 235), Voxel(1, 1,0, 5)]
    @test first(chunk.palette) == first(DEFAULT_PALETTE)
    @test chunk.materials[1] == Material(0, (_type = "_diffuse", _weight = "1", _rough = "0.1", _spec = "0.5", _ior = "0.3"))
    @test chunk == placeholder(collect(DEFAULT_PALETTE), collect(DEFAULT_MATERIALS))

    seekstart(f)
    data = chunk_to_data(chunk)
    modeldata = chunk_to_data(chunk.models)
    @test data[1:16] == read(f, 20)[1:16] # VOX  version MAIN content_size
    @test modeldata == data[21:20+56] == read(f, 56)
end

resource("placeholder-with-materials.vox") do f
    chunk = parse_vox_file(f)
    materials = collect(DEFAULT_MATERIALS)
    materials[216+1] = Material(216, (_type = "_metal", _weight = "0.694737", _plastic = "1", _rough = "0.389474", _spec = "0.821053", _ior = "0.3"))
    @test chunk == placeholder(collect(DEFAULT_PALETTE), materials)
end

end # module test_voxelplots_magicavoxel_parser


text:
export DemoAgent, action

using ReinforcementLearning
import Functors: functor
using Setfield: @set

#This is an extended version of JuliaReinforcementLearning's Agent type. This one can take two trajectories in order to handle expert demonstrations.

Base.@kwdef struct DemoAgent{P<:AbstractPolicy, T<:AbstractTrajectory} <: AbstractPolicy
    policy::P
    trajectory::T
    demo_trajectory::T
end

functor(x::DemoAgent) = (policy = x.policy,), y -> @set x.policy = y.policy

(agent::DemoAgent)(env) = agent.policy(env)

function check(agent::DemoAgent, env::AbstractEnv)
    if ActionStyle(env) === FULL_ACTION_SET &&
        !haskey(agent.trajectory, :legal_actions_mask)
        !haskey(agent.demo_trajectory, :legal_actions_mask)
    end
    check(agent.policy, env)
end

function action(agent::DemoAgent, s; training = true)
    action(agent.policy, s, training=training)
end

function RLBase.update!(agent::DemoAgent)
    if !isnothing(agent.demo_trajectory) && length(agent.demo_trajectory) > 0
        update!(agent.policy, agent.trajectory, agent.demo_trajectory)
    else
        update!(agent.policy, agent.trajectory)
    end
end

Base.nameof(agent::DemoAgent) = nameof(agent.policy)

#####
# Default behaviors
#####

"""
Here we extend the definition of `(p::AbstractPolicy)(::AbstractEnv)` in
`RLBase` to accept an `AbstractStage` as the first argument. Algorithm designers
may customize these behaviors respectively by implementing:
- `(p::YourPolicy)(::AbstractStage, ::AbstractEnv)`
- `(p::YourPolicy)(::PreActStage, ::AbstractEnv, action)`
The default behaviors for `Agent` are:
1. Update the inner `trajectory` given the context of `policy`, `env`, and
   `stage`.
  1. By default we do nothing.
  2. In `PreActStage`, we `push!` the current **state** and the **action** into
     the `trajectory`.
  3. In `PostActStage`, we query the `reward` and `is_terminated` info from
     `env` and push them into `trajectory`.
  4. In the `PostEpisodeStage`, we push the `state` at the end of an episode and
     a dummy action into the `trajectory`.
  5. In the `PreEpisodeStage`, we pop out the lastest `state` and `action` pair
     (which are dummy ones) from `trajectory`.
2. Update the inner `policy` given the context of `trajectory`, `env`, and
   `stage`.
  1. By default, we only `update!` the `policy` in the `PreActStage`. And it's
     despatched to `update!(policy, trajectory)`.
"""

function (agent::DemoAgent)(stage::AbstractStage, env::AbstractEnv)
    update!(agent.trajectory, agent.policy, env, stage)
    update!(agent.policy, agent.trajectory, agent.demo_trajectory, env, stage)
end

function (agent::DemoAgent)(stage::PreActStage, env::AbstractEnv, action)
    update!(agent.trajectory, agent.policy, env, stage, action)
    update!(agent.policy, agent.trajectory, agent.demo_trajectory)
end

function (agent::DemoAgent)(stage::AbstractStage, stage2::InitStage, env::AbstractEnv)
    update!(agent.trajectory, agent.policy, env, stage)
    update!(agent.demo_trajectory, agent.policy, env, stage)
end

function (agent::DemoAgent)(stage::AbstractStage, stage2::PreTrainStage, env::AbstractEnv)
    #update!(agent.policy, agent.demo_trajectory, env, stage)
end

function (agent::DemoAgent)(stage::PreActStage, stage2::InitStage, env::AbstractEnv, action)
    update!(agent.trajectory, agent.policy, env, stage, action)
    update!(agent.demo_trajectory, agent.policy, env, stage, action)

end

function (agent::DemoAgent)(stage::PreTrainStage, env::AbstractEnv)
    update!(agent.policy, agent.demo_trajectory, env, stage)
end

function RLBase.update!(
    ::AbstractPolicy,
    ::AbstractTrajectory,
    ::AbstractTrajectory,
    ::AbstractEnv,
    ::AbstractStage,
    ::AbstractStage,
) end

function RLBase.update!(
    ::AbstractPolicy,
    ::AbstractTrajectory,
    ::AbstractTrajectory,
    ::AbstractEnv,
    ::AbstractStage,
) end

#####
# Default behaviors for known trajectories
#####

function RLBase.update!(
    ::AbstractTrajectory,
    ::AbstractPolicy,
    ::AbstractEnv,
    ::AbstractStage,
    ::AbstractStage,
) end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    demo_trajectory::AbstractTrajectory,
    ::AbstractPolicy,
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    if length(trajectory) > 0
        pop!(trajectory[:state])
        pop!(trajectory[:action])
        if haskey(trajectory, :legal_actions_mask)
            pop!(trajectory[:legal_actions_mask])
        end
    end

    if length(demo_trajectory) > 0
        pop!(demo_trajectory[:state])
        pop!(demo_trajectory[:action])
        if haskey(demo_trajectory, :legal_actions_mask)
            pop!(demo_trajectory[:legal_actions_mask])
        end
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    demo_trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PreActStage,
    action,
)
    s = policy isa NamedPolicy ? state(env, nameof(policy)) : state(env)
    push!(trajectory[:state], s)
    push!(trajectory[:action], action)
    if haskey(trajectory, :legal_actions_mask)
        lasm = policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) : legal_action_space_mask(env)
        push!(trajectory[:legal_actions_mask], lasm)
    end

    push!(demo_trajectory[:state], s)
    push!(demo_trajectory[:action], action)
    if haskey(demo_trajectory, :legal_actions_mask)
        lasm = policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) : legal_action_space_mask(env)
        push!(demo_trajectory[:legal_actions_mask], lasm)
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    demo_trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PostEpisodeStage,
)
    # Note that for trajectories like `CircularArraySARTTrajectory`, data are
    # stored in a SARSA format, which means we still need to generate a dummy
    # action at the end of an episode. Here we simply select a random one using
    # the global rng. In theory it shouldn't affect the performance of specific
    # algorithm.
    # TODO: how to inject a local rng here to avoid polluting the global rng
    action = rand(action_space(env))

    s = policy isa NamedPolicy ? state(env, nameof(policy)) : state(env)
    push!(trajectory[:state], s)
    push!(trajectory[:action], action)
    if haskey(trajectory, :legal_actions_mask)
        lasm = policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) : legal_action_space_mask(env)
        push!(trajectory[:legal_actions_mask], lasm)
    end

    push!(demo_trajectory[:state], s)
    push!(demo_trajectory[:action], action)
    if haskey(demo_trajectory, :legal_actions_mask)
        lasm = policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) : legal_action_space_mask(env)
        push!(demo_trajectory[:legal_actions_mask], lasm)
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    demo_trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PostActStage,
)
    r = policy isa NamedPolicy ? reward(env, nameof(policy)) : reward(env)
    push!(trajectory[:reward], r)
    push!(trajectory[:terminal], is_terminated(env))

    push!(demo_trajectory[:reward], r)
    push!(demo_trajectory[:terminal], is_terminated(env))
end

max_stars_repo_name: polhager/FurutaRL

max_issues_repo_name: polhager/FurutaRL

max_forks_repo_name: polhager/FurutaRL


text:
load("utils/req.jl")
req("dag/transforms.jl")

# -- scattering fusion --------------------------------------------------------

type Scatterer; end
typealias ScatterContext Context{Scatterer}

scattered(node::Node) = evaluate(ScatterContext(), node)
function evaluate(c::ScatterContext, node::SymNode)
    # todo: eliminate duplicate indvars?
    ellipsis = EllipsisNode(SymNode(:indvars, :local))
    node.val.name == :(... ) ? ellipsis : RefNode(node, ellipsis)
end
function evaluate(c::ScatterContext, node::Union(CallNode,RefNode))
    # todo: use first line once ref result type bug is fixed
#    Node(node, { node.args[1]::SymNode, evaluate(c, node.args[2:end])... })
    Node(node, { node.args[1]::SymNode,
                 evaluate(c, Node[node.args[2:end]...])... })
end

# == Front midsection =========================================================

# -- scatter propagation ------------------------------------------------------

type ScatterPropagator; end
typealias ScatterPropContext Context{ScatterPropagator}

scatter_propagated(sink::Node) = evaluate(ScatterPropContext(), sink)
function evaluate(c::ScatterPropContext, node::CallNode)
    if get_op(node).val == SymbolEx(:scatter, :call)
        args = get_callargs(node)
        @expect length(args)==1
        return (@cached scattered(c, args[1]))
    else
        return default_evaluate(c, node)
    end
end

function scattered(c::ScatterPropContext, ns::Vector)#Nodes)
    { (@cached scattered(c, node))|node in ns }
end
function scattered(c::ScatterPropContext, node::CallNode)
    op = get_op(node).val
    if op == SymbolEx(:scatter, :call)
        # let scatter*scatter = scatter
        args = get_callargs[node]
        @expect length(args)==1
        return scattered(c, args[1])
    else
        op = scattered_op(op)
        args = scattered(c, get_callargs(node))
        return CallNode(op, args...)
    end
end
# todo: scattered for other node types: RefNode, ...more?
function scattered(c::ScatterPropContext, node::SymNode)
    @expect node.val.kind == :input "expected kind == :input, got $(node.val.kind)"
    RefNode(node, SymNode(:..., :symbol))
end

function scattered_op(ex::SymbolEx)
    @expect ex.kind == :call
    op = ex.name

    if contains([:*, :/, :\ ], op)
        error("scattered: cannot scatter op = ", op)
    end
    SymNode(op, :call)
end

# == Back midsection ==========================================================

function expand_ellipsis_indexing(sink::Node, indvars::Vector{Symbol})
    indvars = { SymNode(indvar, :local) for indvar in indvars }
    rewrite_dag(sink, (node, args)->expandell_rewrite(node, args, indvars))
end

function expandell_rewrite(oldnode::Node, args::Vector, indvars::Vector)
    Node(oldnode, args)
end
function expandell_rewrite(oldnode::RefNode, args::Vector, indvars::Vector)
    inds = get_inds(oldnode)
    if (length(inds)==1) && (inds[1].val==SymbolEx(:..., :symbol))
        newinds = {args[1], indvars...}
        return Node(oldnode, newinds)
    else
        return Node(oldnode, args)
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
#=
Copyright (c) 2016, 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.
=#

function recreate_parfor_pre(body, linfo)
    @dprintln(3,"DistPass recreate_parfor_pre ast", linfo, body)
    lives = computeLiveness(body, linfo)
    # arrays that their allocations added to prev
    pre_alloc_arrs = LHSVar[]
    # for each basic block:
    # for each parfor, find its write array, add allocation to prestatements
    # remove meta and line numbers
    for bb in collect(values(lives.basic_blocks))
        for i in length(bb.statements):-1:1
            stmt = bb.statements[i].tls.expr
            if isBareParfor(stmt)
                parfor = stmt.args[1]
                rws = CompilerTools.ReadWriteSet.from_exprs(parfor.body, ParallelAccelerator.ParallelIR.pir_rws_cb, linfo, linfo)
                write_arrays = collect(keys(rws.writeSet.arrays))
                # fix to include reduction variable generated by 2nd gemm pattern
                # TODO: generalize rws somehow, or use liveness instead?
                #if length(parfor.preParFor)>=1 && isCall(parfor.preParFor[1]) &&
                #    parfor.preParFor[1].args[1]==GlobalRef(ParallelAccelerator.API,:set_zeros)
                #    push!(write_arrays, parfor.preParFor[1].args[2])
                #end
                @dprintln(3,"DistPass recreate_parfor_pre parfor found: ", parfor, rws, "\nwrite arrs: ",write_arrays)
                # only parfors with one write supported
                # TODO: generalize?
                if length(write_arrays)==1
                    arr = toLHSVar(write_arrays[1])
                    @dprintln(3,"DistPass recreate_parfor_pre write array: ", arr)
                    for j in i-1:-1:1
                        prev_stmt = bb.statements[j].tls.expr
                        if isAllocationAssignment(prev_stmt) && prev_stmt.args[1]==arr
                            @dprintln(3,"DistPass recreate_parfor_pre allocation for parfor found: ", prev_stmt)
                            parfor.preParFor = [prev_stmt; parfor.preParFor]
                            push!(pre_alloc_arrs, arr)
                            break
                        elseif in(arr,bb.statements[j].use) || in(arr,bb.statements[j].def)
                            break
                        end
                    end
                end
                # add initial value assignment of indmin() to preparfor to allow fusion
                if i>1 && isNewMTassign(bb.statements[i-1].tls.expr)
                    @dprintln(3,"DistPass recreate_parfor_pre indmin MT found ")
                    new_mt = bb.statements[i-1].tls.expr
                    parfor.preParFor = [new_mt; parfor.preParFor]
                end
            end
        end
    end
    @dprintln(3,"DistPass recreate_parfor_pre pre_alloc_arrs: ", pre_alloc_arrs)
    out = Any[]
    for i in 1:length(body.args)
        node = body.args[i]
        if !isMeta(node) && !(isAllocationAssignment(node) && in(node.args[1], pre_alloc_arrs))
            push!(out, node)
        end
        if isBareParfor(node)
            parfor = node.args[1]
            # fix parfor generated for indmin() to allow fusion
            if any(isNewMTassign, parfor.preParFor)
                empty!(parfor.first_input.range)
                parfor.first_input.out_dim = 1
                parfor.first_input.dim = 1
                parfor.original_domain_nodes = [ParallelAccelerator.ParallelIR.DomainOperation(:mmap!,Any[])]
                deleteat!(out, length(out)-1)
            end
            fix_parfor_for_fusion(parfor, length(out), linfo)
        end
    end
    body.args = out
end

isNewMTassign(node::Expr) = node.head==:(=) && isNewMT(node.args[2])
isNewMTassign(node::ANY) = false
isNewMT(node::Expr) = node.head==:new && is(node.args[1].name,ParallelAccelerator.API.Lib.MT.name)
isNewMT(node::ANY) = false

function fix_parfor_for_fusion(parfor::PIRParForAst, new_top_level_number, linfo)
    @dprintln(3,"DistPass fix_parfor_for_fusion parfor ", parfor)
    parfor.top_level_number = [new_top_level_number]
    empty!(parfor.array_aliases)
    # remove array assignment from post statements generated by expanding gemm
    if length(parfor.postParFor)==2 && isa(parfor.postParFor[1],Expr) &&
         parfor.postParFor[1].head==:(=) && isa(parfor.postParFor[1].args[2],LHSVar)
        @assert length(parfor.preParFor)>=1 && isAllocationAssignment(parfor.preParFor[1]) "invalid parfor for expanded gemm matched"
        @dprintln(3,"DistPass fix_parfor_for_fusion remove assignment from expanded gemm ")
        lhs_var = parfor.postParFor[1].args[1]
        out_var = parfor.postParFor[1].args[2]
        replaceExprWithDict!(parfor, Dict{LHSVar,Any}(out_var=>lhs_var), linfo, ParallelIR.AstWalk)
        parfor.postParFor = Any[0]
    end
    # use rws to update first_input, which is used for finding correlation in fusion
    # reverse order to match access array
    parfor_indices = [ toLHSVar(parfor.loopNests[i].indexVariable) for i in length(parfor.loopNests):-1:1 ]
    rws = CompilerTools.ReadWriteSet.from_exprs(parfor.body, ParallelAccelerator.ParallelIR.pir_rws_cb, linfo, linfo)
    @dprintln(3,"DistPass fix_parfor_for_fusion parfor indices ", parfor_indices)
    @dprintln(3,"DistPass fix_parfor_for_fusion rws arrays ", union(rws.readSet.arrays, rws.writeSet.arrays))

    for (arr,inds) in union(rws.readSet.arrays, rws.writeSet.arrays)
        # TODO: is this sufficient condition for parfor/array correlation?
        indices = map(x->isa(x,Colon)?x:toLHSVar(x), inds[1])
        if indices==parfor_indices
            @dprintln(3,"DistPass fix_parfor_for_fusion updating first_input.array from ",
                 parfor.first_input.array, " to ", arr)
            parfor.first_input.array = arr
            break
        end
    end
end

"""
Interchange is needed when the parfor is sequential,
but has a parallel (1D) parfor inside.
"""
function parfor_interchange_needed(parfor::PIRParForAst, state)
    @dprintln(3, "parfor_interchange_needed parfor: ", parfor)
    if state.parfor_partitioning[parfor.unique_id]!=SEQ
        return false
    end
    for node in parfor.body
        if isBareParfor(node) && state.parfor_partitioning[node.args[1].unique_id]==ONE_D
            @dprintln(3, "parfor_interchange_needed interchange needed due to: ", node.args[1])
            return true
        end
    end
    return false
end

function dist_optimize(body::Expr, state::DistPassState)
    @assert body.head==:body "invalid body in dist_optimize"
    out_body = Any[]
    for i in 1:length(body.args)
        new_node = dist_optimize_node(body.args[i], i, state)
        if isa(new_node, Array)
            append!(out_body, new_node)
        else
            push!(out_body, new_node)
        end
    end
    body.args = out_body
    recreate_parfor_pre(body, state.LambdaVarInfo)
    @dprintln(3, "dist_optimize after optimizing but before fusion ", body)
    state.LambdaVarInfo, body = ParallelAccelerator.ParallelIR.fusion_pass("dist_opt", state.LambdaVarInfo, body)
    # fusion changes the output of gemm_pp array to the large dimension, we need to change it back
    adjust_pp_alloc(body, state)
    return body
end

function adjust_pp_alloc(body, state)
    for (arr,sizes) in state.adjust_alloc
        @dprintln(3,"DistPass optimize adjust_alloc for ", arr," ", sizes)
        num_dims = length(sizes)
        for node in body.args
            if isBareParfor(node)
                for stmt in node.args[1].preParFor
                    if isAllocationAssignment(stmt)
                        @dprintln(3,"DistPass optimize adjust_alloc allocation found ", stmt)
                        out = stmt.args[1]
                        alloc = stmt.args[2]
                        if out==arr
                            for i in 1:num_dims
                                alloc.args[7+(i-1)*2] = sizes[i]
                            end
                            @dprintln(3,"DistPass optimize adjust_alloc allocation adjusted ", stmt)
                        end
                    end
                end
            end
        end
    end
end

function dist_optimize_node(node::Expr, top_level_number, state)
    @dprintln(3,"DistPass optimize node ", top_level_number, " ", node)
    if isAssignmentNode(node)
        #@dprintln(3,"DistPass optimize assignment: ", node)
        lhs = toLHSVar(node.args[1])
        rhs = node.args[2]
        return dist_optimize_assignment(node, state, top_level_number, lhs, rhs)
    elseif node.head==:parfor
        parfor = node.args[1]
        if parfor_interchange_needed(parfor, state)
            return doParforInterchange(node, state)
        end
        new_body = dist_optimize(Expr(:body, parfor.body...), state)
        parfor.body = new_body.args
    end
    return node
end

function dist_optimize_node(node::ANY, top_level_number, state)
    return node
end

function dist_optimize_assignment(node::Expr, state::DistPassState, top_level_number, lhs::LHSVar, rhs::RHSVar)
    return node
end

function dist_optimize_assignment(node::Expr, state::DistPassState, top_level_number, lhs::LHSVar, rhs::Expr)
    if rhs.head==:call && isBaseFunc(rhs.args[1],:gemm_wrapper!)
        @dprintln(3,"DistPass optimize gemm found: ", node)
        out = toLHSVar(rhs.args[2])
        arr1 = toLHSVar(rhs.args[5])
        t1 = (rhs.args[3]=='T')
        arr2 = toLHSVar(rhs.args[6])
        t2 = (rhs.args[4]=='T')
        # weight multipied by samples (e.g. w*points)
        if isSEQ(arr1,state) && isONE_D(arr2,state) && !t1 && !t2
            @dprintln(3,"DistPass optimize weight times points pattern found")
            return expand_gemm_sp(lhs, out, arr1, arr2, top_level_number, state)
        # labels multipied by samples (e.g. labels*points')
        elseif isSEQ(out,state) && isONE_D(arr1,state) && isONE_D(arr2,state) && !t1 && t2
            @dprintln(3,"DistPass optimize labels times points transpose pattern found")
            # output array is reused for gemm_pp output but its size needs to change
            state.adjust_alloc[lhs] = [state.arrs_dist_info[out].dim_sizes[1], state.arrs_dist_info[out].dim_sizes[end]]
            return expand_gemm_pp(lhs, out, arr1, arr2, top_level_number, state)
        end
    end
    return node
end

function dist_optimize_assignment(node::Expr, state::DistPassState, top_level_number, lhs::ANY, rhs::ANY)
    return node
end

function expand_gemm_sp(lhs, out, arr1, arr2, top_level_number, state)
    size1 = state.arrs_dist_info[arr2].dim_sizes[end]
    size2 = state.arrs_dist_info[arr1].dim_sizes[1]
    size3 = state.arrs_dist_info[arr1].dim_sizes[end]
    # outer loop over samples, inner loop over functions
    parfor_index1 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_parfor_"*string(getDistNewID(state))*"_index1"), Int, ISASSIGNED,state.LambdaVarInfo))
    parfor_index2 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_parfor_"*string(getDistNewID(state))*"_index2"), Int, ISASSIGNED,state.LambdaVarInfo))
    loop_index = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_loop_"*string(getDistNewID(state))*"_index"), Int, ISASSIGNED,state.LambdaVarInfo))
    elem_typ = eltype(CompilerTools.LambdaHandling.getType(arr1, state.LambdaVarInfo))
    temp_var1 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp1_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))
    temp_var2 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp2_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))
    temp_var3 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp3_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))

    loopNests = PIRLoopNest[ PIRLoopNest(parfor_index1, 1, size1, 1), PIRLoopNest(parfor_index2, 1, size2, 1) ]
    parfor_id = getDistNewID(state)
    first_input_info = InputInfo(out)
    first_input_info.dim = 2
    #first_input_info.indexed_dims = ones(Int64, first_input_info.dim)
    first_input_info.indexed_dims = [true,true] # loop over last dimension
    first_input_info.out_dim = 2
    first_input_info.elementTemp = temp_var2
    out_body = Any[]
    pre_statements  = Any[]
    post_statements = Any[ Expr(:(=), lhs, out), 0 ]

    push!(out_body, Expr(:(=), temp_var3, 0))
    # inner loop k dimension
    push!(out_body, Expr(:loophead, loop_index, 1, size3))
    # tmp1 = w[j,k]
    push!(out_body, Expr(:(=), temp_var1, mk_call(GlobalRef(Base,:unsafe_arrayref),[arr1, parfor_index2, loop_index])))
    # tmp2 = points[k,i]
    push!(out_body, Expr(:(=), temp_var2, mk_call(GlobalRef(Base,:unsafe_arrayref),[arr2, loop_index, parfor_index1])))
    push!(out_body, Expr(:(=), temp_var3, mk_add_float_expr(temp_var3, mk_mult_float_expr(temp_var1,temp_var2))))
    push!(out_body, Expr(:loopend, loop_index))
    push!(out_body, mk_call(GlobalRef(Base,:unsafe_arrayset),[out, temp_var3, parfor_index2, parfor_index1]))

    new_parfor = ParallelAccelerator.ParallelIR.PIRParForAst(
        first_input_info,
        out_body,
        pre_statements, Any[],
        loopNests,
        PIRReduction[],
        post_statements,
        [ParallelAccelerator.ParallelIR.DomainOperation(:mmap!,Any[])], # empty domain_oprs
        top_level_number,
        parfor_id,
        Set{LHSVar}(), #arrays_written_past_index
        Set{LHSVar}()) #arrays_read_past_index
    @dprintln(3,"DistPass optimize new_parfor ", new_parfor)
    state.parfor_partitioning[parfor_id] = ONE_D
    state.parfor_arrays[parfor_id] = [lhs,arr2]
    return Expr(:parfor, new_parfor)
end

function expand_gemm_pp(lhs, out, arr1, arr2, top_level_number, state)
    size1 = state.arrs_dist_info[arr2].dim_sizes[end]
    size2 = state.arrs_dist_info[arr1].dim_sizes[1]
    size3 = state.arrs_dist_info[out].dim_sizes[end]
    # outer loop over samples, inner loop over functions
    parfor_index1 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_parfor_"*string(getDistNewID(state))*"_index1"), Int, ISASSIGNED,state.LambdaVarInfo))
    parfor_index2 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_parfor_"*string(getDistNewID(state))*"_index2"), Int, ISASSIGNED,state.LambdaVarInfo))
    loop_index = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_loop_"*string(getDistNewID(state))*"_index"), Int, ISASSIGNED,state.LambdaVarInfo))
    elem_typ = eltype(CompilerTools.LambdaHandling.getType(arr1, state.LambdaVarInfo))
    temp_var1 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp1_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))
    temp_var2 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp2_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))
    temp_var3 = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
        Symbol("_dist_gemm_tmp3_"*string(getDistNewID(state))), elem_typ, ISASSIGNED | ISPRIVATEPARFORLOOP, state.LambdaVarInfo))

    loopNests = PIRLoopNest[ PIRLoopNest(parfor_index1, 1, size1, 1), PIRLoopNest(parfor_index2, 1, size2, 1) ]

    parfor_id = getDistNewID(state)
    first_input_info = InputInfo(arr2)
    first_input_info.dim = 2
    #first_input_info.indexed_dims = ones(Int64, first_input_info.dim)
    first_input_info.indexed_dims = [true,false] # loop over last dimension
    first_input_info.out_dim = 1
    first_input_info.elementTemp = temp_var2
    out_body = Any[]
    pre_statements  = Any[ mk_call(GlobalRef(ParallelAccelerator.API,:set_zeros),[out, mk_mult_int_expr([size3,size2])]) ]
    post_statements = Any[ Expr(:(=), lhs, out), 0 ]

    # tmp2 = labels[j,i]
    push!(out_body, Expr(:(=), temp_var2, mk_call(GlobalRef(Base,:unsafe_arrayref),[arr1, parfor_index2, parfor_index1])))
    # inner loop k dimension
    push!(out_body, Expr(:loophead, loop_index, 1, size3))
    # tmp1 = points[k,i]
    push!(out_body, Expr(:(=), temp_var1, mk_call(GlobalRef(Base,:unsafe_arrayref),[arr2, loop_index, parfor_index1])))
    # tmp3 = w[j,w]
    push!(out_body, Expr(:(=), temp_var3, mk_call(GlobalRef(Base,:unsafe_arrayref),[out, parfor_index2, loop_index])))
    push!(out_body, Expr(:(=), temp_var3, mk_add_float_expr(temp_var3, mk_mult_float_expr(temp_var1,temp_var2))))
    push!(out_body, mk_call(GlobalRef(Base,:unsafe_arrayset),[out, temp_var3, parfor_index2, loop_index]))
    push!(out_body, Expr(:loopend, loop_index))

    new_parfor = ParallelAccelerator.ParallelIR.PIRParForAst(
        first_input_info,
        out_body,
        pre_statements, Any[],
        loopNests,
        [PIRReduction(out, 0, GlobalRef(Base,:(+)))],
        post_statements,
        [ParallelAccelerator.ParallelIR.DomainOperation(:mmap!,Any[])], # empty domain_oprs
        top_level_number,
        parfor_id,
        Set{LHSVar}(), #arrays_written_past_index
        Set{LHSVar}()) #arrays_read_past_index
    @dprintln(3,"DistPass optimize new_parfor ", new_parfor)
    state.parfor_partitioning[parfor_id] = ONE_D
    state.parfor_arrays[parfor_id] = [lhs,arr2]
    return Expr(:parfor, new_parfor)
end

function genLoopHeadFromParfor(parfor)
    ret = Any[]

    for i = 1:length(parfor.loopNests)
        push!(ret, Expr(:loophead, parfor.loopNests[i].indexVariable, parfor.loopNests[i].lower, parfor.loopNests[i].upper))
    end

    return ret
end

function genLoopEndFromParfor(parfor)
    ret = Any[]

    for i = length(parfor.loopNests):-1:1
        push!(ret, Expr(:loopend, parfor.loopNests[i].indexVariable))
    end

    return ret
end

function getParforIndices(parfor)
    ret = Any[]

    for i = 1:length(parfor.loopNests)
        push!(ret, parfor.loopNests[i].indexVariable)
    end

    return ret
end

function genAddInt(x, val)
    return Expr(:call, GlobalRef(Base, :box), Int64, Expr(:call, GlobalRef(Base, :add_int), deepcopy(x), deepcopy(val)))
end
function genSubInt(x, val)
    return Expr(:call, GlobalRef(Base, :box), Int64, Expr(:call, GlobalRef(Base, :sub_int), deepcopy(x), deepcopy(val)))
end

function getParforSizes(parfor)
    ret = Union{RHSVar,Int,Expr}[]

    for i = 1:length(parfor.loopNests)
        if parfor.loopNests[i].step != 1
            throw(string("Skip not yet supported in getParforSizes."))
        end

        if parfor.loopNests[i].lower == 1
            push!(ret, deepcopy(parfor.loopNests[i].upper))
        else
            push!(ret, genAddInt(genSubInt(parfor.loopNests[i].upper, parfor.loopNests[i].lower), 1))
        end
    end

    return ret
end

type InterchangeState
    index_vars
    to_array
end

function interchangeArrayify(node::LHSVar, state::InterchangeState, top_level_number, is_top_level, read)
    if haskey(state.to_array, node)
        if read
            return Expr(:call, GlobalRef(Base, :arrayref), state.to_array[node], state.index_vars...)
        else
            throw(string("Don't handle case of write to arrayified symbol outside lhs of assignment."))
        end
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function interchangeArrayify(node::ANY, state::InterchangeState, top_level_number, is_top_level, read)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function doParforInterchange(parfor_node::Expr, state::DistPassState)
    @dprintln(3, "doParforInterchange ", parfor_node, " ", state)
    ret = Any[]
    new_array_allocs = Any[]
    for_pre = Any[]
    num_inner = 0

    outer_parfor = parfor_node.args[1]

    body_lives = CompilerTools.LivenessAnalysis.from_lambda(state.LambdaVarInfo, outer_parfor.body, ParallelIR.pir_live_cb, state.LambdaVarInfo)
    @dprintln(3, "body_lives = ", body_lives)
    non_nested_region = true

    index_vars = reverse(getParforIndices(outer_parfor))
    sizes = reverse(getParforSizes(outer_parfor))
    num_dims = length(index_vars)
    @dprintln(3, "index_vars = ", index_vars, " num_dims = ", num_dims, " sizes = ", sizes)

    to_array = Dict{LHSVar,LHSVar}()

    append!(for_pre, genLoopHeadFromParfor(outer_parfor))
    new_outer_parfor = nothing

    for i = 1:length(outer_parfor.body)
        node = outer_parfor.body[i]
        @dprintln(3, "Processing ", node, " ", non_nested_region)

        if ParallelAccelerator.ParallelIR.isParforAssignmentNode(node)
            throw(string("Parfor assignment nodes not yet supported in doParforInterchange."))
        elseif ParallelAccelerator.ParallelIR.isBareParfor(node)
            num_inner += 1
            if num_inner > 1
                @dprintln(1, "Multiple nested inner parfors not supported so reverting to original parfor.")
                return parfor_node
            end

            @dprintln(3, "isBareParfor ", non_nested_region)
            if non_nested_region
                append!(for_pre, genLoopEndFromParfor(outer_parfor))
                non_nested_region = false
            end
            inner_parfor = node.args[1]
            @dprintln(3, "isBareParfor ", non_nested_region, " ", inner_parfor)
            new_outer_parfor = deepcopy(inner_parfor)
            new_inner_parfor = deepcopy(outer_parfor)
            new_inner_parfor.body = deepcopy(inner_parfor.body)
            new_outer_parfor.body = Any[Expr(:parfor, new_inner_parfor)]
            new_outer_parfor.top_level_number = deepcopy(outer_parfor.top_level_number)
            for j = 1:length(new_outer_parfor.reductions)
                @dprintln(3, "Processing reduction ",  new_outer_parfor.reductions[j])
                rdsvar = toLHSVar(new_outer_parfor.reductions[j].reductionVar)
                if haskey(to_array, rdsvar)
                    new_outer_parfor.reductions[j].reductionVar = to_array[rdsvar]
                    @dprintln(3, "Changed to ",  new_outer_parfor.reductions[j])
                else
                    throw(string("During parfor interchange, reduction variable ", rdsvar, " on inner parfor did not become an array during interchange."))
                end
            end
            for j = 1:length(new_inner_parfor.body)
                inner_node = new_inner_parfor.body[j]
                @dprintln(3, "Processing inner parfor body node ", inner_node)

                if isAssignmentNode(inner_node)
                    lhs = toLHSVar(inner_node.args[1])
                    rhs = inner_node.args[2]
                    @dprintln(3, "Assignment node: ", lhs)
                    if haskey(to_array, lhs)
                        new_rhs = deepcopy(rhs)
                        new_inner_parfor.body[j] = Expr(:call, GlobalRef(Base, :arrayset), to_array[lhs], ParallelIR.AstWalk(new_rhs, interchangeArrayify, InterchangeState(index_vars, to_array)), index_vars...)
                    else
                        acopy = deepcopy(inner_node)
                        new_inner_parfor.body[j] = ParallelIR.AstWalk(acopy, interchangeArrayify, InterchangeState(index_vars, to_array))
                    end
                else
                    acopy = deepcopy(inner_node)
                    new_inner_parfor.body[j] = ParallelIR.AstWalk(acopy, interchangeArrayify, InterchangeState(index_vars, to_array))
                end
            end
            push!(ret, Expr(:parfor, new_outer_parfor))
        else
            @dprintln(3, "Node is not a parfor ", non_nested_region)
            if !non_nested_region
                append!(ret, genLoopHeadFromParfor(outer_parfor))
                non_nested_region = true
            end

            if isAssignmentNode(node)
                lhs = toLHSVar(node.args[1])
                rhs = node.args[2]
                @dprintln(3, "Assignment node: ", lhs)
                if !haskey(to_array, lhs)
                    lhs_type = CompilerTools.LambdaHandling.getType(lhs, state.LambdaVarInfo)
                    atype = Array{lhs_type, num_dims}
                    new_array_name = Symbol(string("HPAT_",lhs,"_",outer_parfor.unique_id))
                    CompilerTools.LambdaHandling.addLocalVariable(new_array_name, atype, ISASSIGNED, state.LambdaVarInfo)
                    new_array_lhsvar = toLHSVar(new_array_name, state.LambdaVarInfo)
                    @dprintln(3, "New array needed: ", lhs, " ", lhs_type, " ", new_array_name)
                    push!(new_array_allocs, ParallelAccelerator.ParallelIR.mk_assignment_expr(new_array_lhsvar, ParallelAccelerator.ParallelIR.mk_alloc_array_expr(lhs_type, atype, sizes...), state.LambdaVarInfo))
                    @dprintln(3, "new_array_allocs = ", new_array_allocs)
                    to_array[lhs] = new_array_lhsvar
                    @dprintln(3, "to_array = ", to_array)
                    state.arrs_dist_info[new_array_lhsvar] = ArrDistInfo(num_dims)
                    state.arrs_dist_info[new_array_lhsvar].partitioning = SEQ
                    state.arrs_dist_info[new_array_lhsvar].dim_sizes = sizes
                end
                new_rhs = deepcopy(rhs)
                push!(num_inner == 0 ? for_pre : ret, Expr(:call, GlobalRef(Base, :arrayset), new_array_lhsvar, ParallelIR.AstWalk(new_rhs, interchangeArrayify, InterchangeState(index_vars, to_array)), index_vars...))
            else
                acopy = deepcopy(node)
                push!(num_inner == 0 ? for_pre : ret, ParallelIR.AstWalk(acopy, interchangeArrayify, InterchangeState(index_vars, to_array)))
            end
        end
    end

    if non_nested_region
        append!(ret, genLoopEndFromParfor(outer_parfor))
    end

    assert(isa(new_outer_parfor, PIRParForAst))
    append!(new_outer_parfor.hoisted, new_array_allocs)
    append!(new_outer_parfor.hoisted, for_pre)

    @dprintln(3, "Output new_array_allocs = ", new_array_allocs)
    @dprintln(3, "Output regular code = ", ret)
    @dprintln(3, "new outer parfor = ", new_outer_parfor)

    return ret
    #return [new_array_allocs..., ret...]
end


text:
# # Mauna Loa time series example
#
# In this notebook, we apply Gaussian process regression to the Mauna Loa CO₂
# dataset. This showcases a rich combination of kernels, and how to handle and
# optimize all their parameters.

# ## Setup
#
# We make use of the following packages:

using CSV, DataFrames  # data loading
using AbstractGPs  # exact GP regression
using ParameterHandling  # for nested and constrained parameters
using Optim  # optimization
using Zygote  # auto-diff gradient computation
using Plots  # visualisation

# Let's load and visualize the dataset.

# !!! tip
#     The `let` block [creates a new
#     scope](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#scope-of-variables),
#     so any utility variables we define in here won't leak outside. This is
#     particularly helpful to keep notebooks tidy! The return value of the
#     block is given by its last expression.

(xtrain, ytrain), (xtest, ytest) = let
    data = CSV.read(joinpath(@__DIR__, "CO2_data.csv"), Tables.matrix; header=0)
    year = data[:, 1]
    co2 = data[:, 2]

    ## We split the data into training and testing set:
    idx_train = year .< 2004
    idx_test = .!idx_train

    (year[idx_train], co2[idx_train]), (year[idx_test], co2[idx_test])  # block's return value
end
## The utility variables such as `idx_train` and `idx_test` are not available outside the `let` scope

function plotdata()
    plot(; xlabel="year", ylabel="CO₂ [ppm]", legend=:bottomright)
    scatter!(xtrain, ytrain; label="training data", ms=2, markerstrokewidth=0)
    return scatter!(xtest, ytest; label="test data", ms=2, markerstrokewidth=0)
end

plotdata()

# ## Prior
#
# We will model this dataset using a sum of several kernels which describe
#
# - smooth trend: squared exponential kernel with long lengthscale;
# - seasonal component: periodic covariance function with period of one year,
#   multiplied with a squared exponential kernel to allow decay away from exact
#   periodicity;
# - medium-term irregularities: rational quadratic kernel;
# - noise terms: squared exponential kernel with short lengthscale
#   and uncorrelated observation noise.
#
# For more details, see [Rasmussen & Williams (2005), chapter 5](http://www.gaussianprocess.org/gpml/chapters/RW5.pdf).

# We will use
# [ParameterHandling.jl](https://invenia.github.io/ParameterHandling.jl/) for
# handling the (hyper)parameters of our model. It provides functions such as
# `positive` with which we can put constraints on the hyperparameters, and
# allows us to represent all required parameters as a nested NamedTuple:

#! format: off
## initial values to match http://stor-i.github.io/GaussianProcesses.jl/latest/mauna_loa/
θ_init = (;
    se_long = (;
        σ = positive(exp(4.0)),
        ℓ = positive(exp(4.0)),
    ),
    seasonal = (;
        ## product kernels only need a single overall signal variance
        per = (;
            ℓ = positive(exp(0.0)),  # relative to period!
            p = fixed(1.0),  # 1 year, do not optimize over
        ),
        se = (;
            σ = positive(exp(1.0)),
            ℓ = positive(exp(4.0)),
        ),
    ),
    rq = (;
        σ = positive(exp(0.0)),
        ℓ = positive(exp(0.0)),
        α = positive(exp(-1.0)),
    ),
    se_short = (;
        σ = positive(exp(-2.0)),
        ℓ = positive(exp(-2.0)),
    ),
    noise_scale = positive(exp(-2.0)),
)
#! format: on
#md nothing #hide

# We define a couple of helper functions to simplify the kernel construction:

SE(θ) = θ.σ^2 * with_lengthscale(SqExponentialKernel(), θ.ℓ)
## PeriodicKernel is broken, see https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/issues/389
##Per(θ) = with_lengthscale(PeriodicKernel(; r=[θ.ℓ/2]), θ.p)  # NOTE- discrepancy with GaussianProcesses.jl
Per(θ) = with_lengthscale(SqExponentialKernel(), θ.ℓ) ∘ PeriodicTransform(1 / θ.p)
RQ(θ) = θ.σ^2 * with_lengthscale(RationalQuadraticKernel(; α=θ.α), θ.ℓ)
#md nothing #hide

# This allows us to write a function that, given the nested tuple of parameter values, constructs the GP prior:

function build_gp_prior(θ)
    k_smooth_trend = SE(θ.se_long)
    k_seasonality = Per(θ.seasonal.per) * SE(θ.seasonal.se)
    k_medium_term_irregularities = RQ(θ.rq)
    k_noise_terms = SE(θ.se_short) + θ.noise_scale^2 * WhiteKernel()
    kernel = k_smooth_trend + k_seasonality + k_medium_term_irregularities + k_noise_terms
    return GP(kernel)  # [`ZeroMean`](@ref) mean function by default
end
#md nothing #hide

# ## Posterior
#
# To construct the posterior, we need to first build a [`FiniteGP`](@ref),
# which represents the infinite-dimensional GP at a finite number of input
# features:

function build_finite_gp(θ)
    f = build_gp_prior(θ)
    return f(xtrain)
end
#md nothing #hide

# !!! info "`WhiteKernel` vs `FiniteGP` observation noise"
#     In this notebook, we already included observation noise through the
#     `WhiteKernel` as part of the GP prior covariance in `build_gp_prior`. We
#     therefore call `f(xtrain)` which implies zero (additional) observation
#     noise.
#
#     Alternatively, we could have omitted the `θ.noise_scale^2 *
#     WhiteKernel()` term and instead passed the noise variance as a second
#     argument to the GP call in `build_finite_gp`, `f(xtrain,
#     θ.noise_scale^2)`.
#
#     These two approaches have slightly different semantics: In the first one,
#     the `WhiteKernel` contributes non-zero variance to the `[i, j]` element
#     of the covariance matrix of the `FiniteGP` if `xtrain[i] == xtrain[j]`
#     (based on the values of the features). In the second one, the observation
#     noise variance passed to `FiniteGP` only contributes to the diagonal
#     elements of the covariance matrix, i.e. for `i == j`.
#
#     Moreover, the variance (uncertainty) of the posterior predictions
#     includes the variance from the `WhiteKernel`, but does not include the
#     variance of the observation noise passed to the `FiniteGP`. To include
#     the observation noise in posterior predictions from the second approach,
#     call `fpost_opt(xtest, noise_scale^2)`.
#
# !!! tip
#     For most use-cases and if in any doubt, we recommend that you pass in
#     observation noise to the `FiniteGP`, and omit the explicit `WhiteKernel`.
#     This is slightly faster (no need to check `xtrain[i] == xtrain[j]` for
#     all pairs `i`, `j`), and `WhiteKernel` will not give stable gradients if
#     you wish to compute the gradient of the log marginal likelihood
#     `logpdf(f(x), y)` w.r.t. `x`.

# We obtain the posterior, conditioned on the (finite) observations, by calling
# [`posterior`](@ref):

function build_posterior_gp(θ)
    fx = build_finite_gp(θ)
    return posterior(fx, ytrain)
end
#md nothing #hide

# Now we can put it all together to obtain a [`PosteriorGP`](@ref).
# The call to `ParameterHandling.value` is required to replace the constraints
# (such as `positive` in our case) with concrete numbers:

fpost_init = build_posterior_gp(ParameterHandling.value(θ_init))

# Let's visualize what the GP fitted to the data looks like, for the initial choice of kernel hyperparameters.
#
# We use the following function to plot a GP `f` on a specific range, using the
# AbstractGPs [plotting
# recipes](https://juliagaussianprocesses.github.io/AbstractGPs.jl/dev/concrete_features/#Plotting).
# By setting `ribbon_scale=2` we visualize the uncertainty band with ``\pm 2``
# (instead of the default ``\pm 1``) standard deviations.

plot_gp!(f; label) = plot!(f(1920:0.2:2030); ribbon_scale=2, linewidth=1, label)
#md nothing #hide

plotdata()
plot_gp!(fpost_init; label="posterior f(⋅)")

# A reasonable fit to the data, but poor extrapolation away from the observations!

# ## Hyperparameter Optimization
#
# To improve the fit, we want to maximize the (log) marginal likelihood with
# respect to the hyperparameters.
# [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/) expects to
# minimize a loss, so we define it as the negative log marginal likelihood:

function loss(θ)
    fx = build_finite_gp(θ)
    lml = logpdf(fx, ytrain)  # this computes the log marginal likelihood
    return -lml
end
#md nothing #hide

# !!! note "Work-in-progress"
#     In the future, we are planning to provide the `optimize_loss` utility
#     function as part of JuliaGaussianProcesses -- for now, we just define it
#     inline.
#
# The L-BFGS parameters were chosen because they seem to work well empirically.
# You could also try with the defaults.

default_optimizer = LBFGS(;
    alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
    linesearch=Optim.LineSearches.BackTracking(),
)

function optimize_loss(loss, θ_init; optimizer=default_optimizer, maxiter=1_000)
    options = Optim.Options(; iterations=maxiter, show_trace=true)

    θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init)
    loss_packed = loss ∘ unflatten

    ## https://julianlsolvers.github.io/Optim.jl/stable/#user/tipsandtricks/#avoid-repeating-computations
    function fg!(F, G, x)
        if F !== nothing && G !== nothing
            val, grad = Zygote.withgradient(loss_packed, x)
            G .= only(grad)
            return val
        elseif G !== nothing
            grad = Zygote.gradient(loss_packed, x)
            G .= only(grad)
            return nothing
        elseif F !== nothing
            return loss_packed(x)
        end
    end

    result = optimize(Optim.only_fg!(fg!), θ_flat_init, optimizer, options; inplace=false)

    return unflatten(result.minimizer), result
end
#md nothing #hide

# We now run the optimization:

θ_opt, opt_result = optimize_loss(loss, θ_init)
opt_result

# The final value of the log marginal likelihood is:

-opt_result.minimum

# !!! warning
#     To avoid bad local optima, we could (and should) have carried out several
#     random restarts with different initial values for the hyperparameters,
#     and then picked the result with the highest marginal likelihood. We omit
#     this for simplicity. For more details on how to fit GPs in practice,
#     check out [A Practical Guide to Gaussian
#     Processes](https://tinyurl.com/guide2gp).
#
# Let's construct the posterior GP with the optimized hyperparameters:

fpost_opt = build_posterior_gp(ParameterHandling.value(θ_opt))
#md nothing #hide

# This is the kernel with the point-estimated hyperparameters:

fpost_opt.prior.kernel

# Let's print the optimized values of the hyperparameters in a more helpful format:

# !!! note "Work-in-progress"
#     This is another utility function we would eventually like to move out of this notebook:

using Printf

show_params(nt::Union{Dict,NamedTuple}) = String(take!(show_params(IOBuffer(), nt)))
function show_params(io, nt::Union{Dict,NamedTuple}, indent::Int=0)
    for (s, v) in pairs(nt)
        if v isa Union{Dict,NamedTuple}
            println(io, " "^indent, s, ":")
            show_params(io, v, indent + 4)
        else
            println(io, " "^indent, s, " = ", @sprintf("%.3f", v))
        end
    end
    return io
end

print(show_params(ParameterHandling.value(θ_opt)))

# And, finally, we can visualize our optimized posterior GP:

plotdata()
plot_gp!(fpost_opt; label="optimized posterior f(⋅)")

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__); Pkg.instantiate()

include("synapse_setup.jl")
include("helpers.jl")

simplified = structural_simplify(network)

sim_exp = ModelingToolkit.build_torn_function(simplified; expression=true)
cleaned = clean_expr(sim_exp, simplified, rgf=true, pretty=true, name=Symbol("simple_synapse!"))

clipboard(format_text(string(cleaned), YASStyle()))


text:
using LightGraphs
using SparseArrays

import Base: show, eltype, copy
import LightGraphs: nv, has_edge, has_vertex, add_edge!, rem_edge!, rem_vertex!,
    rem_vertices!, add_vertex!, add_vertices!, outneighbors, inneighbors, vertices, edges,
    adjacency_matrix, src, dst, nv, edgetype

export AbstractMultigraph
# export multype

"""
    AbstractMultigraph{T}<:AbstractGraph{T}

An abstract type representing a multigraph.
"""
abstract type AbstractMultigraph{T<:Integer} <:AbstractGraph{T} end

function copy(mg::AbstractMultigraph{T}) where {T} end

function show(io::IO, mg::AbstractMultigraph{T}) where {T}
    dir = is_directed(mg) ? "directed" : "undirected"
    print(io, "{$(nv(mg)), $(ne(mg))} $(dir) $(T) multigraph")
end

eltype(mg::AbstractMultigraph{T}) where {T<:Integer} = T
multype(mg::AbstractMultigraph) = Int
edgetype(mg::AbstractMultigraph) = MultipleEdge{eltype(mg), multype(mg)}

function nv(mg::AbstractMultigraph) end
function vertices(mg::AbstractMultigraph) end

function adjacency_matrix(mg::AbstractMultigraph) end

function has_edge(mg::AbstractMultigraph, e::AbstractMultipleEdge)
    s = src(e)
    d = dst(e)
    if has_vertex(mg, s) && has_vertex(mg, d)
        return mul(mg, s, d) >= mul(e)
    end
    return false
end

has_edge(mg::AbstractMultigraph, t) = has_edge(mg, MultipleEdge(t))
add_edge!(mg::AbstractMultigraph, t) = add_edge!(mg, MultipleEdge(t))
rem_edge!(mg::AbstractMultigraph, t) = rem_edge!(mg, MultipleEdge(t))

has_edge(mg::AbstractMultigraph, x, y) = has_edge(mg, MultipleEdge(x, y))
add_edge!(mg::AbstractMultigraph, x, y) = add_edge!(mg, MultipleEdge(x, y))
rem_edge!(mg::AbstractMultigraph, x, y) = rem_edge!(mg, MultipleEdge(x, y))

"""
    has_edge(mg::AbstractMultigraph, s, d, mul)

Return `true` if `mg` has a multiple edge from `s` to `d` whose multiplicity
is not less than `mul`.

## Examples
```julia
julia> using LightGraphs, Multigraphs

julia> mg = Multigraph(3);

julia> add_edge!(mg, 1, 2, 2);

julia> has_edge(mg, 1, 2, 3)
false

julia> has_edge(mg, 1, 2, 2)
true
```
"""
has_edge(mg::AbstractMultigraph, x, y, z) = has_edge(mg, MultipleEdge(x, y, z))

"""
    add_edge!(mg::AbstractMultigraph, s, d, mul)

Add a multiple edge from `s` to `d` multiplicity `mul`. If there is a multiple
edge from `s` to `d`, it will increase its multiplicity by `mul`.

Return `true` multiple edge was added successfully, otherwise return `false`.

## Examples
```julia
julia> using LightGraphs, Multigraphs

julia> mg = Multigraph(3);

julia> e = MultipleEdge(1, 2, 1);

julia> add_edge!(mg, e);

julia> ne(mg, true)
1

julia> add_edge!(mg, e);

julia> ne(mg, true)
2
```
"""
add_edge!(mg::AbstractMultigraph, x, y, z) = add_edge!(mg, MultipleEdge(x, y, z))

"""
    rem_edge!(mg::AbstractMultigraph, s, d, mul)

Remove the multiplicity of edge from `s` to `d` by `mul` in `mg`, if `mg` has such
a multiple edge.

## Examples
```julia
julia> using LightGraphs, Multigraphs

julia> mg = Multigraph(3);

julia> add_edge!(mg, 1, 2, 2);

julia> rem_edge!(mg, 1, 2, 3)
false

julia> rem_edge!(mg, 1, 2, 2)
true
```
"""
rem_edge!(mg::AbstractMultigraph, x, y, z) = rem_edge!(mg, MultipleEdge(x, y, z))

has_vertex(mg::AbstractMultigraph, v::Integer) = v in vertices(mg)
rem_vertex!(mg::AbstractMultigraph{T}, v::T) where {T<:Integer} = rem_vertices!(mg, [v])
add_vertex!(mg::AbstractMultigraph{T}) where {T<:Integer} = add_vertices!(mg, one(T))

function outneighbors(mg::AbstractMultigraph, v) end
function inneighbors(mg::AbstractMultigraph, v) end

"""
    edges(mg::AbstractMultigraph)

Return a  `MultipleEdgeIter` for `mg`.

## Examples
```jltestdoc
julia>
julia> using LightGraphs, Multigraphs

julia> mg = Multigraph(path_graph(4));

julia> add_edge!(mg, 1, 3, 2);

julia> collect(edges(mg))
4-element Array{Any,1}:
 Multiple edge 1 => 2 with multiplicity 1
 Multiple edge 1 => 3 with multiplicity 2
 Multiple edge 2 => 3 with multiplicity 1
 Multiple edge 3 => 4 with multiplicity 1

```
"""
edges(mg::AbstractMultigraph) = MultipleEdgeIter(mg)

"""
    mul(mg::AbstractMultigraph, src, dst)

Return the multiplicity of the edge from `src` to `dst`.
"""


text:
using Test
using SafeTestsets

@safetestset "LDU Tests" begin
    include("ldu_test.jl")
end


text:
# convex/concave relaxation (Khan 3.1-3.2) of integer powers of 1/x for positive reals
function cv_negpowpos(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  return x^n,n*x^(n-1)
end
function cc_negpowpos(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
end
# convex/concave relaxation of integer powers of 1/x for negative reals
@inline function cv_negpowneg(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  if (isodd(n))
    return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
  else
    return x^n,n*x^(n-1)
  end
end
@inline function cc_negpowneg(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  if (isodd(n))
    return x^n,n*x^(n-1)
  else
    return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
  end
end
# convex/concave relaxation of even powers greater than or equal to 4
@inline function cv_pow4(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  return x^n,n*x^(n-1)
end
@inline function cc_pow4(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
end
# convex/concave relaxation of odd powers
@inline function cv_powodd(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
          if (xU<=zero(T))
             return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
          elseif (zero(T)<=xL)
            return x^n,n*x^(n-1)
          else
            val::T = (xL^n)*(xU-x)/(xU-xL)+(max(zero(T),x))^n
            dval::T = -(xL^n)/(xU-xL)+n*(max(zero(T),x))^(n-1)
            return val,dval
          end
end
@inline function cc_powodd(x::T,xL::T,xU::T,n::Integer) where {T<:AbstractFloat}
  if (xU<=zero(T))
    return x^n,n*x^(n-1)
  elseif (zero(T)<=xL)
    return line_seg(x,xL,xL^n,xU,xU^n),dline_seg(x,xL,xL^n,xU,xU^n,n*x^(n-1))
  else
    val::T = (xU^n)*(x-xL)/(xU-xL)+(min(zero(T),x))^n
    dval::T = (xU^n)/(xU-xL)+n*(min(zero(T),x))^(n-1)
    return val,dval
  end
end
# power of a generalized McCormick object raised to c
function pow(x::SMCg{N,T},c::Integer) where {N,T<:AbstractFloat}
  mid1::T = zero(T)
  mid2::T = zero(T)
  cc::T = zero(T)
  cv::T = zero(T)
  if (c==0)
    return one(T)
  elseif (c>0)
    if (c==2)
      return sqr(x)
    elseif (c==1)
      return x
    elseif (isodd(c))
      eps_max::T = x.Intv.hi
      eps_min::T = x.Intv.lo
      midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
      midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
      if (MC_param.mu >= 1)
        cc,dcc::T = cc_powodd(midcc,x.Intv.lo,x.Intv.hi,c)
        cv,dcv::T = cv_powodd(midcv,x.Intv.lo,x.Intv.hi,c)
        gcc1::T,gdcc1::T = cc_powodd(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcv1::T,gdcv1::T = cv_powodd(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcc2::T,gdcc2::T = cc_powodd(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  gcv2::T,gdcv2::T = cv_powodd(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
      else
        cc,dcc = cc_powodd(midcc,x.Intv.lo,x.Intv.hi,c)
        cv,dcv = cv_powodd(midcv,x.Intv.lo,x.Intv.hi,c)
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
      end
      intv::Interval{T} = pow(x.Intv,c)
      return SMCg{N,T}(cc, cv, cc_grad, cv_grad, intv,x.cnst,x.IntvBox,x.xref)
    else
      if (x.Intv.hi<zero(x.cc))
        eps_min = x.Intv.hi
        eps_max = x.Intv.lo
      elseif (x.Intv.lo>zero(x.cc))
        eps_min = x.Intv.lo
        eps_max = x.Intv.hi
      else
        eps_min = zero(x.cc)
        eps_max = (abs(x.Intv.lo)>=abs(x.Intv.hi))? x.Intv.lo : x.Intv.hi
      end
      midcc,cc_id = mid3(x.cc,x.cv,eps_max)
      midcv,cv_id = mid3(x.cc,x.cv,eps_min)
      cc,dcc = cc_pow4(midcc,x.Intv.lo,x.Intv.hi,c)
      cv,dcv = cv_pow4(midcv,x.Intv.lo,x.Intv.hi,c)
      if (MC_param.mu >= 1)
        gcc1,gdcc1 = cc_pow4(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcv1,gdcv1 = cv_pow4(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcc2,gdcc2 = cc_pow4(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  gcv2,gdcv2 = cv_pow4(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  cv_grad = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    	  cc_grad = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
      else
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
      end
      intv = x.Intv^c
      return SMCg{N,T}(cc, cv, cc_grad, cv_grad, intv,x.cnst,x.IntvBox,x.xref)
    end
  else
    if (x.Intv.hi<zero(x.cc))
      if (isodd(c))
        eps_min = x.Intv.hi
        eps_max = x.Intv.lo
      else
        eps_min = x.Intv.lo
        eps_max = x.Intv.hi
      end
      midcc,cc_id = mid3(x.cc,x.cv,eps_max)
      midcv,cv_id = mid3(x.cc,x.cv,eps_min)
      cc,dcc = cc_negpowneg(midcc,x.Intv.lo,x.Intv.hi,c)
      cv,dcv = cv_negpowneg(midcv,x.Intv.lo,x.Intv.hi,c)
      if (MC_param.mu >= 1)
        gcc1,gdcc1 = cc_negpowneg(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcv1,gdcv1 = cv_negpowneg(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcc2,gdcc2 = cc_negpowneg(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  gcv2,gdcv2 = cv_negpowneg(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  cv_grad = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    	  cc_grad = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
      else
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
      end
      intv = x.Intv^c
      return SMCg{N,T}(cc, cv, cc_grad, cv_grad, intv,x.cnst,x.IntvBox,x.xref)
    elseif (x.Intv.lo>zero(x.cc))
      eps_min = x.Intv.hi
      eps_max = x.Intv.lo
      midcc,cc_id = mid3(x.cc,x.cv,eps_max)
      midcv,cv_id = mid3(x.cc,x.cv,eps_min)
      cc,dcc = cc_negpowpos(midcc,x.Intv.lo,x.Intv.hi,c)
      cv,dcv = cv_negpowpos(midcv,x.Intv.lo,x.Intv.hi,c)
      if (MC_param.mu >= 1)
        gcc1,gdcc1 = cc_negpowpos(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcv1,gdcv1 = cv_negpowpos(x.cv,x.Intv.lo,x.Intv.hi,c)
    	  gcc2,gdcc2 = cc_negpowpos(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  gcv2,gdcv2 = cv_negpowpos(x.cc,x.Intv.lo,x.Intv.hi,c)
    	  cv_grad = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    	  cc_grad = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
      else
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
      end
      intv = x.Intv^c
      return SMCg{N,T}(cc, cv, cc_grad, cv_grad, intv,x.cnst,x.IntvBox,x.xref)
    else
      error("Function unbounded on domain")
    end
  end
end
########### power of a generalized McCormick object raised to c
function  (^)(x::SMCg{N,T},c::Integer) where {N,T<:AbstractFloat}
  pow(x,c)
end
function  (^)(x::SMCg{N,T},c::Float64) where {N,T<:AbstractFloat}
  if (isinteger(c))
    return pow(x,Int64(c))
  else
    error("Relaxation of non-integer powers are unavailable.")
  end
end
########### Defines inverse
function inv(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  #println("ran inverse")
  x^(-1)
end


text:
using Float8s
using Test

@testset "Float8 -> Float32 ordered?" begin
    @testset for i in 0x00:UInt8(floatmax(Float8))
        f1 = Float32(Float8(i))
        f2 = Float32(Float8(i+0x1))
        @test f2 > f1
    end

    @testset for i in 0x80:UInt8(-floatmax(Float8))
        f1 = Float32(Float8(i))
        f2 = Float32(Float8(i+0x1))
        @test f2 < f1
    end
end

@testset "Float8_4 -> Float32 ordered?" begin
    @testset for i in 0x00:UInt8(floatmax(Float8_4))
        f1 = Float32(Float8_4(i))
        f2 = Float32(Float8_4(i+0x1))
        @test f2 > f1
    end

    @testset for i in 0x80:UInt8(-floatmax(Float8_4))
        f1 = Float32(Float8_4(i))
        f2 = Float32(Float8_4(i+0x1))
        @test f2 < f1
    end
end

@testset "Float32 -> Float8 non-subnormals ordered?" begin

    N = 100

    fs = (Float32(floatmax(Float8))-Float32(floatmin(Float8)))*rand(Float32,N) .+ Float32(floatmin(Float8))
    sort!(fs)

    f8s = Float8.(fs)
    @testset for i in 1:N-1
        @test f8s[i+1] >= f8s[i]
    end

    f8s = Float8.(-fs)
    @testset for i in 1:N-1
        @test f8s[i+1] <= f8s[i]
    end
end

@testset "Float32 -> Float8 subnormals ordered?" begin

    N = 100

    fs = Float32(floatmin(Float8))*rand(Float32,N)
    sort!(fs)

    f8s = Float8.(fs)
    @testset for i in 1:N-1
        @test f8s[i+1] >= f8s[i]
    end

    f8s = Float8.(-fs)
    @testset for i in 1:N-1
        @test f8s[i+1] <= f8s[i]
    end
end

@testset "Conversion Float8 <-> Float32" begin

    @testset for i in 0x00:0xff
        if ~isnan(Float8(i))
            @test i == reinterpret(UInt8,Float8(Float32(Float8(i))))
        end
    end
end

# Currently not implemented
# @testset "Conversion Float8_4 <-> Float32" begin
#
#     @testset for i in 0x00:0xff
#         if ~isnan(Float8_4(i))
#             @test i == reinterpret(UInt8,Float8_4(Float32(Float8_4(i))))
#         end
#     end
# end

@testset "Negation" begin

    @testset for i in 0x00:0xff
        f8 = Float8(i)
        f8_4 = Float8_4(i)

        if ~isnan(f8)
            @test f8 == -(-f8)
        end

        if ~isnan(f8_4)
            @test f8_4 == -(-f8_4)
        end
    end
end

@testset "Rounding" begin

    @testset for i in 0x00:0xff
        f8 = Float8(i)
        f8_4 = Float8_4(i)

        if ~isnan(f8)
            @test f8 >= floor(f8)
            @test f8 <= ceil(f8)
        end

        if ~isnan(f8_4)
            @test f8_4 >= floor(f8_4)
            @test f8_4 <= ceil(f8_4)
        end
    end
end

@testset "Triangle inequality Float8" begin

    @testset for i in 0x00:0xff
        for j in 0x00:0xff

            f1 = Float8(i)
            f2 = Float8(j)

            if ~isnan(f1) && ~isnan(f2) && isfinite(f1) && isfinite(f2)
                @test abs(f1) + abs(f2) >= abs(f1+f2)
                @test abs(f1) - abs(f2) <= abs(f1-f2)
                @test abs(f1) * abs(f2) >= f1*f2
            end
        end
    end
end

@testset "Triangle inequality Float8_4" begin

    @testset for i in 0x00:0xff
        for j in 0x00:0xff

            f1 = Float8_4(i)
            f2 = Float8_4(j)

            if ~isnan(f1) && ~isnan(f2) && isfinite(f1) && isfinite(f2)
                @test abs(f1) + abs(f2) >= abs(f1+f2)
                @test abs(f1) - abs(f2) <= abs(f1-f2)
                @test abs(f1) * abs(f2) >= f1*f2
            end
        end
    end
end

f = Float8(2.)
g = Float8(1.)

@testset "Comparison Float8" begin
    @test f >= g
    @test f > g
    @test g < f
    @test g <= g
    @test all([g g] .< [f f])
    @test all([g g] .<= [f f])
    @test all([f f] .> [g g])
    @test all([f f] .>= [g g])
    @test isless(g, f)
    @test !isless(f, g)

    @test Float8(2.5) == Float8(2.5)
    @test Float8(2.5) != Float8(2.6)
end

f = Float8_4(2.)
g = Float8_4(1.)

@testset "Comparison Float8_4" begin
    @test f >= g
    @test f > g
    @test g < f
    @test g <= g
    @test all([g g] .< [f f])
    @test all([g g] .<= [f f])
    @test all([f f] .> [g g])
    @test all([f f] .>= [g g])
    @test isless(g, f)
    @test !isless(f, g)

    @test Float8_4(2.5) == Float8_4(2.5)
    @test Float8_4(2.5) != Float8_4(2.7)
end

@testset "NaN8 and Inf8" begin
    @test isnan(NaN8)
    @test isnan(-NaN8)
    @test !isnan(Inf8)
    @test !isnan(-Inf8)
    @test !isnan(Float8(2.6))
    @test NaN8 != NaN8

    @test isinf(Inf8)
    @test isinf(-Inf8)
    @test !isinf(NaN8)
    @test !isinf(-NaN8)
    @test !isinf(Float8(2.6))
    @test Inf8 == Inf8
    @test Inf8 != -Inf8
    @test -Inf8 < Inf8
end

@testset "NaN8_4 and Inf8_4" begin
    @test isnan(NaN8_4)
    @test isnan(-NaN8_4)
    @test !isnan(Inf8_4)
    @test !isnan(-Inf8_4)
    @test !isnan(Float8(2.6))
    @test NaN8_4 != NaN8_4

    @test isinf(Inf8_4)
    @test isinf(-Inf8_4)
    @test !isinf(NaN8_4)
    @test !isinf(-NaN8_4)
    @test !isinf(Float8(2.6))
    @test Inf8_4 == Inf8_4
    @test Inf8_4 != -Inf8_4
    @test -Inf8_4 < Inf8_4
end

@testset "Nextfloat" begin
    for i in 0x00:0xff

        f = Float8(i)
        if isfinite(f)
            @test f < nextfloat(f)
            #@test 0x1 == UInt8(nextfloat(f))-i     # this is only true for positive numbers
        end

        f = Float8_4(i)
        if isfinite(f)
            @test f < nextfloat(f)
            #@test 0x1 == UInt8(nextfloat(f))-i     # this is only true for positive numbers
        end
    end

    @test NaN8 != nextfloat(NaN8)
    @test Inf8 == nextfloat(Inf8)
    @test -floatmax(Float8) == nextfloat(-Inf8)

    @test NaN8_4 != nextfloat(NaN8_4)
    @test Inf8_4 == nextfloat(Inf8_4)
    @test -floatmax(Float8_4) == nextfloat(-Inf8_4)
end

@testset "Prevfloat" begin
    for i in 0x00:0xff

        f = Float8(i)
        if isfinite(f)
            @test f > prevfloat(f)
            # @test 0x1 == i-UInt8(prevfloat(f))    # this is only true for positive numbers
        end

        f = Float8_4(i)
        if isfinite(f)
            @test f > prevfloat(f)
            # @test 0x1 == i-UInt8(prevfloat(f))    # this is only true for positive numbers
        end
    end

    @test NaN8 != prevfloat(NaN8)
    @test Inf8 > prevfloat(Inf8)
    @test NaN8_4 != prevfloat(NaN8_4)
    @test Inf8_4 > prevfloat(Inf8_4)

    @test -zero(Float8) > prevfloat(-zero(Float8))
    @test -zero(Float8_4) > prevfloat(-zero(Float8_4))

    @test zero(Float8) > prevfloat(zero(Float8))
    @test zero(Float8_4) > prevfloat(zero(Float8_4))

    @test prevfloat(zero(Float8)) == prevfloat(-zero(Float8))
    @test prevfloat(zero(Float8_4)) == prevfloat(-zero(Float8_4))

    @test -Inf8 == prevfloat(-Inf8)
    @test -Inf8_4 == prevfloat(-Inf8_4)
end

@testset "Signbit" begin

    for i in 0x00:0x7f      # positive numbers
        @test ~signbit(Float8(i))
        @test ~signbit(Float8_4(i))
    end

    for i in 0x80:0xff      # positive numbers
        @test signbit(Float8(i))
        @test signbit(Float8_4(i))
    end

end

@testset "Sign" begin

    for i in 0x01:0x7f

        f = Float8(i)
        if isnan(f)
            @test UInt8(f) == UInt8(sign(f))    # NaN == NaN yields false but the bitpattern are identical
        else
            @test one(f) == sign(f)
        end

        f = Float8_4(i)
        if isnan(f)
            @test UInt8(f) == UInt8(sign(f))    # NaN == NaN yields false but the bitpattern are identical
        else
            @test one(f) == sign(f)
        end
    end

    for i in 0x81:0xff

        f = Float8(i)
        if isnan(f)
            @test UInt8(f) == UInt8(sign(f))    # NaN == NaN yields false but the bitpattern are identical
        else
            @test -one(f) == sign(f)
        end

        f = Float8_4(i)
        if isnan(f)
            @test UInt8(f) == UInt8(sign(f))    # NaN == NaN yields false but the bitpattern are identical
        else
            @test -one(f) == sign(f)
        end
    end

    @test zero(Float8) == sign(zero(Float8))
    @test zero(Float8_4) == sign(zero(Float8_4))

    @test -zero(Float8) == sign(-zero(Float8))
    @test -zero(Float8_4) == sign(-zero(Float8_4))
end


text:
using CUDA
using Debugger
function test_cuda_debug(A,B)
    B .= A
end
A = CUDA.rand(Float32, 16)
B = CUDA.rand(Float32, 16)
#test_cuda_debug(A,B)
@run test_cuda_debug(A,B)

using SpikingNN
using Plots

# simulation parameters
T = 1000

# create three SRM0 neurons
η₀ = 30.0
τᵣ = 1.0
vth = 1.0

# create population
weights = Float32[ 0  5;
                   0  0]
pop = Population(weights; cell = () -> SRM0(η₀, τᵣ),
                          synapse = Synapse.Alpha,
                          threshold = () -> Threshold.Ideal(vth),
                          learner = STDP(0.5, 0.5), size(weights, 1)))
pop = gpu(pop)
# create step input currents
input = gpu(InputPopulation([ConstantRate(0.8)]))

# create network
net = Network(Dict([:input => input, :pop => pop]))
connect!(net, :input, :pop; weights = [1 0], synapse = Synapse.Alpha)

# simulate
w = Float64[]
voltages = Dict([(i, Float64[]) for i in 1:30])
cb = () -> begin
    push!(w, net[:pop].weights[1, 2])
    for id in 1:size(pop)
        push!(voltages[id], getvoltage(pop[id]))
    end
end
@time outputs = simulate!(net, T; cb = cb, dense = true)

weight_plot = plot(1:T, w, label = "")
title!("Synaptic Weights Over Simulation")
xlabel!("Time (sec)")
ylabel!("Weight")

raster_plot = rasterplot(outputs[:pop])
title!("Raster Plot")
xlabel!("Time (sec)")

plot(voltages[1], label = "Input")
voltage_plot = plot!(voltages[2], label = "Neuron")
title!("Membrane Voltages")
xlabel!("Time (sec)")
ylabel!("Voltage (V)")

plot(weight_plot, raster_plot, voltage_plot, layout = grid(3, 1))


text:
module SignalView

using GeometryTypes, Reactive, GLVisualize, GLAbstraction, Colors, GLWindow
import GraphLayout
import GraphLayout: AdjList

include("graphlayout.jl")
include("view.jl")

export view_signal

end