text:
function const_k2_2D(h)
    #Define Medium
    a = -0.4
    s0 = 2.0
    #Set up grid
    X = 0.0:h:8.0
    Z = 0.0:h:4.0
    sx = 4.0
    sz = 0.0
    sxi = findfirst(x->x==sx, X)
    szi = findfirst(z->z==sz, Z)
    x0 = CartesianIndex(szi, sxi)
    dist2 = ((X.-sx)').^2 .+ ((Z.-sz)).^2
    depth = 0.0*(X.-sx)'.+(Z.-sz)
    k2 = s0^2 .+ 2*a*depth
    Sb2 = s0^2 .+ a*depth
    sig2 = 2*dist2./(Sb2.+sqrt.(Sb2.^2 .- a*a*dist2))
    t_exact = Sb2.*sqrt.(sig2) .- a*a*sqrt.(sig2).^3 ./ 6
    (x0, k2, t_exact)
end

function const_v2_2D(h)
    #Define Medium
    a = 1.0
    s0 = 2.0
    #Set up grid
    X = 0.0:h:8.0
    Z = 0.0:h:4.0
    sx = 4.0
    sz = 0.0
    sxi = findfirst(x->x==sx, X)
    szi = findfirst(z->z==sz, Z)
    x0 = CartesianIndex(szi, sxi)
    dist2 = ((X.-sx)').^2 .+ ((Z.-sz)).^2
    depth = 0.0*(X.-sx)'.+(Z.-sz)
    k = 1 ./ (1/s0 .+ a*depth)
    t_exact = @. 1/a*acosh(1+1/2*s0*a*a*k*dist2)
    (x0, k.^2, t_exact)
end

function gaussian_factor_2D(h)
    sigmax = 0.4
    sigmaz = 0.1
    x1x_attempt = 2
    x1z_attempt = 4/3
    sx = 2.0
    sz = 1.0
    X = 0.0:h:8.0
    Z = 0.0:h:4.0
    sxi = findfirst(x->x==sx, X)#-1 # -1 here is for consistency with Treister & Haber paper / implementation
    szi = findfirst(z->z==sz, Z)#-1
    x1x = X[argmin(abs.(X.-x1x_attempt))]#-h # -h here is for consistency with Treister & Haber paper / implementation
    x1z = Z[argmin(abs.(Z.-x1z_attempt))]#-h

    x0 = CartesianIndex(szi, sxi)
    t1_exact = @. exp(-(sigmaz*((Z-x1z)^2)+sigmax*((X'-x1x)^2)))/2+1/2
    t1_dz_exact = @. -2*sigmaz*(Z-x1z)*(t1_exact-1/2)
    t1_dx_exact = @. -2*sigmax*(X'-x1x)*(t1_exact-1/2)

    t0 = ones(size(t1_exact))
    FEFMM.mul_analytic!(t0, h, h, size(t1_exact)..., szi, sxi)
    t0_dz, t0_dx = FEFMM.grad_analytic(t0, h, h, size(t1_exact)..., szi, sxi)

    t_dz = @. t0*t1_dz_exact+t0_dz*t1_exact
    t_dx = @. t0*t1_dx_exact+t0_dx*t1_exact
    k2 = @. t_dz*t_dz+t_dx*t_dx
    (x0, k2, t1_exact.*t0)
end

function const_k2_3D(h)
    #Define Medium
    a = -1.65
    s0 = 2.0
    #Set up grid
    X1 = 0.0:h:0.8
    X2 = 0.0:h:1.6
    X3 = 0.0:h:1.6
    sx1 = 0.0
    sx2 = 0.8
    sx3 = 0.8
    sx1i = findfirst(x->x==sx1, X1)
    sx2i = findfirst(x->x==sx2, X2)
    sx3i = findfirst(x->x==sx3, X3)
    x0 = CartesianIndex(sx1i, sx2i, sx3i)

    dist2 = zeros(length(X1), length(X2), length(X3))
    depth = zeros(length(X1), length(X2), length(X3))
    for k = 1:length(X3)
        for j = 1:length(X2)
            for i = 1:length(X1)
                dist2[i,j,k] = (X1[i]-sx1)^2+(X2[j]-sx2)^2+(X3[k]-sx3)^2
                depth[i,j,k] = (X1[i]-sx1)
            end
        end
    end

    k2 = @. s0^2+2*a*depth
    Sb2 = @. s0^2+a*depth
    sig2 = @. 2*dist2/(Sb2+sqrt(Sb2^2-a*a*dist2))
    t_exact = @. Sb2*sqrt(sig2)-a*a*(sqrt(sig2)^3)/6
    (x0, k2, t_exact)
end

function const_v2_3D(h)
    #Define Medium
    a = 1.0
    s0 = 2.0
    #Set up grid
    X1 = 0.0:h:0.8
    X2 = 0.0:h:1.6
    X3 = 0.0:h:1.6
    sx1 = 0.0
    sx2 = 0.8
    sx3 = 0.8
    sx1i = findfirst(x->x==sx1, X1)
    sx2i = findfirst(x->x==sx2, X2)
    sx3i = findfirst(x->x==sx3, X3)
    x0 = CartesianIndex(sx1i, sx2i, sx3i)

    dist2 = zeros(length(X1), length(X2), length(X3))
    depth = zeros(length(X1), length(X2), length(X3))
    for k = 1:length(X3)
        for j = 1:length(X2)
            for i = 1:length(X1)
                dist2[i,j,k] = (X1[i]-sx1)^2+(X2[j]-sx2)^2+(X3[k]-sx3)^2
                depth[i,j,k] = (X1[i]-sx1)
            end
        end
    end

    #get parameters
    k = @. 1/(1/s0+a*depth)
    t_exact = @. 1/a*acosh(1+1/2*s0*a*a*k*dist2)
    (x0, k.^2, t_exact)
end

function gaussian_factor_3D(h)
    #define medium
    sigmax1 = 0.2
    sigmax2 = 0.4
    sigmax3 = 0.1
    x11_attempt = 0.4
    x12_attempt = 1.6/3
    x13_attempt = 0.4
    sx1 = 0.2
    sx2 = 0.4
    sx3 = 0.4
    #set up grid
    X1 = 0.0:h:0.8
    X2 = 0.0:h:1.6
    X3 = 0.0:h:1.6
    sx1i = findfirst(x->x==sx1, X1)
    sx2i = findfirst(x->x==sx2, X2)
    sx3i = findfirst(x->x==sx3, X3)
    x0 = CartesianIndex(sx1i, sx2i, sx3i)

    x11 = X1[argmin(abs.(X2.-x11_attempt))]#-h # -h here is for consistency with Treister & Haber paper / implementation
    x12 = X2[argmin(abs.(X2.-x12_attempt))]#-h
    x13 = X3[argmin(abs.(X3.-x13_attempt))]

    t1_exact = zeros(length(X1), length(X2), length(X3))
    t1_dx1_exact = zeros(length(X1), length(X2), length(X3))
    t1_dx2_exact = zeros(length(X1), length(X2), length(X3))
    t1_dx3_exact = zeros(length(X1), length(X2), length(X3))

    # get medium
    for k = 1:length(X3)
        for j = 1:length(X2)
            for i = 1:length(X1)
                t1_exact[i,j,k] = exp(-(sigmax1*(X1[i]-x11)^2+sigmax2*(X2[j]-x12)^2+sigmax3*(X3[k]-x13)^2))/2+1/2
                t1_dx1_exact[i,j,k] = -2*sigmax1*(X1[i]-x11)*(t1_exact[i,j,k]-1/2)
                t1_dx2_exact[i,j,k] = -2*sigmax2*(X2[j]-x12)*(t1_exact[i,j,k]-1/2)
                t1_dx3_exact[i,j,k] = -2*sigmax3*(X3[k]-x13)*(t1_exact[i,j,k]-1/2)
            end
        end
    end

    t0 = ones(size(t1_exact))
    FEFMM.mul_analytic!(t0, h, h, h, size(t1_exact)..., sx1i, sx2i, sx3i)
    t0_dx1, t0_dx2, t0_dx3 = FEFMM.grad_analytic(t0, h, h, h, size(t1_exact)...,  sx1i, sx2i, sx3i)

    t_dx1 = @. t0*t1_dx1_exact+t0_dx1*t1_exact
    t_dx2 = @. t0*t1_dx2_exact+t0_dx2*t1_exact
    t_dx3 = @. t0*t1_dx3_exact+t0_dx3*t1_exact

    k2 = @. t_dx1*t_dx1+t_dx2*t_dx2+t_dx3*t_dx3
    (x0, k2, t1_exact.*t0)
end


text:
using Compat

fsm = state_machine(@compat Dict{AbstractString, Any}(
    "initial" => "first",
    "final" => "fourth",
    "events" => @compat( [ Dict{AbstractString, Any}(
            "name" => "hop",
            "from" => "first",
            "to" => "second"
        ),  Dict{AbstractString, Any}(
            "name" => "skip",
            "from" => "second",
            "to" => "third"
        ),  Dict{AbstractString, Any}(
            "name" => "jump",
            "from" => "third",
            "to" => "fourth"
        )
    ]),
    "callbacks" => @compat Dict{AbstractString, Any}(
        "onbeforeevent" => (fsm::StateMachine, args...) -> 1,
        "onafterevent" => (fsm::StateMachine, args...) -> 1 + 1,
        "onbeforejump" => (fsm::StateMachine, args...) -> 1 + 2,
        "onafterskip" => (fsm::StateMachine, args...) -> 2 + 3,
        "onchangestate" => (fsm::StateMachine, args...) -> 3 + 5,
        "onenterstate" => (fsm::StateMachine, args...) -> 5 + 8,
        "onentersecond" => (fsm::StateMachine, args...) -> 8 + 13,
        "onleavestate" => (fsm::StateMachine, args...) -> 13 + 21,
        "onleavefirst" => (fsm::StateMachine, args...) -> 21 + 34
    )
))

@test fsm.current == "first"
@test fsm.terminal == "fourth"
@test fsm.map == @compat( Dict{AbstractString, Any}(
    "startup" => @compat( Dict("none" => "first")),
    "skip" => @compat( Dict("second" => "third")),
    "jump" => @compat( Dict("third" => "fourth")),
    "hop" => @compat( Dict("first" => "second"))
))

dfsm = state_machine(@compat Dict{AbstractString, Any}(
    "initial" => @compat( Dict{AbstractString, Any}(
        "state" => "first",
        "event" => "dance",
        "defer" => true,
    )),
    "terminal" => "fourth",
    "events" => @compat [Dict{AbstractString, Any}(
            "name" => "hop",
            "from" => "first",
            "to" => "second",
        ), Dict{AbstractString, Any}(
            "name" => "skip",
            "from" => "second",
            "to" => "third",
        ), Dict{AbstractString, Any}(
            "name" => "sleep",
            "from" => "third",
        ), Dict{AbstractString, Any}(
            "name" => "coupdegrace",
            "from" => "third",
            "to" => "fourth",
        )
    ]
))

@test dfsm.current == "none"
@test dfsm.terminal == "fourth"

mfsm = state_machine(@compat Dict{AbstractString, Any}(
    "initial" => @compat( Dict{AbstractString, Any}(
        "state" => "first",
        "defer" => false,
    )),
    "events" => @compat( [Dict{AbstractString, Any}(
            "name" => "hop",
            "from" => ["first", "third"],
            "to" => "second",
        ), Dict{AbstractString, Any}(
            "name" => "skip",
            "from" => "second",
            "to" => "third",
        ), Dict{AbstractString, Any}(
            "name" => "sleep",
            "from" => "third",
        )
    ]),
    "callbacks" => @compat( Dict{AbstractString, Any}(
        "onevent" => (fsm::StateMachine, args...) -> 1,
        "onstate" => (fsm::StateMachine, args...) -> -1,
    ))
))

@test mfsm.current == "first"

ext: jl


text:
"""
   The adjoint operator of one_step_forward (in-place).
"""
function one_step_adjoint!(spt2::Snapshot, spt1::Snapshot, params::TdParams,
         tmp::Vector{Tv}, tmp_z1::Vector{Tv}, tmp_z2::Vector{Tv},
         tmp_x1::Vector{Tv}, tmp_x2::Vector{Tv}) where {Tv<:AbstractFloat}

    # update vz column-by-column
    for ix = 1 : params.Nx
        ilower = (ix-1) * params.Nz + 1
        iupper = ilower + params.Nz - 1

        idx = 0
        for iz = ilower : iupper
            idx = idx + 1
            tmp_z1[idx] = params.MpzBvz[iz] * spt1.pz[iz]
        end

        At_mul_b!(tmp_z2, params.dvdz, tmp_z1)
        # mul!(tmp_z2, transpose(params.dvdz), tmp_z1)

        idx = 0
        for iz = ilower : iupper
            idx = idx + 1
            spt2.vz[iz] = spt1.vz[iz] + tmp_z2[idx]
            spt2.pz[iz] = spt1.pz[iz] * params.MpzBpz[idx]
        end
    end

    # update vx row-by-row
    for iz = 1 : params.Nz

        idx = iz
        for ix = 1 : params.Nx
            tmp_x1[ix] = params.MpxBvx[idx] * spt1.px[idx]
            idx= idx + params.Nz
        end

        At_mul_b!(tmp_x2, params.dvdx, tmp_x1)
        # mul!(tmp_x2, transpose(params.dvdx), tmp_x1)

        idx = iz
        for ix = 1 : params.Nx
            spt2.vx[idx] = spt1.vx[idx] + tmp_x2[ix]
            spt2.px[idx] = spt1.px[idx] * params.MpxBpx[ix]
            idx= idx + params.Nz
        end
    end

    # update pz column-by-column
    for ix = 1 : params.Nx
        ilower = (ix-1) * params.Nz + 1
        iupper = ilower + params.Nz - 1

        idx = 0
        for iz = ilower : iupper
            idx = idx + 1
            tmp_z1[idx] = params.MvzBp[iz] * spt2.vz[iz]
        end

        At_mul_b!(tmp_z2, params.dpdz, tmp_z1)
        # mul!(tmp_z2, transpose(params.dpdz), tmp_z1)

        idx = 0
        for iz = ilower : iupper
            idx = idx + 1
            tmp[iz] = tmp_z2[idx]
        end
    end

    # update px row-by-row
    for iz = 1 : params.Nz

        idx = iz
        for ix = 1 : params.Nx
            tmp_x1[ix] = params.MvxBp[idx] * spt2.vx[idx]
            idx= idx + params.Nz
        end

        At_mul_b!(tmp_x2, params.dpdx, tmp_x1)
        # mul!(tmp_x2, transpose(params.dpdx), tmp_x1)

        idx = iz
        for ix = 1 : params.Nx
            tmp[idx] = tmp[idx] + tmp_x2[ix]
            idx= idx + params.Nz
        end
    end

    # finalize the result
    for ix = 1 : params.Nx
        row_idx = (ix-1) * params.Nz

        for iz = 1 : params.Nz
            idx = row_idx + iz

            spt2.vz[idx] = params.MvzBvz[iz] * spt2.vz[idx]
            spt2.vx[idx] = params.MvxBvx[ix] * spt2.vx[idx]
            spt2.pz[idx] = tmp[idx] + spt2.pz[idx]
            spt2.px[idx] = tmp[idx] + spt2.px[idx]
        end
    end

    return nothing
end

"""
   save the adjoint snapshot, wavefield or pressure field
"""
function multi_step_adjoint(rec::Recordings, params::TdParams;
                             path_spt="NULL", path_wfd="NULL" , path_pre="NULL", interval=1)

    # initialize intermediate variables
    spt1 = Snapshot(params)
    spt2 = Snapshot(params)
    tmp    = zeros(params.data_format, params.Nz * params.Nx)
    tmp_z1 = zeros(params.data_format, params.Nz)
    tmp_z2 = zeros(params.data_format, params.Nz)
    tmp_x1 = zeros(params.data_format, params.Nx)
    tmp_x2 = zeros(params.data_format, params.Nx)

    # inject the recordings to the last snapshot
    inject_rec2spt!(spt2, rec, params.nt)

    # the time range of source function
    if path_spt != "NULL"
       hdr_spt = snapshot_header(params, interval)
       fid_spt = write_RSheader(path_spt, hdr_spt)
       append_one_snapshot(fid_spt, spt2)
    end

    # save wavefield
    if path_wfd != "NULL"
       hdr_wfd = wavefield_header(params, interval)
       fid_wfd = write_RSheader(path_wfd, hdr_wfd)
       append_one_wavefield(fid_wfd, spt2, params)
    end

    # save pressure
    if path_pre != "NULL"
       hdr_pre = pressure_header(params, interval)
       fid_pre = write_RSheader(path_pre, hdr_pre)
       append_one_pressure(fid_pre, spt2, params)
    end

    # back propagation
    for it = params.nt-1 : -1 : 1

        one_step_adjoint!(spt1, spt2, params, tmp, tmp_z1, tmp_z2, tmp_x1, tmp_x2)
        inject_rec2spt!(spt1, rec, it)
        copy_snapshot!(spt2, spt1)

        # save wavefield at current step to disk
        if mod(it-1, interval) == 0
           path_spt != "NULL" && append_one_snapshot(fid_spt, spt2)
           path_wfd != "NULL" && append_one_wavefield(fid_wfd, spt2, params)
           path_pre != "NULL" && append_one_pressure(fid_pre, spt2, params)
        end
    end

    # close file and flush the buffer
    path_spt != "NULL" && close(fid_spt)
    path_wfd != "NULL" && close(fid_wfd)
    path_pre != "NULL" && close(fid_pre)

    # reverse the time order of the adjoint wavefield
    if path_spt != "NULL"
       tmp_dir = splitdir(path_spt)
       path_tmp= joinpath(tmp_dir[1], "spt_reverse_file.rsb")
       reverse_order(path_tmp, path_spt; save_flag="snapshot")
       mv(path_tmp, path_spt, force=true)
    end

    if path_wfd != "NULL"
       tmp_dir = splitdir(path_wfd)
       path_tmp= joinpath(tmp_dir[1], "wfd_reverse_file.rsb")
       reverse_order(path_tmp, path_wfd; save_flag="wavefield")
       mv(path_tmp, path_wfd, force=true)
    end

    if path_pre != "NULL"
       tmp_dir = splitdir(path_pre)
       path_tmp= joinpath(tmp_dir[1], "pre_reverse_file.rsb")
       reverse_order(path_tmp, path_pre; save_flag="pressure")
       mv(path_tmp, path_pre, force=true)
    end

    return nothing
end

"""
   one step backward reconstruction of wave field from the boundary values, the boundary
wavefield at only one time step is provided
"""
function one_step_backward!(wfd2::Wavefield, wfd1::Wavefield,
         bnd::WavefieldBound, params::TdParams,
         tmp_z1::Vector{Tv}, tmp_z2::Vector{Tv},
         tmp_x1::Vector{Tv}, tmp_x2::Vector{Tv}) where {Ti<:Int64, Tv<:AbstractFloat}

    # update pressure column-by-column
    for ix = 1 : params.nx
        ilower = (ix-1) * params.nz + 1
        iupper = ilower + params.nz - 1

        idx    = 0
        for iz = ilower : iupper
            idx= idx + 1
            tmp_z1[idx] = wfd1.vz[iz]
        end

        A_mul_b!(tmp_z2, params.rvdz, tmp_z1)        # dvdz

        idx    = 0
        for iz = ilower : iupper
            idx = idx + 1
            wfd2.p[iz] = wfd1.p[iz] - params.RpBv[iz]*tmp_z2[idx]
        end
    end

    # update pressure row-by-row
    for iz = 1 : params.nz

        idx = iz
        for ix = 1 : params.nx
            tmp_x1[ix] = wfd1.vx[idx]
            idx= idx + params.nz
        end

        A_mul_b!(tmp_x2, params.rvdx, tmp_x1)         # dvdx

        idx = iz
        for ix = 1 : params.nx
            wfd2.p[idx] = wfd2.p[idx] - params.RpBv[idx]*tmp_x2[ix]
            idx= idx + params.nz
        end
    end

    # correct for the boundary value
    for i = 1 : length(params.wfd2bnd)
        j = params.wfd2bnd[i]
        wfd2.p[j] = bnd.p[i]
    end

    # update vz column-by-column
    for ix = 1 : params.nx
        ilower = (ix-1) * params.nz + 1
        iupper = ilower + params.nz - 1

        idx    = 0
        for iz = ilower : iupper
            idx= idx + 1
            tmp_z1[idx] = wfd2.p[iz]
        end

        A_mul_b!(tmp_z2, params.rpdz, tmp_z1)        # dpdz

        idx    = 0
        for iz = ilower : iupper
            idx = idx + 1
            wfd2.vz[iz] = wfd1.vz[iz] - params.RvzBp[iz]*tmp_z2[idx]
        end
    end

    # update vx row-by-row
    for iz = 1 : params.nz

        idx = iz
        for ix = 1 : params.nx
            tmp_x1[ix] = wfd2.p[idx]
            idx= idx + params.nz
        end

        A_mul_b!(tmp_x2, params.rpdx, tmp_x1)         # dpdx

        idx = iz
        for ix = 1 : params.nx
            wfd2.vx[idx] = wfd1.vx[idx] - params.RvxBp[idx]*tmp_x2[ix]
            idx= idx + params.nz
        end
    end

    # correct for boundary part
    for i = 1 : length(params.wfd2bnd)
        j = params.wfd2bnd[i]
        wfd2.vz[j] = bnd.vz[i]
        wfd2.vx[j] = bnd.vx[i]
    end

    return nothing
end

"""
   Reconstruct pressure field backward using the boundary wavefield value and
the last wavefield
"""
function pressure_reconstruct_backward(path_bnd::Tp, path_lwfd::Tp,
         src::Ts, params::TdParams) where {Tp<:String, Ts<:Union{Source,Vector{Source}}}

    # make the source as a vector has one element
    if Ts <: Source
       srcs    = Vector{Source}(undef,1)
       srcs[1] = src
    else
       srcs    = src
    end
    ns = length(srcs)

    # length of one-step pressure field
    N = params.nz * params.nx

    # initialize intermediate variables
    wfd1 = Wavefield(params)
    wfd2 = read_one_wavefield(path_lwfd, 1)
    tmp_z1 = zeros(params.data_format, params.nz)
    tmp_z2 = zeros(params.data_format, params.nz)
    tmp_x1 = zeros(params.data_format, params.nx)
    tmp_x2 = zeros(params.data_format, params.nx)

    # initialize boundary value
    fid_bnd = open(path_bnd, "r")
    bnd     = WavefieldBound(params)

    # pre-allocate memory for saving the reconstructed pressure field
    pre = zeros(params.data_format, N * params.nt)
    idx_o = N*params.nt - N + 1 # the pressure field at last time step
    copyto!(pre, idx_o, wfd2.p, 1, N)

    # prepare for backward reconstruction
    subtract_source!(wfd2, srcs, params.nt)

    # backward iterations
    for it = params.nt-1 : -1 : 1

        # read the boundary value
        read_one_boundary!(bnd, fid_bnd, it, params)

        # one step backward reconstruction
        one_step_backward!(wfd1, wfd2, bnd, params,
                           tmp_z1, tmp_z2, tmp_x1, tmp_x2)

        # save the pressure field
        idx_o = idx_o - N
        copyto!(pre, idx_o, wfd1.p, 1, N)

        # prepare for next step
        copy_wavefield!(wfd2, wfd1)
        subtract_source!(wfd2, srcs, it)
    end

    close(fid_bnd)
    return reshape(pre, params.nz, params.nx, params.nt)
end


text:
"""Orbit defined by its minimum (L2-norm) and """
mutable struct Orbit{T<:AbstractFloat}
    N::Int              # number of variables
    length::Int         # length of orbit
    min::Vector{T}      # identify orbit by the point with minimum L2 norm
    norm::T             # that L2 norm
    basin::Int          # basin of attractions, number of random initial conditions that yield that orbit
    basin_norm::Real    # normalised to a fraction
end

# genereator function
Orbit(l::Int,min::Vector{T}) where T = Orbit{T}(length(min),l,min,norm(min),1,0.0)

# two orbits are equivalent when their βs, length and min agree
function Base.:(==)(o1::Orbit{T},o2::Orbit{T}) where T
    # if the norms don't match the orbits can't be identical
    # if the lengths match then the mins still might disagree
    condition1 = (o1.norm == o2.norm) && (o1.length == o2.length)

    if condition1
        o2min = copy(o2.min)
        not_equiv = o1.min != o2min
        nshifts = length(o2min)-1
        ishift = 0
        while not_equiv && ishift < nshifts
            o2min = circshift(o2min,1)
            not_equiv = o1.min != o2min ? true : false
            ishift += 1
        end
        return ~not_equiv
    else
        return false
    end
end

# two equivalent orbits added adds their basins
function Base.:(+)(o1::Orbit{T},o2::Orbit{T}) where T
    @assert o1 == o2 "$o1 != $o2."
    o1.basin += o2.basin
    return o1
end

# addition without uniqueness check
function unsafe_add(o1::Orbit{T},o2::Orbit{T}) where T
    o1.basin += o2.basin
    return o1
end

# return index of first element in x equal to y
function findin(x::Array{T,1},y::T) where T
    @inbounds for i in 1:length(x)
        x[i] == y && return i
    end
end

# orbit array reduction methods for distributed calculation
# check for uniqueness and extend orbit array if so
reduce_orbits(o1::Orbit{T},o2::Orbit{T}) where T = o1 == o2 ? [unsafe_add(o1,o2)] : [o1,o2]

function reduce_orbits(orbit::Orbit{T},orbits::Array{Orbit{T},1}) where T
    if orbit in orbits  # check whether orbit already exists in orbits
        # if so add their basins of attraction
        orbits[findin(orbits,orbit)] += orbit
        return orbits
    else  # append the orbit to the orbits array
        return push!(orbits,orbit)
    end
end

function reduce_orbits(orbits::Array{Orbit{T},1},orbit::Orbit{T}) where T
    if orbit in orbits  # check whether orbit already exists in orbits
        # if so add their basins of attraction
        orbits[findin(orbits,orbit)] += orbit
        return orbits
    else  # append the orbit to the orbits array
        return push!(orbits,orbit)
    end
end

function reduce_orbits(orbits1::Array{Orbit{T},1},orbits2::Array{Orbit{T},1}) where T
    for orbit in orbits1
        orbits2 = reduce_orbits(orbit,orbits2)
    end
    return orbits2
end

function normalise_basins!(orbits::Array{Orbit{T},1}) where T

    s = 0                       # sum of orbits
    for orbit in orbits
        if orbit.length > 0     # check that an orbit was actually found
            s += orbit.basin
        end
    end

    # normalise / set to zero for aborted orbit searches
    for orbit in orbits
        if orbit.length > 0
            orbit.basin_norm = orbit.basin/s
        else
            orbit.basin = 0
        end
    end
end

# define < for sorting
function Base.isless(o1::Orbit{T},o2::Orbit{T}) where T
    o1.length < o2.length
end

function Base.show(io::IO, o::Orbit{T}) where T

    N = o.N
    l = @sprintf("%10d",o.length)
    m = @sprintf("%s",repr(o.min))
    b = @sprintf("%s",repr(o.basin_norm))

    print(io,"Orbit{$T,N=$N}(length=$l, min=$m, basin=$b)")
end


text:
if VERSION < v"1.1"
    fieldtypes(t) = Tuple(fieldtype(t, i) for i = 1:fieldcount(t))
end

function coretype(M)
    if isdefined(M, :name)
        return M.name
    else
        return coretype(M.body)
    end
end


text:
using Plots, Test
using RecipesPipeline

@testset "plot" begin
    pl = plot(1:5)
    pl2 = plot(pl, tex_output_standalone = true)
    @test pl[:tex_output_standalone] == false
    @test pl2[:tex_output_standalone] == true
    plot!(pl, tex_output_standalone = true)
    @test pl[:tex_output_standalone] == true
end

@testset "get_axis_limits" begin
    x = [.1, 5]
    p1 = plot(x, [5, .1], yscale=:log10)
    p2 = plot!(identity)
    @test all(RecipesPipeline.get_axis_limits(p1, :x) .== x)
    @test all(RecipesPipeline.get_axis_limits(p2, :x) .== x)
end


text:
@doc raw"""
```
Z = retrieve_surface(algorithm::Pentland, img::AbstractArray)
```
Attempts to produce a heightmap from a grayscale image using Pentland's algorithm.

Under the assumptions of Lambertian surface, orthographic projections,
the surface being illuminated by distant light sources, the surface is
not self-shadowing and the surface has constant albedo, hence it can be ignored.
The algorithm employs Tayler series expansion and Fourier transforms to compute
a non-iterative solution.
# Output
Returns an M by N array (matching dimensions of original image) of Float `Z`
that represents the reconstructed height at the point.
# Details
In Pentlands algorithm the image irradiance is defined as:
```math
E(x,y)=R(p,q)=\dfrac{\rho(i_xp+i_yq-i_z)}{\sqrt{1+p^2+q^2}}=\dfrac{p\sin\sigma\cos
\tau+q\sin\sigma\sin\tau+\cos\sigma}{\sqrt{1+p^2+q^2}}
```
This can be reduced using the Taylor series expansion about ``p=p_0`` and ``p=p_0``
and ignoring the higher order terms becomes:
```math
E(x,y)=R(p,q)\approx R(p_0,q_0)+(p-p_0)\dfrac{\delta R}{\delta p}(p_0,q_0)+(q-q_0)
\dfrac{\delta R}{\delta q}(p_0,q_0)
```
which for ``p_0=q_0=0`` further reduces to:
```math
E(x,y)\approx\cos\sigma+p\cos\tau\sin\sigma+q\sin\tau\sin\sigma
```
This gives the following transform identities:
```math
\begin{gathered}
p=\dfrac{\delta}{\delta x}Z(x,y)\xleftrightarrow{\Im}(-j\omega_x)F_z(
\omega_x,\omega_y)\\q=\dfrac{\delta}{\delta y}Z(x,y)\xleftrightarrow{\Im}(-j
\omega_y)F_z(\omega_x,\omega_y)
\end{gathered}
```
By taking the Fourier transform of both sides if ``E(x,y)`` yields the following:
```math
F_E=(-j\omega_x)F_z(\omega_x,\omega_y)\cos\tau\sin\sigma+(-j\omega_y)F_z(\omega_x
,\omega_y)\sin\tau\sin\sigma
```
where ``F_z`` is the Fourier transform of ``Z(x,y)``.

These can be rearranged, and the Inverse Fourier transform used to recover the
surface ``Z(x,y)`` as per the following:
```math
\begin{gathered}
F_E=F_z(\omega_x,\omega_y)[-j\omega_x\cos\tau\sin\sigma-j\omega_y\sin\tau\sin
\sigma]\\\Rightarrow F_z(\omega_x,\omega_y)=\dfrac{F_E}{-j\omega_x\cos\tau\sin
\sigma-j\omega_y\sin\tau\sin\sigma}\\Z(x,y)=\Im^{-1}\{F_z(\omega_x,\omega_y)\}
\end{gathered}
```

The `slant` and `tilt` can be manually defined using the function signature:
```
Z = retrieve_surface(algorithm::Pentland, img::AbstractArray, slant::Real, tilt::Real)
```
Note: if `slant` and `tilt` are not defined they will be calculated at runtime
using `estimate_img_properties`.
# Arguments
The function arguments are described in more detail below.
##  `img`
An `AbstractArray` storing the grayscale value of each pixel within
the range [0,1].
## `slant`
A `Real` that specifies the slant value to be used by the algorithm. The `slant`
should be a value in the range [0,π/2]. If `slant` is specified to must the `tilt`.
## `tilt`
A `Real` that specifies the tilt value to be used by the algorithm. The `tilt`
should be a value in the range [0,2π]. If `tilt` is specified to must
the `slant`.
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(), 1, [0.2,0,0.9], radius = 5)

#calculate the heightmap
Z = retrieve_surface(Pentland(), img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. A. Pentland, "Shape Information From Shading: A Theory About Human Perception," [1988 Proceedings] Second International Conference on Computer Vision, Tampa, FL, USA, 1988, pp. 404-413. [doi: 10.1109/CCV.1988.590017](https://doi.org/10.1109/ccv.1988.590017)
"""
function retrieve_surface(algorithm::Pentland, img::AbstractArray)
    ρ, I, σ, τ = estimate_img_properties(img)
    return retrieve_surface(Pentland(),img,σ,τ)
end

function retrieve_surface(algorithm::Pentland, img::AbstractArray, slant::Real, tilt::Real)
    #find illumination and albedo
    σ = slant
    τ = tilt
    E = Complex{Float64}.(reinterpret(Float64, img))
    #take Fourier transform
    fft!(E)
    M, N = size(E)

    #setup wx and wy
    wx, wy = setup_transform_values_pentland(M, N)

    #using the illumination direction calculate the transformed Z
    Zₜ = zeros(Complex{Float64}, size(E))
    for i in CartesianIndices(Zₜ)
        Zₜ[i] = E[i] / (-1im * wx[i] * cos(τ) * sin(σ) - 1im * wy[i] * sin(τ)
            * sin(σ))
    end

    #recover Z
    ifft!(Zₜ)
    Z = zeros(Float64, size(E))
    for i in CartesianIndices(Z)
        Z[i] = abs(Zₜ[i])
    end
    return Z
end


text:
export compile_factor_graph,
    fg_to_cnf,
    wmc_chavira,
    get_varprob,
    zoo_fg_file,
    time_compilation,
    time_compilation_cnf

using Pkg.Artifacts

# Might want to move this to LoadSave later
zoo_fg_file(name) =
    artifact"circuit_model_zoo" * LogicCircuits.LoadSave.zoo_version * "/fgs/$name"

"""
    wmc_chavira(root::LogicCircuit; varprob::Function)::Float

Compute weighted model count in the context of a chavira encoding (default negative literals to 1)
Probability of each variable is given by `varprob` Function which defauls to 1/2 for every variable.
"""
function wmc_chavira(root::LogicCircuit,
                  varprob::Function = v -> 0.5)::Float32
    f_con(n) = istrue(n) ? 1.0 : 0.0
    f_lit(n) = ispositive(n) ? varprob(variable(n)) : 1.0
    f_a(n, call) = mapreduce(call, *, children(n))
    f_o(n, call) = mapreduce(call, +, children(n))
    foldup(root, f_con, f_lit, f_a, f_o, Float64)
end

function get_varprob(fg::FactorGraph, var_lits, fac_lits)
    d = Dict()
    # Just give ever variable indicator literal weight 1 for now
    for kv in collect(var_lits)
        for literal in kv[2]
              d[variable(literal)] = 1.0
        end
    end

    for fac in fg.facs
        for k in keys(fac.factor)
            lit = fac_lits[fac][k]
            weight = fac.factor[k]
            d[variable(lit)] = weight
        end
    end

    x -> d[x]
end

function bijection(a, b)
    (!a | b) & (a | !b)
end

function biject_cnf(a, b)
    # Here we assume that A is a list of literals and b is a literal
    clause1 = vcat(b, map(x -> negate(x), a))
    clauses = map(x -> [negate(b), x], a)
    vcat([clause1], clauses)
end

function fg_to_cnf(fg::FactorGraph)
    curr = 0
    # So we can get ids for each
    function nextSDDVar()
        curr += 1
        return curr
    end

    # Use a balanced vtree for now
    total_lits = sum(v.dim for v in fg.vars) + sum(length(keys(f.factor)) for f in fg.facs)
    vtr = PlainVtree(total_lits, :balanced)
    mgr = SddMgr(vtr)

    # Create lits for varnodes
    var_lits = Dict(v => [PlainLiteralNode(nextSDDVar()) for i = 1:v.dim] for v in fg.vars)

    # lits for factor nodes
    fac_lits = Dict(f => Dict(k => PlainLiteralNode(nextSDDVar()) for k in keys(f.factor)) for f in fg.facs)

    @assert total_lits == curr

    # Our cnf is just going to be a list of lists
    cnf = Vector{Vector{PlainLiteralNode}}()

    # Compile indicator CNFs
    for varnode in fg.vars
        for i = 2:varnode.dim
            for j = 1:i-1
                clause1 = [negate(var_lits[varnode][i]), negate(var_lits[varnode][j])]
                push!(cnf, clause1)
                clause2 = [var_lits[varnode][i], var_lits[varnode][j]]
                push!(cnf, clause2)
            end
        end
    end

    # Compile factor CNFs
    for facnode in fg.facs
        for config in keys(facnode.factor)
            fc = [var_lits[varnode][assignment] for (assignment, varnode) in zip(config, facnode.neighbs)]
            append!(cnf, biject_cnf(fc, fac_lits[facnode][config]))
        end
    end
    conjoin(map(x -> disjoin(x), cnf))
end

function compile_factor_graph(fg::FactorGraph, vtr=nothing)
    curr = 0
    # So we can get ids for each
    function nextSDDVar()
        curr += 1
        return curr
    end

    if vtr === nothing
        # Use a balanced vtree for now
        total_lits = sum(v.dim for v in fg.vars) + sum(length(keys(f.factor)) for f in fg.facs)
        vtr = PlainVtree(total_lits, :balanced)
    end
    mgr = SddMgr(vtr)

    # Create lits for varnodes
    var_lits = Dict(v => [compile(mgr, Lit(nextSDDVar())) for i = 1:v.dim] for v in fg.vars)

    # lits for factor nodes
    fac_lits = Dict(f => Dict(k => compile(mgr, Lit(nextSDDVar())) for k in keys(f.factor)) for f in fg.facs)

    # Begin compilation
    c = compile(mgr, true)

    # Compile indicator CNFs
    for varnode in fg.vars
        for i = 2:varnode.dim
            for j = 1:i-1
                c = c & (!var_lits[varnode][i] | !var_lits[varnode][j])
                c = c & (var_lits[varnode][i] | var_lits[varnode][j])
            end
        end
    end

    # Compile factor CNFs
    for facnode in fg.facs
        for config in keys(facnode.factor)
            fc = compile(mgr, true)
            for (assignment, varnode) in zip(config, facnode.neighbs)
                fc = fc & var_lits[varnode][assignment]
            end
            c = c & bijection(fc, fac_lits[facnode][config])
        end
    end
    c, var_lits, fac_lits
end

function time_compilation(fg, vtr)
    @time compile_factor_graph(fg, vtr)
end

function time_compilation_cnf(fg, vtr)
    mgr = SddMgr(vtr)
    cnf = fg_to_cnf(fg)
    @time compile(mgr, cnf)
end


text:
## matern.jl : implementation of Mat\'ern covariance function

## Matern ##
struct Matern{T} <: IsotropicCovarianceStructure{T}
    λ::T
    ν::T
    σ::T
    p::T
end

"""
    Matern(λ, ν, σ=1, p=2)

Create a Mat\u00E9rn covariance structure with correlation length `λ`, smoothness `ν`, (optional) marginal standard deviation `σ` and (optional) `p`-norm.

# Examples
```jldoctest
julia> m1 = Matern(0.1, 1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> m2 = Matern(0.1, 1.0, σ=2.0)
Matérn (λ=0.1, ν=1.0, σ=2.0, p=2.0)

```
"""
function Matern(λ::T where {T<:Real}, ν::T where {T<:Real}; σ=1.0::T where {T<:Real}, p=2::T where {T<:Real})
    λ > 0 || throw(ArgumentError("correlation length λ of Mat\u00E9rn covariance cannot be negative or zero"))
    ν > 0 || throw(ArgumentError("smoothness ν of Mat\u00E9rn covariance cannot be negative or zero"))
    σ > 0 || throw(ArgumentError("marginal standard deviation σ of Mat\u00E9rn covariance cannot be negative or zero"))
    p >= 1 || throw(ArgumentError("in p-norm, p must be greater than or equal to 1"))
	isinf(p) && throw(ArgumentError("in p-norm, p cannot be infinity"))
    Matern{promote_type(typeof(λ),typeof(ν),typeof(σ),typeof(p))}(promote(λ,ν,σ,p)...)
end

# TODO : Evaluation of the besselk function is extremely expensive; maybe try to catch special cases for ν = 1/2, 3/2 etc.
# evaluate Matern covariance
function apply(m::Matern,x::T) where {T<:Real}
    x == zero(T) ? one(T) : 2^(1-m.ν)/gamma(m.ν)*(sqrt(2*m.ν)*x/m.λ)^m.ν*besselk(m.ν,sqrt(2*m.ν)*x/m.λ)
end

show(io::IO, m::Matern) = print(io, "Mat\u00E9rn (λ=$(m.λ), ν=$(m.ν), σ=$(m.σ), p=$(m.p))")


text:
export JacobiWeight

"""
`JacobiWeight`
weights a basis on `[-1,1]` weighted by `(1+x)^α*(1-x)^β`.
Note the inconsistency of the parameters with `Jacobi`.
when the domain is `[a,b]` the weight is inferred by mapping to `[-1,1]`
"""
immutable JacobiWeight{S,DD} <: WeightSpace{S,RealBasis,DD,1}
    α::Float64
    β::Float64
    space::S
    function JacobiWeight(α::Float64,β::Float64,space::S)
        if isa(space,JacobiWeight)
            JacobiWeight(α+space.α,β+space.β,space.space)
        else
            new(α,β,space)
        end
    end
end

JacobiWeight(a::Number,b::Number,d::RealUnivariateSpace)=JacobiWeight{typeof(d),typeof(domain(d))}(Float64(a),Float64(b),d)
JacobiWeight(a::Number,b::Number,d::IntervalDomain)=JacobiWeight(Float64(a),Float64(b),Space(d))
JacobiWeight(a::Number,b::Number,d::Vector)=JacobiWeight(Float64(a),Float64(b),Space(d))
JacobiWeight(a::Number,b::Number)=JacobiWeight(a,b,Chebyshev())

JacobiWeight(a::Number,b::Number,s::PiecewiseSpace) = PiecewiseSpace(JacobiWeight(a,b,vec(s)))

identity_fun(S::JacobiWeight)=isapproxinteger(S.α)&&isapproxinteger(S.β)?Fun(x->x,S):Fun(identity,domain(S))

spacescompatible(A::JacobiWeight,B::JacobiWeight)=A.α==B.α && A.β == B.β && spacescompatible(A.space,B.space)
spacescompatible{DD<:IntervalDomain}(A::JacobiWeight,B::RealUnivariateSpace{DD})=spacescompatible(A,JacobiWeight(0,0,B))
spacescompatible{DD<:IntervalDomain}(B::RealUnivariateSpace{DD},A::JacobiWeight)=spacescompatible(A,JacobiWeight(0,0,B))

transformtimes{JW1<:JacobiWeight,JW2<:JacobiWeight}(f::Fun{JW1},g::Fun{JW2})=
            Fun(coefficients(transformtimes(Fun(f.coefficients,f.space.space),
                                            Fun(g.coefficients,g.space.space))),
                             JacobiWeight(f.space.α+g.space.α,f.space.β+g.space.β,f.space.space))
transformtimes{JW<:JacobiWeight}(f::Fun{JW},g::Fun) = Fun(coefficients(transformtimes(Fun(f.coefficients,f.space.space),g)),f.space)
transformtimes{JW<:JacobiWeight}(f::Fun,g::Fun{JW}) = Fun(coefficients(transformtimes(Fun(g.coefficients,g.space.space),f)),g.space)

##  α and β are opposite the convention for Jacobi polynomials
# Here, α is the left algebraic singularity and β is the right algebraic singularity.

jacobiweight(α,β,x)=(1+x).^α.*(1-x).^β
jacobiweight(α,β,d::Domain)=Fun([1.],JacobiWeight(α,β,ConstantSpace(d)))
jacobiweight(α,β)=jacobiweight(α,β,Interval())

weight(sp::JacobiWeight,x)=jacobiweight(sp.α,sp.β,tocanonical(sp,x))
dimension(sp::JacobiWeight)=dimension(sp.space)

Base.first{JW<:JacobiWeight}(f::Fun{JW}) = space(f).α>0?zero(eltype(f)):f(first(domain(f)))
Base.last{JW<:JacobiWeight}(f::Fun{JW}) = space(f).β>0?zero(eltype(f)):f(last(domain(f)))

setdomain(sp::JacobiWeight,d::Domain)=JacobiWeight(sp.α,sp.β,setdomain(sp.space,d))

# we assume that points avoids singularities

##TODO: paradigm for same space
function coefficients{SJ1,SJ2,DD<:IntervalDomain}(f::Vector,sp1::JacobiWeight{SJ1,DD},sp2::JacobiWeight{SJ2,DD})
    α,β=sp1.α,sp1.β
    c,d=sp2.α,sp2.β

    if isapprox(c,α) && isapprox(d,β)
        # remove wrapper spaces and then convert
        coefficients(f,sp1.space,sp2.space)
    else
        # go back to default
        defaultcoefficients(f,sp1,sp2)
    end
end
coefficients{SJ,S,IT,DD<:IntervalDomain}(f::Vector,sp::JacobiWeight{SJ,DD},
                                         S2::SubSpace{S,IT,RealBasis,DD,1}) = error("Implement")
coefficients{SJ,S,IT,DD<:IntervalDomain}(f::Vector,
                                         S2::SubSpace{S,IT,RealBasis,DD,1},
                                         sp::JacobiWeight{SJ,DD}) = error("Implement")
#TODO: it could be possible that we want to JacobiWeight a SumSpace....
coefficients{SJ,SV,DD<:IntervalDomain}(f::Vector,sp::JacobiWeight{SJ,DD},S2::SumSpace{SV,RealBasis,DD,1}) =
    sumspacecoefficients(f,sp,S2)
coefficients{SJ,DD<:IntervalDomain}(f::Vector,sp::JacobiWeight{SJ,DD},S2::RealUnivariateSpace{DD}) =
    coefficients(f,sp,JacobiWeight(0,0,S2))
coefficients{SJ,DD<:IntervalDomain}(f::Vector,S2::RealUnivariateSpace{DD},sp::JacobiWeight{SJ,DD}) =
    coefficients(f,JacobiWeight(0,0,S2),sp)

"""
`increase_jacobi_parameter(f)` multiplies by `1-x^2` on the unit interval.
`increase_jacobi_parameter(-1,f)` multiplies by `1+x` on the unit interval.
`increase_jacobi_parameter(+1,f)` multiplies by `1-x` on the unit interval.
On other domains this is accomplished by mapping to the unit interval.
"""
increase_jacobi_parameter(f) = Fun(f,JacobiWeight(f.space.α+1,f.space.β+1,space(f).space))
increase_jacobi_parameter(s,f) = s==-1?Fun(f,JacobiWeight(f.space.α+1,f.space.β,space(f).space)):
                                       Fun(f,JacobiWeight(f.space.α,f.space.β+1,space(f).space))

function canonicalspace(S::JacobiWeight)
    if isapprox(S.α,0) && isapprox(S.β,0)
        canonicalspace(S.space)
    else
        #TODO: promote singularities?
        JacobiWeight(S.α,S.β,canonicalspace(S.space))
    end
end

function union_rule{P<:PolynomialSpace}(A::ConstantSpace,B::JacobiWeight{P})
    # we can convert to a space that contains contants provided
    # that the parameters are integers
    # when the parameters are -1 we keep them
    if isapproxinteger(B.α) && isapproxinteger(B.β)
        JacobiWeight(min(B.α,0.),min(B.β,0.),B.space)
    else
        NoSpace()
    end
end

## Algebra

for op in (:/,:./)
    @eval begin
        function ($op){JW<:JacobiWeight}(c::Number,f::Fun{JW})
            g=($op)(c,Fun(f.coefficients,space(f).space))
            Fun(g.coefficients,JacobiWeight(-f.space.α,-f.space.β,space(g)))
        end
    end
end

function .^{JW<:JacobiWeight}(f::Fun{JW},k::Float64)
    S=space(f)
    g=Fun(coefficients(f),S.space)^k
    Fun(coefficients(g),JacobiWeight(k*S.α,k*S.β,space(g)))
end

function .*{JW1<:JacobiWeight,JW2<:JacobiWeight}(f::Fun{JW1},g::Fun{JW2})
    @assert domainscompatible(f,g)
    fα,fβ=f.space.α,f.space.β
    gα,gβ=g.space.α,g.space.β
    m=(Fun(f.coefficients,space(f).space).*Fun(g.coefficients,space(g).space))
    if isapprox(fα+gα,0)&&isapprox(fβ+gβ,0)
        m
    else
        Fun(m.coefficients,JacobiWeight(fα+gα,fβ+gβ,space(m)))
    end
end

./{JW1<:JacobiWeight,JW2<:JacobiWeight}(f::Fun{JW1},g::Fun{JW2})=f*(1/g)

# O(min(m,n)) Ultraspherical conjugated inner product

function conjugatedinnerproduct{S,V}(sp::Ultraspherical,u::Vector{S},v::Vector{V})
    λ=order(sp)
    if λ==1
        mn = min(length(u),length(v))
        if mn > 0
            return dotu(u[1:mn],v[1:mn])*π/2
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    else
        T,mn = promote_type(S,V),min(length(u),length(v))
        if mn > 1
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            ret = u[1]*wi*v[1]
            for i=2:mn
              wi *= (i-2one(T)+2λ)/(i-one(T)+λ)*(i-2one(T)+λ)/(i-one(T))
              ret += u[i]*wi*v[i]
            end
            return ret
        elseif mn > 0
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            return u[1]*wi*v[1]
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    end
end

function conjugatedinnerproduct(::Chebyshev,u::Vector,v::Vector)
    mn = min(length(u),length(v))
    if mn > 1
        return (2u[1]*v[1]+dotu(u[2:mn],v[2:mn]))*π/2
    elseif mn > 0
        return u[1]*v[1]*π
    else
        return zero(promote_type(eltype(u),eltype(v)))
    end
end

function bilinearform{LT,D}(f::Fun{JacobiWeight{Ultraspherical{LT,D},D}},g::Fun{Ultraspherical{LT,D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.α == f.space.β == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform{LT,D}(f::Fun{Ultraspherical{LT,D}},
                            g::Fun{JacobiWeight{Ultraspherical{LT,D},D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ && g.space.α == g.space.β == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform{LT,D}(f::Fun{JacobiWeight{Ultraspherical{LT,D},D}},
                            g::Fun{JacobiWeight{Ultraspherical{LT,D},D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ && f.space.α+g.space.α == f.space.β+g.space.β == λ-0.5
        return complexlength(domain(f))/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function linebilinearform{LT,D}(f::Fun{JacobiWeight{Ultraspherical{LT,D},D}},g::Fun{Ultraspherical{LT,D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.α == f.space.β == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform{LT,D}(f::Fun{Ultraspherical{LT,D}},g::Fun{JacobiWeight{Ultraspherical{LT,D},D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ &&  g.space.α == g.space.β == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform{LT,D}(f::Fun{JacobiWeight{Ultraspherical{LT,D},D}},g::Fun{JacobiWeight{Ultraspherical{LT,D},D}})
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ &&  f.space.α+g.space.α == f.space.β+g.space.β == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end


text:
export Sdd, SddMgr,
    SddLeafNode, SddInnerNode, SddLiteralNode, SddConstantNode,
    Sdd⋀Node, Sdd⋁Node, SddTrueNode, SddFalseNode,
    sdd_mgr_for

#############
# Trimmed Sdds
#############

"Root of the trimmed Sdd circuit node hierarchy"
abstract type Sdd <: StructLogicCircuit end

#############
# Elements and XY Partitions
#############

"Represents elements that are not yet compiled into conjunctions"
struct Element
    prime::Sdd
    sub::Sdd
end

"Represent an XY-partition that has not yet been compiled into a disjunction"
const XYPartition = Vector{Element}

"Unique nodes cache for decision nodes"
const Unique⋁Cache = Dict{XYPartition,Sdd}

"Broader definition of an XYPartition, also including vectors of nodes"

Base.hash(x::XYPartition) = mapreduce(e -> hash(e.prime, hash(e.sub)), ⊻, x)

function Base.isequal(x::XYPartition, y::XYPartition)
    length(x) != length(y) && return false
    (x == y) && return true
    # note: the rest of this function could be removed if the XYPartitions were sorted, but that is not advantageous: the sorting is too slow
    return all(x) do e1
        any(y) do e2
            e1.prime === e2.prime && e1.sub === e2.sub
        end
    end
end

#############
# Apply cache
#############

"Representation of the arguments of an Apply call"
struct ApplyArgs
    a1::Sdd
    a2::Sdd
end

Base.hash(aa::ApplyArgs) = hash(aa.a1) ⊻ hash(aa.a2)
Base.isequal(x::ApplyArgs,y::ApplyArgs) =
    (x.a1 === y.a1 && x.a2 === y.a2 ) || (x.a1 === y.a2 && x.a2 === y.a1)

"Apply cache for the result of conjunctions and disjunctions"
const ApplyCache = Dict{ApplyArgs,Sdd}

#############
# Trimmed Sdd managers
#############

"Root of the SDD manager node hierarchy"
abstract type SddMgr <: Vtree end

"SDD manager inner vtree node for trimmed SDD nodes"
mutable struct SddMgrInnerNode <: SddMgr

    left::SddMgr
    right::SddMgr

    parent::Union{SddMgrInnerNode, Nothing}

    variables::BitSet
    unique⋁cache::Unique⋁Cache

    conjoin_cache::ApplyCache

    SddMgrInnerNode(left::SddMgr, right::SddMgr) where D = begin
        @assert isdisjoint(variables(left), variables(right))
        this = new(left, right, nothing,
            union(variables(left), variables(right)),
            Unique⋁Cache(), ApplyCache())
        @assert left.parent === nothing
        left.parent = this
        @assert right.parent === nothing
        right.parent = this
        this
    end

end

"SDD manager leaf vtree node for trimmed SDD nodes"
mutable struct SddMgrLeafNode <: SddMgr

    var::Var
    parent::Union{SddMgrInnerNode, Nothing}

    positive_literal::Sdd
    negative_literal::Sdd

    SddMgrLeafNode(v::Var) = begin
        this = new(v, nothing)
        this.positive_literal = construct_new_sdd_literal(this, var2lit(v))
        this.negative_literal = construct_new_sdd_literal(this, -var2lit(v))
        this
    end

end

@inline NodeType(::SddMgrLeafNode) = Leaf()
@inline NodeType(::SddMgrInnerNode) = Inner()

#############
# Sdd nodes
#############

"A SDD logical leaf node"
abstract type SddLeafNode <: Sdd end

"A SDD logical inner node"
abstract type SddInnerNode <: Sdd end

"A SDD logical literal leaf node, representing the positive or negative literal of its variable"
mutable struct SddLiteralNode <: SddLeafNode
    literal::Lit
    vtree::SddMgrLeafNode
    data
    counter::UInt32
    SddLiteralNode(l,v) = new(l,v,nothing,false)
end

"""
A SDD logical constant leaf node, representing true or false.
These are the only structured nodes that don't have an associated vtree node (cf. SDD file format)
"""
abstract type SddConstantNode <: SddLeafNode end

"A SDD logical true constant."
mutable struct SddTrueNode <: SddConstantNode
    counter::UInt32
    data
end

# there is an issue with using this canonical node: if someone breaks the `bit` field in one SDD manager, it will be broken for all SDD managers...
"Canonical true Sdd node"
const true_sdd = SddTrueNode(false, nothing)

"A SDD logical false constant."
mutable struct SddFalseNode <: SddConstantNode
    counter::UInt32
    data
end

# there is an issue with using this canonical node: if someone breaks the `bit` field in one SDD manager, it will be broken for all SDD managers...
"Canonical false Sdd node"
const false_sdd = SddFalseNode(false, nothing)

"A SDD logical conjunction node"
mutable struct Sdd⋀Node <: SddInnerNode
    prime::Sdd
    sub::Sdd
    vtree::SddMgrInnerNode
    counter::UInt32
    data
    Sdd⋀Node(p,s,v) = new(p,s,v,false)
end

"A SDD logical disjunction node"
mutable struct Sdd⋁Node <: SddInnerNode
    children::Vector{Sdd⋀Node}
    vtree::SddMgrInnerNode
    counter::UInt32
    negation::Sdd⋁Node
    data
    Sdd⋁Node(ch,v) = new(ch, v, false) # leave negation uninitialized
    Sdd⋁Node(ch,v,neg) = new(ch, v, false, neg)
end

#####################
# traits
#####################

@inline GateType(::Type{<:SddLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:SddConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Sdd⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Sdd⋁Node}) = ⋁Gate()

#####################
# Constructor
#####################

SddMgr(v::Var) = SddMgrLeafNode(v)
SddMgr(left::SddMgr, right::SddMgr) = SddMgrInnerNode(left, right)

"Obtain an SDD manager that can support compiling the given circuit"
sdd_mgr_for(c::Sdd) = mgr(c)
sdd_mgr_for(c::StructLogicCircuit) =
    (vtree(c) isa SddMgr) ? vtree(c) : SddMgr(vtree(c))
sdd_mgr_for(c::LogicCircuit) =
    SddMgr(SddMgr.(Var.(variables(c))), :balanced) # this could be "smarter" finding a good vtree

"Helper function to construct new SDD literal objects in SddMgrLeafNode constructor"
function construct_new_sdd_literal(n::SddMgrLeafNode, l::Lit)::SddLiteralNode
    @assert n.var == lit2var(l) "Cannot compile literal $l respecting vtree leaf for variable $(n.var)"
    SddLiteralNode(l,n)
end


text:
module QuantumInt

include("operators.jl")
include("hamiltonian.jl")
include("observables.jl")
include("absorbers.jl")
include("calc.jl")

end # module


text:
### CONFIGURATION OF THE OCCLUDED INTERSECTION ENVIRONMENT

const LANE_ID_LIST = ["right_from_left", "straight_from_right", "left_from_right", "straight_from_left"]

"""
Container type for the environment parameters
"""
mutable struct SimpleInterParams
    # geometry
    n_lanes::Int64
    hor_roadway_length::Float64
    ver_roadway_length::Float64
    lane_width::Float64
    turn_radius::Float64
    stop_line::Float64
    # cars
    speed_limit::Float64
    car_rate::Float64
end

function SimpleInterParams(;n_lanes::Int64 = 2,
                    hor_roadway_length::Float64 = 20.,
                    ver_roadway_length::Float64 = 10.,
                    lane_width::Float64 = 3.,
                    turn_radius::Float64 = 5.,
                    stop_line::Float64 = 8.,
                    speed_limit::Float64 = 8.,
                    car_rate::Float64 = 0.5)
    return SimpleInterParams(n_lanes, hor_roadway_length, ver_roadway_length, lane_width, turn_radius,
                     stop_line, speed_limit, car_rate)
end

"""
Define the Intersection environment
"""
mutable struct SimpleInterEnv <: OccludedEnv
    roadway::Roadway
    obstacles::Vector{ConvexPolygon}
    params::SimpleInterParams
    lane_map::Dict{String, Lane}
    end_pos::Float64
end

"""
    build_t_shape(params::SimpleInterParams)
generate a T-shape intersection using AutomotiveDrivingModels
"""
function build_t_shape(params::SimpleInterParams)
    roadway = Roadway()
    BEZIER_SAMPLES = 51
    LANE_ID_MAP = Dict{String, Lane}()
    SPEED_LIMIT = SpeedLimit(0., params.speed_limit)
    # Define coordinates of the entry and exit points to the intersection
    r = params.turn_radius # turn radius
    A = VecSE2(0,params.lane_width,-π)
    B = VecSE2(0,0,0)
    B = convert(VecSE2{Float64}, B)
    C = VecSE2(r,-r,-π/2)
    D = VecSE2(r+params.lane_width,-r,π/2)
    E = VecSE2(2r+params.lane_width, 0, 0)
    F = VecSE2(2r+params.lane_width,params.lane_width,-π)

    # Append right turn coming from the left
    curve = gen_straight_curve(convert(VecE2, B+VecE2(-params.hor_roadway_length/2,0)),
                               convert(VecE2, B), 2)
    append_to_curve!(curve, gen_bezier_curve(B, C, 0.6r, 0.6r, BEZIER_SAMPLES)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, C),
                                               convert(VecE2, C+VecE2(0,-params.ver_roadway_length)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    lane.width = params.lane_width
    lane.speed_limit = SPEED_LIMIT
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["right_from_left"] = roadway.segments[1].lanes[1]

    # Append straight right
    curve = gen_straight_curve(convert(VecE2, F+VecE2(params.hor_roadway_length/2,0)),
                               convert(VecE2, F), 2)
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, F), convert(VecE2, A), 2)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, A), convert(VecE2, A+VecE2(-params.hor_roadway_length/2,0)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    lane.width = params.lane_width
    lane.speed_limit = SPEED_LIMIT
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["straight_from_right"] = roadway.segments[2].lanes[1]

    # Append left turn coming from the right
    curve = gen_straight_curve(convert(VecE2, F+VecE2(params.ver_roadway_length,0)), convert(VecE2, F), 2)
    append_to_curve!(curve, gen_bezier_curve(F, C, 0.9r, 0.9r, BEZIER_SAMPLES)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, C), convert(VecE2, C+VecE2(0,-params.ver_roadway_length)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    lane.width = params.lane_width
    lane.speed_limit = SPEED_LIMIT
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["left_from_right"] = roadway.segments[3].lanes[1]

    # Append straight left
    curve = gen_straight_curve(convert(VecE2, B+VecE2(-params.hor_roadway_length/2,0)),
                               convert(VecE2, B), 2)
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, B), convert(VecE2, E), 2)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, E),
                                               convert(VecE2, E+VecE2(params.hor_roadway_length/2,0)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve, width = params.lane_width,speed_limit = SPEED_LIMIT)
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["straight_from_left"] = roadway.segments[4].lanes[1]

    # Append right turn coming from below
    curve = gen_straight_curve(convert(VecE2, D+VecE2(0,-params.ver_roadway_length)), convert(VecE2, D), 2)
    append_to_curve!(curve, gen_bezier_curve(D, E, 0.6r, 0.6r, BEZIER_SAMPLES)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, E), convert(VecE2, E+VecE2(params.ver_roadway_length,0)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    lane.width = params.lane_width
    lane.speed_limit = SPEED_LIMIT
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["ego_right"] = roadway.segments[5].lanes[1]

    # Append left turn coming from below
    curve = gen_straight_curve(convert(VecE2, D+VecE2(0,-params.ver_roadway_length)), convert(VecE2, D), 2)
    append_to_curve!(curve, gen_bezier_curve(D, A, 0.9r, 0.9r, BEZIER_SAMPLES)[2:end])
    append_to_curve!(curve, gen_straight_curve(convert(VecE2, A), convert(VecE2, A+VecE2(-params.hor_roadway_length/2,0)), 2))
    lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
    lane.width = params.lane_width
    lane.speed_limit = SPEED_LIMIT
    push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))
    LANE_ID_MAP["ego_left"] = roadway.segments[6].lanes[1]
    return roadway, LANE_ID_MAP
end

"""
    append_to_curve!(target::Curve, newstuff::Curve)
add a new lane to the roadway
"""
function append_to_curve!(target::Curve, newstuff::Curve)
    s_end = target[end].s
    for c in newstuff
        push!(target, CurvePt(c.pos, c.s+s_end, c.k, c.kd))
    end
    return target
end

function build_obstacles(params::SimpleInterParams)
    HEIGHT = 6
    WIDTH = 6

    #Left Building
    left_building = ConvexPolygon(4)
    top_right = VecE2(3*params.lane_width/4, -params.lane_width/2 - 1.5)
    left_building.pts = [top_right, top_right + VecE2(-WIDTH, 0.),
                    top_right + VecE2(-WIDTH, -HEIGHT), top_right + VecE2(0., -HEIGHT)]
    left_building.npts = 4; # a rectangle

    #Right Building
    right_building = ConvexPolygon(4)
    top_right = top_right + VecE2(top_right.x + 2*params.lane_width + WIDTH, 0.)
    right_building.pts = [top_right, top_right + VecE2(-WIDTH, 0.),
                    top_right + VecE2(-WIDTH, -HEIGHT), top_right + VecE2(0., -HEIGHT)]
    right_building.npts = 4; # a rectangle
    return [left_building, right_building]
end

"""
    CrosswalkEnv(params::SimpleInterParams)
constructor for the crosswalk environment
"""
function SimpleInterEnv(params::SimpleInterParams = SimpleInterParams())
    roadway, lane_id_map = build_t_shape(params)
    obstacles = build_obstacles(params)
    end_pos = round(lane_id_map["ego_left"].curve[end].s)
    return SimpleInterEnv(roadway, obstacles, params, lane_id_map, end_pos)
end


text:
function leak()
	for i=1:100000
		t = Task(eval(:(function() produce() end)))
		consume(t)
		try
		  Base.throwto(t, nothing)
		end
	end
	gc()
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
include("../../GenEvo.jl")

module MNIST

using GenEvo

include("load.jl")

end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
#####
##### Tools for checking profiling
#####

function check_profile(fex::nGraph.FluxExecutable, frame; only_greater = false)
    # Read the profile data from the function
    perf = nGraph.get_performance(fex.ex)

    data = frame.profile_data

    expected_total = 0.0
    actual_total = 0.0
    for node in nodes(data)
        hasprofile(node) || continue
        if nGraph.Lib.can_select_algo(nGraph.getpointer(unx(node)))
            algo_var = frame.model[:algo_var]
            local algo_enum
            for enum in enums(gettime(node))
                if approx_one(algo_var[node, enum])
                    algo_enum = enum
                    break
                end
            end

            actual = perf[nGraph.name(node)]
            expected = timeat(gettime(node), algo_enum)

            # Get the expected move time
            _async = get(frame.model[:tensor_async], nGraph.name(node), nothing)
            if !isnothing(_async)
                async_time = JuMP.value(_async)
            else
                async_time = 0.0
            end

            expected_total += expected
            actual_total += actual

            # Print out the results for this node.
            if !only_greater || actual > expected
                println("Algorithm selection for $(nGraph.name(node)): $algo_enum")
                println("    Actual Time: $(actual)")
                println("    Expected Time: $(expected)")
                println("    Async Time: $(async_time)")
                println()
            end
        else
            actual = perf[nGraph.name(node)]
            expected = gettime(node, getconfig(nGraph.Node(unx(node))))

            # Get the expected move time
            _async = get(frame.model[:tensor_async], nGraph.name(node), nothing)
            if !isnothing(_async)
                async_time = JuMP.value(_async)
            else
                async_time = 0.0
            end

            expected_total += expected
            actual_total += actual

            if !only_greater || actual > expected
                println("No Algorithm selection for $(nGraph.name(node))")
                println("    Actual Time: $(actual)")
                println("    Expected Time: $(expected)")
                println("    Async Time: $(async_time)")
                println()
            end
        end
    end

    @info """
    Expected Total Time: $expected_total
    Actual Total Time: $actual_total
    """

    return nothing
end

function fastest_time(frame)
    data = frame.profile_data

    time = 0.0
    for node in filter(hasprofile, nodes(data))
        if nGraph.Lib.can_select_algo(nGraph.getpointer(unx(node)))
            time += minimum(times(gettime(node)))
        else
            time += gettime(node)
        end
    end

    return time
end

function fastest_time(fex, frame)
    # Take into account actual runtime.
    data = frame.profile_data
    kernel_times = nGraph.get_performance(fex.ex)

    time = 0.0
    for node in filter(hasprofile, nodes(data))
        if nGraph.Lib.can_select_algo(nGraph.getpointer(unx(node)))
            time += min(
                minimum(times(gettime(node))),
                kernel_times[nGraph.name(node)]
            )
        else
            time += min(
                gettime(node),
                kernel_times[nGraph.name(node)]
            )
        end
    end

    return time
end

function show_algorithm_slowdown(frame)
    data = frame.profile_data
    model = frame.model

    for node in Iterators.filter(hasprofile, nodes(data))
        if nGraph.Lib.can_select_algo(nGraph.getpointer(node))
            printstyled("Checking node $(nGraph.name(node))\n"; color = :green)

            # Get the fastest executing algorithm
            time, ind = findmin(get_times(gettime(data, node)))
            enum = get_enums(gettime(data, node))[ind]
            println("    Fastest Enum: $enum. (time) $time")

            # Get the actual used algorithm
            algo_var = frame.model[:algo_var]
            local algo_enum
            for enum in get_enums(gettime(data, node))
                if approx_one(algo_var[node, enum])
                    algo_enum = enum
                    break
                end
            end

            time = get_time(gettime(data, node), algo_enum)
            println("    Actual Enum: $enum. (time) $(time)")
        end
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Test

function runtests(mod)
    for name in names(mod; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        @testset "$(name)" begin
            getfield(mod, name)()
        end
    end
end

@testset "C" begin
    include("C_wrapper.jl")
end

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end

@testset "MathProgBase" begin
    include("MPB_wrapper.jl")
end


text:
function args_train(s)
  @add_arg_table s begin
    "--model-type", "-m"
      help = "model to be used - one of 'unetlike', 'multiscale3', 'fcrna'"
      arg_type = String
      default = "multiscale3"
    "--patch-number", "-n"
      help = "number of patches used for training"
      arg_type = Int
      default = 10
    "--patch-size", "-s"
      help = "number of patches used for training"
      arg_type = Int
      default = 10
    "--test-size", "--testset-size", "-t"
      help = "number of patches used for testing"
      arg_type = Int
      default = 5
    "--batch-size", "-b"
      help = "number of patches per mini-batch during training"
      arg_type = Int
      default = 1
    "--batch-normalization", "-B"
      help = "whether to activate batch normalization"
      action = :store_true
    "--no-gpu"
      help = "prevent usage of gpu acceleration even if gpu support is detected"
      action = :store_true
    "--epochs", "-e"
      help = "number of training epochs"
      arg_type = Int
      default = 10
    "--learning-rate", "-l"
      help = "initial learning rate for the chosen optimizer [-o]"
      arg_type = Float64
      default = 5e-5
    "--noise", "-N"
      help = "magnitude of the applied random noise on the patches (TODO)"
      arg_type = Float64
      default = 0.2
    "--proximity-kernel-size", "-k"
      help = "size of the peaks in the proximity map"
      arg_type = Int
      default = 7
    "--proximity-kernel-height", "-K"
      help = "height of the peaks in the proximity map"
      arg_type = Float64
      default = 100.
    "--optimizer", "-o"
      help = "optimizer used - one of 'adam', 'rmsprop', 'nesterov'"
      arg_type = String
      default = "adam"
    "modelfile"
      required = true
      arg_type = String
      help = "DCellC model (.dccm) file that the trained model is saved to"
    "images"
      nargs='+'
      required = true
      arg_type = String
      help = "image (and corresponding annotation) files that the model trained on"
  end
end

function cmd_train(args)

  # Load labeled images

  lmgs = [ lmgload(imgfile) for imgfile in args["images"] ]

  # Make sure that all images have the same type

  I = typeof(first(lmgs[1]))

  if !all(typeof.(first.(lmgs)) .== I)
    println("Error: all image files must have the same color format")
    exit()
  end

  # Extract some more options

  mtype = args["model-type"]
  otype = args["optimizer"]
  bn    = args["batch-normalization"]

  # Initialize the model

  if mtype == "unetlike"
    model = UNetLike(I, bn=bn)
  elseif mtype == "multiscale3"
    model = Multiscale3(I, bn=bn)
  elseif mtype == "fcrna"
    model = FCRNA(I, bn=bn)
  end

  if otype == "adam"
    opt = Adam
  elseif otype == "rmsprop"
    opt = Rmsprop
  elseif otype == "nesterov"
    opt = Nesterov
  else
    println(STDERR, "Optimizer $(args["optimizer"]) invalid")
  end

  if args["no-gpu"] || gpu() < 0
    println("# CPU mode")
    at = Array{Float32}
  else
    println("# GPU mode")
    at = KnetArray{Float32}
  end

  train!(model, lmgs,
         epochs = args["epochs"],
         modelpath = args["modelfile"],
         patchmode = args["patch-number"],
         patchsize = args["patch-size"],
         batchsize = args["batch-size"],
         peakheight = args["proximity-kernel-height"],
         kernelsize = args["proximity-kernel-size"],
         record = false, at = at,
         lr = args["learning-rate"])
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
"""
Calculate distances, angles, speeds and points for aviation navigation purposes.

All angles are in degrees, distances are standard in meters, and speeds are
standard in m/s.

These methods are not to be used for operational purposes.

Implemented Functions:
* distance
* angular_distance
* bearing
* final_bearing
* midpoint
* intermediate_point
* destination_point
* intersection_point
* along_track_distance
* cross_track_distance
* Vground
* head_wind
* cross_wind
* normalize
* closest_point_to_pole
* isinside

Implemented Types:
* Point(ϕ, λ)
* Airspace

Implemented constants:
* Rₑ_m    Radius Earth in [m]

Based on:
* source: www.movable-type.co.uk/scripts/latlong.html
* source: edwilliams.org/avform.htm
"""
module Navigation

import Base.-, Base.Math.rad2deg, Base.Math.deg2rad

include("utility.jl")
include("airspace.jl")
include("distances.jl")
include("bearings.jl")
include("points.jl")
include("wind.jl")

#TODO Crossing parallels
#TODO Rhumb lines

end # module


text:
using SpikingNetworks

parseParameters"""
τ = 10
σ=100"""

parseEquations"""
dv/dt = (Iₑ(t)-v+σ*ξ(t))/τ
Iₑ(t)=1"""

net=init"""
v=zeros(N)"""

stateMonitor"[1]"

spike"v.>0.8"

reset"v[i]=0"

@time run!(net,100/net.dt)
using MatlabPlot
mplot(net.vRecord[1,:]')


text:
#####################################################################################
#                   Neighborhood.jl Interface & convenience functions               #
#####################################################################################
using Neighborhood, Distances

export WithinRange, NeighborNumber
export Euclidean, Chebyshev, Cityblock

Neighborhood.KDTree(D::AbstractDataset, metric::Metric = Euclidean(); kwargs...) =
KDTree(D.data, metric; kwargs...)

# Convenience extensions for ::Dataset in bulksearches
for f ∈ (:bulkisearch, :bulksearch)
    for nt ∈ (:NeighborNumber, :WithinRange)
        @eval Neighborhood.$(f)(ss::KDTree, D::AbstractDataset, st::$nt, args...; kwargs...) =
        $(f)(ss, D.data, st, args...; kwargs...)
    end
end

#=
    all_neighbors(vtree, vs, ns, K, w)
Return the `maximum(K)`-th nearest neighbors for all input points `vs`,
with indices `ns` in original data, while respecting the theiler window `w`.

This function is nothing more than a convinience call to `Neighborhood.bulksearch`.

It is internal, convenience function.
=#
function all_neighbors(vtree, vs, ns, K, w)
    w ≥ length(vtree.data)-1 && error("Theiler window larger than the entire data span!")
    k = maximum(K)
    tw = Theiler(w, ns)
    idxs, dists = bulksearch(vtree, vs, NeighborNumber(k), tw)
end

"""
    all_neighbors(A::Dataset, stype, w = 0) → idxs, dists
Find the neighbors of all points in `A` using search type `stype` (either
[`NeighborNumber`](@ref) or [`WithinRange`](@ref)) and `w` the [Theiler window](@ref).

This function is nothing more than a convinience call to `Neighborhood.bulksearch`.
"""
function all_neighbors(A::AbstractDataset, stype, w::Int = 0)
    theiler = Theiler(w)
    tree = KDTree(A)
    idxs, dists = bulksearch(tree, A, stype, theiler)
end


text:
using LightGraphs, JuMP, Cbc, Plots
gr(size=(600,600))

function dollar_game(g, W)
  nv = length(vertices(g))

  model = Model(solver = CbcSolver())
  @variable(model, give[1:nv] >= 0, Int)
  @variable(model, take[1:nv] >= 0, Int)
  @objective(model, Min, sum(give[i] + take[i] for i = 1:nv))
  @expression(model, x[i=1:nv], W[i] +
              (take[i] - give[i]) * length(neighbors(g, i)) +
              sum(give[j] - take[j] for j = neighbors(g, i)))
  @constraint(model, x .>= 0)
  status = solve(model)

  println("status = $status")
  if status != :Optimal
    return zeros(nv), W
  end

  give = Int.(getvalue(give))
  take = Int.(getvalue(take))
  sol = Int.(getvalue(x))
  return give - take, sol
end

function plot_common!(p, g, W, game=[])
  nv = length(vertices(g))
  θ = 2π / nv
  t = (1:nv) * θ - 1.5θ
  r = 2.5
  x, y = r * sin.(t), r * cos.(t)

  for i = 1:nv
    for j = neighbors(g, i)
      plot!([x[i], x[j]], [y[i], y[j]], lw=2, c=:black)
    end
  end
  scatter!(x, y, ann=[(x[i],y[i],W[i]) for i = 1:nv], m=(30,:white,stroke(2,:black)))
  if length(game) > 0
    annotate!([(1.2x[i], 1.2y[i], game[i], :red) for i = 1:nv])
  end
  xlims!(-1.3r, 1.3r)
  ylims!(-1.3r, 1.3r)
  return p
end

function plot_common(g, W, game=[])
  p = plot(leg=false, axis=false, grid=false)
  return plot_common!(p, g, W, game)
end

function plot_graph(g, W, filename, game=[])
  p = plot_common(g, W, game)
  png(filename)
end

function plot_game(g, W, prefix, game)
  p = plot_common(g, W)
  nv = length(vertices(g))
  k = 0
  fname = @sprintf("%s-%03d", prefix, k)
  title!("Initial")
  png(p, fname)

  θ = 2π / nv
  t = (1:nv) * θ - 1.5θ
  r = 2.5
  x, y = r * sin.(t), r * cos.(t)

  for i = 1:nv
    for n = 1:abs(game[i])
      k += 1
      ng = neighbors(g, i)
      s = sign(game[i])
      W[i] -= s * length(ng)
      c = s > 0 ? RGB(0.1, 0.4, 0.9) : :red

      W[ng] .+= s

      p = plot(leg=false, axis=false, grid=false)
      plot_common!(p, g, W)
      for j = ng
        v = [x[j] - x[i]; y[j] - y[i]]
        σ = 0.4 / norm(v)
        X = [x[i] + v[1] * σ; x[j] - v[1] * σ]
        Y = [y[i] + v[2] * σ; y[j] - v[2] * σ]
        if s < 0
          X, Y = X[[2;1]], Y[[2;1]]
        end
        plot!(X, Y, c=c, lw=3, l=:arrow)
      end
      scatter!([x[i]], [y[i]], m=(30,c,stroke(3,:black)))
      scatter!([x[ng]], [y[ng]], m=(30,c,stroke(3,:black)))
      I = setdiff(1:nv, [i; ng])
      scatter!([x[I]], [y[I]], m=(30,:white,stroke(2,:black)))

      fname = @sprintf("%s-%03d", prefix, k)
      title!("Move $k")
      png(fname)
    end
  end
end

function example()
  g = Graph(5)
  W = [1; -2; 3; 2; -1]
  add_edge!(g, 1, 2)
  add_edge!(g, 1, 3)
  add_edge!(g, 2, 3)
  add_edge!(g, 2, 5)
  add_edge!(g, 3, 4)

  plot_graph(g, W, "example")
  game, sol = dollar_game(g, W)
  plot_graph(g, sol, "example-sol", game)
  plot_game(g, W, "example", game)
end

function example2()
  g = Graph(3)
  W = [0; 1; -1]
  add_edge!(g, 1, 2)
  add_edge!(g, 1, 3)
  add_edge!(g, 2, 3)

  game, sol = dollar_game(g, W)
end

function example3()
  g = Graph(3)
  W = [-2; -1; 4]
  add_edge!(g, 1, 2)
  add_edge!(g, 2, 3)

  game, sol = dollar_game(g, W)
end

example()
example2()
example3()


text:
include("sand_04.jl")
using .Sand_04
using Plots
fig = Plots.plot()

NN = 20 # no of cells
tmax = 0.4 # max time
dt = 0.4 # how often to plot
tp = 0.0#5000 # when to start plotting

r= Sand_04.ring(NN; U=5.5)
for t in 0.0:dt:tmax
#for i in 1:10
    t = r.t
    if t ≥ tp
        x0_ = 0.5 * (r.X_[1:end-1] + r.X_[2:end])
#        Plots.plot!(x0_[1:end],r.Σ_[1:end],label="t: $t",legend=false)#:bottomright)#,title="Σ")#;legend=false)
#        display(fig)
        Plots.plot!(x0_[1:end],r.J_[1:end],label=false)#,label="t: $t",legend=:topleft)#,title="J")#;legend=false)
        display(fig)
    end
    Sand_04.tstep!(r,dt)
    #Sand_03.step!(r)
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
using Images, MXNet

### LOADING THE MODEL
const MODEL_NAME = "weights/mobilenet-v2/mobilenet_v2"
const MODEL_CLASS_NAMES = "weights/mobilenet-v2/synset.txt"

nnet = mx.load_checkpoint(MODEL_NAME, 0, mx.FeedForward; context = mx.gpu());
synset = readlines(MODEL_CLASS_NAMES);

### SEARCH FOR A LAYER OF INTERESTS
layers = mx.get_internals(nnet.arch);
layers_flatten = nothing
layers_to_remove = Symbol[]

# We iterate over all layers until we find the one matching our requirements
# and remove the ones to follow after
for i = 1:2000

    layer = layers[i];
    layer_name = mx.get_name(layer)

    if layers_flatten == nothing && layer_name == :pool6
        layers_flatten = layer
    elseif layers_flatten != nothing
        push!(layers_to_remove, layer_name)
        if layer_name in [:softmax, :label, :prob] break end
    end
end

nnet.arch = @mx.chain layers_flatten => Flatten()
map(x -> delete!(nnet.arg_params, x), layers_to_remove);
map(x -> delete!(nnet.aux_params, x), layers_to_remove);

mx.save_checkpoint(nnet, "weights/mobilenet-v2/MobiletNet-FE", mx.OptimizationState(1, 0, 0, 0))


text:
using RecipesBase
using Base.Test

srand(1)
RecipesBase.is_key_supported(k::Symbol) = true

type T end

@recipe function plot{N<:Integer}(t::T, n::N = 1; customcolor = :green)
    :markershape --> :auto, :require
    :markercolor --> customcolor, :force
    :xrotation --> 5
    :zrotation --> 6, :quiet
    rand(10,n)
end

# this is similar to how Plots would call the method
typealias KW Dict{Symbol,Any}
d = KW(:customcolor => :red)
data_list = RecipesBase.apply_recipe(d, T(), 2)

# make sure the attribute dictionary was populated correctly, and the returned arguments are as expected
@test data_list[1].args == (srand(1); (rand(10,2),))
@test d == KW(
	:customcolor => :red,
    :markershape => :auto,
    :markercolor => :red,
    :xrotation => 5,
    :zrotation => 6
)


text:
include("../../../src/Sparso.jl")
include("../../../src/simple-show.jl")
using Sparso
using CompilerTools.OptFramework
using MatrixMarket
using MatrixMarket2

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_REPLACE_CALLS)

include("dss.jl")

function speye_int32(m::Integer)
  rowval = [Int32(x) for x in [1:m;]]
  colptr = [Int32(x) for x in [rowval; m + 1]]
  nzval = ones(Float64, m)
  return SparseMatrixCSC(m, m, colptr, rowval, nzval)
end

# A*D*B, where D is a diagonal matrix
# Since it's implemented in CSR, we pass tranposes of A and B
# in CSC
function adb_inspect(AT::SparseMatrixCSC, BT::SparseMatrixCSC)
  csrA = Sparso.create_CSR(AT)
  csrB = Sparso.create_CSR(BT)

  # A*D*A' inspection
  csrADAT = ccall((:CSR_ADBInspect, LIB_CSR_PATH), Ptr{Void},
                  (Ptr{Void}, Ptr{Void}),
                  csrA, csrB)

  m = size(A, 1)
  rowptr = pointer_to_array(
    ccall((:CSR_GetRowPtr, LIB_CSR_PATH), Ptr{Cint},
          (Ptr{Void},),
          csrADAT),
    (m + 1,))
  nnz = rowptr[m + 1] - 1
  colidx = pointer_to_array(
    ccall((:CSR_GetColIdx, LIB_CSR_PATH), Ptr{Cint},
          (Ptr{Void},),
          csrADAT),
    (nnz,))
  values = pointer_to_array(
    ccall((:CSR_GetValues, LIB_CSR_PATH), Ptr{Cdouble},
          (Ptr{Void},),
          csrADAT),
    (nnz,))

  Sparso.destroy_CSR(csrA)
  Sparso.destroy_CSR(csrB)

  ADB = SparseMatrixCSC{Cdouble, Cint}(
    m, m, rowptr, colidx, values)
end

# A*D*B, where D is a diagonal matrix
function adb_execute!(ADB::SparseMatrixCSC, AT::SparseMatrixCSC, BT::SparseMatrixCSC, d::Vector)
  csrA = Sparso.create_CSR(AT)
  csrB = Sparso.create_CSR(BT)
  csrADB = Sparso.create_CSR(ADB)

  ccall((:CSR_ADB, LIB_CSR_PATH), Void,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Cdouble}),
        csrADB, csrA, csrB, d)

  Sparso.destroy_CSR(csrA)
  Sparso.destroy_CSR(csrB)
  Sparso.destroy_CSR(csrADB)
end

# optimized implementation of interior-point method with inspector hoisted
function ipm(A, b, p, print = true) # A: constraint coefficients, b: constraint rhs, p: objective
  t0 = time()
  (m,n) = size(A)

# set initial point, based on largest element in (A,b,p)
  bigM = maximum(A)
  bigM = maximum([norm(b, Inf) norm(p, Inf) bigM])
  x = 100*bigM*ones(n)
  s = x
  y = zeros(m)

  bc = 1 + maximum([norm(b) norm(p)])

  relResidual = NaN
  iter=1

  blas1_time = 0.
  spgemm_time = 0.
  fact_time = 0.
  trslv_time = 0.
  spmv_time = 0.

  spmv_inspector_time = 0.
  spgemm_inspector_time = 0.
  fact_inspector_time = 0.

  D = speye_int32(n)

  spmv_inspector_time -= time()
  AT = A'
  spmv_inspector_time += time()

  spgemm_inspector_time -= time()
  ADAT = adb_inspect(AT, A)
  spgemm_inspector_time += time()

  fact_inspector_time -= time()

  dss_handle = dss_analyze(ADAT)
  dy = Array(Cdouble, m)

  fact_inspector_time += time()

  for iter=1:200

    # compute residuals
    spmv_time -= time()
    Rd = AT*y + s - p
    Rp = A*x - b
    spmv_time += time()

    blas1_time -= time()
    Rc = x.*s
    mu = mean(Rc)
    relResidual = norm([Rd; Rp; Rc])/bc
    blas1_time += time()

    if print
      #@printf "iter %2i: log10(mu) = %5.2f, resid = %9.2e, obj = %9.2e\n" iter log10(mu) relResidual (p'*x)[1][1]
    end
    if (relResidual <= 1e-7 && mu <= 1e-7) break; end

    blas1_time -= time()
    Rc = Rc - min(0.1, 100*mu)*mu

    # set up the scaling matrix, and form the coefficient matrix for
    # the linear system
    d = min(5.e+15, x./s)
    blas1_time += time()

    spgemm_time -= time()
    adb_execute!(ADAT, AT, A, d)
    spgemm_time += time()

    # use the form of the Cholesky routine "cholinc" that's best
    # suited to interior-point methods
    fact_time -= time()
    dss_factor(dss_handle, ADAT)
    fact_time += time()

    # set up the right-hand side
    blas1_time -= time()
    t1 = x.*Rd - Rc;
    blas1_time += time()

    spmv_time -= time()
    t2 = -(Rp + A*(t1./s));
    spmv_time += time()

    # solve it and recover the other step components
    trslv_time -= time()
    opt = MKL_DSS_DEFAULTS
    dss_solve!(dss_handle, t2, dy)
    trslv_time += time()

    spmv_time -= time()
    temp = AT*dy
    spmv_time += time()

    blas1_time -= time()
    dx = (x.*temp + t1)./s
    ds = -(s.*dx + Rc)./x

    tau = max(.9995, 1 - mu)
    ap = -1/minimum([dx./x; -1])
    ad = -1/minimum([ds./s; -1])
    ap = tau*ap
    ad = tau*ad
    x = x + ap*dx
    s = s + ad*ds
    y = y + ad*dy
    blas1_time += time()
  end

  blas1_time -= time()
  f = p'*x
  blas1_time += time()

  if print
    @printf "\nwith_inspector_hoisting_total_time = %f\n" time() - t0
    @printf "spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
    @printf "spgemm_inspector = %f fact_inspector_time = %f spmv_inspector = %f\n" spgemm_inspector_time fact_inspector_time spmv_inspector_time
    @printf "iter %2i, resid = %9.2e, objval = %e\n" iter relResidual f[1][1]
  end

  x
end

# unoptimized implementation of interior-point method w/o inspector hoisting
function ipm_unopt(A, b, p, print = true) # A: constraint coefficients, b: constraint rhs, p: objective
  t0 = time()
  (m,n) = size(A)

# set initial point, based on largest element in (A,b,p)
  bigM = maximum(A)
  bigM = maximum([norm(b, Inf) norm(p, Inf) bigM])
  x = 100*bigM*ones(n)
  s = x
  y = zeros(m)

  bc = 1 + maximum([norm(b) norm(p)])

  relResidual = NaN
  iter=1

  blas1_time = 0.
  spgemm_time = 0.
  fact_time = 0.
  trslv_time = 0.
  spmv_time = 0.

  D = speye_int32(n)
  dy = Array(Cdouble, m)

  for iter=1:200

    # compute residuals
    spmv_time -= time()
    Rd = A'*y + s - p
    Rp = A*x - b
    spmv_time += time()

    blas1_time -= time()
    Rc = x.*s
    mu = mean(Rc)
    relResidual = norm([Rd; Rp; Rc])/bc
    blas1_time += time()

    if print
      #@printf "iter %2i: log10(mu) = %5.2f, resid = %9.2e, obj = %9.2e\n" iter log10(mu) relResidual (p'*x)[1][1]
    end
    if (relResidual <= 1e-7 && mu <= 1e-7) break; end

    blas1_time -= time()
    Rc = Rc - min(0.1, 100*mu)*mu

    # set up the scaling matrix, and form the coefficient matrix for
    # the linear system
    d = min(5.e+15, x./s)
    blas1_time += time()

    spgemm_time -= time()
    ADAT = adb_inspect(A', A)
    adb_execute!(ADAT, A', A, d)
    spgemm_time += time()

    # use the form of the Cholesky routine "cholinc" that's best
    # suited to interior-point methods
    fact_time -= time()
    dss_handle = dss_analyze(ADAT)
    dss_factor(dss_handle, ADAT)
    fact_time += time()

    # set up the right-hand side
    blas1_time -= time()
    t1 = x.*Rd - Rc;
    blas1_time += time()

    spmv_time -= time()
    t2 = -(Rp + A*(t1./s));
    spmv_time += time()

    # solve it and recover the other step components
    trslv_time -= time()
    opt = MKL_DSS_DEFAULTS
    dss_solve!(dss_handle, t2, dy)
    trslv_time += time()

    spmv_time -= time()
    temp = A'*dy
    spmv_time += time()

    blas1_time -= time()
    dx = (x.*temp + t1)./s
    ds = -(s.*dx + Rc)./x

    tau = max(.9995, 1 - mu)
    ap = -1/minimum([dx./x; -1])
    ad = -1/minimum([ds./s; -1])
    ap = tau*ap
    ad = tau*ad
    x = x + ap*dx
    s = s + ad*ds
    y = y + ad*dy
    blas1_time += time()
  end

  blas1_time -= time()
  f = p'*x
  blas1_time += time()

  if print
    @printf "\nwo_inspector_hoisting_total_time = %f\n" time() - t0
    @printf "spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
    @printf "iter %2i, resid = %9.2e, objval = %e\n" iter relResidual f[1][1]
  end

  x
end

# a reference implementation of interior-point method
# that uses Julia's default SpGEMM and Cholesky
function ipm_ref(A, b, p, print = true) # A: constraint coefficients, b: constraint rhs, p: objective
  t0 = time()
  (m,n) = size(A)

# set initial point, based on largest element in (A,b,p)
  bigM = maximum(A)
  bigM = maximum([norm(b, Inf) norm(p, Inf) bigM])
  x = 100*bigM*ones(n)
  s = x
  y = zeros(m)

  bc = 1 + maximum([norm(b) norm(p)])

  relResidual = NaN
  iter=1

  blas1_time = 0.
  spgemm_time = 0.
  fact_time = 0.
  trslv_time = 0.
  spmv_time = 0.

  D = speye_int32(n)

  for iter=1:200

    # compute residuals
    spmv_time -= time()
    Rd = A'*y + s - p
    Rp = A*x - b
    spmv_time += time()

    blas1_time -= time()
    Rc = x.*s
    mu = mean(Rc)
    relResidual = norm([Rd; Rp; Rc])/bc
    blas1_time += time()

    if print
      #@printf "iter %2i: log10(mu) = %5.2f, resid = %9.2e, obj = %9.2e\n" iter log10(mu) relResidual (p'*x)[1][1]
    end
    if (relResidual <= 1e-7 && mu <= 1e-7) break; end

    blas1_time -= time()
    Rc = Rc - min(0.1, 100*mu)*mu

    # set up the scaling matrix, and form the coefficient matrix for
    # the linear system
    d = min(5.e+15, x./s)
    blas1_time += time()

    spgemm_time -= time()
    D.nzval = d
    B = A*D*A'
    spgemm_time += time()

    # use the form of the Cholesky routine "cholinc" that's best
    # suited to interior-point methods
    fact_time -= time()
    R = cholfact(SparseMatrixCSC{Float64, Int64}(B))
    fact_time += time()

    # set up the right-hand side
    blas1_time -= time()
    t1 = x.*Rd - Rc;
    blas1_time += time()

    spmv_time -= time()
    t2 = -(Rp + A*(t1./s));
    spmv_time += time()

    # solve it and recover the other step components
    trslv_time -= time()
    dy = R\t2
    trslv_time += time()

    spmv_time -= time()
    temp = A'*dy
    spmv_time += time()

    blas1_time -= time()
    dx = (x.*temp + t1)./s
    ds = -(s.*dx + Rc)./x

    tau = max(.9995, 1 - mu)
    ap = -1/minimum([dx./x; -1])
    ad = -1/minimum([ds./s; -1])
    ap = tau*ap
    ad = tau*ad
    x = x + ap*dx
    s = s + ad*ds
    y = y + ad*dy
    blas1_time += time()
  end

  blas1_time -= time()
  f = p'*x
  blas1_time += time()

  if print
    @printf "\nref_total_time = %f\n" time() - t0
    @printf "spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
    @printf "iter %2i, resid = %9.2e, objval = %e\n" iter relResidual f[1][1]
  end

  x
end

# min 2*x1 + x2 subject to x1 + x2 = 1, x1 >= 0, x2 >= 0
# expected solution: x1 = 0, x2 = 1, obj = 1
#A = sparse([1 1])
#b = [ 1 ]'
#p = [ 2 1 ]'
A = MatrixMarket2.mmread(string(ARGS[1], "-A.mtx"))'
b = vec(MatrixMarket.mmread(string(ARGS[1], "-b.mtx")))
p = vec(MatrixMarket.mmread(string(ARGS[1], "-p.mtx")))

m = size(A, 1)
n = size(A, 2)
println("Problem size = [$m $n]")

#print_time = false
#@acc ipm_ref(A, b, p, print_time) # ignore timing of the first run
#ipm_ref(A, b, p)

ipm_unopt(A, b, p, false)
ipm_unopt(A, b, p)

#ipm(A, b, p, false)
#ipm(A, b, p)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
export   functionalf!,
         functionalg!,
         functionalg1D!,
         inverse!,
         inversehybrid!,
         inverse1D!

# function functionalf!(F, xk, cache, cache_vander, ψoff, output::Array{Float64,1}, R::IntegratedFunction)

function functionalf!(F, xk, cache, cache_vander, ψoff, output, R::IntegratedFunction)
    function integrand!(v::Vector{Float64}, t::Float64)
        repeated_grad_xk_basis!(cache, cache_vander, R.f,  t*xk)
        # cache .= repeated_grad_xk_basis(R.f,  t*xk)
        @avx @. v = (cache .* ψoff) *ˡ R.f.coeff
        evaluate!(v, R.g, v)
        # v .= R.g((repeated_grad_xk_basis(R.f,  t*xk) .* ψoff)*R.f.coeff)
    end

    # quadgk!(integrand!, F, 0.0, 1.0)
    quadgk!(integrand!, F, 0.0, 1.0)
    F .*= xk
    F .-= output
    nothing
end

# functionalf!(cache, cache_vander, ψoff, output::Array{Float64,1}, R::IntegratedFunction) =

functionalf!(cache, cache_vander, ψoff, output, R::IntegratedFunction) =
    (F, xk) -> functionalf!(F, xk, cache, cache_vander, ψoff, output, R)

function functionalg!(J, xk, cache, cache_vander, ψoff, output, R::IntegratedFunction)
    # cache .= repeated_grad_xk_basis(R.f,  xk)
    repeated_grad_xk_basis!(cache, cache_vander, R.f,  xk)
    @avx @. J.diag = (cache .* ψoff) *ˡ R.f.coeff
    evaluate!(J.diag, R.g, J.diag)
    nothing
end

functionalg!(cache, cache_vander, ψoff, output, R::IntegratedFunction) =
    (J, xk) -> functionalg!(J, xk, cache, cache_vander, ψoff, output, R)

# In this version, the Jacobian is stored as a vector (it is diagonal in this case).
function functionalg1D!(J::AbstractVector{Float64}, xk, cache, cache_vander, ψoff, output, R::IntegratedFunction)
    # cache .= repeated_grad_xk_basis(R.f,  xk)
    repeated_grad_xk_basis!(cache, cache_vander, R.f,  xk)
    @avx @. J = (cache .* ψoff) *ˡ R.f.coeff
    evaluate!(J, R.g, J)
    nothing
end

functionalg1D!(cache, cache_vander, ψoff, output, R::IntegratedFunction) =
    (J, xk) -> functionalg1D!(J, xk, cache, cache_vander, ψoff, output, R)

# The state is modified in-place
# function inverse!(X::Array{Float64,2}, F, R::IntegratedFunction, S::Storage)
"""
$(TYPEDSIGNATURES)

Inverts in-place the ensemble matrix `X` such that `C(X) = F`.
"""
function inverse!(X, F, R::IntegratedFunction, S::Storage)
    Nψ = R.Nψ
    Nx = R.Nx
    NxX, Ne = size(X)

    @assert NxX == R.Nx "Wrong dimension of the sample X"

    cache  = zeros(Ne, Nψ)
    cache_vander = zeros(Ne, maximum(R.f.idx[:,Nx])+1)
    f0 = zeros(Ne)

    # Remove f(x_{1:k-1},0) from the output F
    @avx for i=1:Ne
        f0i = zero(Float64)
        for j=1:Nψ
            f0i += (S.ψoffψd0[i,j])*R.f.coeff[j]
        end
        F[i] -= f0i
    end

    df_inverse = OnceDifferentiable(functionalf!(cache, cache_vander, S.ψoff, F, R),
                                    functionalg!(cache, cache_vander, S.ψoff, F, R),
                                    f0, f0, Diagonal(f0))

    # Start minimization from the prior value
    result = nlsolve(df_inverse, X[end,:]; method = :newton, linesearch = LineSearches.MoreThuente());

    # Check convergence
    # @show converged(result)
    if converged(result) == true
    # @assert converged(result) "Optimization hasn't converged"
        X[end,:] .= result.zero
    else
        println("Optimization hasn't converged")
    end
end

"""
$(TYPEDSIGNATURES)

Inverts in-place the ensemble matrix `X` such that `C(X) = F`.
"""
inverse!(X, F, C::HermiteMapComponent, S::Storage) = inverse!(X, F, C.I, S)

function inverse!(X::Array{Float64,2}, F, L::LinHermiteMapComponent, S::Storage)
    # Pay attention that S is computed in the renormalized space for improved performance !!!
    transform!(L.L, X)
    inverse!(X, F, L.C, S)
    itransform!(L.L, X)
end

# function functionalfg!(F, J, xk, cache, ψoff, output::Array{Float64,1}, R::IntegratedFunction{m, Nψ, Nx}) where {m, Nψ, Nx}
#
#     if !(F == nothing)
#         function integrand!(v::Vector{Float64}, t::Float64)
#             v .= R.g((repeated_grad_xk_basis(R.f,  t*xk) .* ψoff)*R.f.coeff)
#         end
#         F .= xk .* quadgk!(integrand!, cache, 0, 1)[1]
#         F .-= output
#     end
#
#     if !(J == nothing)
#         J .= Diagonal(R.g((repeated_grad_xk_basis(R.f,  xk) .* ψoff)*R.f.coeff))
#     end
# end
#
#
# functionalfg!(cache, ψoff, output::Array{Float64,1}, R::IntegratedFunction{m, Nψ, Nx}) where {m, Nψ, Nx} =
#     (F, J, xk) -> functionalfg!(F, J, xk, cache, ψoff, output, R)


text:
# return info comment strings for UF sparse matrix
if isdefined(:SparseArrays)
    Sparse = Base.SparseArrays.CHOLMOD.Sparse
else
    Sparse = Base.SparseMatrix.CHOLMOD.Sparse
end

function ufinfo(filename::AbstractString)
    mmfile = open(filename,"r")
    info = "\n"
    ll = readline(mmfile)
    while length(ll) > 0 && ll[1] == '%'
        info = string(info, ll)
        ll = readline(mmfile)
    end
    return info
end

# read Matrix Market data
"""
`mmreader(dir, name, info)`

dir: directory of the file
name: file name
info: whether to return infomation
"""
function mmreader(dir::AbstractString, name::AbstractString; info::Bool = true)
    pathfilename = string(dir, '/', name, ".mtx")
    if info
        println(ufinfo(pathfilename))
        println("use matrixdepot(\"$name\", :read) to read the data")
    else
        sparse(Sparse(pathfilename))
    end
end

# read UF sparse matrix data
"""
`ufreader(dir, name, info, meta)`

dir: directory of the file
name: file name
info: whether to return information
meta: wehterh to return metadata
"""
function ufreader(dir::AbstractString, name::AbstractString;
                  info::Bool = true, meta::Bool = false)
    dirname = string(dir, '/', name)
    files = filenames(dirname)
    if info
        println(ufinfo(string(dirname, '/', name, ".mtx")))
        if length(files) > 1
            println("metadata:")
            display(files)
        end
        #println("use matrixdepot(\"$name\", :read) to read the data")
    else
        A = sparse(Sparse(string(dirname, '/', name, ".mtx")))
        if meta
            metadict = Dict{AbstractString, Any}()
            datafiles = readdir(dirname)
            for data in datafiles
                dataname = split(data, '.')[1]
                if endswith(data, "mtx")
                    try
                        metadict[dataname] =  sparse(Sparse(string(dirname, '/', data)))
                    catch
                        metadict[dataname] = denseread(string(dirname,'/', data))
                    end
                else
                    metadict[dataname] = readall(string(dirname,'/', data))
                end
            end
            metadict
        else
            A
        end
    end
end


text:
using Crayons
Crayons.force_color(true)

const FORM_C = string(Crayon(foreground = :light_blue))
const HEAD_C = string(Crayon(foreground = :light_yellow))
const NODE_C = string(Crayon(foreground = :light_green))
const RES_C = string(Crayon(reset = true))

abstract type View end

struct FormView <: View
    size::Tuple{Int, Int}
    components::Vector{View}
end

function paint(form_view::FormView; state=:normal)
    h, w = form_view.size

    T.buffered() do buffer
        T.print(buffer, FORM_C)
        T.cmove(buffer, 1, 1)
        T.print(buffer, Char(0x250C)) # ┌
        T.print(buffer, Char(0x2500)^(w-2))
        T.print(buffer, Char(0x2510)) # ┐
        T.cmove_line_down(buffer)

        for i=2:(h-1)
            T.print(buffer, Char(0x2502))
            T.print(buffer, ' '^(w-2))
            T.print(buffer, Char(0x2502))
            T.cmove_line_down(buffer)
        end

        T.print(buffer, Char(0x2514)) # └
        T.print(buffer, Char(0x2500)^(w-2))
        T.print(buffer, Char(0x2518)) # ┘

        T.cmove(buffer, 1, 2)
        T.print(buffer, " Snake Game ")

        T.print(buffer, RES_C)
    end

    if state === :lose
        str = "You Lose ~~"
        T.cmove(trunc(Int, h/2), trunc(Int, (w-textwidth(str))/2))
        T.print(str)
        return
    end

    for component in form_view.components
        paint(component)
    end
end

struct SnakeView <: View
    snake::SnakeModel
end

function paint(snake_view::SnakeView)
    T.buffered() do buffer
        T.print(HEAD_C)
        for (i, node_pos) in enumerate(snake_view.snake.node_anchors)
            y, x = node_pos
            T.cmove(buffer, y, x)
            T.csave(buffer)
            T.print(buffer, "┌┐")
            T.crestore(buffer)
            T.cmove_down(buffer)
            T.print(buffer, "└┘")

            (i == 1) && (T.print(buffer, NODE_C))
        end
        T.print(buffer, RES_C)
    end
end


text:
#Run this in Escher to see a rotating cube in the browser.

import ThreeJS
using FileIO

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    x = Signal(1.0)
    y = Signal(1.0)
    z = Signal(1.0)
    mesh_geom = load("cat.obj")
    vbox(
        "ThreeJS Example",
        vskip(2em),
        vbox(
            "Scale",
            hbox("x",slider(1.0:5.0) >>> x),
            hbox("y",slider(1.0:5.0) >>> y),
            hbox("z",slider(1.0:5.0) >>> z),
        ),
        vskip(2em),
        map(x,y,z) do x,y,z
        ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                    ThreeJS.mesh(0.0, 0.0, 0.0) <<
                    [
                        ThreeJS.geometry(mesh_geom, scale=(x,y,z)), ThreeJS.material(Dict(:kind=>"lambert",:color=>"red"))
                    ],
                    ThreeJS.pointlight(10.0, 10.0, 10.0),
                    ThreeJS.camera(2.0, 2.0, 2.0)
                ]
            )
        end
    ) |> pad(2em)
end