text:
"""
    SparseMatrixCLIL{T, Ti}

The SparseMatrixCLIL represents a sparse matrix in two distinct ways:

1. As a sparse (in both row and column) n x m matrix
2. As a row-dense, column-sparse k x m matrix

The data structure keeps a permutation between the row order of the two representations.
Swapping the rows in one does not affect the other.

On construction, the second representation is equivalent to the first with fully-sparse
rows removed, though this may cease being true as row permutations are being applied
to the matrix.

The default structure of the `SparseMatrixCLIL` type is the second structure, while
the first is available via the thin `AsSubMatrix` wrapper.
"""
struct SparseMatrixCLIL{T, Ti<:Integer} <: AbstractSparseMatrix{T, Ti}
    nparentrows::Int
    ncols::Int
    nzrows::Vector{Ti}
    row_cols::Vector{Vector{Ti}} # issorted
    row_vals::Vector{Vector{T}}
end
Base.size(S::SparseMatrixCLIL) = (length(S.nzrows), S.ncols)
Base.copy(S::SparseMatrixCLIL{T, Ti}) where {T, Ti} =
    SparseMatrixCLIL(S.nparentrows, S.ncols, copy(S.nzrows), map(copy, S.row_cols), map(copy, S.row_vals))
function swaprows!(S::SparseMatrixCLIL, i, j)
    swap!(S.nzrows, i, j)
    swap!(S.row_cols, i, j)
    swap!(S.row_vals, i, j)
end

function SparseMatrixCLIL(mm::AbstractMatrix)
    nrows, ncols = size(mm)
    row_cols = [findall(!iszero, row) for row in eachrow(mm)]
    row_vals = [row[cols] for (row, cols) in zip(eachrow(mm), row_cols)]
    SparseMatrixCLIL(nrows, ncols, Int[1:length(row_cols);], row_cols, row_vals)
end

struct CLILVector{T, Ti} <: AbstractSparseVector{T, Ti}
    vec::SparseVector{T, Ti}
end
Base.size(v::CLILVector) = Base.size(v.vec)
Base.getindex(v::CLILVector, idx::Integer...) = Base.getindex(v.vec, idx...)
Base.setindex!(vec::CLILVector, v, idx::Integer...) = Base.setindex!(vec.vec, v, idx...)
Base.view(a::SparseMatrixCLIL, i::Integer, ::Colon) =
    CLILVector(SparseVector(a.ncols, a.row_cols[i], a.row_vals[i]))
SparseArrays.nonzeroinds(a::CLILVector) = SparseArrays.nonzeroinds(a.vec)
SparseArrays.nonzeros(a::CLILVector) = SparseArrays.nonzeros(a.vec)

function Base.setindex!(S::SparseMatrixCLIL, v::CLILVector, i::Integer, c::Colon)
    if v.vec.n != S.ncols
        throw(BoundsError(v, 1:S.ncols))
    end
    S.row_cols[i] = copy(v.vec.nzind)
    S.row_vals[i] = copy(v.vec.nzval)
    return v
end

zero!(a::AbstractArray{T}) where {T} = a[:] .= zero(T)
zero!(a::SparseVector) = (empty!(a.nzind); empty!(a.nzval))
zero!(a::CLILVector) = zero!(a.vec)

struct NonZeros{T <: AbstractArray}
    v::T
end
Base.pairs(nz::NonZeros{<:CLILVector}) = NonZerosPairs(nz.v)

struct NonZerosPairs{T <: AbstractArray}
    v::T
end

# N.B.: Because of how we're using this, this must be robust to modification of
# the underlying vector. As such, we treat this as an iteration over indices
# that happens to short cut using the sparse structure and sortedness of the
# array.
function Base.iterate(nzp::NonZerosPairs{<:CLILVector}, (idx, col))
    v = nzp.v.vec
    nzind = v.nzind
    nzval = v.nzval
    if idx > length(nzind)
        idx = length(col)
    end
    oldcol = nzind[idx]
    if col !== oldcol
        # The vector was changed since the last iteration. Find our
        # place in the vector again.
        tail = col > oldcol ? (@view nzind[idx+1:end]) : (@view nzind[1:idx])
        tail_i = searchsortedfirst(tail, col + 1)
        # No remaining indices.
        tail_i > length(tail) && return nothing
        new_idx = col > oldcol ? idx + tail_i : tail_i
        new_col = nzind[new_idx]
        return (new_col=>nzval[new_idx], (new_idx, new_col))
    end
    idx == length(nzind) && return nothing
    new_col = nzind[idx+1]
    return (new_col=>nzval[idx+1], (idx+1, new_col))
end

function Base.iterate(nzp::NonZerosPairs{<:CLILVector})
    v = nzp.v.vec
    nzind = v.nzind
    nzval = v.nzval
    isempty(nzind) && return nothing
    return nzind[1]=>nzval[1], (1, nzind[1])
end

# Arguably this is how nonzeros should behave in the first place, but let's
# build something that works for us here and worry about it later.
nonzerosmap(a::CLILVector) = NonZeros(a)

function bareiss_update_virtual_colswap_mtk!(zero!, M::SparseMatrixCLIL, k, swapto, pivot, last_pivot; pivot_equal_optimization=true)
    # for ei in nzrows(>= k)
    eadj = M.row_cols
    old_cadj = M.row_vals
    vpivot = swapto[2]

    ## N.B.: Micro-optimization
    #
    # For rows that do not have an entry in the eliminated column, all this
    # update does is multiply the row in question by `pivot/last_pivot` (the
    # result of which is guaranteed to be integer by general properties of the
    # bareiss algorithm, even if `pivot/last_pivot` is not).
    #
    # Thus, when `pivot == last pivot`, we can skip the update for any rows that
    # do not have an entry in the eliminated column (because we'd simply be
    # multiplying by 1).
    #
    # As an additional MTK-specific enhancement, we further allow the case
    # when the absolute values are equal, i.e. effectively multiplying the row
    # by `-1`. To ensure this is legal, we need to show two things.
    # 1. The multiplication does not change the answer and
    # 2. The multiplication does not affect the fraction-freeness of the Bareiss
    #    algorithm.
    #
    # For point 1, remember that we're working on a system of linear equations,
    # so it is always legal for us to multiply any row by a sclar without changing
    # the underlying system of equations.
    #
    # For point 2, note that the factorization we're now computing is the same
    # as if we had multiplied the corresponding row (accounting for row swaps)
    # in the original matrix by `last_pivot/pivot`, ensuring that the matrix
    # itself has integral entries when `last_pivot/pivot` is integral (here we
    # have -1, which counts). We could use the more general integrality
    # condition, but that would in turn disturb the growth bounds on the
    # factorization matrix entries that the bareiss algorithm guarantees. To be
    # conservative, we leave it at this, as this captures the most important
    # case for MTK (where most pivots are `1` or `-1`).
    pivot_equal = pivot_equal_optimization && abs(pivot) == abs(last_pivot)

    for ei in k+1:size(M, 1)
        # elimate `v`
        coeff = 0
        ivars = eadj[ei]
        vj = findfirst(isequal(vpivot), ivars)
        if vj !== nothing
            coeff = old_cadj[ei][vj]
            deleteat!(old_cadj[ei], vj)
            deleteat!(eadj[ei], vj)
        elseif pivot_equal
            continue
        end

        # the pivot row
        kvars = eadj[k]
        kcoeffs = old_cadj[k]
        # the elimination target
        ivars = eadj[ei]
        icoeffs = old_cadj[ei]

        tmp_incidence = similar(eadj[ei], 0)
        tmp_coeffs = similar(old_cadj[ei], 0)
        # TODO: We know both ivars and kvars are sorted, we could just write
        # a quick iterator here that does this without allocation/faster.
        vars = sort(union(ivars, kvars))

        for v in vars
            v == vpivot && continue
            ck = getcoeff(kvars, kcoeffs, v)
            ci = getcoeff(ivars, icoeffs, v)
            ci = (pivot*ci - coeff*ck) ÷ last_pivot
            if ci !== 0
                push!(tmp_incidence, v)
                push!(tmp_coeffs, ci)
            end
        end

        eadj[ei] = tmp_incidence
        old_cadj[ei] = tmp_coeffs
    end
end

function bareiss_update_virtual_colswap_mtk!(zero!, M::AbstractMatrix, k, swapto, pivot, last_pivot; pivot_equal_optimization=true)
    if pivot_equal_optimization
        error("MTK pivot micro-optimization not implemented for `$(typeof(M))`.
            Turn off the optimization for debugging or use a different matrix type.")
    end
    bareiss_update_virtual_colswap!(zero!, M, k, swapto, pivot, last_pivot)
end

struct AsSubMatrix{T, Ti<:Integer} <: AbstractSparseMatrix{T, Ti}
    M::SparseMatrixCLIL{T, Ti}
end
Base.size(S::AsSubMatrix) = (S.M.nparentrows, S.M.ncols)

function Base.getindex(S::SparseMatrixCLIL{T}, i1::Integer, i2::Integer) where {T}
    checkbounds(S, i1, i2)

    col = S.row_cols[i1]
    nncol = searchsortedfirst(col, i2)
    (nncol > length(col) || col[nncol] != i2) && return zero(T)

    return S.row_vals[i1][nncol]
end

function Base.getindex(S::AsSubMatrix{T}, i1::Integer, i2::Integer) where {T}
    checkbounds(S, i1, i2)
    S = S.M

    nnrow = findfirst(==(i1), S.nzrows)
    isnothing(nnrow) && return zero(T)

    col = S.row_cols[nnrow]
    nncol = searchsortedfirst(col, i2)
    (nncol > length(col) || col[nncol] != i2) && return zero(T)

    return S.row_vals[nnrow][nncol]
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using LinearAlgebra, FFTW
import BSON, Statistics, PyPlot, Random

import NMRSpectraSimulator

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

# for loading something with Interpolations.jl
import OffsetArrays
import Interpolations

import MultistartOptimization, NLopt

#import Clustering

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

# run prep_full.jl first.

r = 1
T = Float64
N_starts = 500
local_optim_algorithm = NLopt.LN_BOBYQA
xtol_rel = 1e-9
maxeval = 200 # 2 # 50
maxtime = Inf
β_optim_algorithm = :GN_DIRECT_L
w_lb_default = 1e-1
w_ub_default = 100.0
β_max_iters = 500 # 2 # 500
β_xtol_rel = 1e-9
β_ftol_rel = 1e-9
β_maxtime = Inf

N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
#N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

shift_lb = -ones(N_d)
shift_ub = ones(N_d)

Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )

a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)
#

y_cost = y[cost_inds_set[r]]
U_cost = U_y[cost_inds_set[r]]
P_cost = P_y[cost_inds_set[r]]
LS_inds = 1:length(U_cost)

# prepare.
#N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
@assert length(shift_ub) == length(shift_lb) == N_d

U_rad_cost = U_cost .* (2*π)

# setup inner optim over β.
q0, updatedfunc, getshiftfunc, N_vars_set,
run_optim, obj_func_β, E_BLS, w_BLS, b_BLS, updateβfunc, updatewfunc,
q_β = NMRCalibrate.setupcostnesteddwarpw(Es, Bs, As, fs, SW, LS_inds, U_rad_cost,
    y_cost, Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings;
    β_optim_algorithm = β_optim_algorithm,
    w_lb_default = w_lb_default,
    w_ub_default = w_ub_default,
    β_max_iters = β_max_iters,
    β_xtol_rel = β_xtol_rel,
    β_ftol_rel = β_ftol_rel,
    β_maxtime = β_maxtime)

# set up outer optim over shifts.
#N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )
p_β = zeros(T, N_β) # persistant buffer.

obj_func = pp->NMRCalibrate.costnesteddw(U_rad_cost, y_cost, updatedfunc,
updatewfunc, pp, Es, Bs, κs_β_orderings, κs_β_DOFs, run_optim,
E_BLS, w_BLS, b_BLS, p_β)

#p = zeros(N_d)
p = [0.02;] # Serine 700 MHz BMRB.
#p = [0.13513301735885624;]
updatedfunc(p)

fill!(p_β, zero(T)) # always start from 0-phase?
minf, minx, ret, N_evals = run_optim(p_β)
p_β[:] = minx # take out?

#p_β[2] = p_β[1]

#NMRCalibrate.updateβ!(Bs, κs_β_orderings, κs_β_DOFs, p_β, 1)
NMRCalibrate.updateβ!(Bs, p_β, 1)
updatewfunc(1.0)

cost = norm(reinterpret(T, E_BLS)*w_BLS - b_BLS)^2

#w_test = copy(ws)
w_test = copy(w_BLS) #[5.5;]

#### visualize.
N_viz = 50000

U = LinRange(u_min, u_max, N_viz)
P = hz2ppmfunc.(U)
U_rad = U .* (2*π)

q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = w_test )
q_U = q2.(U_rad)

plotregion(P, U, q_U, P_y, y, P_cost, y_cost, 0.0, 1,
    "", "test", "";
    save_plot_flag = false,
    display_plot_flag = true,
    canvas_size = (1000,400))

# revisit cost again.

function evalcost(q2, U_cost_rad, y_cost)
    cost2 = 0.0
    for m = 1:length(U_cost_rad)
        cost2 += abs2(q2(U_cost_rad[m]) - y_cost[m])
    end

    return cost2
end

U_cost_rad = U_cost .* (2*π)
cost2 = evalcost(q2, U_cost_rad, y_cost)

#norm(b_BLS - reinterpret(T, y_cost))
@assert 12==3

q3, updateβfunc3, updatewfunc3, E_BLS3, w_BLS3, b_BLS3,
getβfunc3 = NMRCalibrate.setupcostβLSw(Es, Bs, As, LS_inds, U_rad_cost, y_cost,
    κs_β_DOFs, κs_β_orderings;
    w_lb_default = w_lb_default,
    w_ub_default = w_ub_default)
#
updatewfunc3(1.0)
E_BLS3*w_BLS3 - E_BLS*w_BLS
norm(E_BLS3- E_BLS)

(E_BLS3*w_BLS3)[1]
u_rad = U_cost[1]*2*pi
println("q0(u_rad) = ", q0(u_rad))
println("q3(u_rad) = ", q3(u_rad))

ext: jl


text:
module PatchTests

using Fixtures
using FactCheck

# This module has an example of some kinds of thing we can patch and is used by the tests
module TestModuleOne
  export test1variable, test1constmethod, test1nonconstmethod, test1function, test1lambda
  test1variable = 7
  test1constmethod(x) = x + test1variable
  test1constmethod(x, y) = x + y + test1variable
  test1nonconstmethod = test1constmethod
  test1function = function(x)
    x + test1variable
  end
  test1lambda = x -> x + test1variable
end

# This module has an example of most kinds of thing we can patch plus it uses TestModuleOne
# it is used by the tests
module TestModuleTwo
  using PatchTests.TestModuleOne
  export test2variable, test2constmethod, test2nonconstmethod, test2function, test2lambda
  test2variable = 7
  test2constmethod(x) = x + test2variable
  test2constmethod(x, y) = x + y + test2variable
  test2nonconstmethod = test2constmethod
  test2function = function(x)
    x + test2variable
  end
  test2lambda = x -> x + test2variable

  module nestedmodule
    value = 12
    func() = 200
  end

  type List
    value::Any
    next::List
    List(v::Any) = new(v)
    List(v::Any, nxt::List) = new(v, nxt)
    List{T}(vs::Vector{T}) = length(vs) > 1 ? List(vs[1], List(vs[2:end])) : List(vs[1])
  end

  lst = List([40,30,20,10])
end

# A global (outside of a module) we can use for testing.
mylist = TestModuleTwo.List(1, TestModuleTwo.List(2, TestModuleTwo.List(3, TestModuleTwo.List(4))))

# This function checks that TestModuleOne and TestModuleOne work as defined above, we run it
# before and after each test context to show that the patching has worked
function checkmodules1and2()
  @fact TestModuleOne.test1variable --> 7
  @fact TestModuleOne.test1constmethod(1) --> 8
  @fact TestModuleOne.test1constmethod(1,10) --> 18
  @fact TestModuleOne.test1nonconstmethod(1) --> 8
  @fact TestModuleOne.test1nonconstmethod(1,10) --> 18
  @fact TestModuleOne.test1function(1) --> 8
  @fact TestModuleOne.test1lambda(1) --> 8

  @fact TestModuleTwo.test1variable --> 7
  @fact TestModuleTwo.test1constmethod(1) --> 8
  @fact TestModuleTwo.test1constmethod(1,10) --> 18
  @fact TestModuleTwo.test1nonconstmethod(1) --> 8
  @fact TestModuleTwo.test1nonconstmethod(1,10) --> 18
  @fact TestModuleTwo.test1function(1) --> 8
  @fact TestModuleTwo.test1lambda(1) --> 8

  @fact TestModuleTwo.test2variable --> 7
  @fact TestModuleTwo.test2constmethod(1) --> 8
  @fact TestModuleTwo.test2constmethod(1,10) --> 18
  @fact TestModuleTwo.test2nonconstmethod(1) --> 8
  @fact TestModuleTwo.test2nonconstmethod(1,10) --> 18
  @fact TestModuleTwo.test2function(1) --> 8
  @fact TestModuleTwo.test2lambda(1) --> 8
  @fact TestModuleTwo.nestedmodule.value --> 12
  @fact TestModuleTwo.nestedmodule.func() --> 200

  @fact TestModuleTwo.lst.value --> 40
  @fact TestModuleTwo.lst.next.value --> 30
  @fact TestModuleTwo.lst.next.next.value --> 20
  @fact TestModuleTwo.lst.next.next.next.value --> 10

  @fact typeof(TestModuleTwo.List(100)) --> TestModuleTwo.List
  @fact typeof(TestModuleTwo.List(100, TestModuleTwo.List(200))) --> TestModuleTwo.List
  @fact typeof(TestModuleTwo.List([1, 2, 3])) --> TestModuleTwo.List
  @fact_throws TestModuleTwo.List(1,2,3)
end

# Checks that mylist is still as it was defined to be above, we run it before and after tests
# to show that the patching is working
function checkmylist()
  @fact mylist.value --> 1
  @fact mylist.next.value --> 2
  @fact mylist.next.next.value --> 3
  @fact mylist.next.next.next.value --> 4
end

# Our Patch test suite
facts("Patch tests") do

  context("patch test2variable with another value") do
    patch(TestModuleTwo, :test2variable, 100) do
      @fact TestModuleTwo.test2variable --> 100
      @fact TestModuleTwo.test2constmethod(1) --> 101
      @fact TestModuleTwo.test2constmethod(1,10) --> 111
      @fact TestModuleTwo.test2nonconstmethod(1) --> 101
      @fact TestModuleTwo.test2nonconstmethod(1,10) --> 111
      @fact TestModuleTwo.test2function(1) --> 101
      @fact TestModuleTwo.test2lambda(1) --> 101
    end
  end

  context("patch test1variable with another value") do
    patch(TestModuleOne, :test1variable, 100) do
      @fact TestModuleTwo.test1variable --> 100
      @fact TestModuleTwo.test1constmethod(1) --> 101
      @fact TestModuleTwo.test1constmethod(1,10) --> 111
      @fact TestModuleTwo.test1nonconstmethod(1) --> 101
      @fact TestModuleTwo.test1nonconstmethod(1,10) --> 111
      @fact TestModuleTwo.test1function(1) --> 101
      @fact TestModuleTwo.test1lambda(1) --> 101
    end
  end

  # This test shows the importance of patching globals in the module in which
  # they are defined. At least we get a warning.
  context("patch test1variable with another value in wrong module") do
    patch(TestModuleTwo, :test1variable, 100) do
      @fact TestModuleTwo.test1variable --> 100 # This is all we have acheived
      @fact TestModuleTwo.test1constmethod(1) --> 8
      @fact TestModuleTwo.test1constmethod(1,10) --> 18
      @fact TestModuleTwo.test1nonconstmethod(1) --> 8
      @fact TestModuleTwo.test1nonconstmethod(1,10) --> 18
      @fact TestModuleTwo.test1function(1) --> 8
      @fact TestModuleTwo.test1lambda(1) --> 8
    end
  end

  context("patch test2lambda with a lambda") do
    patch(TestModuleTwo, :test2lambda, x->200) do
      @fact TestModuleTwo.test2lambda(1) --> 200
      @fact TestModuleTwo.test2lambda(2) --> 200
      @fact TestModuleTwo.test2lambda(3) --> 200
    end
  end

  context("patch test1lambda with a lambda") do
    patch(TestModuleOne, :test1lambda, x->200) do
      @fact TestModuleTwo.test1lambda(1) --> 200
      @fact TestModuleTwo.test1lambda(2) --> 200
      @fact TestModuleTwo.test1lambda(3) --> 200
    end
  end

  context("patch test2function with a lambda") do
    patch(TestModuleTwo, :test2function, x->200) do
      @fact TestModuleTwo.test2function(1) --> 200
      @fact TestModuleTwo.test2function(2) --> 200
      @fact TestModuleTwo.test2function(3) --> 200
    end
  end

  context("patch test1function with a lambda") do
    patch(TestModuleOne, :test1function, x->200) do
      @fact TestModuleTwo.test1function(1) --> 200
      @fact TestModuleTwo.test1function(2) --> 200
      @fact TestModuleTwo.test1function(3) --> 200
    end
  end

  context("patch test2nonconstmethod with a lambda") do
    patch(TestModuleTwo, :test2nonconstmethod, x->200) do
      @fact TestModuleTwo.test2nonconstmethod(1) --> 200
      @fact TestModuleTwo.test2nonconstmethod(2) --> 200
      @fact TestModuleTwo.test2nonconstmethod(3) --> 200
      @fact_throws TestModuleTwo.test2nonconstmethod(1,10)
    end
  end

  context("patch test1nonconstmethod with a lambda") do
    patch(TestModuleOne, :test1nonconstmethod, x->200) do
      @fact TestModuleTwo.test1nonconstmethod(1) --> 200
      @fact TestModuleTwo.test1nonconstmethod(2) --> 200
      @fact TestModuleTwo.test1nonconstmethod(3) --> 200
      @fact_throws TestModuleTwo.test1nonconstmethod(1,10)
    end
  end

  context("patch test2constmethod with a lambda") do
    patch(TestModuleTwo, :test2constmethod, x->200) do
      @fact TestModuleTwo.test2constmethod(1) --> 200
      @fact TestModuleTwo.test2constmethod(2) --> 200
      @fact TestModuleTwo.test2constmethod(3) --> 200
      @fact_throws TestModuleTwo.test2constmethod(1,10)
    end
  end

  context("patch test1constmethod with a lambda") do
    patch(TestModuleOne, :test1constmethod, x->200) do
      @fact TestModuleTwo.test1constmethod(1) --> 200
      @fact TestModuleTwo.test1constmethod(2) --> 200
      @fact TestModuleTwo.test1constmethod(3) --> 200
      @fact_throws TestModuleTwo.test1constmethod(1,10)
    end
  end

  # a method we will patch in over other methods for testing
  testmethod(x) = 200
  testmethod(x,y) = 400

  context("patch test2lambda with a method") do
    patch(TestModuleTwo, :test2lambda, testmethod) do
      @fact TestModuleTwo.test2lambda(1) --> 200
      @fact TestModuleTwo.test2lambda(2) --> 200
      @fact TestModuleTwo.test2lambda(1,10) --> 400
      @fact TestModuleTwo.test2lambda(2,10) --> 400
    end
  end

  context("patch test1lambda with a method") do
    patch(TestModuleTwo, :test1lambda, testmethod) do
      @fact TestModuleTwo.test1lambda(1) --> 200
      @fact TestModuleTwo.test1lambda(2) --> 200
      @fact TestModuleTwo.test1lambda(1,10) --> 400
      @fact TestModuleTwo.test1lambda(2,10) --> 400
    end
  end

  context("patch test2function with a method") do
    patch(TestModuleTwo, :test2function, testmethod) do
      @fact TestModuleTwo.test2function(1) --> 200
      @fact TestModuleTwo.test2function(2) --> 200
      @fact TestModuleTwo.test2function(1,10) --> 400
      @fact TestModuleTwo.test2function(2,10) --> 400
    end
  end

  context("patch test2nonconstmethod with a method") do
    patch(TestModuleTwo, :test2nonconstmethod, testmethod) do
      @fact TestModuleTwo.test2nonconstmethod(1) --> 200
      @fact TestModuleTwo.test2nonconstmethod(2) --> 200
      @fact TestModuleTwo.test2nonconstmethod(1,10) --> 400
      @fact TestModuleTwo.test2nonconstmethod(2,10) --> 400
    end
  end

  context("patch test1nonconstmethod with a method") do
    patch(TestModuleTwo, :test1nonconstmethod, testmethod) do
      @fact TestModuleTwo.test1nonconstmethod(1) --> 200
      @fact TestModuleTwo.test1nonconstmethod(2) --> 200
      @fact TestModuleTwo.test1nonconstmethod(1,10) --> 400
      @fact TestModuleTwo.test1nonconstmethod(2,10) --> 400
    end
  end

  context("patch test2constmethod with a method") do
    patch(TestModuleTwo, :test2constmethod, testmethod) do
      @fact TestModuleTwo.test2constmethod(1) --> 200
      @fact TestModuleTwo.test2constmethod(2) --> 200
      @fact TestModuleTwo.test2constmethod(1,10) --> 400
      @fact TestModuleTwo.test2constmethod(2,10) --> 400
    end
  end

  context("patch test1constmethod with a method") do
    patch(TestModuleOne, :test1constmethod, testmethod) do
      @fact TestModuleTwo.test1constmethod(1) --> 200
      @fact TestModuleTwo.test1constmethod(2) --> 200
      @fact TestModuleTwo.test1constmethod(1,10) --> 400
      @fact TestModuleTwo.test1constmethod(2,10) --> 400
    end
  end

  context("patch test2lambda with a value") do
    patch(TestModuleTwo, :test2lambda, 200) do
      @fact TestModuleTwo.test2lambda --> 200
    end
  end

  context("patch test2function with a value") do
    patch(TestModuleTwo, :test2function, 200) do
      @fact TestModuleTwo.test2function --> 200
    end
  end

  context("patch test2nonconstmethod with a value") do
    patch(TestModuleTwo, :test2nonconstmethod, 200) do
      @fact TestModuleTwo.test2nonconstmethod --> 200
    end
  end

  context("patching test2constmethod with a value throws an error") do
    @fact_throws patch(TestModuleTwo, :test2constmethod, 200) do
    end
  end

  context("Ensure that the return value of the wrapped fn is always returned") do
    @fact patch(()->123, TestModuleTwo, :test2variable, 100) --> 123
    @fact patch(()->456, TestModuleTwo, :test2constmethod, testmethod) --> 456
    @fact patch(()->789, TestModuleTwo, :test2nonconstmethod, testmethod) --> 789
    @fact patch(()->123, TestModuleTwo, :test2function, testmethod) --> 123
    @fact patch(()->456, TestModuleTwo, :test2lambda, testmethod) --> 456

    @fact patch(()->789, TestModuleTwo, :test2nonconstmethod, ()->200) --> 789
    @fact patch(()->123, TestModuleTwo, :test2function, ()->200) --> 123
    @fact patch(()->456, TestModuleTwo, :test2lambda, ()->200) --> 456

    @fact patch(()->789, TestModuleTwo, :test2nonconstmethod, 200) --> 789
    @fact patch(()->123, TestModuleTwo, :test2function, 200) --> 123
    @fact patch(()->456, TestModuleTwo, :test2lambda, 200) --> 456
  end

  context("patch a value in a nested module") do
    patch(TestModuleTwo.nestedmodule, :value, 4000) do
      @fact TestModuleTwo.nestedmodule.value --> 4000
    end
  end

  context("patch a method in a nested module") do
    patch(TestModuleTwo.nestedmodule, :func, ()->1234) do
      @fact TestModuleTwo.nestedmodule.func() --> 1234
    end
  end

  context("patch first value") do
    patch(TestModuleTwo.lst, :value, 100) do
      @fact TestModuleTwo.lst.value --> 100
      @fact TestModuleTwo.lst.next.value --> 30
      @fact TestModuleTwo.lst.next.next.value --> 20
      @fact TestModuleTwo.lst.next.next.next.value --> 10
    end
  end

  context("patch second value") do
    patch(TestModuleTwo.lst.next, :value, 100) do
      @fact TestModuleTwo.lst.value --> 40
      @fact TestModuleTwo.lst.next.value --> 100
      @fact TestModuleTwo.lst.next.next.value --> 20
      @fact TestModuleTwo.lst.next.next.next.value --> 10
    end
  end

  context("patch third value") do
    patch(TestModuleTwo.lst.next.next, :value, 100) do
      @fact TestModuleTwo.lst.value --> 40
      @fact TestModuleTwo.lst.next.value --> 30
      @fact TestModuleTwo.lst.next.next.value --> 100
      @fact TestModuleTwo.lst.next.next.next.value --> 10
    end
  end

  context("patch fourth value") do
    patch(TestModuleTwo.lst.next.next.next, :value, 100) do
      @fact TestModuleTwo.lst.value --> 40
      @fact TestModuleTwo.lst.next.value --> 30
      @fact TestModuleTwo.lst.next.next.value --> 20
      @fact TestModuleTwo.lst.next.next.next.value --> 100
    end
  end

#  context("patch constructor with a lambda") do
#    patch(TestModuleTwo, :List, (x,y,z) -> 100) do
#      @fact_throws TestModuleTwo.List(100)
#      @fact_throws TestModuleTwo.List(100, TestModuleTwo.List(200))
#      @fact_throws TestModuleTwo.List([1, 2, 3])
#      @fact TestModuleTwo.List(1,2,3) --> 100
#      @fact typeof(TestModuleTwo.List(1,2,3)) --> not(TestModuleTwo.List)
#    end
#  end
#
#  # A method to patch in as the List constructor
#  mynewconstructor(x,y,z) = 100
#
#  context("patch constructor with a method") do
#    patch(TestModuleTwo, :List, mynewconstructor) do
#      @fact_throws TestModuleTwo.List(100)
#      @fact_throws TestModuleTwo.List(100, TestModuleTwo.List(200))
#      @fact_throws TestModuleTwo.List([1, 2, 3])
#      @fact TestModuleTwo.List(1,2,3) --> 100
#    end
#  end

end

# Our Patch test suite
facts("Patch tests") do

  context("patch first value of mylist") do
    patch(mylist, :value, 100) do
      @fact mylist.value --> 100
      @fact mylist.next.value --> 2
      @fact mylist.next.next.value --> 3
      @fact mylist.next.next.next.value --> 4
    end
  end

  context("patch second value of mylist") do
    patch(mylist.next, :value, 100) do
      @fact mylist.value --> 1
      @fact mylist.next.value --> 100
      @fact mylist.next.next.value --> 3
      @fact mylist.next.next.next.value --> 4
    end
  end

  context("patch third value of mylist") do
    patch(mylist.next.next, :value, 100) do
      @fact mylist.value --> 1
      @fact mylist.next.value --> 2
      @fact mylist.next.next.value --> 100
      @fact mylist.next.next.next.value --> 4
    end
  end

  context("patch fourth value of mylist") do
    patch(mylist.next.next.next, :value, 100) do
      @fact mylist.value --> 1
      @fact mylist.next.value --> 2
      @fact mylist.next.next.value --> 3
      @fact mylist.next.next.next.value --> 100
    end
  end

  # A different list to use in tests
  const myotherlist = TestModuleTwo.List([4,3,2,1])

  context("patch global variable") do
    @patch(mylist, myotherlist) do
      @fact mylist.value --> 4
      @fact mylist.next.value --> 3
      @fact mylist.next.next.value --> 2
      @fact mylist.next.next.next.value --> 1
    end
  end

  context("patch global variable and check the return value") do
    const result = @patch(mylist, myotherlist) do
      mylist
    end
    @fact result --> myotherlist
  end
end

module BrokenExample
  inner_method(x) = 12345
  outer_method(x) = inner_method(x)
end

facts("Broken stuff") do
  context("patch is nobbled by method caching/inlining GH265") do
    @fact BrokenExample.outer_method(100) --> 12345

    # because outer_method(int64) has already been compiled
    # we can't patch it to use our implementation of inner_method
    # See https://github.com/JuliaLang/julia/issues/265
    patch(BrokenExample, :inner_method, x->54321) do
      @fact BrokenExample.outer_method(100) --> not(54321) # :(
    end
    @fact BrokenExample.outer_method(100) --> 12345

    # But if we change type then we are able to patch it
    patch(BrokenExample, :inner_method, x->987654) do
      @fact BrokenExample.outer_method(100.0) --> 987654
    end

    # and it is all ok afterwards
    @fact BrokenExample.outer_method(100) --> 12345
    @fact BrokenExample.outer_method(100.0) --> 12345
  end
end

facts("Multiple patches") do
  context("Trying to do multiple pathces") do

    testmethod(x) = 200
    testmethod(x,y) = 400

    patchers = [
        Patcher(TestModuleTwo, :test2variable, 100),
        Patcher(TestModuleOne, :test1variable, 100),
        Patcher(TestModuleTwo, :test1variable, 100),
        Patcher(TestModuleOne, :test1lambda, x->200),
        Patcher(TestModuleTwo, :test2function, x->200),
        Patcher(TestModuleOne, :test1function, x->200),
        Patcher(TestModuleOne, :test1nonconstmethod, x->200),
        Patcher(TestModuleTwo, :test2constmethod, x->200),
        Patcher(TestModuleTwo, :test2lambda, testmethod),
        Patcher(TestModuleTwo, :test1lambda, testmethod),
        Patcher(TestModuleTwo, :test1nonconstmethod, testmethod),
        Patcher(TestModuleOne, :test1constmethod, testmethod),
        Patcher(TestModuleTwo, :test2nonconstmethod, 200),
        Patcher(TestModuleTwo.nestedmodule, :value, 4000),
        Patcher(TestModuleTwo.nestedmodule, :func, ()->1234),
        Patcher(TestModuleTwo.lst, :value, 101),
        Patcher(TestModuleTwo.lst.next, :value, 102),
        Patcher(TestModuleTwo.lst.next.next, :value, 103),
        Patcher(TestModuleTwo.lst.next.next.next, :value, 104),
        Patcher(mylist, :value, 101),
        Patcher(mylist.next, :value, 102),
        Patcher(mylist.next.next, :value, 103),
        Patcher(mylist.next.next.next, :value, 104),
      ]

    @fact TestModuleOne.test1variable --> 7
    @fact TestModuleOne.test1constmethod(1) --> 8
    @fact TestModuleOne.test1constmethod(1,10) --> 18
    @fact TestModuleOne.test1nonconstmethod(1) --> 8
    @fact TestModuleOne.test1nonconstmethod(1,10) --> 18
    @fact TestModuleOne.test1function(1) --> 8
    @fact TestModuleOne.test1lambda(1) --> 8
    @fact TestModuleTwo.test2variable --> 7
    @fact TestModuleTwo.test2constmethod(1) --> 8
    @fact TestModuleTwo.test2constmethod(1,10) --> 18
    @fact TestModuleTwo.test2nonconstmethod(1) --> 8
    @fact TestModuleTwo.test2nonconstmethod(1,10) --> 18
    @fact TestModuleTwo.test2function(1) --> 8
    @fact TestModuleTwo.test2lambda(1) --> 8
    @fact TestModuleTwo.nestedmodule.value --> 12
    @fact TestModuleTwo.nestedmodule.func() --> 200
    @fact TestModuleTwo.lst.value --> 40
    @fact TestModuleTwo.lst.next.value --> 30
    @fact TestModuleTwo.lst.next.next.value --> 20
    @fact TestModuleTwo.lst.next.next.next.value --> 10
    @fact typeof(TestModuleTwo.List(100)) --> TestModuleTwo.List
    @fact typeof(TestModuleTwo.List(100, TestModuleTwo.List(200))) --> TestModuleTwo.List
    @fact typeof(TestModuleTwo.List([1, 2, 3])) --> TestModuleTwo.List
    @fact_throws TestModuleTwo.List(1,2,3)
    @fact mylist.value --> 1
    @fact mylist.next.value --> 2
    @fact mylist.next.next.value --> 3
    @fact mylist.next.next.next.value --> 4

    random_number() = rand(1:1000)

    patch(patchers) do
      @fact TestModuleTwo.test2variable --> 100
      @fact TestModuleTwo.test2function(1) --> 200
      @fact TestModuleTwo.test1function(1) --> 200
      @fact TestModuleTwo.test1variable --> 100 # This is all we have acheived
      @fact TestModuleTwo.test2constmethod(random_number()) --> 200
      @fact_throws TestModuleTwo.test2constmethod(random_number(),random_number())
      @fact TestModuleTwo.test2lambda(random_number()) --> 200
      @fact TestModuleTwo.test2lambda(random_number(),random_number()) --> 400
      @fact TestModuleTwo.test1lambda(random_number()) --> 200
      @fact TestModuleTwo.test1lambda(random_number(),random_number()) --> 400
      @fact TestModuleTwo.test1nonconstmethod(random_number()) --> 200
      @fact TestModuleTwo.test1nonconstmethod(random_number(),random_number()) --> 400
      @fact TestModuleTwo.test1constmethod(random_number()) --> 200
      @fact TestModuleTwo.test1constmethod(random_number(),random_number()) --> 400
      @fact TestModuleTwo.test2nonconstmethod --> 200
      @fact TestModuleTwo.nestedmodule.value --> 4000
      @fact TestModuleTwo.nestedmodule.func() --> 1234
      @fact TestModuleTwo.lst.value --> 101
      @fact TestModuleTwo.lst.next.value --> 102
      @fact TestModuleTwo.lst.next.next.value --> 103
      @fact TestModuleTwo.lst.next.next.next.value --> 104
      @fact mylist.value --> 101
      @fact mylist.next.value --> 102
      @fact mylist.next.next.value --> 103
      @fact mylist.next.next.next.value --> 104
    end

    @fact TestModuleOne.test1variable --> 7
    @fact TestModuleOne.test1constmethod(1) --> 8
    @fact TestModuleOne.test1constmethod(1,10) --> 18
    @fact TestModuleOne.test1nonconstmethod(1) --> 8
    @fact TestModuleOne.test1nonconstmethod(1,10) --> 18
    @fact TestModuleOne.test1function(1) --> 8
    @fact TestModuleOne.test1lambda(1) --> 8
    @fact TestModuleTwo.test2variable --> 7
    @fact TestModuleTwo.test2constmethod(1) --> 8
    @fact TestModuleTwo.test2constmethod(1,10) --> 18
    @fact TestModuleTwo.test2nonconstmethod(1) --> 8
    @fact TestModuleTwo.test2nonconstmethod(1,10) --> 18
    @fact TestModuleTwo.test2function(1) --> 8
    @fact TestModuleTwo.test2lambda(1) --> 8
    @fact TestModuleTwo.nestedmodule.value --> 12
    @fact TestModuleTwo.nestedmodule.func() --> 200
    @fact TestModuleTwo.lst.value --> 40
    @fact TestModuleTwo.lst.next.value --> 30
    @fact TestModuleTwo.lst.next.next.value --> 20
    @fact TestModuleTwo.lst.next.next.next.value --> 10
    @fact typeof(TestModuleTwo.List(100)) --> TestModuleTwo.List
    @fact typeof(TestModuleTwo.List(100, TestModuleTwo.List(200))) --> TestModuleTwo.List
    @fact typeof(TestModuleTwo.List([1, 2, 3])) --> TestModuleTwo.List
    @fact_throws TestModuleTwo.List(1,2,3)
    @fact mylist.value --> 1
    @fact mylist.next.value --> 2
    @fact mylist.next.next.value --> 3
    @fact mylist.next.next.next.value --> 4

  end
end

end


text:
# Base.eltype{T}(::Type{Polynomial{T}}) = T
# Base.eltype{T}(::Type{TransferFunction{T}}) = T
# Base.eltype{T}(::Type{SisoZpk{T}}) = T
# Base.eltype{T}(::Type{SisoRational{T}}) = T
# Base.eltype{T}(::Type{StateSpace{T}}) = T
#

# """
#     responsetype(T)
# Return the numerical type of the time-response for composite LTISystem type `T`, e.g., `responsetype(ss(1.)) == Float64`
# If `T` is not a type, return `responsetype(typeof(T))`
# """
# responsetype{T<:Real}(::Type{T}) = T
# responsetype{T<:Complex}(::Type{T}) = promote_type(T.types...)
# responsetype(x) = typeof(x) # Fallback method
# responsetype{T}(::Type{T}) = T # Catch all

#matrix_type{T}(::Type{StateSpace{T, M{T}}}) = T

# Abstract type pyramid =============================================================

# NOTE: a lot that doesnt seem to be used...
#Base.promote_rule(::Type{StateSpace{T1}}, ::Type{StateSpace{T2}}) where {T1,T2} = StateSpace{promote_type(T1, T2)}

# NOTE: Is the below thing correct always?
function Base.promote_rule(::Type{StateSpace{T1,MT1}}, ::Type{StateSpace{T2,MT2}}) where {T1,T2,MT1,MT2}
    MT = promote_type(MT1, MT2)
    if isconcretetype(MT) # This sometimes fails and produces a Unionall. This condition can be checked at compile time, checked by f(x,y) = (MT = promote_type(x, y); isconcretetype(MT)); @code_llvm f(Int,Float64)
        StateSpace{promote_type(T1, T2),MT}
    else # If fail, fall back to an ordinary matrix
        FT = promote_type(T1, T2)
        StateSpace{FT,Matrix{FT}}
    end
end

Base.promote_rule(::Type{StateSpace{T,MT}}, ::Type{HeteroStateSpace{AT,BT,CT,DT}}) where {T,MT,AT,BT,CT,DT} = HeteroStateSpace{promote_type(MT,AT),promote_type(MT,BT),promote_type(MT,CT),promote_type(MT,DT)}

function Base.promote_rule(::Type{TransferFunction{S1}}, ::Type{StateSpace{T2,MT}}) where {T1,S1<:SisoTf{T1},T2,MT}
    #promote_type(Base.promote_op(/, T1, T1), T2)
    StateSpace{promote_type(T1,T2), promote_type(Matrix{T1},MT)}
end
# NOTE: Perhaps should try to keep matrix structure?

function Base.promote_rule(::Type{TransferFunction{S1}}, ::Type{DelayLtiSystem{T2}}) where {T1,S1<:SisoTf{T1},T2}
    DelayLtiSystem{promote_type(T1,T2)}
end
function Base.promote_rule(::Type{StateSpace{T1,MT}}, ::Type{DelayLtiSystem{T2}}) where {T1,MT,T2}
    DelayLtiSystem{promote_type(T1,T2)}
end

Base.promote_rule(::Type{TransferFunction{S1}}, ::Type{TransferFunction{S2}}) where {S1, S2} = TransferFunction{promote_type(S1, S2)}
#Base.promote_rule(::Type{SisoTf}, ::Type{TransferFunction}) = TransferFunction
#Base.promote_rule(::Type{SisoZpk}, ::Type{TransferFunction}) = TransferFunction
#Base.promote_rule(::Type{SisoRational}, ::Type{TransferFunction}) = TransferFunction

#function Base.promote_rule{T<:StateSpace,P<:TransferFunction}(::Type{T}, ::Type{P})  #where T <: StateSpace where P <: TransferFunction
#    S = promote_type(primitivereal(P), primitivereal(T)) # TODO: this is not correct for P <: zpk, in which case the StateSpace system gets complex matrices
#    StateSpace{Matrix{S}}
#end

# Promotion of Number ==========================================================

#Base.promote_rule{T1<:Number,T2<:Number}(::Type{TransferFunction{S}}, ::T2) = promote_type(T1, T2)

Base.promote_rule(::Type{TransferFunction{S}}, ::Type{T2}) where {S<:SisoTf, T2<:Number} =
    TransferFunction{promote_type(S, T2)}

# TODO Figure out a better way
function Base.promote_rule(::Type{StateSpace{T1, MT}}, ::Type{T2}) where {T1, MT, T2<:Number}
    NewMT = Base.promote_op(*, MT, T2)
    StateSpace{eltype(NewMT), NewMT}
end

Base.promote_rule(::Type{TransferFunction{SisoZpk{T1,TR1}}}, ::Type{M2}) where {T1, TR1, T2, M2<:AbstractMatrix{T2}} = TransferFunction{SisoZpk{T1, promote_type(TR1, T2)}}

Base.promote_rule(::Type{TransferFunction{SisoRational{T1}}}, ::Type{M2}) where {T1, T2, M2<:AbstractMatrix{T2}} = TransferFunction{SisoRational{promote_type(T1, T2)}}

function Base.promote_rule(::Type{StateSpace{T1, MT1}}, ::Type{MT2}) where {T1, MT1, MT2<:AbstractMatrix}
    MT = promote_type(MT1, MT2)
    StateSpace{eltype(MT), MT}
end

Base.promote_rule(::Type{DelayLtiSystem{T1,S}}, ::Type{MT1}) where {T1, S, MT1<:AbstractMatrix} =
    DelayLtiSystem{promote_type(T1, eltype(MT1)),S}

#Base.promote_rule{S<:TransferFunction{<:SisoTf}}(::Type{S}, ::Type{<:Real}) = S

# We want this, but not possible, so hardcode for SisoTypes
#Base.promote_rule(::Type{S}, ::Type{T2}) where {T1, S<:SisoTf{T1}, T2<:Number} = S{promote_type(T1, T2)}

function Base.promote_rule(::Type{SisoZpk{T1,C2}}, ::Type{T2}) where {T1, C2, T2<:Number}
    GainType = promote_type(T1, T2)
    return SisoZpk{GainType, complex(GainType)}
end
Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1, T2<:Number} = SisoRational{promote_type(T1,T2)}
#Base.promote_rule(::Type{StateSpace{T1, MT}}, ::Type{T2}) where {T1, T2<:Number} = SisoRational{promote_type(T1,T2)}

#Base.promote_rule{T<:SisoZpk}(::Type{T}, ::Type{<:Real}) = T
#Base.promote_rule{T}(::Type{SisoRational{T}}, ::Type{<:Real}) = SisoRational{T}
#Base.promote_rule{T}(::Type{StateSpace{T}}, ::Type{<:Real}) = StateSpace{T}
#Base.promote_rule{S<:TransferFunction{<:SisoTf},T<:Real}(::Type{S}, ::Union{Type{Array{T,2}},Type{Array{T,1}}}) = S

# Less abstract promotions
#Base.promote_rule(::Type{TransferFunction{SisoRational}}, ::Type{TransferFunction{SisoZpk}}) = TransferFunction{SisoZpk} # NOTE: Is this what we want?

# function Base.promote_rule{T1c<:SisoZpk,T2c<:SisoRational}(::Type{T1c}, ::Type{T2c})
#     upper_type = SisoZpk
#     inner_type = promote_type(arraytype(T1c), arraytype(T2c))
#     upper_type{inner_type}
# end

# @show siso_type = promote_type(eltype.(systems)...)
# if isleaftype(siso_type)
#     siso_type = siso_type.name.wrapper
# end
# # array_type = promote_type(arraytype.(systems)...)
# # mat = hcat([(e->convert(siso_type{array_type}, e)).(s.matrix) for s in systems]...)
# @show promote_type([s.matrix for s in systems]...)

#

#
#
# # Promote_op types # NOTE: Needed?
# Base.promote_op{T<:SisoTf}(::Any, ::Type{T}, ::Type{T}) = T
#
# #Just default SisoTf to SisoRational
# SisoTf(args...) = SisoRational(args...)
#
# Base.zero(::Type{<:SisoTf}) = zero(SisoRational)
# Base.zero(::SisoTf) = zero(SisoRational)
# Base.zero(::Type{<:SisoZpk}) = SisoZpk(Float64[],Float64[],0.0)
# Base.zero(::SisoZpk) = Base.zero(SisoZpk)
# Base.zero{T}(::Type{SisoRational{T}}) = SisoRational(zero(Polynomial{T}), one(Polynomial{T})) # FIXME: Is the in analogy with how zero for SisoRational is createdzer?
# Base.zero{T}(::SisoRational{T}) = Base.zero(SisoRational{T})


text:
# `shiftby` is equal to the number of bits required to represent index information
# for one level of the BitmappedTrie.
#
# Here, `shiftby` is 5, which means that the BitmappedTrie Arrays will be length 32.
const shiftby = 5
const trielen = 2^shiftby

abstract BitmappedTrie{T}

# Copy elements from one Array to another, up to `n` elements.
#
function copy_to{T}(from::Array{T}, to::Array{T}, n::Int)
    for i=1:n
        to[i] = from[i]
    end
    to
end

# Copies elements from one Array to another of size `len`.
#
copy_to_len{T}(from::Array{T}, len::Int) =
    copy_to(from, Array(T, len), min(len, length(from)))

mask(t::BitmappedTrie, i::Int) = (((i - 1) >>> shift(t)) & (trielen - 1)) + 1

Base.endof(t::BitmappedTrie) = length(t)

Base.length(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement Base.length")
shift(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.shift")
maxlength(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.maxlength")
arrayof(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.arrayof")

function Base.isequal(t1::BitmappedTrie, t2::BitmappedTrie)
    length(t1)    == length(t2)    &&
    shift(t1)     == shift(t2)     &&
    maxlength(t1) == maxlength(t2) &&
    arrayof(t1)   == arrayof(t2)
end

==(t1::BitmappedTrie, t2::BitmappedTrie) = isequal(t1, t2)


text:
using Test, PolynomialBases
import SymPy

x, α = SymPy.symbols("x, alpha")

@test 0 == SymPy.simplify( gegenbauer(x, 0, α) - 1 )
@test 0 == SymPy.simplify( gegenbauer(x, 1, α) - ( 2α*x ) )
@test 0 == SymPy.simplify( gegenbauer(x, 2, α) - ( -α + 2α*(1+α)*x^2 ) )
@test 0 == SymPy.simplify( gegenbauer(x, 3, α) - ( -2α*(1+α)*x + 4*α*(1+α)*(2+α)*x^3/3 ) )

@inferred gegenbauer(x, 4, α)
@inferred gegenbauer(10., 4, 3)


text:
module feastUtil
    using LinearAlgebra
    using Random
    export linquad,trapezoidal,getMaxResInside,Quadrrsolve,beyn,nseig,getInsideIndex
    include("util.jl")
end

module feastCore
    using LinearAlgebra
    using Random
    using ..feastUtil
    export feast_core,feastNS_core
    include("feast_core.jl")
end

module feastLinear
    using SparseArrays
    using LinearAlgebra
    using Random
    using ..feastCore,..feastUtil,IterativeSolvers
    export feast_linear,ifeast_linear,feastNS_linear,ifeastNS_linear
    include("feast_linear.jl")
end

module feastNonlinear
    using SparseArrays
    using LinearAlgebra
    using Random
    using IterativeSolvers
    using ..feastCore,..feastUtil
    export nlfeast_quad,nlfeast_beyn,inlfeast_beyn
    include("feast_nonlinear.jl")
end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
using HDF5, JLD, Test

# names is deprecated to keys v0.12 JLD
if !isdefined(JLD, :keys)
   keys = names
end

# Define variables of different types
x = 3.7
A = reshape(collect(1:15), 3, 5)
str = "Hello"
stringsA = String["It", "was", "a", "dark", "and", "stormy", "night"]
stringsU = String["It", "was", "a", "dark", "and", "stormy", "night"]
empty_string = ""
empty_string_array = String[]
empty_array_of_strings = String[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = Float32(3)+Float32(7)im
cint = 1+im  # issue 108
C = reinterpret(ComplexF64, vec(B))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    mutable struct MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1], "aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = UInt8
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]] # issue #123
typevar = Array{Int}[[1]]
typevar_lb = Vector{<:Integer}[[1]]
typevar_ub = (Vector{U} where U>:Int)[[1]]
typevar_lb_ub = (Vector{U} where Int<:U<:Real)[[1]]
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
objwithpointer = r"julia"
# Custom BitsType (#99)
primitive type MyBT 64 end
bt = reinterpret(MyBT, Int64(55))
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = view([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = Tuple{}[]

iseq(x,y) = isequal(x,y)
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in Base.Sys.CPUinfo.names
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                @warn("Error reading ", $(string(sym)))
                rethrow(e)
            end
            if !iseq(tmp, $sym)
                error("For ", $(string(sym)), ", read value does not agree with written value")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @test a == b
function checkexpr(a::Expr, b::Expr)
    @test a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if (isa(a.args[i], Expr) && a.args[i].head == :line) || isa(a.args[i], LineNumberNode)
            i += 1
            continue
        end
        if (isa(b.args[j], Expr) && b.args[j].head == :line) || isa(b.args[j], LineNumberNode)
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @test i >= length(a.args) && j >= length(b.args)
end

for fn in ("v0.2.26.jld", "v0.2.28.jld")
    for mmap = (true, false)
        undefv = Array{Any}(undef, 1)
        undefm = Array{Any}(undef, 2, 2)
        ms_undef = MyStruct(0)

        fidr = jldopen(joinpath(splitdir(@__FILE__)[1], fn), "r"; mmaparrays=mmap)
        @check fidr x
        !mmap && @check fidr A # FIXME: fails on 0.7 due to a data alignment issue
        @check fidr str
        @check fidr stringsA
        @check fidr stringsU
        @check fidr empty_string
        @check fidr empty_string_array
        @check fidr empty_array_of_strings
        @check fidr tf
        !mmap && @check fidr TF # FIXME: fails on 0.7 due to a data alignment issue
        !mmap && @check fidr AB # FIXME: fails on 0.7 due to a data alignment issue
        @check fidr t
        @check fidr c
        @check fidr cint
        !mmap && @check fidr C # FIXME: fails on 0.7 due to a data alignment issue
        @check fidr emptyA
        @check fidr emptyB
        !mmap && @check fidr ms # FIXME: fails on 0.7 due to a data alignment issue
        @check fidr msempty
        @check fidr sym
        @check fidr syms
        @check fidr d
        exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
        checkexpr(ex, exr)
        @check fidr T
        @check fidr char
        @check fidr unicode_char
        @check fidr α
        !mmap && @check fidr β # FIXME: fails on 0.7 due to a data alignment issue
        !mmap && @check fidr vv # FIXME: fails on 0.7 due to a data alignment issue
        @check fidr rng
        !mmap && @check fidr typevar # FIXME: fails on 0.7 due to a data alignment issue
        !mmap && @check fidr typevar_lb # FIXME: fails on 0.7 due to a data alignment issue
        !mmap && @check fidr typevar_ub # FIXME: fails on 0.7 due to a data alignment issue
        !mmap && @check fidr typevar_lb_ub # FIXME: fails on 0.7 due to a data alignment issue

        # Special cases for reading undefs
        undefv = read(fidr, "undef")
        if !isa(undefv, Array{Any, 1}) || length(undefv) != 1 || isassigned(undefv, 1)
            error("For undef, read value does not agree with written value")
        end
        undefm = read(fidr, "undefs")
        if !isa(undefm, Array{Any, 2}) || length(undefm) != 4 || any(map(i->isassigned(undefm, i), 1:4))
            error("For undefs, read value does not agree with written value")
        end
        ms_undef = read(fidr, "ms_undef")
        if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
            error("For ms_undef, read value does not agree with written value")
        end

        @test !in("objwithpointer", keys(fidr))
        @check fidr bt
        @check fidr sa_asc
        @check fidr sa_utf8
        @check fidr arr_empty_tuple

        x1 = read(fidr, "group1/x")
        @test x1 == Any[1]
        x2 = read(fidr, "group2/x")
        @test x2 == Any[2]

        # load but don't check
        read(fidr, "cpus")

        # subarray (changed representation in julia 0.4)
        if !mmap # FIXME: fails on 0.7 due to a data alignment issue
            subarray_safe = JLD.JLD00.readsafely(fidr, "subarray")
            @test subarray_safe["indexes"] == (1:5,)
            @test subarray_safe["parent"] == [1:5;]
        end

        close(fidr)
    end
end


text:
# using Elliptic
"""
This function calculates the inductance for an N array of coaxial wire loops

The input is an N x 2 array of coordinates [R₁ Θ₁;...;Rₙ Θₙ]
The output is the overall inductance including any mutual inductances

For strange arrangements of loops, this may be useful, but for simple solenoids, other methods are desireable
Empirical formulas (``SimpleInduct'') will be more accurate and less cumbersome in those cases.
"""
function CalcInductance(Coords)
    #Thus function takes in two vectors listing the centroid [in polar coords]
    #of series-wound strands of wire to form a coil. It also takes in the index
    #(of the vector) that different coil "segments" are wound from. An example
    #being if you wanted to break up a Helmholtz pair into the two halves to
    #look at the self-inductances & mutual inductances of each independently
    #you would mark the middle index value in Segment Index List

    for i = 1:(length(Coords[:, 1])-1)
        println(i)
        for j = (i+1):length(Coords[:, 1])
            if Coords[j, :] == Coords[i, :]
                Coords[j, :] = Coords[j, :] .+ 1e-8
            end
        end
    end

    R_List = Coords[:, 1]
    Z_List = Coords[:, 2]
    L = length(R_List)

    μ₀ = 4 * π * 1e-7

    RMesh, AMesh = MPI_Tools.meshgrid(R_List,dims=2)

    KVal = similar(RMesh)
    EValue = similar(RMesh)

    KVal, EValue = (
        [ellipke.(fFunk(R_List, Z_List) .^ 2)[i, j][1] for i = 1:L, j = 1:L],
        [ellipke.(fFunk(R_List, Z_List) .^ 2)[i, j][2] for i = 1:L, j = 1:L],
    )

    M =
        μ₀ .* sqrt.(RMesh .* AMesh) .* (
            (2 .* (KVal .- EValue)) ./ fFunk(R_List, Z_List) .-
            fFunk(R_List, Z_List) .* KVal
        ) #https://thompsonrd.com/OSEE-inductance.pdf
    TotalL = sum(M)

    return TotalL
end

function fFunk(r,z)

    RMesh, AMesh = MPI_Tools.meshgrid(r,dims=2)
    ZMesh_A, ZMesh_B = MPI_Tools.meshgrid(z,dims=2)
    ZDelta = ZMesh_A .- ZMesh_B
    out = sqrt.( (4 .* RMesh.*AMesh ) ./ (ZDelta.^2 .+ (RMesh .+ AMesh).^2) )
    replace!(x -> x>=1 ? 1-eps() : x,out)

    return out
end

"""
SIMPLE INDUCTANCE FORMULAS FOR RADIO COILS*
BY HAROLD A. WHEELER
This is an empirical formula for inductance
    Takes in turns(N), radius (a), and length (b)
    the formulation was designed for inches, but by default
    this function assumes meters.
"""
function SimpleInduct(N,a,b;units = "meters")
    if units=="meters"
        a = a/25.4e-3
        b = b/25.4e-3
    elseif units == "inches"
    else
    error("unknown units")
    end

    L = (a^2*N^2 / (9*a + 10*b))*1e-6
end


text:
include("PRM.jl")
#using Distributions
glvisualize()

####################################
# Randomly Generate Clutter (obstacles)
####################################

# implement function to create clutter
# Let's say on average we want to create 4 obstacles... (otherwise, *most* of
# the time we will product one large rectangle that is falling out of the room

#targetNumObs = 2
function genClutter(roomperimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    obstacles = Vector{HyperRectangle}()
    roomArea = roomWidth * roomHeight
    targetSumObsArea= roomArea * clutterPercentage
    sumObsArea = 0
    while sumObsArea < targetSumObsArea
        x,y = rand(Uniform(0, roomWidth),2)

        randWidth, randHeight = rand(Uniform(1, roomWidth/targetNumObs),2)
        #randWidth, randHeight = rand(1:roomWidth/targetNumObs,2)
        protoObstacle = HyperRectangle(Vec(x, y), Vec(randWidth, randHeight)) #Todo
        if contains(roomperimeter, protoObstacle)
            push!(obstacles, protoObstacle)
            sumObsArea += randWidth*randHeight
        end
    end
    actualSumObsArea = sumObsArea
    #print("obstacles generated")
    return actualSumObsArea, obstacles
end
####

function testGenClutter()

    ####################################
    ## PARAMETERS
    ####################################
    numSamples = 100
    connectRadius = 5
    param = algT.AlgParameters(numSamples, connectRadius)

    targetNumObs = 3
    clutterPercentage = 0.25
    roomWidth,roomHeight  = 20,20

    startstate = Point(1.,1)
    goalstate = Point(19.,19)

    ####################################
    ## INIT
    ####################################
    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()
    for l in wallsPerimeter
        push!(walls, l)
    end
    #plot!(walls, color =:black)
    #plot!(obstacles, fillalpha=0.5)

####################################
# Run PRM once and display config space plot
####################################

    obsArea, obstacles = genClutter(perimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    r = algT.Room(roomWidth,roomHeight,walls,obstacles)
    roomPlot = plotfxn.plotRoom(r)

    ## Run preprocessing
    nodeslist, edgeslist = preprocessPRM(r, param)
    print("\n ---- pre-processed --- \n")

    ## Query created RM

    pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)
    print("\n ---- queried ---- \n")

    ## Plot path found
    timestamp = Base.Dates.now()
    title = "PRM with # samples=$numSamples, maxDist=$connectRadius, \npathcost = $pathcost,
            timestamp=$timestamp)\n\n"
    roadmap = algT.roadmap(startstate, goalstate, nodeslist, edgeslist)

    #prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String, afont)
    prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String)

    print("\n --- Time --- \n")
    @show timestamp
    #print("\n --- Obstacles --- \n")
    #@show obstacles
    #print("\n --- Nodes --- \n")
    #@show nodeslist
    #print("\n --- Edges --- \n")
    #@show edgeslist
    # print("\n -------- \n")

    f = Plots.Font(05,"Liberation Serif")
    plot!(prmPlot,titlefont = f)

    gui(prmPlot)
        print("\n --- Press ENTER to continue --- \n")
        #readline(STDIN)
        # todo: save one "representative" figure from each run
        # print("\n --- continue --- \n")
    @show  isPathFound

end

function clutterExp()
    ####################################
    # Run multiple trials and scatterplot all cost vs area for all runs
    # due to varying absolute area
    ####################################
    ####################################
    ## PARAMETERS
    ####################################
    connectRadius = 8.
    numSamples = 100.
    param = algT.AlgParameters(numSamples, connectRadius)
    targetNumObs = 5
    #clutterPercentage = 0.15
    roomWidth,roomHeight  = 20,20

    startstate = Point(1.,1)
    goalstate = Point(19.,19)

    ####################################
    ## PARAMETERS PART TWO
    ####################################
    nTrials = 10

    #clutterPercentageList = [0 0.01 0.05 0.1 0.12 0.15 0.2 0.25 0.3 0.4 0.5]
    clutterPercentageList = [0 0.1 0.2 0.5]
    #clutterPercentageList = 0 : 0.05 :0.3

    ####################################
    ## Init
    ####################################
    cost = 0
    listCosts = Vector{Float32}()
    listpSucc= Vector{Float32}()

    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()
    for l in wallsPerimeter
        push!(walls, l)
    end

    ####################################
    # Run Trials
    ####################################

#@show listpSucc
    plot()
    tic()
    for pClutter in clutterPercentageList
        totalcost = 0.
        idx = 0
        nSuccess = 0
        pathcost = 0.
        sumSqdError = 0.
        param = algT.AlgParameters(numSamples, connectRadius)
        print("\n ------ \n")
        #@show pClutter
        x_nTrial = Vector{Float32}() #this is unused
        while idx < nTrials
        #    @show idx
            obsAreaList = Vector{Float32}() #this is unused
            idx += 1
            obstacles = Vector{HyperRectangle}()
            obsArea, obstacles = genClutter(perimeter, param, targetNumObs, pClutter, roomWidth, roomHeight)
            r = algT.Room(roomWidth,roomHeight,walls,obstacles)
            nodeslist, edgeslist = preprocessPRM(r, param)
            pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)
            push!(obsAreaList, obsArea) #i am not using this
            if isPathFound
                nSuccess += 1
        #        @show nSuccess
                push!(x_nTrial, pathcost)
            end
            if !(pathcost == Void)
                totalcost += pathcost
            end
        end

        avgCost = totalcost /  nSuccess
        pSucc = nSuccess / nTrials
        #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
        push!(listCosts, avgCost)
        @show pSucc
        push!(listpSucc, pSucc)
    end

    timeExperiment = toc();
    print("Time to run experiment: $timeExperiment\n")
    timestamp = Base.Dates.now()

    ####################################
    # Plot Results
    ####################################
    # clutter (area) on X
    # pathcost on Y
    sizeplot = (800,600)

    @show clutterPercentageList
    @show listpSucc

    supTitle= ""
    # why does glvisualize() have terrible title() issues? oh well comment out for now
    # supTitle="\nPRM with maxDist=$connectRadius,
            # # samples = $numSamples, nTrials=$nTrials.
            # (time to run:$timeExperiment, timestamp=$timestamp)\n\n"
    costTitle= "clutter % vs pathcost\n"

    pPRMcost = scatter(clutterPercentageList, listCosts',
        color=:black,
        title = supTitle * costTitle, ylabel = "euclidean path cost", xlabel = "clutter %",
        yaxis=((0,40), 0:20:40))

    successTitle = "\nclutter % vs pSuccess\n"
    pPRMsuccess = scatter(clutterPercentageList, listpSucc',
        color = :orange, markersize= 6,
        title = successTitle, ylabel = ("P(success)=numSucc/$nTrials trials"), xlabel = "clutter %",
        yaxis=((0, 1.2), 0:0.1:1))

    plot(pPRMcost, pPRMsuccess, layout=(2,1), legend=false,
        xaxis=((0, 0.5), 0:0.05:0.5),
        size = sizeplot)
end
#testGenClutter()

#testGenClutter()
clutterExp()


text:
function sym_mul(a, b)
    @match (a, b) begin
        ( 0, B) => 0
        ( 1, B) => B
        (-1, B::Number) => -B
        (-1, B) => :(-$B)

        (A,  0) => 0
        (A,  1) => A
        (A::Number, -1) => -A
        (A, -1) => :(-$A)

        _ => :($a * $b)
    end
end

function sym_plus(a, b)
    @match (a, b) begin
        (0, 0) => 0
        (0, B) => B
        (A, 0) => A
        _ => :($a+$b)
    end
end

function nexprs(f, n::Int)
    ex = Expr(:block)
    for k in 1:n
        push!(ex.args, f(k))
    end
    return ex
end

function isperm(x::Matrix)
    all(count(x->!isa(x, Number) || !iszero(x), col) == 1 for col in eachcol(x))
end

function isdiag(x::Matrix)
    for i in 1:size(x, 1), j in 1:size(x, 2)
        e = x[i, j]
        if i != j
            e isa Number && iszero(e) && continue
            return false
        end
    end
    return true
end