text:
using DataFrames
using CSV

include("definitions.jl")

#= GLOBAL VARIABLES =#
CURRENT_RECORDS = Ref{Records}()

#= RECORD =#
function record!(run::Run; records = CURRENT_RECORDS)
    @assert (run.name != "") "Not a valid name"
    records.size += 1
    push!(records.head.children, run)
    run.parent = records.head
    records.head = run
end

#= GET =#
function get(run::Run, name::AbstractString)
    @assert (name != "") "Not a valid name"
    nodes = [run]
    while !isempty(nodes)
        item = popfirst!(nodes)
        if item.name == name
            return item
        else
            nodes = vcat(nodes, item.children)
            println(nodes)
        end
    end
    return missing
end

function get(name::AbstractString; records = CURRENT_RECORDS)
    return get(records.sentinel, name)
end

#= CHECKOUT =#
function checkout!(name::AbstractString; records = CURRENT_RECORDS)
    records.head = get(records, name)
    return nothing
end

function checkout!(run::Run; records = CURRENT_RECORDS)
    records.head = run
    return nothing
end

#= Population =#
function population(run::Run; index = 1)
    @assert index <= length(run.population) "Index is larger than population array size"
    return run.population[index]
end

#= INIT =#
function init(;dir::AbstractString = pwd(), runs = "runs/", models = "models/")
    # if !isdir(dir)
    #     dir = joinpath(pwd(), dir)
    # end
    # @assert isdir(dir) "Directory creation failure"
    # @assert !isdir(joinpath(dir, runs)) "Records already exist"
    # @assert !isdir(joinpath(dir, models)) "Models already exist"
    # try
    #     mkdir(joinpath(dir, runs))
    #     mkdir(joinpath(dir, models))
    # catch
    #     "Failure creating directory"
    # end
    records = Records(dir, runs, models)
    records!(records)
    return records
end

function current()
    return CURRENT_RECORDS
end

function records!(records::Records)
    global CURRENT_RECORDS = records
    return nothing
end

#= REINIT =#
function reinit!(;dir = wd)
    records = deserialize(dir)
    global CURRENT_RECORDS = records
    return records
end


text:
using MarketData
using Test

@testset "utils" begin
  @testset "cleanup_colname!" begin  # issue #57
    for (old, new) ∈ [Symbol("Adj. Close")  => :AdjClose,
                      Symbol("Ex-Dividend") => :ExDividend,
                      Symbol("Split Ratio") => :SplitRatio,
                      Symbol("Adj. Open")   => :AdjOpen,
                      Symbol("Adj. High")   => :AdjHigh,
                      Symbol("Adj. Low")    => :AdjLow,
                      Symbol("Adj. Close")  => :AdjClose,
                      Symbol("Adj. Volume") => :AdjVolume]
      ta = TimeArray(cl, colnames = [old])
      MarketData.cleanup_colname!(ta)
      @test colnames(ta) == [new]
    end
  end # @testset "cleanup_colname!"
end # @testset "utils"

ext: jl


text:
using BandedMatrices, LinearAlgebra, LazyArrays, Random, Test
    import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedColumnMajor

@testset "Symmetric" begin
    A = Symmetric(brand(10,10,1,2))
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  2
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (2,2)
    @test MemoryLayout(A) == SymmetricLayout(BandedColumnMajor(), 'U')

    @test A[1,2] == A[2,1]
    @test A[1,4] == 0
    x=rand(10)
    @test A*x ≈ Matrix(A)*x
    @test all(A*x .=== (similar(x) .= Mul(A,x)) .=== (similar(x) .= 1.0.*Mul(A,x) .+ 0.0.*similar(x)) .===
                BLAS.sbmv!('U', 2, 1.0, parent(A).data, x, 0.0, similar(x)))

    A = Symmetric(brand(10,10,1,2),:L)
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  1
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (1,1)
    @test MemoryLayout(A) == SymmetricLayout(BandedColumnMajor(), 'L')

    @test A[1,2] == A[2,1]
    @test A[1,3] == 0
    x=rand(10)
    @test A*x ≈ Matrix(A)*x

    @test all(A*x .=== (similar(x) .= Mul(A,x)) .=== (similar(x) .= 1.0.*Mul(A,x) .+ 0.0.*similar(x)) .===
                BLAS.sbmv!('L', 1, 1.0, view(parent(A).data,3:4,:), x, 0.0, similar(x)))

    # eigvals
    Random.seed!(0)

    A = Symmetric(brand(Float64, 100, 100, 0, 4))
    @test eigvals(A) ≈ eigvals(Symmetric(Matrix(A)))

    # generalized eigvals

    function An(::Type{T}, N::Int) where {T}
        A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
        for n = 0:N-1
            parent(A).data[3,n+1] = T((n+1)*(n+2))
        end
        A
    end

    function Bn(::Type{T}, N::Int) where {T}
        B = Symmetric(BandedMatrix(Zeros{T}(N,N), (0,2)))
        for n = 0:N-1
            parent(B).data[3,n+1] = T(2*(n+1)*(n+2))/T((2n+1)*(2n+5))
        end
        for n = 0:N-3
            parent(B).data[1,n+3] = -sqrt(T((n+1)*(n+2)*(n+3)*(n+4))/T((2n+3)*(2n+5)*(2n+5)*(2n+7)))
        end
        B
    end

    A = An(Float64, 100)
    B = Bn(Float64, 100)

    λ = eigvals(A, B)
    @test λ ≈ eigvals(Symmetric(Matrix(A)), Symmetric(Matrix(B)))

    err = λ*(2/π)^2 ./ (1:length(λ)).^2 .- 1

    @test norm(err[1:40]) < 100eps(Float64)

    A = An(Float32, 100)
    B = Bn(Float32, 100)

    λ = eigvals(A, B)

    @test λ ≈ eigvals(Symmetric(Matrix(A)), Symmetric(Matrix(B)))

    err = λ*(2.f0/π)^2 ./ (1:length(λ)).^2 .- 1

    @test norm(err[1:40]) < 100eps(Float32)
end

@testset "Hermitian" begin
    T = ComplexF64
    A = Hermitian(brand(T,10,10,1,2))
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  2
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (2,2)
    @test MemoryLayout(A) == HermitianLayout(BandedColumnMajor(), 'U')

    @test A[1,2] == conj(A[2,1])
    @test A[1,4] == 0
    x=rand(T,10)
    @test A*x ≈ Matrix(A)*x
    @test all(A*x .=== (similar(x) .= Mul(A,x)) .=== (similar(x) .= T(1.0).*Mul(A,x) .+ T(0.0).*similar(x)) .===
                BLAS.hbmv!('U', 2, T(1.0), parent(A).data, x, T(0.0), similar(x)))

    A = Hermitian(brand(T,10,10,1,2),:L)
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  1
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (1,1)
    @test MemoryLayout(A) == HermitianLayout(BandedColumnMajor(), 'L')

    @test A[1,2] == conj(A[2,1])
    @test A[1,3] == 0
    x=rand(T, 10)
    @test A*x ≈ Matrix(A)*x

    @test all(A*x .=== (similar(x) .= Mul(A,x)) .=== (similar(x) .= one(T).*Mul(A,x) .+ zero(T).*similar(x)) .===
                BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x)))
end


text:
# This file was generated, do not modify it. # hide
tm = TunedModel(model=ensemble_model,
                tuning=Grid(resolution=10), # 10x10 grid
                resampling=Holdout(fraction_train=0.8, rng=StableRNG(42)),
                ranges=[B_range, K_range])

tuned_ensemble = machine(tm, X, y)
fit!(tuned_ensemble, rows=train);


text:
using Test
function testShattered(X::Array{T,N},layer::stParallel,m::Int64) where {T <:Real,N}
    thing = Scattered(layer, X)
    @test typeof(thing.data[1]) == Array{T, N+1}
    @test typeof(thing.output[1]) == Array{T, N+1}
    @test thing.m ==layer.m
    # test the subsampling of the array is correct
    n, q, dataSizes, outputSizes, resultingSize =
        ScatteringTransform.calculateSizes(layer, (-1,-1), size(X))
    @test [size(thing.data[i])[end-2] for i=1:m+1] == sizes(bsplineType(),
                                                       layer.subsampling,
                                                       size(X)[end-1])[1:3]
    @test [size(thing.data[i])[end-1] for i=1:m+1] == sizes(bsplineType(),
                                                       layer.subsampling,
                                                       size(X)[end])[1:3]
    @test [size(thing.output[i])[end-2] for i=1:m+1] == [outputSizes[i][end-2] for
                                                    i=1:m+1]
    @test [size(thing.output[i])[end-1] for i=1:m+1] == [outputSizes[i][end-1] for
                                                   i=1:m+1]
    # check the number of scales
    @test [size(thing.data[i])[end] for i=1:m+1] == [prod([1; q][1:i]) for i=1:m+1]
    @test [size(thing.output[i])[end] for i=1:m+1] == [prod([1; q][1:i]) for i=1:m+1]
end

function testShattering(X::Array{Float64}, nonlinear::S, layer::stParallel, m::Int64) where S <: nonlinearity
    @time outputFull = st(X,layer,nonlinear,thin=false)

    # check that the data isn't NaN
    @test minimum([isfinite(maximum(abs.(outputFull.data[i]))) for i=1:m+1])
    @test minimum([isfinite(maximum(abs.(outputFull.output[i]))) for i=1:m+1])

    # check that the data isn't actually zero when the input is nonzero TODO:
    # when data is functional
    @test minimum(minimum(abs.(outputFull.output[i])) for i=1:m+1) != 0.0

    # see how the thin version is doing
    @time outputThin = st(X,layer,nonlinear,thin=true)
    @test isfinite(maximum(abs.(outputThin)))
    @test minimum(minimum(abs.(outputFull.output[i])) for i=1:m+1) != 0.0

    # compare the two
    @test maximum(abs.(flatten(outputFull, layer)-outputThin))==0.0
    return outputFull
end

@testset "actually shattering" begin
    testShattering(X1, absType(), layer1, 2)
    testShattering(X2, absType(), layer2, 2)
    testShattering(X1, ReLUType(), layer1, 2)
    testShattering(X2, ReLUType(), layer2, 2)
    testShattering(X1, tanhType(), layer1, 2)
    testShattering(X2, tanhType(), layer2, 2)
    testShattering(X1, softplusType(), layer1, 2)
    testShattering(X2, softplusType(), layer2, 2)
end


text:
module AndorSIF

using Images, FileIO

# SIF.jl, adds an imread function for Andor .sif images
# 2013 Ronald S. Rock, Jr.

function FileIO.load(f::File{format"AndorSIF"})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function FileIO.load(fs::Stream{format"AndorSIF"})
    # line 1
    io = stream(fs)
    # line 2
    l = ""
    while !eof(io) && isempty(l)
        l = strip(readline(io))
    end
    l == "65538 1" || error("Unknown Andor version number at line 2: " * l)

    # line 3 TInstaImage thru "Head model"
    l = strip(readline(io))
    fields = split(l)
    fields[1] == "65547" || fields[1] == "65558" || fields[1] == "65567" ||
        error("Unknown TInstaImage version number at line 3: " * fields[1])

    ixon = Dict{Any,Any}()
    ixon["data_type"] = parse(Int,fields[2])
    ixon["active"] = parse(Int,fields[3])
    ixon["structure_vers"] = parse(Int,fields[4]) # (== 1)
    # date is recored as seconds counted from 1970.1.1 00:00:00
    ixon["date"] = parse(Int,fields[5]) # need to convert to actual date
    ixon["temperature"] = max(parse(Float64,fields[6]), parse(Float64,fields[48]))
    ixon["temperature_stable"] = parse(Float64,fields[6]) != -999
    ixon["head"] = fields[7]
    ixon["store_type"] = fields[8]
    ixon["data_type"] = fields[9]
    ixon["mode"] = fields[10]
    ixon["trigger_source"] = fields[11]
    ixon["trigger_level"] = fields[12]
    ixon["exposure_time"] = parse(Float64,fields[13])
    ixon["frame_delay"] = parse(Float64,fields[14])
    ixon["integration_cycle_time"] = parse(Float64,fields[15])
    ixon["no_integrations"] = parse(Int,fields[16])
    ixon["sync"] = fields[17]
    ixon["kin_cycle_time"] = parse(Float64,fields[18])
    ixon["pixel_readout_time"] = parse(Float64,fields[19])
    ixon["no_points"] = parse(Int,fields[20])
    ixon["fast_track_height"] = parse(Int,fields[21])
    ixon["gain"] = parse(Int,fields[22])
    ixon["gate_delay"] = parse(Float64,fields[23])
    ixon["gate_width"] = parse(Float64,fields[24])
    ixon["gate_step"] = parse(Float64,fields[25])
    ixon["track_height"] = parse(Int,fields[26])
    ixon["series_length"] = parse(Int,fields[27])
    ixon["read_pattern"] = fields[28]
    ixon["shutter_delay"] = fields[29]
    ixon["st_center_row"] = parse(Int,fields[30])
    ixon["mt_offset"] = parse(Int,fields[31])
    ixon["operation_mode"] = fields[32]
    ixon["flipx"] = fields[33]
    ixon["flipy"] = fields[34]
    ixon["clock"] = fields[35]
    ixon["aclock"] = fields[36]
    ixon["MCP"] = fields[37]
    ixon["prop"] = fields[38]
    ixon["IOC"] = fields[39]
    ixon["freq"] = fields[40]
    ixon["vert_clock_amp"] = fields[41]
    ixon["data_v_shift_speed"] = parse(Float64,fields[42])
    ixon["output_amp"] = fields[43]
    ixon["pre_amp_gain"] = parse(Float64,fields[44])
    ixon["serial"] = parse(Int,fields[45])
    ixon["num_pulses"] = parse(Int,fields[46])
    ixon["m_frame_transfer_acq_mode"] = parse(Int,fields[47])
    ixon["unstabilized_temperature"] = parse(Float64,fields[48])
    ixon["m_baseline_clamp"] = parse(Int,fields[49])
    ixon["m_pre_scan"] = parse(Int,fields[50])
    ixon["m_em_real_gain"] = parse(Int,fields[51])
    ixon["m_baseline_offset"] = parse(Int,fields[52])
    _ = fields[53]
    _ = fields[54]
    ixon["sw_vers1"] = parse(Int,fields[55])
    ixon["sw_vers2"] = parse(Int,fields[56])
    ixon["sw_vers3"] = parse(Int,fields[57])
    ixon["sw_vers4"] = parse(Int,fields[58])

    # line 4
    ixon["camera_model"] = strip(readline(io))

    # line 5 something like camera dimensions??
    _ = readline(io)

    # line 6
    ixon["original_filename"] = strip(readline(io))

    # line 7
    l = strip(readline(io))
    fields = split(l)
    fields[1] == "65538" || error("Unknown TUserText version number in line 7: $fields[1]")
    usertextlen = parse(Int,fields[2]) # don't need?

    # line 8
    usertext = strip(readline(io))
    # ixon["usertext"] = usertext # Not useful

    # line 9 TShutter
    l = strip(readline(io)) # Weird!

    # line 10 TCalibImage
    l = strip(readline(io))

    # Skip a number of lines until we get to this part:
    # -1
    # 65538 0 512 512 512 512 1 1 1 1
    # 65540
    # 0 1 0 0
    # 0 1 0 0
    # 0 1 0 0
    l = strip(readline(io))
    while l != "-1"
        # repeat until we get there
        l = strip(readline(io))
    end
    # Then skip some lines
    for i = 1:9
        _ = strip(readline(io))
    end

    # what a bizarre file format here
    # length of the next string is in this line
    next_str_len = parse(Int,strip(readline(io)))
    # and here is the next string, followed by the length
    # of the following string, with no delimeter in between!
    l = strip(readline(io))
    next_str_len = parse(Int,l[(next_str_len + 1):end])
    # lather, rinse, repeat...
    l = strip(readline(io))
    next_str_len = parse(Int,l[(next_str_len + 1):end])

    l = strip(readline(io))
    l = l[(next_str_len + 1):end]
    fields = split(l)
    fields[1] == "65538" || fields[1] == "65541" ||
        error("Unknown version number at image dims record")
    ixon["image_format_left"] = parse(Int,fields[2])
    ixon["image_format_top"] = parse(Int,fields[3])
    ixon["image_format_right"] = parse(Int,fields[4])
    ixon["image_format_bottom"] = parse(Int,fields[5])
    frames = parse(Int,fields[6])
    ixon["frames"] = frames
    ixon["num_subimages"] = parse(Int,fields[7])
    ixon["total_length"] = parse(Int,fields[8]) # in pixels across all frames
    ixon["single_frame_length"] = parse(Int,fields[9])

    # Now at the first (and only) subimage
    l = strip(readline(io))
    fields = split(l)
    fields[1] == "65538" || error("unknown TSubImage version number: " * fields[1])
    left = parse(Int,fields[2])
    top = parse(Int,fields[3])
    right = parse(Int,fields[4])
    bottom = parse(Int,fields[5])
    vertical_bin = parse(Int,fields[6])
    horizontal_bin = parse(Int,fields[7])
    subimage_offset = parse(Int,fields[8])

    # calculate frame width, height, with binning
    width = right - left + 1
    mod = width%horizontal_bin
    width = (width - mod) ÷ horizontal_bin
    height = top - bottom + 1
    mod = height%vertical_bin
    height = (height - mod) ÷ vertical_bin

    ixon["left"] = left
    ixon["top"] = top
    ixon["right"] = right
    ixon["bottom"] = bottom
    ixon["vertical_bin"] = vertical_bin
    ixon["horizontal_bin"] = horizontal_bin
    ixon["subimage_offset"] = subimage_offset

    # rest of the header is a timestamp for each frame
    # (actually, just a bunch of zeros). Skip
    for i = 1:frames
        _ = readline(io)
    end

    # two additional lines?
    _ = readline(io)
    _ = readline(io)

    offset = position(io) # start of the actual pixel data, 32-bit float, little-endian

    pixelmatrix = Array{Gray{Float32}}(undef, height, width, frames)
    read!(io, pixelmatrix)
    properties = Dict(
        :ixon => ixon,
    )
    ImageMeta(pixelmatrix, properties)
end

end # module


text:
using MPI
using MPIArrays
using Trilinos
using BenchmarkTools
using Test

include("cartesian.jl")
using .Cartesian

const IndexT = Int32

"""
Determine if a point is part of the Dirichlet boundary conditions
"""
function isdirichlet(grid, I)
  if(I[1] == 1 || I[1] == grid.nx || I[2] == 1 || I[2] == grid.ny)
    return true
  end
  return false
end

"""
Get the dirichlet boundary condition value for the given point
"""
function dirichlet(grid, I)
  return 0.0
end

function buildgraph(comm, grid)
  mpicomm = MPI.Comm(comm)
  cartinds = cartesianindices(grid)

  # Map system GID to grid GID
  systogrid = MPIArray{Int}(mpicomm, (nb_procs,), length(cartinds))
  forlocalpart!(lp -> lp .= localindices(systogrid)[1], systogrid)
  sync(systogrid)
  filter!(i -> !isdirichlet(grid, cartinds[i]), systogrid)
  redistribute!(systogrid)

  # Mapping from grid GID to system GID (or -1 if not part of system)
  gridtosys = MPIArray{Int}(mpicomm, (nb_procs,), length(cartinds))
  forlocalpart!(lp -> fill!(lp, -1), gridtosys)
  sync(gridtosys)

  forlocalpart(systogrid) do sys_lp
    sysgids = localindices(systogrid)[1]
    gridblock = gridtosys[minimum(sys_lp):maximum(sys_lp)]
    gridarr = getblock(gridblock)
    gridglb = GlobalBlock(gridarr, gridblock)
    for (i,gid) in enumerate(sys_lp)
      gridglb[gid] = sysgids[i]
    end
    putblock!(gridarr, gridblock)
  end

  sync(gridtosys)

  # Count the number of non-zero entries in the graph
  my_nnz = forlocalpart(systogrid) do sys_lp
    result = 0
    for gridgid in sys_lp
      result += nb_neighbors(grid, gridgid, !isdirichlet)
    end
    return result
  end

  # Fill the graph with mesh global IDs
  my_nrows = length(localindices(systogrid)[1])
  staticgraph = Kokkos.StaticCrsGraph(my_nrows,my_nnz)
  forlocalpart(systogrid) do sys_lp
    counter = 0
    for (i,gridgid) in enumerate(sys_lp)
      staticgraph.rowmap[i] = counter
      counter += setneighbors!(staticgraph.entries, staticgraph.rowmap[i]+1, grid, gridgid, !isdirichlet)
    end
    @assert counter == staticgraph.rowmap[end]
  end

  # Sync missing off-processor gridnodes
  gridghosts = GhostedBlock(gridtosys, staticgraph.entries)

  # Convert the graph from mesh global IDs to system global IDs
  map!(x -> getglobal(gridghosts,x), staticgraph.entries, staticgraph.entries)

  # sync missing off-processor linear system entries
  sysghosts = GhostedBlock(systogrid, staticgraph.entries)

  # Convert the graph from system global IDs to system local IDs
  map!(x -> globaltolocal(sysghosts,x)-1, staticgraph.entries, staticgraph.entries)

  # Construct the Tpetra maps
  rowmap = Tpetra.Map(length(systogrid), Teuchos.ArrayView(collect(localindices(systogrid)...)), 1, comm)
  col_glb_ids = globalids(sysghosts)
  nb_glb_col_ids = MPI.Allreduce(length(col_glb_ids), +, mpicomm)
  colmap = Tpetra.Map(nb_glb_col_ids, Teuchos.ArrayView(col_glb_ids), 1, comm)

  graph = Tpetra.CrsGraph(rowmap, colmap, staticgraph)
  @assert Tpetra.isFillComplete(graph)

  return graph, sysghosts
end

function fillmatrix!(A, grid)
  nb_rows = Int(Tpetra.getNodeNumRows(A))
  @assert Tpetra.isLocallyIndexed(A)

  row_indices = zeros(IndexT,5)
  row_values_dummy = zeros(Float64,5)
  row_values = fill(-1.0, 5)
  row_values[1] = 4.0

  indices_view = Teuchos.ArrayView(row_indices)
  values_view = Teuchos.ArrayView(row_values_dummy)
  rowlength_ref = Ref(UInt(0))

  Tpetra.resumeFill(A)
  for row in 1:nb_rows
    Tpetra.getLocalRowCopy(A, row-1, indices_view, values_view, rowlength_ref)
    rowlength = Int(rowlength_ref[])
    Tpetra.replaceLocalValues(A, row-1, Teuchos.ArrayView(row_indices,rowlength), Teuchos.ArrayView(row_values,rowlength))
  end
  Tpetra.fillComplete(A)
end

function fillrhs!(b, grid, systogrid)
  b_view = Tpetra.device_view(b)
  gcoords = coordinates(grid)
  for i in LinearIndices(b_view)
    gid = systogrid[i+1]
    x = gcoords[gid,1]
    y = gcoords[gid,2]
    b_view[i] = 2*grid.h^2*((1-x^2)+(1-y^2)) + dirichletrhs(grid, gid, isdirichlet, dirichlet)
  end
end

function laplace2d(comm, grid)
  println("Graph time:")
  @time (graph, systogrid) = buildgraph(comm, grid)

  A = Tpetra.CrsMatrix(graph)
  println("Matrix fill time:")
  @time fillmatrix!(A, grid)

  rangemap = Tpetra.getRangeMap(graph)
  b = Tpetra.Vector(rangemap)
  println("RHS fill time:")
  @time fillrhs!(b, grid, systogrid)

  maxiter = 1000

  params = Trilinos.default_parameters("CG")
  solver_params = params["Linear Solver Types"]["Belos"]["Solver Types"]["CG"]
  solver_params["Convergence Tolerance"] = 1e-12
  solver_params["Verbosity"] = Belos.StatusTestDetails + Belos.FinalSummary + Belos.TimingDetails
  solver_params["Maximum Iterations"] = Int32(maxiter)
  params["Preconditioner Type"] = "MueLu"

  # params = Trilinos.default_parameters()
  # solver_params = params["Linear Solver Types"]["Belos"]["Solver Types"]["BLOCK GMRES"]
  # solver_params["Convergence Tolerance"] = 1e-12
  # solver_params["Verbosity"] = Belos.StatusTestDetails + Belos.FinalSummary + Belos.TimingDetails
  # solver_params["Num Blocks"] = Int32(maxiter)
  # solver_params["Maximum Iterations"] = Int32(maxiter)
  # params["Preconditioner Type"] = "Ifpack2"

  solver = TpetraSolver(A, params)

  return (solver,b,systogrid)
end

# MPI setup
if !MPI.Initialized()
  MPI.Init()
  MPI.finalize_atexit()
end

const comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))
const my_rank = Teuchos.getRank(comm)
const nb_procs = Int(Teuchos.getSize(comm))

if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

function check_solution(sol, g::CartesianGrid)
  coords = coordinates(g)
  gids = LinearIndices(sol)
  cartinds = CartesianIndices(sol)[localindices(sol)...]

  return forlocalpart(sol) do lp
    errorcount = 0
    for (φ,I) in zip(lp,cartinds)
      gid = gids[I]
      x = coords[gid,1]
      y = coords[gid,2]
      errorcount +=  Int(abs(φ - (1-x^2)*(1-y^2)) > 1e-10)
    end
    return errorcount
  end
end

const nx = 1001

function solve_laplace2d(comm)
  h = 2/(nx-1)
  grid = CartesianGrid(nx,nx,h)

  (lows,b,systogrid) = laplace2d(comm, grid)
  (lows,b,systogrid) = laplace2d(comm, grid)
  (lows,b,systogrid) = laplace2d(comm, grid)

  # Compute the solution
  solvec = lows \ b

  cartinds = cartesianindices(grid)
  # Solution in grid coordinates
  gridsol = MPIArray{Float64}(MPI.Comm(comm), (nb_procs,1), grid.nx, grid.ny)

  myrange = CartesianIndices(localindices(gridsol))

  solview = Tpetra.device_view(solvec)
  ghosts = Dict{Int,CartesianIndex{2}}()
  forlocalpart!(gridsol) do lp
    myblock = GlobalBlock(lp, gridsol[localindices(gridsol)...])
    for (i,x) in enumerate(solview)
      gid = systogrid[i]
      I = cartinds[gid]
      if I ∈ myrange
        myblock[I] = x
      else
        ghosts[i-1] = I
      end
    end
  end

  for (i,I) in ghosts
    gridsol[I] = solview[i]
  end

  for gid in myrange
    I = cartinds[gid]
    if isdirichlet(grid,I)
      gridsol[I] = dirichlet(grid, I)
    end
  end

  sync(gridsol)

  @test check_solution(gridsol, grid) == 0
  println("-----------------------------------")
  println("test time:")
  @time check_solution(gridsol, grid)
  @time check_solution(gridsol, grid)
  @time check_solution(gridsol, grid)

  return (grid,gridsol)
end

(grid,sol) = solve_laplace2d(comm)

if my_rank == 0
  using Plots
  heatmap(xrange(grid),yrange(grid),getblock(sol[:,:]),aspect_ratio=1)
  gui()
  readline()
  savefig("laplace2d-result.png")
end


text:
#md# This animation shows the process of fitting the layout.
using CSV
using DataFrames
using WordCloud

df = CSV.File(pkgdir(WordCloud) * "/res/guxiang_frequency.txt", header=false) |> DataFrame;
words = df[!, "Column2"]
weights = df[!, "Column3"]

wc = wordcloud(words, weights, density=0.6)
gifdirectory = "animation2"
@record gifdirectory overwrite=true generate!(wc, 100, optimiser=WordCloud.Momentum())
println("results are saved in animation2")
wc
#eval# runexample(:animation2)
#md# ![](animation2/animation.gif)


text:
"""
    γ = growthrate(τ,μ,ω,ϵ)

Evaluate the number damping rate for the stochastic projected Gross-Pitaevskii equation,
in the simple-growth formulation.

Arguments are in units of `ħω`, where `ω` is the chosen reference trap frequency.

`τ`: reservoir temperature.

`μ`: reservoir chemical potential.

`ω`: frequency unit of simulation (usually one of the trap frequencies).

`ϵ`: cutoff energy, above which the particles form an incoherent reservoir.

External links

[Bose-Einstein condensation from a rotating thermal cloud: Vortex nucleation
and lattice formation, Bradley et al, Physical Review A 77, 033616, 2008](https://arxiv.org/abs/1507.02023)

"""

function growthrate(τ,μ,ω,ϵ)
## bare rate and full cutoff dependent rates.
# requires lerch.m

#opts = F;
MHz     = 1e6
kHz     = 1e3
cm      = 1e-2
mm      = 1e-3
mum     = 1e-6
nm      = 1e-9
muK     = 1e-6
nK      = 1e-9
kB      = 1.38e-23
ħ       = 1.034e-34

amu     = 1.6605402e-27
m       = 87*amu #Rb87
bohR    = 5.29177249e-11
as      = 100*bohR
T       = τ*ħ*ω/kB         #T in Kelvin

γ0      = 4*m/π/ħ^3*(as*kB*T)^2   #All with dimensions

μ       = μ*ħ*ω            #make μ, Ecut have dimensions
ϵ       = ϵ*hbar*w0

μ       = μ/kB/T          #in units of kB*T (for exp(beta...))
ϵ       = ϵ/kB/T
G1i     = log(1-exp(μ-ϵ))^2

Nmax2 = 1000
NmaxLerch = 1000
R = 1:Nmax2
temp =  exp.(2*(μ-ϵ))*exp.(R*(μ-2*ϵ)).*lerch.(exp.(μ-ϵ),1,R+1,NmaxLerch).^2
G2i = sum(temp)

Gi = G1i + G2i

g = γ0*Gi

#Rates in computational units:
gbareComp  = γ0*ħ/kB/T
gComp      = g*ħ/kB/T

#change sinfo.gamma here:

#sinfo.gamma =
return gComp
end


text:
#
# Structure to contain force vector and partial vector for
# parallel computations
#

struct Forces

  # Force array
  f :: Array{Float64}

  # Partial vectors for parallel computation
  fpartial :: Array{Array{Float64}}
  upartial :: Vector{Float64}
  nenc_partial :: Vector{Int64}

  # Constants that have to be computed every time
  sig2 :: Float64
  cutoff2 :: Float64
  sig6 :: Float64
  sig12 :: Float64
  sig66 :: Float64
  sig1212 :: Float64
  eps4 :: Float64

  # Parameters and arrays for linked lists
  nboxes :: Int64
  boxside :: Float64
  ifirst :: Array{Int64}
  inext :: Vector{Int64}

  # Array to check for previous and new encounters
  encij :: Vector{Bool}

end

# Initializer

function Forces(input :: MDInput)

  boxside = input.side/input.cutoff
  nboxes = round(Int64,input.side/boxside)
  ifirst = Array{Int64}(undef,nboxes,nboxes)
  inext = Vector{Int64}(undef,input.n)

  return Forces(zeros(input.n,2), # actual forces
                [ zeros(input.n,2) for i in 1:Threads.nthreads() ], # partials for parallel computations
                zeros(Float64,Threads.nthreads()), # partial energies for parallel computation
                zeros(Int64,Threads.nthreads()), # partial number of encounters
                input.sig^2, input.cutoff^2, input.sig^6, input.sig^12, -6*input.sig^6, -12*input.sig^12, 4*input.eps,
                nboxes, boxside, ifirst, inext,
                zeros(round(Int64,(input.n*(input.n-1))/2)) )

end


text:
const IUP_BUFFER = "BUFFER"
const IUP_STEREO = "STEREO"
const IUP_BUFFER_SIZE = "BUFFER_SIZE"
const IUP_RED_SIZE = "RED_SIZE"
const IUP_GREEN_SIZE = "GREEN_SIZE"
const IUP_BLUE_SIZE = "BLUE_SIZE"
const IUP_ALPHA_SIZE = "ALPHA_SIZE"
const IUP_DEPTH_SIZE = "DEPTH_SIZE"
const IUP_STENCIL_SIZE = "STENCIL_SIZE"
const IUP_ACCUM_RED_SIZE = "ACCUM_RED_SIZE"
const IUP_ACCUM_GREEN_SIZE = "ACCUM_GREEN_SIZE"
const IUP_ACCUM_BLUE_SIZE = "ACCUM_BLUE_SIZE"
const IUP_ACCUM_ALPHA_SIZE = "ACCUM_ALPHA_SIZE"
const IUP_DOUBLE = "DOUBLE"
const IUP_SINGLE = "SINGLE"
const IUP_INDEX = "INDEX"
const IUP_RGBA = "RGBA"
const IUP_YES = "YES"
const IUP_NO = "NO"


text: include("Singletons.jl")


text:
# geothermal energy simulation; see http://fvt.simulkade.com/posts/2017-04-07-heat-transport-in-porous-media.html
# written by Ali A. Eftekhari, while on paternity leave!
# 14 June 2017; Lyngby public library
# output:
# pressure gradient (pressure difference between the injection and the production well)
# outlet temperature: temperature at the middle of the right boundary
# vs simulation time
#

# You can install CoolProp by:
# Pkg.clone("https://github.com/vimalaad/CoolProp.jl")
# Pkg.build("CoolProp")
using CoolProp, DataFrames, PyPlot, Polynomials, JFVM
import GR
# a table of density values
function run_geothermal_sim()
g_gravity = 9.81 # [m/s^2]
depth_reservoir = 2500         # [m]
p_surf= 101325.0               # [Pa]
T_res = 80 + 273.15            # [K]
T_inj = 40 + 273.15            # [K]
p_res = PropsSI("D", "T", T_res, "P", p_surf, "water")*g_gravity*depth_reservoir # [Pa]
T     = collect(linspace(T_inj, T_res, 10))
rho_water = zeros(length(T))
mu_water  = zeros(length(T))
for i in eachindex(T)
    rho_water[i] = PropsSI("D", "T", T[i], "P", p_res, "water")           # [kg/m^3]
    mu_water[i]  = PropsSI("viscosity", "T", T[i], "P", p_res, "water")   # [Pa.s]
end

mu_fit  = polyfit(T, mu_water, 2)
rho_fit = polyfit(T, rho_water, 2)
# subplot(2,1,1)
# scatter(T, mu_water, label = "μ data")
# plot(T, polyval(mu_fit, T))
# xlabel("T [K]")
# ylabel("Viscosity [Pa.s]")
# legend()
# axis([minimum(T), maximum(T), minimum(mu_water), maximum(mu_water)])
#
# subplot(2,1,2)
# scatter(T, rho_water, label = "ρ data")
# plot(T, polyval(rho_fit, T))
# xlabel("T [K]")
# ylabel("Density [kg/m^3]")
# legend()
# run it in a loop:
df_res = DataFrame(q_m3_h = [], k_D = [], L_m = [], csv_name =[])
flow_range = [50.0, 100.0, 150.0, 200.0, 250.0] # [m3/h]
perm_range = [20.0, 50.0, 100.0, 200.0, 500.0, 1000.0]/1000.0 # [D]
L_range    = [500.0, 650.0, 800.0, 900.0, 1000.0] # [m]
for flow_m3_h in flow_range
  for perm_D in perm_range
    for L_res in L_range
      # flow_m3_h = 100        # [m3/h]; change in loop
      flow_prod  = flow_m3_h/3600 # m3/s
      flow_inj   = flow_m3_h/3600*PropsSI("D", "T", T_res, "P", p_surf, "H2O")/PropsSI("D", "T", T_inj, "P", p_surf, "H2O")
      # perm_D  = 0.02        # Darcy; change in loop
      perm    = perm_D*1e-12 # [m^2]
      # L_res    = 500            # [m]; change in loop
      L_left = 500           # [m] distance between the left boundary and injection well
      L_right = 500          # [m] distance between the right boundary and injection well
      L = L_res + L_left + L_right # [m] domain length
      W     = 500            # [m]
      thick_res = 100        # [m]
      poros = 0.20           # [-]
      T0    = 25.0 + 273.15  # [K] reference temperature
      V_dp  = 0.2            # Dykstra-Parsons coef
      clx   = 0.1            # Correlation length in x direction
      cly   = 0.1            # Correlation length in y direction
      T_inj = 40 + 273.15    # [K]
      theta_inj = T_inj - T0 # [K]
      T_res = 80 + 273.15    # [K]
      theta_res = T_res - T0 # [K]
      λ_water  = PropsSI("conductivity", "T", T_res, "P", p_res, "water") # W/m/K
      λ_rock   = 0.16        # [W/m/K] rock conductivity (for sandstone from my thesis)
      λ_eff    = λ_water^poros*λ_rock^(1-poros) # effective conductivity
      cp_water = PropsSI("CPMASS", "T", T_res, "P", p_res, "water") # J/kg/K
      cp_rock  = 837.0       # [J/kg/K] rock heat capacity (for sandstone)
      rho_rock = 2650.0      # [kg/m^3] rock density (for sandstone)

      well_diameter = 7*0.0254 # [m]; converted from a well diameter of 7 inch
      nx_left = 20
      nx_right = 20
      nx_well           = 10
      ny_well           = 10
      nx                = 40
      ny                = 15
      x_range =
        collect([linspace(0, L_left-(nx_well+1)*well_diameter, nx_left);
                 linspace(L_left-nx_well*well_diameter, L_left+nx_well*well_diameter, 2*nx_well);
                 linspace(L_left+(nx_well+1)*well_diameter, L_left+L_res-(nx_well+1)*well_diameter, nx);
                 linspace(L_left+L_res-nx_well*well_diameter, L_left+L_res+nx_well*well_diameter, 2*nx_well);
                 linspace(L_left+L_res+(nx_well+1)*well_diameter, L, nx_right)])
      y_range = collect([linspace(0, 0.5*W-(ny_well+1)*well_diameter, ny);
                         linspace(0.5*W-ny_well*well_diameter, 0.5*W+ny_well*well_diameter, 2*ny_well);
                         linspace(0.5*W+(ny_well+1)*well_diameter, W, ny)])
      ind_inj  = [nx_left+nx_well+1, ny+ny_well]
      ind_prod = [nx_left+nx_well+nx_well+nx+nx_well, ny+ny_well]
      m   = createMesh2D(x_range, y_range)  # 2D Cartesian grid
      v_cell = thick_res*cellVolume(m)
      Nx  = length(x_range)-1
      Ny  = length(y_range)-1

      perm_val   = perm # permfieldlogrnde(Nx, Ny, perm, V_dp, clx, cly)
      perm_field = createCellVariable(m, perm_val)

      left_range  = ny+ny_well:ny+ny_well
      right_range = ny+ny_well:ny+ny_well
      u_inj = flow_inj/(m.cellsize.y[ny+ny_well+1]*thick_res)
      BCp = createBC(m)                 # pressure boundary
      BCp.left.a[:] = 0.0
      BCp.left.b[:] = 1.0
      BCp.left.c[:] = p_res # this makes the convergence easier and does not affect the physics much
      # BCp.right.a[:] = 0.0
      # BCp.right.b[:] = 1.0
      # BCp.right.c[:] = p_res

      # BCp.left.a[left_range]   =
      #     perm_field.value[2, left_range]/polyval(mu_fit, T_inj)
      # BCp.left.b[left_range]   = 0.0
      # BCp.left.c[left_range]   = -u_inj
      # BCp.right.a[right_range] = 0.0
      # BCp.right.b[right_range] = 1.0
      # BCp.right.c[right_range] = p_res

      BCt = createBC(m)                 # temperature boundary
      # Danckwertz with assumptions (density of water is almost constant)
      rho_water_inj  = polyval(rho_fit, T_inj)
      # BCt.left.a[left_range]  = λ_eff
      # BCt.left.b[left_range]  = -rho_water_inj*cp_water*u_inj
      # BCt.left.c[left_range]  = -rho_water_inj*cp_water*u_inj*theta_inj

      # figure(figsize=(8,2))
      # visualizeCells(perm_field)
      # title("Permeability [m^2]", fontsize = 10)
      # colorbar()

      theta_init = createCellVariable(m, theta_res, BCt)
      theta_val  = createCellVariable(m, theta_res, BCt)
      p_init     = createCellVariable(m, p_res, BCp)
      p_val      = createCellVariable(m, p_res, BCp)
      rho_init   = createCellVariable(m,
          polyval(rho_fit, theta_init.value+T0))
      rho_val    = copyCell(rho_init)
      mu_init    = createCellVariable(m,
          polyval(mu_fit, theta_init.value+T0))
      mu_val     = copyCell(mu_init)

      dt_init    = (L*W*thick_res*poros)/flow_inj/50 # (L/Nx)/(u_inj/poros)   # [s] time step
      final_time = 100*dt_init              # [s]
      t_step     = collect([linspace(0, 5*dt_init, 30);
                            linspace(5.5*dt_init, final_time, 70)])
      n_steps = length(t_step) # number of time steps
      # variables to be stored:
      dp_res = zeros(n_steps)
      T_out  = zeros(n_steps)
      p_inj  = zeros(n_steps)
      p_out  = zeros(n_steps)
      T_out[1] = T_res

      # discretization
      M_BCp, RHS_BCp = boundaryConditionTerm(BCp)
      M_BCt, RHS_BCt = boundaryConditionTerm(BCt)
      M_conductivity = diffusionTerm(-harmonicMean(
              createCellVariable(m, λ_eff)))
      q_inj = createCellVariable(m, 0)
      q_inj.value[ind_inj+1 ...] = flow_inj/v_cell.value[ind_inj+1 ...]
      RHSinjt = constantSourceTerm(rho_water_inj*cp_water*(T_inj-T0)*q_inj)

      q_prod = createCellVariable(m, 0)
      q_prod.value[ind_prod+1 ...] = -flow_prod/v_cell.value[ind_prod+1 ...]

      q_cont = createCellVariable(m, 0)

      for t_ind in 2:length(t_step)
        dt = t_step[t_ind]-t_step[t_ind-1]
        t  = t_step[t_ind]
        println(t_ind)
        for i in 1:3 # internal loop
          # solve pressure equation
          RHS_ddt_p   = constantSourceTerm(poros/dt*(rho_val - rho_init))
          water_mobil = harmonicMean(perm_field./mu_val)
          rho_face = arithmeticMean(rho_val)
          M_diff_p    = diffusionTerm(-rho_face.*water_mobil)
          q_cont.value[ind_inj+1 ...] = rho_val.value[ind_inj+1 ...]*flow_inj/v_cell.value[ind_inj+1 ...]
          q_cont.value[ind_prod+1 ...] = -rho_val.value[ind_prod+1 ...]*flow_prod/v_cell.value[ind_prod+1 ...]
          RHScont = constantSourceTerm(q_cont)
          p_val       = solveLinearPDE(m, M_diff_p + M_BCp,
              RHS_BCp - RHS_ddt_p + RHScont)

          # velocity vector
          u = -water_mobil.*gradientTerm(p_val)

          # solve heat equation
          α                  = poros*rho_val*cp_water+(1-poros)*rho_rock*cp_rock
          M_trans, RHS_trans = transientTerm(theta_init, dt, α)
          RHS_ddt_t          = constantSourceTerm(cp_water*poros*theta_val.*(rho_val - rho_init)/dt)
          M_conv             = convectionUpwindTerm(cp_water*rho_face.*u)
          # RHSprodt = constantSourceTerm(-rho_val.value[ind_prod+1 ...]*cp_water*T0*q_prod)
          q_inj.value[ind_inj+1 ...] = flow_inj/v_cell.value[ind_inj+1 ...]
          RHSinjt = constantSourceTerm(rho_water_inj*cp_water*(T_inj-T0)*q_inj)
          Mprodt   = linearSourceTerm(rho_val.value[ind_prod+1 ...]*cp_water*q_prod)
          theta_val          = solveLinearPDE(m,
              M_BCt + M_conv + M_trans + M_conductivity - Mprodt,
              RHS_BCt + RHS_trans + RHSinjt - RHS_ddt_t)

          # update density and viscosity values
          rho_val.value[:] = polyval(rho_fit, theta_val.value[:] + T0)
          mu_val.value[:]  = polyval(mu_fit, theta_val.value[:] + T0)
          flow_inj   = flow_prod*rho_val.value[ind_prod+1 ...]/rho_water_inj
        end # end of inner loop
        # GR.imshow(internalCells(theta_val)+T0)
        T_out[t_ind]  = theta_val.value[ind_prod+1 ...]+T0
        p_inj[t_ind]  = p_val.value[ind_inj+1 ...]
        p_out[t_ind]  = p_val.value[ind_prod+1 ...]
        dp_res[t_ind] = p_val.value[ind_inj+1 ...]-p_val.value[ind_prod+1 ...]
                      #  0.5*(p_val.value[end, 1+ny+ny_well]+p_val.value[end-1, 1+ny+ny_well])
        rho_init   = copyCell(rho_val)
        p_init     = copyCell(p_val) # not necessary
        theta_init = copyCell(theta_val)
      end
      dp_res[1] = dp_res[2] # pressure grad at time zero; assume that process has already started
      p_inj[1] = p_inj[2]
      p_out[1] = p_out[2]
      df = DataFrame(T_s = t_step, dp_Pa = dp_res, T_K = T_out, p_inj_Pa = p_inj, p_out_Pa = p_out)
      f_name = "results/Q-$flow_m3_h-L-$L_res-k$perm_D.csv"
      DataFrames.writetable(f_name, df)
      push!(df_res, [flow_m3_h, perm_D, L_res, f_name])
      # figure(figsize=(8,2))
      # visualizeCells(theta_val+T0)
      # title("temperature profile")
      # colorbar()
      # figure()
      # plot(t_step/(365*24*3600), T_out, ".")
      # xlabel("time [year]")
      # ylabel("Water temperature [K]")
      # figure()
      # plot(t_step/(365*24*3600), dp_res/1e5)
      # xlabel("time [year]")
      # ylabel("Pressure difference [bar]")
      # figure()
      # plot(t_step/(365*24*3600), p_inj/1e5)
      # plot(t_step/(365*24*3600), p_out/1e5)
      # xlabel("time [year]")
      # ylabel("Pressure [bar]")
    end
  end
end
DataFrames.writetable("results/res_guide.csv", df_res)
end

max_stars_repo_name: simulkade/geothermal

max_issues_repo_name: simulkade/geothermal

max_forks_repo_name: simulkade/geothermal


text:
# load my files -
include("Include.jl")

function update_model_dictionary(parameter_array,default_model_dictionary)

    # what is the host_type?
    host_type = :cell_free

    # make a deepcopy -
    model_data_dictionary = deepcopy(default_model_dictionary)

    # Phase 1: parameter update =========================================================================== #
    # update the paramaters in the model data dictionary -
    # for now - lets only search over dG's -
    R = model_data_dictionary["R"]
    T_K = model_data_dictionary["T_K"]

    # compute W -
    tmp_W_array = Float64[]
    for index = 1:2
        parameter_guess = parameter_array[index]
        value = exp(-1*(parameter_guess/1.0)/(R*T_K))
        push!(tmp_W_array,value)
    end

    # update the control W's -
    control_parameter_dictionary = model_data_dictionary["control_parameter_dictionary"]
	control_parameter_dictionary["W_deGFP_RNAP"] = tmp_W_array[1]       # 3
	control_parameter_dictionary["W_deGFP_sigma_70"] = tmp_W_array[2]   # 4
    model_data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary

    binding_parameter_dictionary = model_data_dictionary["binding_parameter_dictionary"]
	binding_parameter_dictionary["n_deGFP_sigma_70"] = parameter_array[3]
	binding_parameter_dictionary["K_deGFP_sigma_70"] = parameter_array[4]
    model_data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary

    # time constant modifier -
    # time constant modifiers -
	time_constant_modifier_array = [
		0.0	                        ;	# 1	deGFP
		0.0	                        ;	# 2	sigma_70
		parameter_array[5] 	        ;	# 3	mRNA_deGFP
		1.0	                        ;	# 4	mRNA_sigma_70
		parameter_array[6] 	        ;	# 5	protein_deGFP
		1.0	                        ;	# 6	protein_sigma_70
    ]
    model_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

    # setup degradation_modifier_array -
    degradation_modifier_array = [
		0.0	;	# 1	deGFP
		0.0	;	# 2	sigma_70
		parameter_array[7]	;	# 3	mRNA_deGFP
		1.0	;	# 4	mRNA_sigma_70
		parameter_array[8]	;	# 5	protein_deGFP
		parameter_array[9]	;	# 6	protein_sigma_70
    ]

    # update the translation time -
    model_data_dictionary["half_life_translation_capacity"] = parameter_array[10]

    # lastly, update KL -
    biophysical_constants_dictionary = model_data_dictionary["biophysical_constants_dictionary"]
    biophysical_constants_dictionary["translation_saturation_constant"] = parameter_array[11]
    model_data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary

    # grab defaults -
    species_symbol_type_array = model_data_dictionary["species_symbol_type_array"]
    protein_coding_length_array = model_data_dictionary["protein_coding_length_array"]
    gene_coding_length_array = model_data_dictionary["gene_coding_length_array"]
    time_constant_modifier_array = model_data_dictionary["time_constant_modifier_array"]
    initial_condition_array = model_data_dictionary["initial_condition_array"]

    # # get gene IC -
    idx_gene = findall(x->x==:gene,species_symbol_type_array)
    gene_abundance_array = initial_condition_array[idx_gene]

    # Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
    model_data_dictionary["translation_parameter_array"] = translation_parameter_array

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
    model_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

    # Dilution degrdation matrix -
    dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary, species_symbol_type_array, degradation_modifier_array)
    model_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix

    # return -
    return model_data_dictionary
end

function main(path_to_ensemble_file::String, path_to_sim_dir::String; sample_fraction::Float64=1.0)

    # load the default data_dictionary -
    time_start = 0.0
    time_stop = 16.0
    time_step_size = 0.01

    # what is the host_type?
    host_type = :cell_free

    # path to parameters -
    path_to_biophysical_constants_file = "./CellFree.json"

    # Load the data dictionary (uses the default biophysical_constants file)
    default_data_dictionary = build_data_dictionary((time_start,time_stop,time_step_size), path_to_biophysical_constants_file, host_type)

    # customize -
    customized_data_dictionary = customize_data_dictionary(default_data_dictionary, host_type)

    # load the ensemble file -
    parameter_ensemble_array = readdlm(path_to_ensemble_file)

    # sort the array -
    parameter_ensemble_array = sort_parameter_ensemble_array(parameter_ensemble_array)

    # what is the size of the ensemble -
    (number_of_parameters, number_of_samples) = size(parameter_ensemble_array)

    # take the top sample fraction -
    number_of_samples = round(Int,sample_fraction*number_of_samples)
    for sample_index = 1:number_of_samples

        # grab the parameter array -
        local_parameter_array = parameter_ensemble_array[:,sample_index]

        # update the default data_dictionary to reflect the new parameters -
        model_data_dictionary = update_model_dictionary(local_parameter_array, customized_data_dictionary)

        # solve the model equations -
        (T,X) = SolveBalances(time_start,time_stop,time_step_size,model_data_dictionary)

        # dump -
        data_array = [T X]
        filename = "$(path_to_sim_dir)/simulation-P$(sample_index).dat"
        writedlm(filename, data_array)

        # give user some notification -
        @show sample_index
    end
end

# setup paths -
path_to_sim_file = "$(pwd())/simulations"
path_to_ensemble_file = "$(pwd())/Ensemble-poets-MVP-T20.dat"

# call -
main(path_to_ensemble_file, path_to_sim_file; sample_fraction = 1.0)

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
# Uniform Stein Distribution
#
# Represents a uniform distribution with independent components

type SteinUniform <: SteinDistribution
    range::Array{Float64, 2}
    c1::Float64
    c2::Float64
    c3::Float64
end

# Constructor for p independent uniform([0,1]) variables with non-uniform constants
SteinUniform(p::Int64, c1::Float64, c2::Float64, c3::Float64) =
SteinUniform(repmat([0.0; 1.0], 1, p), c1, c2, c3)

# Constructor with best known Stein factors setting (c1,c2,c3) = (1,1,1)
SteinUniform(range::Array{Float64, 2}) = SteinUniform(range, 1.0, 1.0, 1.0)

# Constructor for p independent uniform([0,1]) variables
SteinUniform(p::Int64) = SteinUniform(repmat([0.0; 1.0], 1, p))

# Univariate uniform([0,1]) constructor
SteinUniform() = SteinUniform(1)

# Draw n independent samples from distribution
function rand(d::SteinUniform, n::Int64)
    p = size(d.range, 2)
    # Sample uniform [0,1] variables
    x = zeros(n, p)
    rand!(x)
    # Shift and scale as needed
    shift = d.range[1,:]
    scale = (d.range[2,:] - d.range[1,:])
    scaled_x = broadcast(*, x, scale')
    broadcast(+, scaled_x, shift')
end

function supportlowerbound(d::SteinUniform, j::Int64)
    d.range[1, j];
end

function supportupperbound(d::SteinUniform, j::Int64)
    d.range[2, j];
end

# Gradient of log density of distribution evaluated at each row of
# n x p matrix X
function gradlogdensity(d::SteinUniform, X::Array{Float64,1})
    zeros(size(X));
end

function gradlogdensity(d::SteinUniform, X::Array{Float64,2})
    zeros(size(X));
end

# Cumulative distribution function (only valid when X is univariate)
function cdf(d::SteinUniform,t)
    @assert numdimensions(d) == 1
    lower = supportlowerbound(d,1);
    upper = supportupperbound(d,1);
    min(max((t-lower)/(upper-lower),0),1);
end

function numdimensions(d::SteinUniform)
    size(d.range, 2)
end


text:
__precompile__()

@doc raw"""
Define a Parsing Expression Grammar via a macro and abuse of Julia syntax.

* Rules: `@rule name = expression`
* Choice: infix `,`
* Sequence: infix `&`
* Positive lookahead: prefix `+`
* Negative lookahead: prefix `-`
* Option (zero or one time): postfix `[:?]` (≡ `[0:1]`)
  * (note that [?] won't work in Julia >= 1.0 per JuliaLang/julia#22712)
* Zero or more times: postfix `[*]` (≡ `[0:end]`)
* One or more times: postfix `[+]` (≡ `[1:end]`)
* Exactly `m` times: postfix `[m]` (≡ `[m:m]`) (where m is an integer)
* Between `m` and `n` times inclusive: postfix `[m:n]`
* At most `n` times: postfix `[0:n]`
* At least `m` times: postfix `[m:end]`
* Terminals: `r"regex"`, `"string"`
  * Extra regex flags: `p` is for punctuation, and eats whitespace (`\s*`)
    after the match; `w` is for word, and implies `p`, but also makes sure
    match boundaries are word boundaries (`\b`); `h` modifies `p` and `w` to
    eat only horizontal whitespace (`\h`). Values passed to semantics functions
    exclude eaten whitespace.
* Semantics: `expression |> unary_function` (like ParserCombinator)
  * or `expression > nary_function` to interpolate args.
  * Returning the special singleton value `PEG.Failure()` from a semantics
    function causes the parsing expression it's attached to to fail (return
    `nothing` instead of a tuple). Returning `nothing` from a semantics
    function is not special; it just makes the first part of the tuple
    `nothing`. See the parsing function signature below.

Put another way:

```julia
using PEG
@rule grammar = "using PEG\n" & rule[*]
@rule rule = r"@rule"p & nonterminal & r"="p & choice
@rule choice = seq & (r","p & seq)[*]
@rule seq = item & (r"&"p & item)[*] & (r"\|?>"p & julia_function)[:?]
@rule item = lookahead , counted
@rule lookahead = r"\("p & (r"[+-]"p) & seq & r"\)"p
@rule counted = single & (count)[:?]
@rule count = range , r"\["p & (":?" , r"[\*\+]"p) & r"]"p
@rule range = r"\["p & integer & (r":"p & (integer , r"end"w))[:?] & r"]"p
@rule integer = r"\d+"w
@rule single = parens , terminal , nonterminal
@rule parens = r"\("p & choice & r"\)"p
@rule nonterminal = r"\pL\w+"w
@rule terminal = regex , string & r"\s*"
@rule regex = r"\br" & string & r"[himpswx]*\s*"
@rule string = r"\"(\\.|[^\"])*\""
@rule julia_function = # left as an exercise ;)
```

Each rule defines a parsing function with the following signature:

```julia
nonterminal(input::T, cache=PEG.Cache()) where T <: AbstractString
  ::Union{Nothing,Tuple{Any,SubString}}
```

The `Any` part of the return value is the abstract syntax tree, while the
`SubString` is the remaining input after the parsed portion. If parsing fails,
`nothing` is returned.

While you can use rules defined in this way directly, it might be more
convenient to use the functions `parse_next(rule, input; whole=false)` or
`parse_whole(rule, input)`. See their documentation for more information.

Call `PEG.setdebug!()` to have debugging information printed during parsing.
Call `PEG.setdebug!(false)` to turn it off again.
"""
module PEG
export @rule, parse_next, parse_whole, squash_ast

global debug = false
"""
    setdebug!(val=true)

Enable or disable printing of debugging information.
"""
function setdebug!(val::Bool=true)
  global debug = val
end

"""
A Cache maps from (rule name or gensym, input string length remaining) to
either nothing or (parsed value, remaining input substring) (which is what rule
functions return).
"""
const Cache = Dict{Tuple{Symbol,Int},Union{Nothing,Tuple{Any,SubString}}}

function cache_rule(sym::Symbol, fn::Function, input::SubString, cache::Cache)
  local key = (sym, lastindex(input))
  haskey(cache, key) && return cache[key]
  if debug
    if !occursin(r"^##", string(sym))
      println(" "^length(stacktrace()) * string(sym))
    end
    cache[key] = fn(input, cache)
    if cache[key] != nothing
      println(" "^length(stacktrace()) * "$sym matched " * string(lastindex(input)) * ":" * string(nextind(cache[key][2],lastindex(cache[key][2]))) * " bytes from end of input, returning " * string(cache[key][1]))
#    else
#      println(" "^length(stacktrace()) * string(sym) * " failed to match " * string(lastindex(input)) * " bytes from end of input")
    end
    return cache[key]
  else
    return (cache[key] = fn(input, cache))
  end
end

# terminal ""
function to_rule(str::AbstractString)
  # not cached because I'm guessing startswith is probably faster
  (input, cache)->begin
    if startswith(input, str)
      (str, input[nextind(str,lastindex(str)):end])
    end
  end
end

# nonterminal
function to_rule(sym::Symbol)
  # not cached here because it's already cached in the callee rule
  :((input, cache) -> $(esc(sym))(input, cache))
  #                   ^^^^^^^^^^^
  # This is why some to_rule methods return expressions instead of anonymous
  # functions. Anything that could contain a nonterminal needs to be an
  # expression, so that this esc call appears literally in the expression the
  # @rule macro returns. If we used eval instead of esc here, and raw anonymous
  # functions everywhere, then sym would be evaluated in the context of the PEG
  # module instead of in the caller's module.
end

to_rule(expr::Expr) = to_rule(Val(expr.head), expr.args...)
to_rule(head::Val{:call}, name::Symbol, args...) =
  to_rule(head, Val(name), args...)
to_rule(head::Val{:macrocall}, name::Symbol, source::LineNumberNode, args...) =
  to_rule(head, Val(name), args...)
to_rule(head, args...) = error("can't convert $head expression to PEG rule; args=$args")

# variant of to_rule(::Symbol) for qualified rule references
to_rule(::Val{:.}, args...) =
  :((input, cache) -> $(esc(Expr(:., args...)))(input, cache))

"""
    remove_re_flag(flags, flag)

Return (a copy of) regex flags with flag removed, and a boolean indicating whether the flag was there in the first place.
"""
function remove_re_flag(flags::T, flag::Char) where T <: AbstractString
  local i = findfirst(isequal(flag), flags)
  if i == nothing
    (flags, false)
  else
    (flags[1:prevind(flags, i)] * flags[nextind(flags, i):end], true)
  end
end

# terminal r""
function to_rule(::Val{:macrocall}, ::Val{Symbol("@r_str")}, str::String, flags...)
  if length(flags) == 0
    flags = ""
  else
    flags = flags[1]
  end
  # p = punctuation; allows whitespace after the match
  local p
  flags, p = remove_re_flag(flags, 'p')
  # w = word; ensures both sides of the match are word boundaries; implies p
  local w
  flags, w = remove_re_flag(flags, 'w')
  p = p || w
  local h
  flags, h = remove_re_flag(flags, 'h')
  str = "^" * (w ? "\\b" : "") * "($str)" * (w ? "\\b" : "") * (p ? (h ? "\\h*" : "\\s*") : "")
  re = Regex(str, flags)
  local sym = gensym(string(re))
  (input, cache)->cache_rule(sym, (input, cache)->begin
    local m = match(re, input)
    m == nothing && return
    (m.captures[1], input[nextind(m.match,lastindex(m.match)):end])
  end, input, cache)
end

# fn[m:n]
function to_rule(::Val{:ref}, fn, range::Expr)
  fn = to_rule(fn)
  (range.head == :call && range.args[1] == :(:)) ||
    error("can't convert non-range ref expression to PEG rule: $range")
  local m
  local n
  m, n = range.args[2:3]
  n == :end && (n = Inf)
  local sym = gensym("[$m:$n]")
  :((input, cache)->cache_rule($(Meta.quot(sym)), (input, cache)->begin
    local results = []
    local i = 0
    while i < $n
      local m = $fn(input, cache)
      m == nothing && break
      local result
      result, input = m
      push!(results, result)
      i += 1
    end
    i < $m && return
    (results, input)
  end, input, cache))
end

# defined in terms of the above method
to_rule(::Val{:ref}, fn, n::Int) = to_rule(Val(:ref), fn,:($n:$n))
to_rule(::Val{:ref}, fn, n::Symbol) = to_rule(Val(:ref), fn, Val(n))
to_rule(::Val{:ref}, fn, n::QuoteNode) = to_rule(Val(:ref), fn, Val(n.value))
to_rule(::Val{:ref}, fn, ::Val{:?}) = to_rule(Val(:ref), fn,:(0:1))
to_rule(::Val{:ref}, fn, ::Val{:*}) = to_rule(Val(:ref), fn,:(foo[0:end]).args[2])
to_rule(::Val{:ref}, fn, ::Val{:+}) = to_rule(Val(:ref), fn,:(foo[1:end]).args[2])

# positive lookahead
function to_rule(::Val{:call}, ::Val{:+}, fn)
  fn = to_rule(fn)
  # not cached here because it's already cached in fn
  :((input, cache)->begin
    local m = $fn(input, cache)
    m == nothing ? nothing : (nothing, input) # NOTE: original input
  end)
end

# negative lookahead
function to_rule(::Val{:call}, ::Val{:-}, fn)
  fn = to_rule(fn)
  # not cached here because it's already cached in fn
  :((input, cache)->begin
    local m = $fn(input, cache)
    m == nothing ? (nothing, input) : nothing
  end)
end

"Given e.g. x=:(((a op b) op c) op d), return [a,b,c,d]."
flatten_op(x, op::Symbol) = Any[x]
function flatten_op(x::Expr, op::Symbol)
  if x.head == :call && x.args[1] == op
    push!(flatten_op(x.args[2], op), x.args[3])
  else
    Any[x]
  end
end

# sequence
function to_rule(::Val{:call}, ::Val{:&}, a, b)
  local seq = map(to_rule, [flatten_op(a, :&); b])
  local sym = gensym(:&)
  :((input, cache)->cache_rule($(Meta.quot(sym)), (input, cache)->begin
    local results = []
    local item
    for item ∈ [$(seq...)]
      local m = item(input, cache)
      m == nothing && return
      local result
      result, input = m
      push!(results, result)
    end
    (results, input)
  end, input, cache))
end

"Singleton value that semantics functions may return to cause the parsing
expression they were called from to fail, e.g. `return Failure()`."
struct Failure
end

# semantics
function to_rule(::Val{:call}, ::Val{:|>}, pe, fn)
  pe = to_rule(pe)
  local sym = gensym(:|>)
  :((input, cache)->cache_rule($(Meta.quot(sym)), (input, cache)->begin
    local m = ($pe)(input, cache)
    m == nothing && return nothing
    local pe_result
    pe_result, input = m
    fn_result = $(esc(fn))(pe_result)
    fn_result == Failure() && return nothing
    (fn_result, input)
  end, input, cache))
end

function to_rule(::Val{:call}, ::Val{:>}, pe, fn)
  pe = to_rule(pe)
  local sym = gensym(:>)
  :((input, cache)->cache_rule($(Meta.quot(sym)), (input, cache)->begin
    local m = ($pe)(input, cache)
    m == nothing && return nothing
    local pe_result
    pe_result, input = m
    fn_result = $(esc(fn))(pe_result...) # "..." is the only difference from ^^^
    fn_result == Failure() && return nothing
    (fn_result, input)
  end, input, cache))
end

# choice
function to_rule(::Val{:tuple}, args...)
  local choice = map(to_rule, args)
  local sym = gensym(",")
  :((input, cache)->cache_rule($(Meta.quot(sym)), (input, cache)->begin
    local item
    for item ∈ [$(choice...)]
      local m = item(input, cache)
      m == nothing || return m
    end
    return nothing
  end, input, cache))
end

"""
    @rule name = parsing_expression...

Defines a single rule in a parsing expression grammar as a function. See the
documentation of the PEG module for a description of parsing expression
grammars, and the type and usage of the rule functions defined using `@rule`.
"""
macro rule(assignment::Expr)
  assignment.head == :(=) || error("expected = Expr, but got $assignment")
  local name
  local value_expr
  name, value_expr = assignment.args
  local value_fn = to_rule(value_expr)
  :(@Base.__doc__ function $(esc(name))(input::T, cache::Cache=Cache()) where T <: AbstractString
      if !isa(input, SubString)
	input = SubString(input, 1, lastindex(input))
      end
      local sym = Symbol($(string(name)))
      cache_rule(sym, $value_fn, input, cache)
    end)
end

"""
    parse_next(rule, input; whole=false)

Parse a prefix of the `input` string as one instance of the given `rule`. This
differs from just calling the `rule` itself on the `input` in a few important
ways:

* When parsing succeeds, `parse_next` returns only the parsed value, while
  `rule` returns a Tuple of the parsed value and the remaining unparsed input.
* When parsing fails, `parse_next` throws an exception, with information on
  what failed to match starting at the latest point in the string that any
  parsing expression matched up to. `rule` just returns `nothing`.
* When `whole=true`, parsing only succeeds if the whole input is consumed, i.e.
  the second value in the tuple that `rule` returns is `""`.
"""
function parse_next(rule::Function, input::T; whole=false) where T <: AbstractString
  local cache = Cache()
  local m = rule(input, cache)
  if m == nothing || (whole && m[2] != "") # failed to parse
    # find the last index we tried to parse anything starting at, and all the
    # cache keys for the expressions we tried to parse there
    # FIXME plain strings aren't cached, so we'll miss them as keys
    local last_index = nextind(input,lastindex(input)) # start off the end of input
    local last_keys = Symbol[]
    for pair ∈ cache
      if pair[1][2] < last_index
	last_index = pair[1][2]
	last_keys = [pair[1][1]]
      elseif pair[1][2] == last_index
        push!(last_keys, pair[1][1])
      end
    end
    last_index = nextind(input,lastindex(input)) - last_index # ugh, backwards 1-based indexing
    # convert regex gensym keys to something more readable, remove other
    # gensyms, and convert everything to strings
    last_keys = map(x->begin
      local s = string(x)
      local m = match(r"^##(r[^#]+)#\d+$", s)
      if m != nothing
	m.captures[1]
      elseif occursin(r"^##", s)
	nothing
      else
	s
      end
    end, last_keys)
    filter!(x->(x != nothing), last_keys)
    # translate the index into line and column and get the text on the line
    local before = split(input[1:prevind(input,last_index)], r"\r\n|\n\r|\n|\r")
    local after = replace(input[last_index:end], r"[\r\n].*" => "")
    local line_num = length(before)
    local column_num = length(before[end]) + 1
    local line = before[end] * after
    local message = "On line $line_num, at column $column_num (byte $last_index):\n$line\n" * " "^Int(clamp(column_num-1, 0, Inf)) * "^ here\nexpected one of the following: " * join(last_keys, ", ") * "\n"
    debug && print(message)
    throw(Meta.ParseError(message))
  else # parse succeeded
    m[1]
  end
end

"    parse_whole(rule, input) = parse_next(rule, input; whole=true)"
parse_whole(rule, input) = parse_next(rule, input; whole=true)

"""
    squash_ast(x)

Given an ast returned by a grammar with no semantics, return a simplified
version, such that there are no `nothing`s, and no `Any[]`s with exactly zero
or one element.
"""
squash_ast(x) = x
function squash_ast(children::Array{Any,1})
  squashed_children = map(squash_ast, children)
  filter!(x->(x != nothing), squashed_children)
  if length(squashed_children) == 0
    nothing
  elseif length(squashed_children) == 1
    squashed_children[1]
  else
    squashed_children
  end
end

end


text:
#########
# setk! #
#########

function setk!(x, k_requested, v)
    k = setk!_inner(x, k_requested, v)
    if k > 0
        error("Index ", k_requested, " is out of bounds.")
    end
end
@generated function setk!_inner(x::T, k, v) where {T}
    if T <: Array
        if eltype(T) <: Real
            quote
                if k <= length(x)
                    x[k] = v
                    return 0
                end
                return k-length(x)
            end
        else
            quote
                for c = 1:length(x)
                    k = setk!_inner(x[c], k, v)
                    if k <= 0
                        return 0
                    end
                end
            end
        end
    elseif T <: Nothing
        quote
            # We shouldn't be trying to set a Nothing. Skip it.
            # if k <= 1
            #     # x[k] = v
            #     return 0
            # end
            return k
        end
    elseif !isempty(fieldnames(T))
        expr = Vector{Expr}()
        for f in fieldnames(x)
            push!(expr, quote
                if typeof(x.$f) <: Real # For scalars, we must overwrite the field directly.
                    if k == 1
                        x.$f = v
                        return 0
                    end
                    k -= 1
                elseif typeof(x.$f) <: Nothing # For nothings, we ignore the write.
                    # if k == 1
                    #     x.$f = nothing
                    #     return 0
                    # end
                    k # We didn't use any index for this.
                else # For anything else, it's assumed that we can recur and set something inside the field.
                    k = setk!_inner(x.$f, k, v)
                    if k <= 0
                        return 0
                    end
                end
            end)
        end
        push!(expr, :( return k ))
        Expr(:block, expr...)
    else
        error("I don't know how to set an index of a ", T, ".")
    end
end

########
# getk #
########

# Get index k of x, which might be a struct of struct of vectors or something.
function getk(x, k_requested)
    r, k = getk_inner(x, k_requested)
    if k > 0
        error("Index ", k_requested, " is out of bounds.")
    end
    return r
end
@generated function getk_inner(x::T, k) where {T}
    if T <: Real
        quote
            if k == 1
                return (x, 0)
            end
            return (0., k-1)
        end
    elseif T <: Nothing
            quote
                # if k == 1
                #     return (0., 0)
                # end
                return (0., k)
            end
    elseif T <: Array
        if eltype(T) <: Real
            quote
                if k <= length(x)
                    return (x[k], 0)
                end
                return (0., k-length(x))
            end
        else
            quote
                for c = 1:length(x)
                    r, k = getk_inner(x[c], k)
                    if k <= 0
                        return (r, 0)
                    end
                end
            end
        end
    elseif !isempty(fieldnames(T))
        expr = Vector{Expr}()
        push!(expr, :( r = 0. ))
        for f in fieldnames(x)
            push!(expr, quote
                r, k = getk_inner(x.$f, k)
                if k <= 0
                    return (r, 0)
                end
            end)
        end
        push!(expr, :( return (0., k) ))
        Expr(:block, expr...)
    else
        error("I don't know how to get an index of a ", T, ".")
    end
end

###########
# stackem #
###########

# Stacks reals, arrays, and anything with fieldnames into a new vector, returning the vector.
stackem(y::T) where {T} = stackem!(Vector{Float64}(), y)

# Stacks reals, arrays, and anything with fieldnames into an existing vector. Super fast!
# We write this as an if on the type, because this allows us to do more than multiple dispatch.
# E.g., we can detect a struct by asking for its field names.
@generated function stackem!(x, y::T) where {T}
    if T <: Real
        :( push!(x, y) )
    elseif T <: Enum
        :( push!(x, Base.Enums.basetype(T)(y)) )
    elseif T <: Array
        if eltype(T) <: Real
            :( append!(x, y) )
        else
            quote
                for c = 1:length(y)
                    stackem!(x, y[c])
                end
                x
            end
        end
    elseif !isempty(fieldnames(T)) # This is nice; we can't do this with multiple dispatch. Do this for + and * and remove the ModelStates thing?
        Expr(:block, ( :( stackem!(x, y.$f) ) for f in fieldnames(T) )...)
    elseif T <: Nothing # There's nothing to stack.
        :(x) # :o
    else
        error("Sorry, I don't know how to stack a ", T, ".")
    end
end

# Stack into a preallocated vector, starting at k+1. Returns k, the number elements written.
@generated function stackem!(x, y::T, k) where {T}
    if T <: Real || T <: Enum
        quote
            x[k+1] = y
            k+1
        end
    elseif T <: Array
        if eltype(T) <: Real
            quote
                x[k+1:k+length(y)] = y[:]
                k+length(y)
            end
        else
            quote
                for c = 1:length(y)
                    k = stackem!(x, y[c], k)
                end
            end
        end
    elseif !isempty(fieldnames(T)) # This is nice; we can't do this with multiple dispatch. Do this for + and * and remove the ModelStates thing?
        Expr(:block, ( :( k = stackem!(x, y.$f, k) ) for f in fieldnames(T) )...)
    elseif T <: Nothing
        :(k) # Don't choke on that k, my friend.
    else
        error("Sorry, I don't know how to stack a ", T, ".")
    end
end


text:
using Test, LocalStore

struct A
    x::Int
end

Base.hash(a::A, h::UInt) = hash(A.name, hash(a.x, h))

LocalStore.save(a::A, path::String) = open("$path/file.txt", "w") do f
    write(f, string(a.x^2))
end

LocalStore.load(a::A, path::String) = open("$path/file.txt", "r") do f
    parse(Int, readline(f))
end

a = LocalStore.load(A(5))
@test LocalStore.issaved(A(5))
@test a == 25
@test LocalStore.load(A(5)) == 25

LocalStore.clean(A(5))
@test !LocalStore.issaved(A(5))


text:
using JW3126PkgDev
using Test

@testset "JW3126PkgDev.jl" begin
    # Write your own tests here.
end


text:
# https://github.com/JuliaCloud/AWSCore.jl/commit/b719614b812fbec385833e607a2f8772d6009b59
import Base: copy!
Base.@deprecate copy!(dest::AWSCredentials, src::AWSCredentials) copyto!(dest, src) false


text:
@testset "AWSConfig" begin
    @testset "default profile assumes role" begin
        access_key_id = "assumed_access_key_id"
        config_dir = joinpath(@__DIR__, "configs", "default-role")

        # Avoid calling out to STS with invalid credentials
        patch = Patches._assume_role_patch("AssumeRole"; access_key=access_key_id)

        config = withenv(
            "AWS_CONFIG_FILE" => joinpath(config_dir, "config"),
            "AWS_SHARED_CREDENTIALS_FILE" => joinpath(config_dir, "credentials"),
            "AWS_ACCESS_KEY_ID" => nothing,
            "AWS_SECRET_ACCESS_KEY" => nothing,
        ) do
            apply(patch) do
                AWSConfig(; profile="default")
            end
        end

        @test config.credentials.access_key_id == access_key_id
    end

    @testset "default profile section names" begin
        allowed_default_sections = ["default", "profile default"]
        mktemp() do config_path, _
            for default_section_str in allowed_default_sections
                config = """
                [$default_section_str]
                region = xx-yy-1
                """
                write(config_path, config)
                @test aws_get_region(; profile="default", config=config_path) == "xx-yy-1"
            end
        end
    end
end


text:
using SimpleElastixRegistration
using Test

@testset "SimpleElastixRegistration.jl" begin
    # Write your tests here.
end


text:
# Type Hierarchy Exercise
#
# Demonstration of multiple dispatch as opposed to conditional formatting

#module LJtypeHierarchy

abstract type AbstractPerson end

struct Person <: AbstractPerson
    name::String
end

struct Student <: AbstractPerson
    name::String
    grade::Int
end

struct GroupLeader <: AbstractPerson
    name::String
    group::String
end

#end

println("Hello World")

adam = Person("adam")
println("$(adam.name) is a Person")
cameron = Student("cameron", 95)
john = GroupLeader("John", "CDAS")

person_info(p::Person) = println("$(p.name) is a Person")
person_info(p::Student) = println("$(p.name) is a Student with grade $(p.grade)")
person_info(p::GroupLeader) = println("$(p.name) is a GroupLeader of $(p.group)")

people = (adam, cameron, john)

for person in people
    person_info(person)
end

max_stars_repo_name: AdamVStephen/learning-julia

max_issues_repo_name: AdamVStephen/learning-julia

max_forks_repo_name: AdamVStephen/learning-julia


text:
"""
Beta Uniform Mixture (BUM) Model

# Arguments

- π0 : Contributing fraction of the uniform distribution to the full model
- α, β : Parameters of the Beta distribution, Float64, default: 0.5, 3.0

# Return values

`MixtureModel`, as defined in the `Distributions` package, composed of

- a uniform distribution in the interval [0, 1], with weight/prior π₀
- a Beta distribution with parameters α and β, with weight/prior 1-π₀

# Examples

```jldoctest
julia> bum = BetaUniformMixtureModel(0.2, 0.5, 1.0);

julia> using Distributions

julia> pdf.(bum, 0.2:0.2:1.0)
5-element Array{Float64,1}:
 1.09443
 0.832456
 0.716398
 0.647214
 0.6

```

# References

Pounds, S., and Morris, S.W. (2003). Estimating the occurrence of false
positives and false negatives in microarray studies by approximating and
partitioning the empirical distribution of p-values. Bioinformatics 19,
1236–1242.

"""
function BetaUniformMixtureModel(π0::AbstractFloat, α::AbstractFloat = 0.5, β::AbstractFloat = 1.0)
    isin(π0, 0, 1) || throw(DomainError("π0 must be in [0, 1]"))
    bum = MixtureModel([Beta(α, β), Uniform()], [1-π0, π0])
    return bum
end


text:
using LsqFit, Calculus

# a simple wrapper for Levenberg-Marquardt with finite differencing for the Jacobian
lmdif(f, x0; show_trace::Bool=false) =
 LsqFit.levenberg_marquardt(f, jacobian(f), x0, show_trace=show_trace)

# custom MEX-like wrapper
using MATLAB: jvalue
lmdif(args::Vector{MATLAB.MxArray}) =
 lmdif(jvalue(args[1]), jvalue(args[2]), show_trace=(length(args)>2 ? jvalue(args[3]) : false))

max_stars_repo_name: taylormcd/mexjulia


text:
#using Pkg
#pkg"activate .."
#push!(LOAD_PATH,"../src/")

using MPT, Documenter

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

makedocs(;
    modules=[MPT],
    authors="Kiar Fatah",
    repo="https://github.com/Xiar-fatah/MPT.jl/blob/{commit}{path}#{line}",
    sitename="MPT.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Xiar-fatah.github.io/MPT.jl",
        assets=String[],
    ),
    pages=[
        #"Home" => "index.md",
            "index.md",
            "portfolio.md",
            "backtesting.md"
        ]
)

deploydocs(;
    repo="github.com/Xiar-fatah/MPT.jl",
    versions = ["stable" => "master", "dev" => "master"]
    #devbranch="master",
)


text:
function initialize_outer!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]

    Vr_cnv = get_inner_vars(dynamic_device)[Vr_cnv_var]
    Vi_cnv = get_inner_vars(dynamic_device)[Vi_cnv_var]
    θ0_oc = angle(Vr_cnv + 1im * Vi_cnv)

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Update inner_vars
    get_inner_vars(dynamic_device)[P_ES_var] = p_elec_out
    #Update states
    outer_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )
    outer_states = @view device_states[outer_ix]
    outer_states[1] = PSY.get_ω_ref(dynamic_device) #ω
    outer_states[2] = θ0_oc #θ_oc
    outer_states[3] = q_elec_out #qm

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ0_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = PSY.get_ω_ref(dynamic_device)
    #Update Q_ref. Initialization assumes q_ref = q_elec_out of PF solution
    PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = p_elec_out
    PSY.set_P_ref!(PSY.get_active_power(PSY.get_outer_control(dynamic_device)), p_elec_out)
    PSY.get_ext(dynamic_device)[CONTROL_REFS][Q_ref_index] = q_elec_out
end

function initialize_outer!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]

    Vr_cnv = get_inner_vars(dynamic_device)[Vr_cnv_var]
    Vi_cnv = get_inner_vars(dynamic_device)[Vi_cnv_var]
    θ0_oc = angle(Vr_cnv + 1im * Vi_cnv)

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Update inner_vars
    get_inner_vars(dynamic_device)[P_ES_var] = p_elec_out
    #Update states
    outer_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerDroop, PSY.ReactivePowerDroop},
    )
    outer_states = @view device_states[outer_ix]
    outer_states[1] = θ0_oc #θ_oc
    outer_states[2] = p_elec_out #pm
    outer_states[3] = q_elec_out #qm

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ0_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = PSY.get_ω_ref(dynamic_device)
    #Update Q_ref. Initialization assumes q_ref = q_elec_out of PF solution
    PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = p_elec_out
    PSY.set_P_ref!(PSY.get_active_power(PSY.get_outer_control(dynamic_device)), p_elec_out)
    PSY.get_ext(dynamic_device)[CONTROL_REFS][Q_ref_index] = q_elec_out
end

function initialize_outer!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
    )
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]
    Ir_filter = device_states[external_ix[3]]
    Ii_filter = device_states[external_ix[4]]
    Ir_cnv = device_states[external_ix[5]]
    Ii_cnv = device_states[external_ix[6]]

    #Obtain additional expressions
    θ0_oc = get_inner_vars(dynamic_device)[θ_freq_estimator_var]
    I_dq_cnv = ri_dq(θ0_oc + pi / 2) * [Ir_cnv; Ii_cnv]
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Get Outer Controller parameters
    outer_control = PSY.get_outer_control(dynamic_device)
    active_power_control = PSY.get_active_power(outer_control)
    Ki_p = PSY.get_Ki_p(active_power_control) #Integral Gain
    reactive_power_control = PSY.get_reactive_power(outer_control)
    Ki_q = PSY.get_Ki_q(reactive_power_control) #Integral Gain

    #Update inner_vars
    get_inner_vars(dynamic_device)[P_ES_var] = p_elec_out
    #Update states
    outer_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.ActivePowerPI, PSY.ReactivePowerPI},
    )
    outer_states = @view device_states[outer_ix]
    outer_states[1] = I_dq_cnv[q] / Ki_p #σp_oc
    outer_states[2] = p_elec_out #p_oc
    outer_states[3] = I_dq_cnv[d] / Ki_q #σq_oc
    outer_states[4] = q_elec_out #q_oc

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ0_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = PSY.get_ω_ref(dynamic_device)
    get_inner_vars(dynamic_device)[Id_oc_var] = I_dq_cnv[d]
    get_inner_vars(dynamic_device)[Iq_oc_var] = I_dq_cnv[q]
    #Update Q_ref. Initialization assumes q_ref = q_elec_out from PF solution
    PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = p_elec_out
    PSY.set_P_ref!(PSY.get_active_power(PSY.get_outer_control(dynamic_device)), p_elec_out)
    PSY.get_ext(dynamic_device)[CONTROL_REFS][Q_ref_index] = q_elec_out
    PSY.set_Q_ref!(
        PSY.get_reactive_power(PSY.get_outer_control(dynamic_device)),
        q_elec_out,
    )
end


text:
using HDF5

function generate_files(path,N)
    hdf5_file = path*"unique.hdf5"
    csv_file = path*"unique.csv"
    if isfile(hdf5_file)
        rm(hdf5_file)
    end
    if isfile(csv_file)
        rm(csv_file)
    end
    # uniform ids in a range for load balance among nodes
    id = rand(1:256, N)
    x = rand(1:256*1000,N)
    h5write(hdf5_file,"/id", id)
    h5write(hdf5_file,"/x", x)
    writecsv(csv_file, zip(id,x))
end

N = 2*10^6
path = ENV["HOME"]*"/tmp/"
generate_files(path, N)

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
module Test_gsklf

using LinearAlgebra
using MatrixPencils
using Test

@testset "gsklf" begin

fast = false; Ty = Float64; #Ty = Complex{Float64}

for fast in (true, false)
# test example for SLICOT subroutine AB08ND
A = [
   1.0   0.0   0.0   0.0   0.0   0.0
   0.0   1.0   0.0   0.0   0.0   0.0
   0.0   0.0   3.0   0.0   0.0   0.0
   0.0   0.0   0.0  -4.0   0.0   0.0
   0.0   0.0   0.0   0.0  -1.0   0.0
   0.0   0.0   0.0   0.0   0.0   3.0];
E = I;
B = [
   0.0  -1.0
  -1.0   0.0
   1.0  -1.0
   0.0   0.0
   0.0   1.0
  -1.0  -1.0];
C = [
   1.0   0.0   0.0   1.0   0.0   0.0
   0.0   1.0   0.0   1.0   0.0   1.0
   0.0   0.0   1.0   0.0   0.0   1.0];
D = [
   0.0   0.0
   0.0   0.0
   0.0   0.0];
n = 6; m = 2; p = 3;

M2 = [A B; C D]; N2 = [E zeros(n,m); zeros(p,n+m)];
zerref, ziref, = pzeros(M2,N2)

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "none", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (4, 2, 2, 0) && ismissing(nmsz) && niz == 2 &&
      zer == Ty[] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (3, 3, 2, 0) && nmsz == 0 && niz == 2 &&
      zer ≈ [2.0] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (1, 5, 2, 0) && nmsz == 1 && niz == 0 &&
      zer ≈ [2.0, Inf, Inf] && zi == [1, 1]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "all", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (0, 6, 2, 0) && ismissing(nmsz) && ismissing(niz) &&
      sort(zer) ≈ sort(zerref) && zi == ziref

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "infinite", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (2, 4, 2, 0) && ismissing(nmsz) && niz == 2 &&
      zer == [Inf, Inf] && zi == [1,1]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "stable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (3, 3, 2, 0) && nmsz == 0 && niz == 2 &&
      zer ≈ [-1] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "stable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (4, 2, 2, 0) && nmsz == 1 && niz == 0 &&
      zer == Ty[] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "s-unstable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (1, 5, 2, 0) && nmsz == 0 && niz == 2 &&
      zer ≈ [2, Inf, Inf] && zi == [1, 1]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "s-unstable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (0, 6, 2, 0) && nmsz == 1 && niz == 0 &&
      sort(zer) ≈ sort(zerref) && zi == ziref

# test example for SLICOT subroutine AB08ND
A = [
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E = [
    0     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B = [
    -1     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C = [
    0     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];
D = [
    1     2    -2
    0    -1    -2
    0     0     0];
n = 9; m = 3; p = 3;
M2 = [A B; C D]; N2 = [E zeros(n,m); zeros(p,n+m)];
zerref, ziref, = pzeros(M2,N2)

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "none", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (6, 1, 2, 3) && ismissing(nmsz) && niz == 2 &&
      zer == Ty[] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (5, 2, 2, 3) && nmsz == 0 && niz == 2 &&
      zer ≈ [1] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (4, 3, 2, 3) && nmsz == 1 && niz == 0 &&
      zer ≈ [Inf, Inf] && zi == [2]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "all", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (3, 4, 2, 3) && ismissing(nmsz) && ismissing(niz) &&
      sort(zer) ≈ sort(zerref) && zi == ziref

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "infinite", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (4, 3, 2, 3) && ismissing(nmsz) && niz == 2 &&
      zer ≈ [Inf, Inf] && zi == [2]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "finite", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (5, 2, 2, 3) && ismissing(nmsz) && ismissing(niz) &&
      zer ≈ [1] && zi == []

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "stable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (6, 1, 2, 3) && nmsz == 0 && niz == 2 &&
      zer ≈ Ty[] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "stable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (6, 1, 2, 3) && nmsz == 1 && niz == 0 &&
      zer == Ty[] && zi == Int[]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "s-unstable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (3, 4, 2, 3) && nmsz == 0 && niz == 2 &&
      zer ≈ [1, Inf, Inf] && zi == [2]

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "s-unstable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (3, 4, 2, 3) && nmsz == 1 && niz == 0 &&
      sort(zer) ≈ sort(zerref) && zi == ziref

Ty = Complex{Float64};  fast = false;
for Ty in (Float64, Complex{Float64})
n = 5; m = 4; p = 3;
A = rand(Ty,n,n); B = rand(Ty,n,m); C = rand(Ty,p,n); D = rand(Ty,p,m);
E = I;
M2 = [A B; C D]; N2 = [E zeros(n,m); zeros(p,n+m)];
zerref, ziref, = pzeros(M2,N2)

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "none", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)

@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (6, 0, 3, 0) && ismissing(nmsz) && niz == 0 &&
      zer ≈ Ty[] && zi == []

n = 9; m = 3; p = 3; rE = 6;
A = rand(Ty,n,n); B = rand(Ty,n,m); C = rand(Ty,p,n); D = rand(Ty,p,m);
E = [rand(Ty,n,rE)  zeros(Ty,n,n-rE)]*qr(rand(Ty,n,n)).Q
M2 = [A B; C D]; N2 = [E zeros(n,m); zeros(p,n+m)];
zerref, ziref, = pzeros(M2,N2)

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "none", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)

@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.)) &&
      dimsc == (6, 0, 3, 3) && ismissing(nmsz) && niz == 0 &&
      zer ≈ Ty[] && zi == []

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl,atol1=1.e-7,atol2=1.e-7)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.))

M, N, Q, Z, dimsc, nmsz, niz = gsklf(A, E, B, C, D, jobopt = "unstable", disc = true, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
nr = dimsc[1]; nl = dimsc[2]; ml = dimsc[3]; nsinf = dimsc[4];
il = n-nsinf-nl+1:n-nsinf; jal = nr+1:nr+nl; jbl = nr+nl+1:nr+nl+ml; icl = n+1:n+p
Abl = M[il,jal]; Ebl = N[il,jal]; Bbl = M[il,jbl]; Cbl = M[icl,jal]; Dbl = M[icl,jbl];
zer, zi, = spzeros(Abl,Ebl,Bbl,Cbl,Dbl)
@test norm(Q'*[A B]*Z-M[1:n,:]) < sqrt(eps(1.)) &&
      norm([C D]*Z-M[n+1:n+p,:]) < sqrt(eps(1.)) &&
      norm(Q'*[E zeros(n,m)]*Z-N[1:n,:]) < sqrt(eps(1.))
end

end
end #fast

end #module