text:
#===============================================================================

    verify_request.jl

    Author: Tom Price
    Date: Dec 2018

    This Julia script tests the JSON data structures
    that are supplied in the body of GET requests

    *     based on documentation of REST API from juliaapi_new.txt

    *     Julia listens on http://localhost:8081 via HttpServer.jl and dispatch.jl
          and sends responses to the Rails app on http://localhost:3000

    *     the relevant controller in the Rails app is
          chaipcr/web/app/controllers/experiments_controller.rb

    *     currently calls are GET-ed to
          http://127.0.0.1:8081/experiments/#{experiment.id}/standard_curve
          http://127.0.0.1:8081/experiments/#{experiment.id}/amplification
          http://127.0.0.1:8081/experiments/#{experiment.id}/meltcurve
          http://127.0.0.1:8081/experiments/#{experiment.id}/analyze

===============================================================================#

import JSON, DataStructures.OrderedDict
import FactCheck: facts, context, @fact, clear_results, exitstatus, less_than_or_equal

FactCheck.setstyle(:default)

#===============================================================================
    Here are the REST APIs using HTTP GET
===============================================================================#

#===============================================================================

    call: experiments/:experiment_id/standard_curve

===============================================================================#

function verify_request(
    ::Type{Val{QpcrAnalysis.standard_curve}},
    request ::AbstractArray
)
    facts("Standard curve requested") do
        context("Verifying request body") do
            @fact (isa(request,Vector)) --> true
            n_targets=-99
            for i in range(1,length(request))
                well=request[i]
                @fact (isa(well,OrderedDict)) --> true
                if (length(well)>0)
                    @fact (length(well)) --> 1
                    @fact (haskey(well,"well")) --> true
                    array=well["well"]
                    for j in range(1,length(array))
                        dict=array[j]
                        @fact (isa(dict,OrderedDict)) --> true
                        if (length(dict)>0)
                            # fact (length(dict)) --> 3
                            @fact (haskey(dict,"target")) --> true
                            @fact (isa(dict["target"],Integer)) --> true
                            @fact (haskey(dict,"cq")) --> true
                            @fact (isa(dict["cq"],Number)) --> true
                            @fact (haskey(dict,"quantity")) --> true
                            subdict=dict["quantity"]
                            @fact (isa(subdict,OrderedDict)) --> true
                            @fact (length(subdict)) --> 2
                            @fact (haskey(subdict,"m")) --> true
                            @fact (haskey(subdict,"b")) --> true
                            @fact (isa(subdict["m"],Number)) --> true
                            @fact (isa(subdict["b"],Number)) --> true
                        end
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/amplification

===============================================================================#

## according to juliaapi_new.txt:
#
## each set of calibration data (water, channel_1, channel_2) comes from the following SQL query:
#
# SELECT fluorescence_value, well_num, channel
#     FROM fluorescence_data
#     WHERE experiment_id = $calib_id AND step_id = $step_id
#     ORDER BY channel, well_num
# ;
#
# channel_2 will be NULL for single channel:
#     "calibration_info": {
#         "water": {
#             "fluorescence_value": [
#                 [water_1__well_01, water_1__well_02, …, water_1__well_16],
#                 null
#             ]
#         },
#         "channel_1": {
#             "fluorescence_value": [
#                 [signal_1__well_01, signal_1__well_02, …, signal_1__well_16],
#                 null
#             ]
#         }
#     }

function calibration_test(
    calib ::Associative,
    n_channels ::Integer =length(CHANNELS),
    conditions ::AbstractArray =["water","channel_1","channel_2"][range(1, n_channels+1)]
)
    n_conditions=length(conditions)
    @fact (isa(calib,OrderedDict)) --> true
    # @fact (length(calib)) --> n_conditions
    @fact (isa(calib[conditions[1]],OrderedDict))--> true
    @fact (haskey(calib[conditions[1]],"fluorescence_value")) --> true
    @fact (isa(calib[conditions[1]]["fluorescence_value"],Vector)) --> true
    n_wells=length(calib[conditions[1]]["fluorescence_value"][1])
    for condition in conditions
        @fact (haskey(calib,condition)) --> true
        @fact (isa(calib[condition],OrderedDict)) --> true
        @fact (length(calib[condition])) --> 1
        @fact (haskey(calib[condition],"fluorescence_value")) --> true
        @fact (isa(calib[condition]["fluorescence_value"],Vector)) --> true
        @fact (length(calib[condition]["fluorescence_value"])) --> less_than_or_equal(2)
        for channel in range(1,n_channels)
            @fact (isa(calib[condition]["fluorescence_value"][channel],Vector)) --> true
            @fact (length(calib[condition]["fluorescence_value"][channel])) --> n_wells
            for i in range(1,n_wells)
                @fact (isa(calib[condition]["fluorescence_value"][channel][i],Number)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

## Raw Data comes from the following SQL query:
# SELECT fluorescence_value, well_num, cycle_num, channel
#     FROM fluorescence_data
#     WHERE experiment_id = $exp_id AND step_id = $step_id
#   ORDER BY channel, well_num, cycle_num
# ;

function raw_test(raw ::Associative)
    @fact (isa(raw,OrderedDict)) --> true
    variables=["fluorescence_value", "channel", "well_num"]
    if (haskey(raw,"temperature"))
        push!(variables,"temperature")
    else
        push!(variables,"cycle_num")
    end
    n_raw=length(raw["fluorescence_value"])
    for v in variables
        @fact (haskey(raw,v)) --> true
        @fact (isa(raw[v], Vector)) --> true
        @fact abs(length(raw[v]) - n_raw) --> less_than_or_equal(1)
        for i in range(1,length(raw[v]))
            if (v=="fluorescence_value" || v=="temperature")
                @fact (isa(raw[v][i], Number)) --> true
            else
                @fact (isa(raw["well_num"][i], Integer)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

function verify_request(
    ::Type{Val{amplification}},
    request ::Associative
)
    facts("Amplification requested") do
        context("Verifying request body") do
            @fact (isa(request,OrderedDict)) --> true
            @fact (haskey(request,"experiment_id")) --> true
            @fact (isa(request["experiment_id"],Integer)) --> true
            if (haskey(request,"step_id"))
                id="step_id"
            else
                id="ramp_id"
            end
            @fact (haskey(request,id)) --> true
            @fact (isa(request[id],Integer)) --> true
            if (haskey(request,"min_reliable_cyc"))
                @fact (isa(request["min_reliable_cyc"],Integer)) --> true
            end
            if (haskey(request,"baseline_cyc_bounds"))
                @fact (isa(request["baseline_cyc_bounds"],Vector)) --> true
                if (length(request["baseline_cyc_bounds"]) > 0)
                    @fact (length(request["baseline_cyc_bounds"])) --> 2
                    @fact (isa(request["baseline_cyc_bounds"][1],Integer)) --> true
                    @fact (isa(request["baseline_cyc_bounds"][2],Integer)) --> true
                end
            end
            if (haskey(request,"baseline_method"))
                @fact (isa(request["baseline_method"],String)) --> true
                @fact (
                    request["baseline_method"] == "sigmoid" ||
                    request["baseline_method"] == "linear"  ||
                    request["baseline_method"] == "median"
                ) --> true
            end
            if (haskey(request,"cq_method"))
                @fact (isa(request["cq_method"],String)) --> true
                @fact (request["cq_method"] in ["cp_dr1","cp_dr2","Cy0","ct"]) --> true
            end
            if (haskey(request,"min_fluomax"))
                @fact (isa(request["min_fluomax"],Number)) --> true
            end
            if (haskey(request,"min_D1max"))
                @fact (isa(request["min_D1max"],Number)) --> true
            end
            if (haskey(request,"min_D2max"))
                @fact (isa(request["min_D2max"],Number)) --> true
            end
            @fact (haskey(request,"calibration_info")) --> true
            calib=request["calibration_info"]
            @fact (isa(calib,OrderedDict)) --> true
            @fact (haskey(calib,"water")) --> true
            @fact (isa(calib["water"],OrderedDict)) --> true
            @fact (haskey(calib["water"],"fluorescence_value")) --> true
            @fact (isa(calib["water"]["fluorescence_value"],Vector)) --> true
            if length(calib["water"]["fluorescence_value"]) < 2 ||
                calib["water"]["fluorescence_value"][2]==nothing
                n_channels=1
            else
                n_channels=2
            end
            @fact (haskey(request,"raw_data")) --> true
            raw=request["raw_data"]
            @fact (isa(raw,OrderedDict)) --> true
            calibration_test(calib,n_channels)
            raw_test(raw)
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/meltcurve

===============================================================================#

#=
    Notes:
    channel_nums = [1] for 1 channel, [1,2] for 2 channels, etc.
    top_N = number of Tm peaks to report
=#

# Calibration (water, channel_1, channel_2) data comes from the following SQL query:
# SELECT fluorescence_value, well_num, channel
#     FROM fluorescence_data
#     WHERE experiment_id = $calib_id AND step_id = $step_id
#     ORDER BY channel, well_num
# ;

# Raw Data comes from the following SQL query:
# SELECT fluorescence_value, temperature, well_num, channel
#     FROM melt_curve_data
#     WHERE
#         experiment_id = $exp_id AND
#         stage_id = $stage_id
#     ORDER BY channel, well_num
# ;

function verify_request(
    ::Type{Val{QpcrAnalysis.meltcurve}},
    request ::Any
)
    facts("Melting curve requested") do
        context("Verifying request body") do
            @fact (isa(request,OrderedDict)) --> true
            @fact (haskey(request,"experiment_id")) --> true
            @fact (isa(request["experiment_id"],Integer)) --> true
            @fact (haskey(request,"stage_id")) --> true
            @fact (isa(request["stage_id"],Integer)) --> true
            @fact (haskey(request,"calibration_info")) --> true
            @fact (haskey(request,"channel_nums")) --> true
            @fact (isa(request["channel_nums"],Vector)) --> true
            calib=request["calibration_info"]
            @fact (isa(calib,OrderedDict)) --> true
            @fact (haskey(calib,"water")) --> true
            @fact (isa(calib["water"],OrderedDict)) --> true
            @fact (haskey(calib["water"],"fluorescence_value")) --> true
            @fact (isa(calib["water"]["fluorescence_value"],Vector)) --> true
            if length(calib["water"]["fluorescence_value"])<2 ||
                calib["water"]["fluorescence_value"][2]==nothing
                @fact (request["channel_nums"]) --> [1]
                n_channels=1
            else
                @fact (request["channel_nums"]) --> [1,2]
                n_channels=2
            end
            if (haskey(request,"qt_prob"))
                @fact (isa(request["qt_prob"],Number)) --> true
            end
            if (haskey(request,"max_normd_qtv"))
                @fact (isa(request["max_normd_qtv"],Number)) --> true
            end
            if (haskey(request,"top_N"))
                @fact (isa(request["top_N"],Integer)) --> true
            end
            @fact (haskey(request,"raw_data")) --> true
            raw=request["raw_data"]
            calibration_test(calib,n_channels)
            raw_test(raw)
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/thermal_performance_diagnostic

===============================================================================#

# SQL query:
# SELECT *
#     FROM temperature_logs
#     WHERE experiment_id = $exp_id
#     ORDER BY id
# ;

function verify_request(
    ::Type{Val{QpcrAnalysis.thermal_performance_diagnostic}},
    request ::Any
)
    facts("Thermal performance diagnostic requested") do
        context("Verifying request body") do
            variables=["lid_temp","heat_block_zone_1_temp","heat_block_zone_2_temp","elapsed_time"] # ,"cycle_num"]
            @fact (isa(request,OrderedDict)) --> true
            @fact (length(request)) --> length(variables)
            n_cycles=length(request["elapsed_time"])
            for v in variables
                @fact (haskey(request,v)) --> true
                @fact (length(request[v])) --> n_cycles
                for i in range(1,n_cycles)
                    if (v=="elapsed_time")
                        @fact (isa(request[v][i],Integer)) --> true
                    else
                        @fact (isa(request[v][i],Number)) --> true
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/thermal_consistency

===============================================================================#

# SQL query:
# SELECT fluorescence_value, temperature, well_num, channel
#     FROM melt_curve_data
#     WHERE
#         experiment_id = $exp_id AND
#         stage_id = $stage_id
#   ORDER BY channel, well_num
# ;

function verify_request(
    ::Type{Val{QpcrAnalysis.thermal_consistency}},
    request ::Any
)
    facts("Thermal consistency requested") do
        context("Verifying request body") do
            @fact (isa(request,OrderedDict)) --> true
            # @fact (length(request)) --> 2
            @fact (haskey(request,"raw_data")) --> true
            @fact (isa(request["raw_data"],OrderedDict)) --> true
            @fact (haskey(request["raw_data"],"channel")) --> true
            @fact (haskey(request,"calibration_info")) --> true
            calib=request["calibration_info"]
            @fact (isa(calib,OrderedDict)) --> true
            @fact (haskey(calib,"water")) --> true
            @fact (isa(calib["water"],OrderedDict)) --> true
            @fact (haskey(calib["water"],"fluorescence_value")) --> true
            @fact (isa(calib["water"]["fluorescence_value"],Vector)) --> true
            if length(calib["water"]["fluorescence_value"]) < 2 ||
                calib["water"]["fluorescence_value"][2]==nothing
                n_channels=1
            else
                n_channels=2
            end
            calibration_test(calib,n_channels)
            raw_test(request["raw_data"])
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_cal

===============================================================================#

function verify_request(
    ::Type{Val{QpcrAnalysis.optical_cal}},
    request ::Any
)
    facts("Optical calibration requested") do
        context("Verifying request body") do
            @fact (isa(request,OrderedDict)) --> true
            @fact (haskey(request,"calibration_info")) --> true
            calib=request["calibration_info"]
            @fact (isa(calib,OrderedDict)) --> true
            @fact (haskey(calib,"water")) --> true
            @fact (isa(calib["water"],OrderedDict)) --> true
            @fact (haskey(calib["water"],"fluorescence_value")) --> true
            @fact (isa(calib["water"]["fluorescence_value"],Vector)) --> true
            if length(calib["water"]["fluorescence_value"]) < 2 ||
                calib["water"]["fluorescence_value"][2]==nothing
                n_channels=1
            else
                n_channels=2
            end
            calibration_test(calib,n_channels)
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_test_single_channel

===============================================================================#

## SQL query:
# SELECT fluorescence_value, well_num, cycle_num
#     FROM fluorescence_data
#     WHERE experiment_id = $exp_id AND step_id = $step_id
#     ORDER BY well_num, cycle_num
# ;

function verify_request(
    ::Type{Val{QpcrAnalysis.optical_test_single_channel}},
    request ::Any
)
    facts("Single channel optical test requested") do
        context("Verifying request body") do
            conditions=["baseline","excitation"]
            @fact (isa(request,OrderedDict)) --> true
            # @fact (length(request)) --> length(conditions)
            @fact (isa(request["baseline"],OrderedDict))--> true
            @fact (haskey(request["baseline"],"fluorescence_value")) --> true
            @fact (isa(request["baseline"]["fluorescence_value"],Vector)) --> true
            n_wells=length(request["baseline"]["fluorescence_value"])
            for condition in conditions
                @fact (haskey(request,condition)) --> true
                @fact (isa(request[condition],OrderedDict)) --> true
                @fact (length(request[condition])) --> 1
                @fact (haskey(request[condition],"fluorescence_value")) --> true
                @fact (isa(request[condition]["fluorescence_value"],Vector)) --> true
                @fact (isa(request[condition]["fluorescence_value"],Vector)) --> true
                @fact (length(request[condition]["fluorescence_value"])) --> n_wells
                for i in range(1,n_wells)
                    @fact (isa(request[condition]["fluorescence_value"][i],Number)) --> true
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_test_dual_channel

===============================================================================#

## SQL query:
# SELECT fluorescence_value, well_num, cycle_num
#     FROM fluorescence_data
#     WHERE experiment_id = $exp_id AND step_id = $step_id
#     ORDER BY well_num, cycle_num
# ;

function verify_request(
    ::Type{Val{QpcrAnalysis.optical_test_dual_channel}},
    request ::Any
)
    facts("Dual channel optical test requested") do
        context("Verifying request body") do
            @fact (isa(request,OrderedDict)) --> true
            if (haskey(request,"channel_1"))
                calibration_test(request,2,["baseline","water","channel_1","channel_2"])
            else
                calibration_test(request,2,["baseline","water","FAM","HEX"])
            end
        end
    end
    FactCheck.exitstatus()
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#=
findpeaks.jl
=#

using Sound: findpeaks
using Test: @testset, @test

@testset "findpeaks" begin
    N = 32
    i = [7, 12, 25]
    M = length(i)
    x = zeros(N)
    x[i] .= [2, 5, 3]
    p = findpeaks(x)
    @test p[:,1] == [2; 2 .+ i[1:(M-1)]]
    @test p[:,2] == i
    @test p[:,3] == [1 .+ i[1:(M-1)]; N]
end

ext: jl


text:
function _raw_setindex(A::nmod_mat, i::Int, j::Int, x::UInt)
  ccall((:nmod_mat_set_entry, libflint), Nothing, (Ref{nmod_mat}, Int, Int, UInt), A, i - 1, j - 1, x)
end

###############################################################################
#
#  Howell form for Generic.Mat{fmpz}
#
###############################################################################

if Nemo.version() > v"0.15.1"
  function howell_form(A::Generic.Mat{Nemo.fmpz_mod})
    local B::fmpz_mat
    if nrows(A) < ncols(A)
      B = vcat(lift(A), zero_matrix(FlintZZ, ncols(A)-nrows(A), ncols(A)))
    else
      B = lift(A)
    end
    R = base_ring(A)
    ccall((:fmpz_mat_howell_form_mod, libflint), Nothing,
                (Ref{fmpz_mat}, Ref{fmpz}), B, modulus(R))
    return change_base_ring(B, R)
  end

  #
  #  for the in-place function, the number of rows must be at least equal to the number of columns
  #
  function howell_form!(A::Generic.Mat{Nemo.fmpz_mod})

    R = base_ring(A)
    A1 = lift(A)
    ccall((:fmpz_mat_howell_form_mod, libflint), Nothing,
                  (Ref{fmpz_mat}, Ref{fmpz}), A1, modulus(R))
    for i in 1:nrows(A)
      for j in 1:ncols(A)
        A[i, j] = A1[i, j]
      end
    end
    return A
  end

  function triangularize!(A::Generic.Mat{Nemo.fmpz_mod})
    R=base_ring(A)
    n=R.modulus

    #
    #  Get an upper triangular matrix
    #

    for j=1:ncols(A)
      for i=j+1:ncols(A)
        g,s,t,u,v = _xxgcd(A[j,j].data,A[i,j].data,n)
        for k in 1:ncols(A)
          t1 = s* A[j,k] + t* A[i,k]
          t2 = u* A[j,k] + v* A[i,k]
          A[j,k] = t1
          A[i,k] = t2
        end
      end
    end
  end

  function triangularize(A::Generic.Mat{Nemo.fmpz_mod})
    B= triangularize!(deepcopy(A))
    return B
  end
else
  function howell_form(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
    local B::fmpz_mat
    if nrows(A) < ncols(A)
      B = vcat(lift(A), zero_matrix(FlintZZ, ncols(A)-nrows(A), ncols(A)))
    else
      B = lift(A)
    end
    R = base_ring(A)
    ccall((:fmpz_mat_howell_form_mod, libflint), Nothing,
          (Ref{fmpz_mat}, Ref{fmpz}), B, modulus(R))
    return change_base_ring(B, R)
  end

  #
  #  for the in-place function, the number of rows must be at least equal to the number of columns
  #
  function howell_form!(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})

    R = base_ring(A)
    A1 = lift(A)
    ccall((:fmpz_mat_howell_form_mod, libflint), Nothing,
          (Ref{fmpz_mat}, Ref{fmpz}), A1, modulus(R))
    for i in 1:nrows(A)
      for j in 1:ncols(A)
        A[i, j] = A1[i, j]
      end
    end
    return A
  end

  function triangularize!(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
    R=base_ring(A)
    n=R.modulus

    #
    #  Get an upper triangular matrix
    #

    for j=1:ncols(A)
      for i=j+1:ncols(A)
        g,s,t,u,v = _xxgcd(A[j,j].data,A[i,j].data,n)
        for k in 1:ncols(A)
          t1 = s* A[j,k] + t* A[i,k]
          t2 = u* A[j,k] + v* A[i,k]
          A[j,k] = t1
          A[i,k] = t2
        end
      end
    end
  end

  function triangularize(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
    B= triangularize!(deepcopy(A))
    return B
  end
end


text:
# functions to read an MP3 file, uses mpg123 under the hood
include("mpg123.jl")

function initialize_readers()
    # initialize mpg123; this needs to be done only once
    mpg123_init()
end

mutable struct MP3FileSource{T} <: SampleSource
    path::AbstractString
    mpg123::MPG123
    info::MP3INFO
    pos::Int64
    readbuf::Array{T, 2}
end

function MP3FileSource(path::AbstractString, mpg123::MPG123, info::MP3INFO, bufsize::Integer)
    readbuf = Array{info.datatype}(undef, info.nchannels, bufsize)
    MP3FileSource(path, mpg123, info, Int64(0), readbuf)
end

@inline nchannels(source::MP3FileSource) = Int(source.info.nchannels)
@inline samplerate(source::MP3FileSource) = source.info.samplerate
@inline nframes(source::MP3FileSource) = source.info.nframes
@inline Base.eltype(source::MP3FileSource{T}) where {T}= T

"""convert mpg123 encoding to julia datatype"""
function encoding_to_type(encoding)
    mapping = Dict(
       MPG123_ENC_SIGNED_16 => PCM16Sample,
       # TODO: support more
    )

    encoding in keys(mapping) || error("Unsupported encoding $encoding")
    mapping[encoding]
end

"""
loads an MP3 file as SampledSignals.SampleBuf.

# Arguments
* `file::File{format"MP3"}`: the MP3 file to open
* `blocksize::Int`: the size of block to read from the disk at one time.
                   defaults to the outblock size of the MP3 file.
"""
function load(file::File{format"MP3"}; blocksize = -1)
    source = loadstream(file; blocksize = blocksize)
    buffer = try
        read(source)
    finally
        close(source)
    end
    buffer
end

@inline loadstream(path::AbstractString, args...; kwargs...) =
    loadstream(query(path), args...; kwargs...)

function loadstream(f::Function, args...; kwargs...)
    str = loadstream(args...; kwargs...)
    try
        f(str)
    finally
        close(str)
    end
end

function loadstream(path::File{format"MP3"}; blocksize = -1)
    mpg123 = mpg123_new()
    mpg123_open(mpg123, path.filename)
    nframes = mpg123_length(mpg123)
    samplerate, nchannels, encoding = mpg123_getformat(mpg123)
    if blocksize < 0
        blocksize = mpg123_outblock(mpg123)
    end
    datatype = encoding_to_type(encoding)
    encsize = sizeof(datatype)

    info = MP3INFO(nframes, nchannels, samplerate, datatype)
    bufsize = div(blocksize, encsize * nchannels)
    MP3FileSource(filename(path), mpg123, info, bufsize)
end

function unsafe_read!(source::MP3FileSource, buf::Array, frameoffset, framecount)
    total = min(framecount, nframes(source) - source.pos)
    nread = 0

    mpg123 = source.mpg123
    encsize = sizeof(source.info.datatype)
    readbuf = source.readbuf
    nchans = nchannels(source)

    while nread < total
        n = min(size(readbuf, 2), total - nread)
        nr = mpg123_read!(mpg123, readbuf, n * encsize * nchans)
        nr = div(nr, encsize * nchans)

        start_ind = 1 + nread+frameoffset
        buf[range(start_ind, stop=start_ind+nr-1), :] = transpose(readbuf[:, 1:nr])

        source.pos += nr
        nread += nr
        nr == n || break
    end

    nread
end

# @inline function Base.read(source::MP3FileSource)
#     read(source, nframes(source) - source.pos)
# end

@inline function Base.close(source::MP3FileSource)
    mpg123_close(source.mpg123)
    mpg123_delete(source.mpg123)
end


text:
module SVM
    using Pluto

    run_notebook() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/svm.jl"))
end


text:
Base.iterate(ic::Core.Compiler.IncrementalCompact) = Core.Compiler.iterate(ic)
Base.iterate(ic::Core.Compiler.IncrementalCompact, st) = Core.Compiler.iterate(ic, st)
Base.getindex(ic::Core.Compiler.IncrementalCompact, idx) = Core.Compiler.getindex(ic, idx)
Base.setindex!(ic::Core.Compiler.IncrementalCompact, v, idx) = Core.Compiler.setindex!(ic, v, idx)

Base.getindex(ic::Core.Compiler.Instruction, idx) = Core.Compiler.getindex(ic, idx)
Base.setindex!(ic::Core.Compiler.Instruction, v, idx) = Core.Compiler.setindex!(ic, v, idx)

Base.getindex(ir::Core.Compiler.IRCode, idx) = Core.Compiler.getindex(ir, idx)
Base.setindex!(ir::Core.Compiler.IRCode, v, idx) = Core.Compiler.setindex!(ir, v, idx)

Base.getindex(ref::UseRef) = Core.Compiler.getindex(ref)
Base.iterate(uses::UseRefIterator) = Core.Compiler.iterate(uses)
Base.iterate(uses::UseRefIterator, st) = Core.Compiler.iterate(uses, st)

Base.iterate(p::Core.Compiler.Pair) = Core.Compiler.iterate(p)
Base.iterate(p::Core.Compiler.Pair, st) = Core.Compiler.iterate(p, st)

Base.getindex(m::Core.Compiler.MethodLookupResult, idx::Int) = Core.Compiler.getindex(m, idx)

function Core.Compiler.optimize(interp::YaoInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    nargs = Int(opt.nargs) - 1
    @timeit "optimizer" ir = Core.Compiler.run_passes(opt.src, nargs, opt)

    # make sure all const are inlined
    # Julia itself may not inline all
    # the const values we want, e.g gates
    ir = inline_const!(ir)
    ir = compact!(ir)
    ir = elim_map_check!(ir)
    ir = compact!(ir)
    # group quantum statements so we can work on
    # larger quantum circuits before we start optimizations
    ir = group_quantum_stmts!(ir)

    # run quantum passes
    if !isempty(interp.passes)
        ir = convert_to_yaoir(ir)

        if :zx in interp.passes
            ir = run_zx_passes(ir)::YaoIR
        end

        ir = ir.ir
    end

    ir = compact!(ir)
    Core.Compiler.finish(opt, params, ir, result)
end

function group_quantum_stmts_perm(ir::IRCode)
    perms = Int[]
    cstmts_tape = Int[]
    qstmts_tape = Int[]

    for b in ir.cfg.blocks
        for v in b.stmts
            e = ir.stmts[v][:inst]
            if is_quantum_statement(e)
                if quantum_stmt_type(e) in [:measure, :barrier]
                    exit_block!(perms, cstmts_tape, qstmts_tape)
                    push!(perms, v)
                else
                    push!(qstmts_tape, v)
                end
            elseif e isa Core.ReturnNode || e isa Core.GotoIfNot || e isa Core.GotoNode
                exit_block!(perms, cstmts_tape, qstmts_tape)
                push!(cstmts_tape, v)
            elseif e isa Expr && e.head === :enter
                exit_block!(perms, cstmts_tape, qstmts_tape)
                push!(cstmts_tape, v)
            else
                push!(cstmts_tape, v)
            end
        end

        exit_block!(perms, cstmts_tape, qstmts_tape)
    end

    append!(perms, cstmts_tape)
    append!(perms, qstmts_tape)

    return perms # permute_stmts(ci, perms)
end

function permute_stmts!(stmt::InstructionStream, perm::Vector{Int})
    inst = []

    for v in perm
        e = stmt.inst[v]

        if e isa Expr
            ex = replace_from_perm(e, perm)
            push!(inst, ex)
        elseif e isa Core.GotoIfNot
            if e.cond isa Core.SSAValue
                cond = Core.SSAValue(findfirst(isequal(e.cond.id), perm))
            else
                # TODO: figure out which case is this
                # and maybe apply permute to this
                cond = e.cond
            end

            dest = findfirst(isequal(e.dest), perm)
            push!(inst, Core.GotoIfNot(cond, dest))
        elseif e isa Core.GotoNode
            push!(inst, Core.GotoNode(findfirst(isequal(e.label), perm)))
        elseif e isa Core.ReturnNode
            if isdefined(e, :val) && e.val isa Core.SSAValue
                push!(inst, Core.ReturnNode(Core.SSAValue(findfirst(isequal(e.val.id), perm))))
            else
                push!(inst, e)
            end
        else
            # RL: I think
            # other nodes won't contain SSAValue
            # let's just ignore them, but if we
            # find any we can add them here
            push!(inst, e)
            # if e isa Core.SlotNumber
            #     push!(inst, e)
            # elseif e isa Core.NewvarNode
            #     push!(inst, e)
            # else
            # end
            # error("unrecognized statement $e :: ($(typeof(e)))")
        end
    end

    copyto!(stmt.inst, inst)
    permute!(stmt.flag, perm)
    permute!(stmt.line, perm)
    permute!(stmt.type, perm)
    permute!(stmt.flag, perm)
    return stmt
end

function replace_from_perm(stmt, perm)
    stmt isa Core.SSAValue && return Core.SSAValue(findfirst(isequal(stmt.id), perm))

    if stmt isa Expr
        return Expr(stmt.head, map(x -> replace_from_perm(x, perm), stmt.args)...)
    else
        return stmt
    end
end

function exit_block!(perms::Vector, cstmts_tape::Vector, qstmts_tape::Vector)
    append!(perms, cstmts_tape)
    append!(perms, qstmts_tape)
    empty!(cstmts_tape)
    empty!(qstmts_tape)
    return perms
end

function group_quantum_stmts!(ir::IRCode)
    perm = group_quantum_stmts_perm(ir)
    permute_stmts!(ir.stmts, perm)
    return ir
end

# NOTE: this perform simple constant propagation
# inside basic blocks to get better format of
# the quantum statements

# foce inline all constants
function is_arg_allconst(args)
    for arg in args
        if arg isa SSAValue || arg isa Argument
            return false
        elseif !is_inlineable_constant(arg) && !isa(arg, QuoteNode)
            return false
        end
    end
    return true
end

function inline_const!(ir::IRCode)
    ssa_const = Any[nothing for _ in 1:length(ir.stmts)]

    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        stmt isa Expr || continue

        stmt.args = map(stmt.args) do x
            if x isa SSAValue && !(ssa_const[x.id] === nothing)
                return quoted(ssa_const[x.id])
            else
                return x
            end
        end

        if stmt.head === :call
            sig = Core.Compiler.call_sig(ir, stmt)
            f, ft, atypes = sig.f, sig.ft, sig.atypes
            allconst = true
            for atype in sig.atypes
                if !isa(atype, Const)
                    allconst = false
                    break
                end
            end

            if allconst &&
               isa(f, Core.IntrinsicFunction) &&
               is_pure_intrinsic_infer(f) &&
               intrinsic_nothrow(f, atypes[2:end])
                fargs = anymap(x::Const -> x.val, atypes[2:end])
                val = f(fargs...)
                ir.stmts[i][:inst] = quoted(val)
                ssa_const[i] = val
            elseif allconst && isa(f, Core.Builtin) && f === Core.tuple
                fargs = anymap(x::Const -> x.val, atypes[2:end])
                val = (fargs...,)
                ir.stmts[i][:inst] = quoted(val)
                ssa_const[i] = val
            end
        elseif stmt.head === :new
            exargs = stmt.args[2:end]
            allconst = is_arg_allconst(exargs)
            t = stmt.args[1]
            if allconst && isconcretetype(t) && !t.mutable
                args = anymap(exargs) do x
                    if x isa QuoteNode
                        return x.value
                    else
                        return x
                    end
                end
                val = ccall(:jl_new_structv, Any, (Any, Ptr{Cvoid}, UInt32), t, args, length(args))

                ir.stmts[i][:inst] = quoted(val)
                ssa_const[i] = val
            end
        end
    end
    return ir
end

function elim_map_check!(ir::IRCode)
    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        stmt isa Expr || continue

        if stmt.head === :invoke && stmt.args[2] === GlobalRef(YaoCompiler, :map_check)
            exargs = stmt.args[3:end]
            allconst = is_arg_allconst(exargs)

            if allconst
                args = anymap(exargs) do x
                    if x isa QuoteNode
                        return x.value
                    else
                        return x
                    end
                end

                if map_check_nothrow(args[1], args[2])
                    ir.stmts[i][:inst] = nothing
                end
            end
        end
    end
    return ir
end

struct YaoIR
    ir::IRCode
    qb::Vector{UnitRange{Int}}
end

function compute_quantum_blocks(ir::IRCode)
    quantum_blocks = UnitRange{Int}[]
    last_stmt_is_measure_or_barrier = false

    for b in ir.cfg.blocks
        start, stop = 0, 0
        for v in b.stmts
            st = ir.stmts[v][:inst]
            if is_quantum_statement(st)
                if start > 0
                    stop += 1
                else
                    start = stop = v
                end
            else
                if start > 0
                    push!(quantum_blocks, start:stop)
                    start = stop = 0
                end
            end
        end

        if start > 0
            push!(quantum_blocks, start:stop)
        end
    end
    return quantum_blocks
end

function convert_to_yaoir(ir::IRCode)
    quantum_blocks = compute_quantum_blocks(ir)
    return YaoIR(ir, quantum_blocks)
end

function map_virtual_location_expr(@nospecialize(e), regmap)
    qt = quantum_stmt_type(e)
    if qt === :gate
        loc = e.args[2]
    elseif qt === :ctrl
        loc = e.args[2]
        ctrl = e.args[3]
    elseif qt === :barrier
        e.args[2]
    elseif qt === :measure
    else
        error("unknown quantum statement: $tt")
    end
end

function map_virtual_location(@nospecialize(loc), regmap)
    if loc isa AbstractLocations
        return loc
    else
        # TODO: use a custom type
        regmap[length(keys(regmap))+1] = loc
        return Locations
    end
end

function run_zx_passes(ir::YaoIR)
    n = count_qubits(ir)
    # NOTE: we can't optimize
    # non-constant location program
    isnothing(n) && return ir
    iszero(n) && return ir

    compact = Core.Compiler.IncrementalCompact(ir.ir, true)
    for b in ir.qb
        qc = QCircuit(n)
        # if there is no quantum terminator
        # this will be the first classical
        # terminator
        first_terminator = nothing
        for v in b
            e = ir.ir.stmts[v][:inst]
            qt = quantum_stmt_type(e)
            if qt === :gate
                e.args[3] isa IntrinsicRoutine || break
                e.args[4] isa Locations || break
                # if can't convert to ZXDiagram, stop
                zx_push_gate!(qc, e.args[3], e.args[4]) || break
                # set old stmts to nothing
                compact[v] = nothing
            elseif qt === :ctrl
                e.args[3] isa IntrinsicRoutine || break
                e.args[4] isa Locations || break
                e.args[5] isa CtrlLocations || break
                zx_push_gate!(qc, e.args[3], e.args[4], e.args[5]) || break
                compact[v] = nothing
            elseif qt === :measure || qt === :barrier
                # we attach terminator after optimization
                if isnothing(first_terminator)
                    first_terminator = v
                end
            end
        end

        if isnothing(first_terminator)
            first_terminator = last(b)
        end

        zxd = ZXDiagram(qc)
        # ZX passes
        phase_teleportation(zxd)
        clifford_simplification(zxd)

        # TODO: we might want to check if
        # the result circuit is simpler indeed
        qc = QCircuit(zxd)
        for g in ZXCalculus.gates(qc)
            if g.name in (:H, :Z, :X, :S, :T, :Sdag, :Tdag)
                spec = getfield(Intrinsics, g.name)
                mi = specialize_gate(typeof(spec), Locations{Int})
                e = Expr(:invoke, mi, Semantic.gate, spec, Locations(g.loc))
            elseif g.name in (:shift, :Rz, :Rx)
                spec = getfield(Intrinsics, g.name)
                mi = specialize_gate(typeof(spec), Locations{Int})
                e = Expr(:invoke, mi, Semantic.gate, spec, Locations(g.loc))
            elseif g.name === :CNOT
                mi = specialize_ctrl(typeof(YaoCompiler.X), Locations{Int}, CtrlLocations{Int})
                e = Expr(
                    :invoke,
                    mi,
                    Semantic.ctrl,
                    YaoCompiler.X,
                    Locations(g.loc),
                    CtrlLocations(g.ctrl),
                )
            elseif g.name === :CZ
                mi = specialize_ctrl(typeof(YaoCompiler.Z), Locations{Int}, CtrlLocations{Int})
                e = Expr(
                    :invoke,
                    mi,
                    Semantic.ctrl,
                    YaoCompiler.Z,
                    Locations(g.loc),
                    CtrlLocations(g.ctrl),
                )
            else
                error("unknown gate $g")
            end
            Core.Compiler.insert_node!(compact, SSAValue(first_terminator), Nothing, e)
        end
    end
    for _ in compact
    end
    new = Core.Compiler.finish(compact)
    return YaoIR(new, compute_quantum_blocks(new))
end

function specialize_gate(spec, loc)
    atypes = Tuple{typeof(Semantic.gate),spec,loc}
    if spec <: IntrinsicRoutine
        method = first(methods(Semantic.gate, Tuple{IntrinsicRoutine,Locations}))
    elseif spec <: RoutineSpec
        method = first(methods(Semantic.gate, Tuple{RoutineSpec,Locations}))
    end

    return Core.Compiler.specialize_method(method, atypes, Core.svec())
end

function specialize_ctrl(spec, loc, ctrl)
    atypes = Tuple{typeof(Semantic.ctrl),spec,loc,ctrl}
    if spec <: IntrinsicRoutine
        method = first(methods(Semantic.ctrl, Tuple{IntrinsicRoutine,Locations,CtrlLocations}))
    elseif spec <: RoutineSpec
        method = first(methods(Semantic.ctrl, Tuple{RoutineSpec,Locations,CtrlLocations}))
    end

    return Core.Compiler.specialize_method(method, atypes, Core.svec())
end

function zx_push_gate!(qc::QCircuit, gate::IntrinsicRoutine, locs::Locations)
    name = routine_name(gate)

    # NOTE: locs can be UnitRange etc. for instruction
    # remember to expand it
    for each in locs
        if name in (:H, :X, :Z, :S, :Sdag, :T, :Tdag)
            push_gate!(qc, QGate(name, each))
        elseif name === :Y
            push_gate!(qc, QGate(:X, each))
            push_gate!(qc, QGate(:Z, each))
        elseif name in (:Rz, :Rx, :shift)
            push_gate!(qc, QGate(name, each; param = gate.variables[1]))
        elseif name === :Ry
            push_gate!(qc, QGate(:Sdag, each))
            push_gate!(qc, QGate(:Rx, each; param = gate.variables[1]))
            push_gate!(qc, QGate(:S, each))
        else
            return false
        end
    end
    return true
end

function zx_push_gate!(qc::QCircuit, gate::IntrinsicRoutine, locs::Locations, ctrl::CtrlLocations)
    name = routine_name(gate)
    ctrl = ctrl.storage.storage

    for each in locs
        if name === :Z
            push_gate!(qc, QGate(:CZ, each; ctrl = ctrl))
        elseif name === :X
            if ctrl isa Tuple && length(ctrl) == 2
                a, b = ctrl
                c = each
                push_gate!(qc, Val(:H), c)
                push_gate!(qc, Val(:CNOT), c, b)
                push_gate!(qc, Val(:Tdag), c)
                push_gate!(qc, Val(:CNOT), c, a)
                push_gate!(qc, Val(:T), c)
                push_gate!(qc, Val(:CNOT), c, b)
                push_gate!(qc, Val(:Tdag), c)
                push_gate!(qc, Val(:CNOT), c, a)
                push_gate!(qc, Val(:T), b)
                push_gate!(qc, Val(:T), c)
                push_gate!(qc, Val(:H), c)
                push_gate!(qc, Val(:CNOT), b, a)
                push_gate!(qc, Val(:T), a)
                push_gate!(qc, Val(:Tdag), b)
                push_gate!(qc, Val(:CNOT), b, a)
            elseif ctrl isa Int
                push_gate!(qc, Val(:CNOT), each, ctrl)
            end
        else
            return false
        end
    end

    return true
end


text:
using Revise; using BenchmarkTools; using DSP

#= test modulus =#
function mod2test(t, params)
    return exp.(-t)
end

#= convolution convenience function =#
function boltzconvolve(modulus, time_series, dt, prescribed_dot)
    Modulus = modulus(time_series)
    β = conv(Modulus, prescribed_dot)
    β = β[1:length(time_series)]*dt
end

#= to test =#
function noIndex_objective(params, modulus, time_series, dt, prescribed_dot, measured)
    mod = (t->modulus(t,params))
    convolved = boltzconvolve(mod, time_series, dt, prescribed_dot)
    cost = sum((measured - convolved).^2)
    return cost
end

function Index_objective(params, modulus, time_series, dt, prescribed_dot, measured, indices)
    mod = (t->modulus(t,params))
    convolved = boltzconvolve(mod, time_series, dt, prescribed_dot)
    cost = sum((measured - convolved[indices]).^2)
    return cost
end

# get signal
dt = 0.01
t = Vector{Float64}(0.0:dt:1000.0)
loading_derivative = -2*(t .- 50)
exact_response = 102 .- 102*exp.(-t) .- 2t
params = ones(5)

# to bench whether function needs to be specialised, use full range of indices
indicestouse = collect(1:length(t))

# in real code, can pass selected indices only
exact_response_indices_only = exact_response[indicestouse]

@btime noIndex_objective($params, $mod2test, $t, $dt, $loading_derivative, $exact_response)
@btime Index_objective($params, $mod2test, $t, $dt, $loading_derivative, $exact_response_indices_only, $indicestouse)

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
import Random.randperm
import RandomV06.randperm_ver

"""
    ijv = empty_graph_ijv(n)
"""
empty_graph_ijv(n::Integer) = IJV{Float64,Int64}(Int64(n),0,[],[],[])

"""
    ijv = empty_graph(n)
"""
empty_graph(n::Integer) = sparse(empty_graph_ijv(n))

"""
    graph = path_graph(n)
"""
path_graph(n) = sparse(path_graph_ijv(n))

"""
    ijv = path_graph_ijv(n::Int64)
"""
function path_graph_ijv(n::Integer)
    IJV(n, 2*(n-1),
        [collect(1:(n-1)) ; collect(2:n)],
        [collect(2:n); collect(1:(n-1))],
        ones(2*(n-1)))
end

"""
    graph = star_graph(n)
"""
star_graph(n) = sparse(star_graph_ijv(n))

"""
    ijv = star_graph_ijv(n::Int64)
"""
function star_graph_ijv(n::Integer)
    IJV(n, 2*(n-1),
        [ones(Int, n-1) ; collect(2:n)],
        [collect(2:n); ones(Int, n-1)],
        ones(2*(n-1)))
end

"""
    graph = complete_graph(n)
"""
function complete_graph(n::Integer)
  return sparse(ones(n,n) - Matrix(I,n,n))
end

"""
    ijv = complete_graph_ijv(n)
"""
complete_graph_ijv(n::Integer) = IJV(complete_graph(n))

"""
    graph = complete_bipartite_graph(n)
"""
function complete_bipartite_graph(n::Integer)
    return sparse([zeros(n,n) ones(n,n) ; ones(n,n) zeros(n,n)])
end

"""
    ijv = complete_bipartite_graph_ijv(n)
"""
complete_bipartite_graph_ijv(n::Integer) = IJV(complete_bipartite_graph(n))

"""
    graph = ring_graph(n)
"""
ring_graph(n) = sparse(ring_graph_ijv(n))

"""
    ijv = ring_graph_ijv(n)
"""
function ring_graph_ijv(n::Integer)
    if n == 1
        return empty_graph_ijv(1)
    else
        return IJV(n, 2*n,
        [collect(1:(n-1)) ; collect(2:n); 1; n],
        [collect(2:n); collect(1:(n-1)); n; 1],
        ones(2*n))
    end
end

"""
    graph = generalized_ring(n, gens)

A generalization of a ring graph.
The vertices are integers modulo n.
Two are connected if their difference is in gens.
For example,

```
generalized_ring(17, [1 5])
```
"""
generalized_ring(n::T, gens::Array{T}) where T <: Integer =
    sparse(generalized_ring_ijv(n, gens))

function generalized_ring_ijv(n::T, gens::Array{T}) where T <: Integer
    gens = gens[mod.(gens, n) .> 0]
    if isempty(gens)
        return empty_graph_ijv(n)
    end

    k = length(gens)
    m = 2*n*k
    ai = zeros(Int64,m)
    aj = zeros(Int64,m)
    ind = 1
    for i in 0:(n-1)
        for j in 1:k
            ai[ind] = i
            aj[ind] = mod(i+gens[j],n)
            ind = ind + 1
            ai[ind] = i
            aj[ind] = mod(i-gens[j],n)
            ind = ind + 1
        end
    end
    return IJV(n, m, 1 .+ ai,1 .+ aj,ones(m))
end

"""
    graph = rand_gen_ring(n, k; verbose = false, ver=Vcur)

A random generalized ring graph of degree k.
Gens always contains 1, and the other k-1 edge types
are chosen from an exponential distribution
"""
rand_gen_ring(n::Integer, k::Integer; verbose=false, ver=Vcur) =
    sparse(rand_gen_ring_ijv(n,k,verbose=verbose, ver=ver))

function rand_gen_ring_ijv(n::Integer, k::Integer; verbose=false, ver=Vcur)

    # if any of n, 2n, 3n etc. is in gens we will have self loops
    gens = [0]
    while 0 in (gens .% n)
        gens = [1; 1 .+ ceil.(Integer,exp.(rand_ver(ver, k-1)*log(n-1)))]
    end

    if verbose
        println("gens: $(gens)")
    end
    return generalized_ring_ijv(n, gens)
end

"""
    graph = hyperCube(d::Int64)

The d dimensional hypercube.  Has 2^d vertices and d*2^(d-1) edges.
"""
hypercube(d::Integer) = sparse(hypercube_ijv(d))

function hypercube_ijv(d::Integer)
    @assert d >= 0

    if d == 0
        return empty_graph_ijv(1)
    end

    ijvm = hypercube_ijv(d-1)

    ijv = disjoin(ijvm, ijvm)
    append!(ijv.i, [collect(1:2^(d-1)); 2^(d-1) .+ collect(1:2^(d-1))]  )
    append!(ijv.j, [2^(d-1) .+ collect(1:2^(d-1)); collect(1:2^(d-1))] )
    append!(ijv.v, ones(2^d))
    ijv.nnz += 2^d

    return ijv

end

"""
    graph = completeBinaryTree(n::Int64)

The complete binary tree on n vertices
"""
complete_binary_tree(n::Integer) = sparse(cbt_ijv(n))

function cbt_ijv(n::Integer)

    k = div(n-1,2)

    if 2*k+1 < n
        ii0 = Int[n-1]
        jj0 = Int[n]
    else
        ii0 = Int[]
        jj0 = Int[]
    end

    ii = [collect(1:k); collect(1:k); ii0]
    jj = [2*collect(1:k); 2*collect(1:k) .+ 1; jj0]

    return IJV(n, 2*(n-1),
        [ii;jj], [jj;ii], ones(2*length(ii)))

end

"""
    graph = grid2(n::Int64, m::Int64; isotropy=1)

An n-by-m grid graph.  iostropy is the weighting on edges in one direction.
"""
grid2(n::Integer, m::Integer; isotropy=1.0) =
    sparse(grid2_ijv(n, m; isotropy=isotropy))

grid2_ijv(n::Integer, m::Integer; isotropy=1.0) =
    product_graph(isotropy*path_graph_ijv(n), path_graph_ijv(m))

grid2(n::Integer) = grid2(n,n)
grid2_ijv(n::Integer) = grid2_ijv(n,n)

"""
    graph = grid3(n1, n2, n3)
    graph = grid3(n)

An n1-by-n2-by-n3 grid graph.
"""
grid3(n1::Integer, n2::Integer, n3::Integer) =
    sparse(grid3_ijv(n1, n2, n3))

grid3_ijv(n1::Integer, n2::Integer, n3::Integer) =
    product_graph(path_graph(n1), product_graph(path_graph(n2), path_graph(n3)))

grid3(n) = grid3(n,n,n)
grid3_ijv(n) = grid3_ijv(n,n,n)

#=
"""
    graph = wGrid2(n::Integer; weightGen::Function=rand)

An n by n grid with random weights. User can specify the weighting scheme.
"""
wgrid2(n::Integer; weightGen::Function=r(x->rand(x)) =
    sparse(wgrid2_ijv(n, weightGen = weightGen))

function wgrid2_ijv(n::Integer; weightGen::Function=rand)
    gr2 = compress(grid2_ijv(n))

    # inefficient for backwards compatibility
    for i in 1:gr2.nnz
        gr2.v[i] = weightGen()
        if gr2.i[i] < gr2.j[i]
            gr2.v[i] = 0
        end
    end

    return compress(gr2 + gr2')

end
=#

"""
    graph = grid2coords(n::Int64, m::Int64)
    graph = grid2coords(n::Int64)

Coordinates for plotting the vertices of the n-by-m grid graph
"""
function grid2coords(n::Int64, m::Int64)
  x = kron(collect(1:n),ones(m))
  y = kron(ones(n),collect(1:m))
  return x, y
end # grid2coords

grid2coords(n) = grid2coords(n, n)

"""
    graph = rand_matching(n::Integer; ver=Vcur)

A random matching on n vertices
"""
rand_matching(n::Integer; ver=Vcur) = sparse(rand_matching_ijv(n, ver=ver))

function rand_matching_ijv(n::Integer; ver=Vcur)

  p = randperm_ver(ver,n)
  n1 = convert(Int64,floor(n/2))
  n2 = 2*n1

  ii = p[1:n1]
  jj = p[(n1+1):n2]

  return IJV(n, n2,
    [ii;jj], [jj; ii], ones(n2))

end

"""
    graph = rand_regular(n, k; ver=Vcur)

A sum of k random matchings on n vertices
"""
rand_regular(n::Integer, k::Integer; ver=Vcur) = sparse(rand_regular_ijv(n, k, ver=ver))

function rand_regular_ijv(n::Integer, k::Integer; ver=Vcur)

    n1 = convert(Int64,floor(n/2))
    n2 = 2*n1

    ii = Array{Int64}(undef, n1*k)
    jj = Array{Int64}(undef, n1*k)

    ind = 0
    for i in 1:k
        p = randperm_ver(ver, n)
        for j in 1:n1
            ind += 1
            ii[ind] = p[j]
            jj[ind] = p[n1+j]
        end
    end

    return IJV(n, k*n2,
        [ii;jj], [jj;ii], ones(k*n2))

end

"""
    graph = grown_graph(n, k; ver=Vcur)

Create a graph on n vertices.
For each vertex, give it k edges to randomly chosen prior
vertices.
This is a variety of a preferential attachment graph.
"""
grown_graph(n::Integer, k::Integer; ver=Vcur) = sparse(grown_graph_ijv(n,k, ver=ver))

function grown_graph_ijv(n::Integer, k::Integer; ver=Vcur)
    ii = Int[]
    jj = Int[]

    for i = 1:k
        append!(ii, collect(2:n))
        append!(jj, ceil.(Integer,collect(1:n-1).*rand_ver(ver, n-1)))
    end

    return IJV(n, 2*k*(n-1),
        [ii;jj], [jj;ii], ones(2*k*(n-1)))

end # grownGraph

# used in grownGraphD
function randSet(n::Integer,k::Integer; ver=Vcur)
    if n == k
        return collect(1:n)
    elseif n < k
        error("n must be at least k")
    else

        s = sort(ceil.(Integer,n*rand_ver(ver, k)))
        good = (minimum(s[2:end]-s[1:(end-1)]) > 0)
        while good == false
            s = sort(ceil.(Integer,n*rand_ver(ver, k)))
            good = (minimum(s[2:end]-s[1:(end-1)]) > 0)
        end

        return s

    end
end

"""
    graph = grown_graph_d(n::Integer, k::Integer; ver=Vcur)

Like a grownGraph, but it forces the edges to all be distinct.
It starts out with a k+1 clique on the first k vertices
"""
grown_graph_d(n::Integer, k::Integer; ver=Vcur) =
    sparse(grown_graph_d_ijv(n::Integer, k::Integer, ver=ver))

function grown_graph_d_ijv(n::Integer, k::Integer; ver=Vcur)
    @assert n > k > 1

    u = zeros(Int64, k*(n-k-1))
    v = zeros(Int64, k*(n-k-1))

    for i in (k+2):n
        nb = randSet(i-1,k; ver=ver)
        u[(i-k-2)*k .+ collect(1:k)] .= i
        v[(i-k-2)*k .+ collect(1:k)] .= nb
    end

    ijv = IJV(n, 2*length(u),
        [u;v], [v;u], ones(2*length(u)))

    clique = complete_graph_ijv(k+1)
    clique.n = n

    return ijv + clique
end

"""
    graph = pref_attach(n::Int64, k::Int64, p::Float64; ver=Vcur)

A preferential attachment graph in which each vertex has k edges to those
that come before.  These are chosen with probability p to be from a random vertex,
and with probability 1-p to come from the endpoint of a random edge.
It begins with a k-clique on the first k+1 vertices.
"""
pref_attach(n::Integer, k::Integer, p::Float64; ver=Vcur) = sparse(pref_attach_ijv(n,k,p,ver=ver))

function pref_attach_ijv(n::Integer, k::Integer, p::Float64; ver=Vcur)
    @assert n > k
    if n == (k+1)
        return complete_graph_ijv(n)
    end

    u = zeros(Int64,n*k)
    v = zeros(Int64,n*k)

    # fill in the initial clique
    # this will accidentally double every edge in the clique
    # we clean it up at the end
    ind = 1
    for i in 1:(k+1)
        for j in 1:(k+1)
            if i != j
                u[ind] = i
                v[ind] = j
                ind += 1
            end
        end
    end

    s = zeros(Int64,k)
    for i in (k+2):n
        distinct = false
        while distinct == false
            for j in 1:k
                if rand_ver(ver) < p
                    s[j] = rand_ver(ver,1:(i-1))
                else
                    s[j] = v[rand_ver(ver,1:(k*(i-1)))]
                end
            end
            s = sort(s)
            distinct = true
            for ii in 1:(k-1)
                if s[ii] == s[ii+1]
                    distinct = false
                end
            end
            # distinct = (minimum(s[2:end]-s[1:(end-1)]) > 0)

        end

        for j in 1:k
            u[ind] = i
            v[ind] = s[j]
            ind += 1
        end

    end # for i

    w = ones(Float64,n*k)

    w[1:(k*(k+1))] .= 0.5

    return IJV(n, length(w),
        [u;v], [v;u], [w;w])

end

"""
    graph = randperm(mat::AbstractMatrix)
            randperm(f::Expr)

Randomly permutes the vertex indices
"""
randperm(mat::AbstractMatrix) = randperm_ver(Vcur, mat)

function randperm_ver(::Type{V06}, mat::AbstractMatrix)
    perm = randperm_ver(V06, mat.n)
    return mat[perm,perm]
end

function randperm_ver(::Type{Vcur}, mat::AbstractMatrix)
    perm = randperm_ver(Vcur, mat.n)
    return mat[perm,perm]
end

randperm(f::Expr; ver=Vcur) = randperm(eval(f), ver=ver)

function randperm_ver(::Type{V06}, a::IJV)
    perm = invperm(randperm_ver(V06, a.n))  # invperm is for V06 complete_binary_tree
    return IJV(a.n, a.nnz,
        perm[a.i], perm[a.j], a.v)
end

function randperm_ver(::Type{Vcur}, a::IJV)
    perm = randperm_ver(Vcur, a.n)
    return IJV(a.n, a.nnz,
        perm[a.i], perm[a.j], a.v)
end

function randperm_ver!(::Type{Vcur}, a::IJV)
    perm = randperm_ver(Vcur, a.n)
    a.i = perm[a.i];
    a.j = perm[a.j];

end

function randperm_ver!(::Type{V06}, a::IJV)
    perm = randperm_ver(V06, a.n)
    a.i = perm[a.i];
    a.j = perm[a.j];

end

randperm(a::IJV) = randperm_ver(Vcur, a)

"""
    graph = ErdosRenyi(n::Integer, m::Integer; ver=Vcur)

Generate a random graph on n vertices with m edges.
The actual number of edges will probably be smaller, as we sample
with replacement
"""
function ErdosRenyi(n::Integer, m::Integer; ver=Vcur)
    ai = rand_ver(ver, 1:n, m)
    aj = rand_ver(ver, 1:n, m)
    ind = (ai .!= aj)

    mat = sparse(ai[ind],aj[ind],1.0,n,n)
    mat = mat + mat'
    unweight!(mat)
    return mat
end

ErdosRenyi_ijv(n::Integer, m::Integer; ver=Vcur) = IJV(ErdosRenyi(n::Integer, m::Integer, ver=ver))

"""
    graph = ErdosRenyiCluster(n::Integer, k::Integer; ver=Vcur)

Generate an ER graph with average degree k,
and then return the largest component.
Will probably have fewer than n vertices.
If you want to add a tree to bring it back to n,
try ErdosRenyiClusterFix.
"""
function ErdosRenyiCluster(n::Integer, k::Integer; ver=Vcur)
    m = ceil(Integer,n*k/2)
    ai = rand_ver(ver, 1:n, m)
    aj = rand_ver(ver, 1:n, m)
    ind = (ai .!= aj)
    mat = sparse(ai[ind],aj[ind],1.0,n,n)
    mat = mat + mat'

    return biggestComp(mat)
end

ErdosRenyiCluster_ijv(n::Integer, k::Integer; ver=Vcur) = IJV(ErdosRenyiCluster(n, k, ver=ver))

"""
    graph = ErdosRenyiClusterFix(n::Integer, k::Integer; ver=Vcur)

Like an Erdos-Renyi cluster, but add back a tree so
it has n vertices
"""
function ErdosRenyiClusterFix(n::Integer, k::Integer; ver=Vcur)
    m1 = ErdosRenyiCluster(n, k)
    n2 = n - size(m1)[1]
    if (n2 > 0)
        m2 = complete_binary_tree(n2)
        return join_graphs(m1,m2,1,ver=ver)
    else
        return m1
    end
end

function ErdosRenyiClusterFix_ijv(n::Integer, k::Integer; ver=Vcur)
    m1 = ErdosRenyiCluster_ijv(n, k, ver=ver)
    n2 = n - m1.n
    if (n2 > 0)
        m2 = cbt_ijv(n2)
        join_graphs!(m1,m2,1,ver=ver)
    end

    return m1

end

"""
    graph = pure_random_graph(n::Integer; verbose=false, ver=Vcur)

Generate a random graph with n vertices from one of our natural distributions
"""
pure_random_graph(n::Integer; verbose=false, prefix="", ver=Vcur) =
    sparse(pure_random_ijv(n; verbose=verbose, prefix=prefix, ver=ver))

function pure_random_ijv_v6(n::Integer; verbose=false, prefix="")

    gr = []
    wt = []

    push!(gr,:(path_graph_ijv($n)))
    push!(wt,1)

    push!(gr,:(ring_graph_ijv($n)))
    push!(wt,3)

    push!(gr,:(cbt_ijv($n)))
    push!(wt,3)

    push!(gr,:(grown_graph_ijv($n,2, ver=V06)))
    push!(wt,6)

    push!(gr,:(firstn(grid2_ijv(ceil(Integer,sqrt($n))), $n)))
    push!(wt,6)

    push!(gr,:(rand_regular_ijv($n,3, ver=V06)))
    push!(wt,6)

    push!(gr,:(ErdosRenyiClusterFix_ijv($n,2, ver=V06)))
    push!(wt,6)

    if n >= 4
        push!(gr,:(rand_gen_ring_ijv($n,4,ver=V06, verbose=$(verbose))))
        push!(wt,6)
    end

    i = sampleByWeight(wt, ver=V06)

    # make sure get a connected graph - will want to remove.
    if verbose
        println(prefix, gr[i])
    end
    ijv = eval(gr[i])

    its = 0
    while (~isConnected(sparse(ijv))) && (its < 100)
        i = sampleByWeight(wt, ver=V06)

        ijv = eval(gr[i])
        its += 1
    end
    if its == 100
        error("Getting a disconnected graph from $(gr[i])")
    end

    return ijv

end

"""
    a = pure_random_ijv(n::Integer; verbose=false, prefix="", ver=Vcur)

Chooses among path_graph, ring_graph, grid_graph, complete_binary_tree, rand_gen_ring, grown_graph and ErdosRenyiClusterFix.
It can produce a disconnected graph.
For code that always produces a connected graph (and is the same as with Julia v0.6, use pure_random_ijv_v6)
"""
function pure_random_ijv(n::Integer; verbose=false, prefix="", ver=Vcur)
    if ver == Vcur
        pure_random_ijv_v7(n, verbose=verbose, prefix=prefix)
    else
        pure_random_ijv_v6(n, verbose=verbose, prefix=prefix)
    end
end

function pure_random_ijv_v7(n::Integer; verbose=false, prefix="")

    ver = Vcur

    n >= 4 ? rmax = 37 : rmax = 31

    r = rmax*rand_ver(ver)

    if r <= 1
        ijv = path_graph_ijv(n)
        verbose && println("$(prefix) path_graph($(n))")

    elseif r <= 4
        ijv = ring_graph_ijv(n)
        verbose && println("$(prefix) ring_graph($(n))")

    elseif r <= 7
        ijv = cbt_ijv(n)
        verbose && println("$(prefix) complete_binary_tree($(n))")

    elseif r <= 13
        ijv = grown_graph_ijv(n,2, ver=ver)
        verbose && println("$(prefix) grown_graph($(n))")

    elseif r <= 19
        ijv = firstn(grid2_ijv(ceil(Integer,sqrt(n))), n)
        verbose && println("$(prefix) firstn_grid2($(n))")

    elseif r <= 25
        ijv = rand_regular_ijv(n,3, ver=ver)
        verbose && println("$(prefix) rand_regular_ijv($(n),3)")

    elseif r <= 31
        ijv = ErdosRenyiClusterFix_ijv(n,2, ver=ver)
        verbose && println("$(prefix) ErdosRenyiClusterFix_ijv($n,2)")

    else
        ijv = rand_gen_ring_ijv(n,4, verbose=verbose, ver=ver)
        verbose && println("$(prefix) rand_gen_ring($(n), 4)")

    end

    return ijv

end

"""
    ind = sampleByWeight(wt; ver=Vcur)

sample an index with probability proportional to its weight given here
"""
function sampleByWeight(wt; ver=Vcur)
    r = rand_ver(ver, 1)*sum(wt)
    findall(cumsum(wt) .> r)[1]
end

"""
    graph = semiwted_chimera(n::Integer; verbose=false, ver=Vcur)

A Chimera graph with some weights.  The weights just appear when graphs are combined.
For more interesting weights, use `wted_chimera`
"""
function semiwted_chimera(n::Integer; verbose=false, prefix="", ver=Vcur)
    if ver == V06
        return sparse(semiwted_chimera_ijv_v6(n; verbose=verbose, prefix=prefix, ver=ver))
    else
        return sparse(semiwted_chimera_ijv(n; verbose=verbose, prefix=prefix, ver=ver))
    end
end

function semiwted_chimera(n::Integer, k::Integer; verbose=false, prefix="", ver=Vcur)
    if ver == V06
        seed_ver!(ver, 100*n+k)
        return sparse(semiwted_chimera_ijv_v6(n; verbose=verbose, prefix=prefix, ver=ver))
    end

    seed_ver!(ver, hash(n, hash(k)))

    return semiwted_chimera(n; verbose=verbose, prefix=prefix, ver=ver)

end

function semiwted_chimera_ijv(n::Integer; verbose=false, prefix="", ver=Vcur)
    if ver == V06
        return semiwted_chimera_ijv_v6(n, verbose=verbose, prefix=prefix, ver=ver)
    else
        return semiwted_chimera_ijv_v7(n, verbose=verbose, prefix=prefix, ver=ver)
    end
end

function semiwted_chimera_ijv_v6(n::Integer; verbose=false, prefix="", ver=Vcur)

    if (n < 2)
        return empty_graph_ijv(1)
    end

    r = rand_ver(ver)^2

    if (n < 30) || (rand_ver(ver) < .2)

        gr = pure_random_ijv(n, verbose=verbose, prefix=prefix, ver=ver)

        return randperm_ver(ver, gr)
    end

    if (n < 200)
        # just join disjoint copies of graphs

        n1 = 10 + floor(Integer,(n-20)*rand_ver(ver))
        n2 = n - n1
        k = ceil(Integer,exp(rand_ver(ver)*log(min(n1,n2)/2)))

        if verbose
            println(prefix,"joinGraphs($(r)*chimera($(n1)),chimera($(n2)),$(k))")
        end

        pr = string(" ",prefix)
        gr = join_graphs(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
          chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
          k, ver=ver)

        return randperm_ver(ver, gr)
    end

    # split with probability .7

    if (rand_ver(ver) < .7)
        n1 = ceil(Integer,10*exp(rand_ver(ver)*log(n/20)))

        n2 = n - n1
        k = floor(Integer,1+exp(rand_ver(ver)*log(min(n1,n2)/2)))

        if verbose
            println(prefix,"joinGraphs($(r)*chimera($(n1)),chimera($(n2)),$(k))")
        end

        pr = string(" ",prefix)
        gr = join_graphs(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
          chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
          k, ver=ver)

        return randperm_ver(ver, gr)

    else
        n1 = floor(Integer,10*exp(rand_ver(ver)*log(n/100)))

        n2 = floor(Integer, n / n1)

        if (rand_ver(ver) < .5)

            if verbose
                println(prefix,"productGraph($(r)*chimera($(n1)),chimera($(n2)))")
            end
            pr = string(" ",prefix)
            gr = product_graph(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
              chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver))
        else

            k = floor(Integer,1+exp(rand_ver(ver)*log(min(n1,n2)/10)))

            if verbose
                println(prefix, "generalizedNecklace($(r)*chimera($(n1)),chimera($(n2)),$(k))")
            end
            pr = string(" ",prefix)
            gr = generalized_necklace(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
              chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
              k, ver=ver)

        end

        n3 = n - gr.n
        if (n3 > 0)

            if verbose
                println(prefix, "join_graphs!(gr($(gr.n)),chimera($(n3)),2), $(ver)")
            end

            pr = string(" ",prefix)
            join_graphs!(gr,chimera_ijv(n3;verbose=verbose,prefix=pr, ver=ver),
            2, ver=ver)

        end

        return randperm_ver(ver, gr)

    end
end

function semiwted_chimera_ijv_v7(n::Integer; verbose=false, prefix="", ver=Vcur)

    if (n < 2)
        return empty_graph_ijv(1)
    end

    gr = empty_graph(n)

    r = rand_ver(ver)^2

    pr = string(" ",prefix)

    if (n < 30) || (rand_ver(ver) < .2)

        gr = pure_random_ijv(n, verbose=verbose, prefix=prefix, ver=ver)

    elseif (n < 200)
        # just join disjoint copies of graphs

        n1 = 10 + floor(Integer,(n-20)*rand_ver(ver))
        n2 = n - n1
        k = ceil(Integer,exp(rand_ver(ver)*log(min(n1,n2)/2)))

        if verbose
            println(prefix,"joinGraphs($(r)*chimera($(n1)),chimera($(n2)),$(k))")
        end

        gr = join_graphs(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
          chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
          k, ver=ver)

    elseif rand_ver(ver) < .7

            # split with probability .7

            n2 = ceil(Integer,10*exp(rand_ver(ver)*log(n/20)))

            if n2 < n/2
                n1 = n - n2
            else
                n1 = n2
                n2 = n - n2
            end

            k = floor(Integer,1+exp(rand_ver(ver)*log(min(n1,n2)/2)))

            if verbose
                println(prefix,"joinGraphs($(r)*chimera($(n1)),chimera($(n2)),$(k))")
            end

            gr = r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver)
            join_graphs!(gr,
                chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
                k, ver=ver)

        else
            n1 = floor(Integer,10*exp(rand_ver(ver)*log(n/100)))

            n2 = floor(Integer, n / n1)

            choice = rand_ver(ver)

            if choice < .45

                if verbose
                    println(prefix,"productGraph($(r)*chimera($(n1)),chimera($(n2)))")
                end

                gr = product_graph(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
                chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver))

            elseif choice < 0.9

                k = floor(Integer,1+exp(rand_ver(ver)*log(min(n1,n2)/10)))

                if verbose
                    println(prefix, "generalizedNecklace($(r)*chimera($(n1)),chimera($(n2)),$(k))")
                end

                gr = generalized_necklace(r*chimera_ijv(n1;verbose=verbose,prefix=pr, ver=ver),
                chimera_ijv(n2;verbose=verbose,prefix=pr, ver=ver),
                k, ver=ver)

            else

                n1 = div(n,2)
                k = ceil(Integer,exp(rand_ver(ver)*log(n/2)))

                if verbose
                    println(prefix, "two_lift(chimera($(n1), $(k)")
                end

                gr1 = chimera_ijv(n1;verbose=verbose, prefix=pr, ver=ver)
                gr = two_lift(gr1, k)

            end

        end

    n3 = n - gr.n
    if (n3 > 0)

        if verbose
            println(prefix, "join_graphs!(gr($(gr.n)),chimera($(n3)),2), $(ver)")
        end

        join_graphs!(gr,chimera_ijv(n3;verbose=verbose,prefix=pr, ver=ver),
        2, ver=ver)
    end

    randperm_ver!(ver, gr)
    return gr

end

"""
    graph = chimera(n::Integer; verbose=false, ver=Vcur)

Builds a chimeric graph on n vertices.
The components come from pureRandomGraph,
connected by joinGraphs, productGraph and generalizedNecklace
"""
function chimera(n::Integer; verbose=false, prefix="", ver=Vcur)

    gr = sparse(chimera_ijv(n, verbose=verbose, prefix=prefix, ver=ver))

    if ver != V06
        k = 1+floor(Integer,-log.(rand_ver(ver))/2)

        if k > 1
            if verbose
                println(prefix, "thicken($(k))")
            end

            gr = thicken(gr, k)
        end
    end

    return gr
end

function chimera_ijv(n::Integer; verbose=false, prefix="", ver=Vcur)

    ijv = semiwted_chimera_ijv(n; verbose=verbose, prefix=prefix, ver=ver)

    if ver == V06
        gr = sparse(ijv)
        unweight!(gr)
        ijv = IJV(gr)
    end

    return ijv

end

"""
    graph = chimera(n::Integer, k::Integer; verbose=false, ver=Vcur)

Builds the kth chimeric graph on n vertices.
It does this by resetting the random number generator seed.
It should captute the state of the generator before that and then
return it, but it does not yet.
"""
function chimera(n::Integer, k::Integer; verbose=false, prefix="", ver=Vcur)
    if ver == V06
        seed_ver!(ver, 100*n+k)
    else
        seed_ver!(ver, hash(n, hash(k)))
    end

    g = chimera(n; verbose=verbose, prefix=prefix, ver=ver)
    return g
end

"""
    graph = randWeight(graph; ver=Vcur)

Applies one of a number of random weighting schemes to the edges of the graph
"""
function rand_weight(a; ver=Vcur)

    if (rand_ver(ver) < .2)
        return a
    else
        return rand_weight_sub(a, ver=ver)
    end
end

function rand_weight_sub(a; ver=Vcur)

    n = a.n
    (ai,aj) = findnz(a)
    m = length(ai)

    # potentials or edge-based

    if (rand_ver(ver) < .3)
        w = rand_ver(ver, m)

    else
        v = randn_ver(ver, a.n)

        # mult by matrix ?
        if (rand_ver(ver) < .5)

            invdeg = sparse(Diagonal(1 ./(a*ones(size(a)[1]))))
            if (rand_ver(ver) < .5)
                for i in 1:10
                    v = a * (invdeg * v)
                    v = v .- mean(v)
                end
            else

                for i in 1:10
                    v = v - a * (invdeg * v)
                    v = v .- mean(v)
                end
            end
        end

        w = abs.(v[ai] - v[aj])

    end

    # reciprocate or not?

    w[w.==0] .= 1
    w[isnan.(w)] .= 1

    if (rand_ver(ver) < .5)

        w = 1 ./w
    end

    w = w / mean(w)

    ar = sparse(ai,aj,w,n,n)
    ar = ar + ar';
    return ar
end

"""
    graph = wted_chimera(n::Integer, k::Integer; verbose=false, ver=Vcur)

Builds the kth wted chimeric graph on n vertices.
It does this by resetting the random number generator seed.
It should captute the state of the generator before that and then
return it, but it does not yet.
"""
function wted_chimera(n::Integer, k::Integer; verbose=false, ver=Vcur)
    if ver == V06
        seed_ver!(ver, 100*n+k)
    else
        seed_ver!(ver, hash(n, hash(k)))
    end
    g = wted_chimera(n; verbose=verbose, ver=ver)
    return g
end

"""
    graph = wted_chimera(n::Integer)

Generate a chimera, and then apply a random weighting scheme
"""
function wted_chimera(n::Integer; verbose=false, ver=Vcur)

    gr = semiwted_chimera(n; verbose=verbose, ver=ver)

    if ver != V06
        k = 1+floor(Integer,-log.(rand_ver(ver))/2)

        if k > 1
            if verbose
                println("thicken($(k))")
            end

            gr = thicken(gr, k)
        end
    end

    return rand_weight(gr, ver=ver)
end


text:
using Atum
using Atum.ShallowWater

using StaticArrays: SVector
using WriteVTK

import Atum: boundarystate
function boundarystate(law::ShallowWaterLaw, n⃗, q⁻, aux⁻, _)
  ρ⁻, ρu⃗⁻, ρθ⁻ = ShallowWater.unpackstate(law, q⁻)
  ρ⁺, ρθ⁺ = ρ⁻, ρθ⁻
  ρu⃗⁺ = ρu⃗⁻ - 2 * (n⃗' * ρu⃗⁻) * n⃗
  SVector(ρ⁺, ρu⃗⁺..., ρθ⁺), aux⁻
end

function bickleyjet(law, x⃗)
  FT = eltype(law)
  x, y = x⃗

  ϵ = FT(1 / 10)
  l = FT(1 / 2)
  k = FT(1 / 2)

  U = cosh(y)^(-2)

  Ψ = exp(-(y + l / 10)^2 / (2 * (l^2))) * cos(k * x) * cos(k * y)

  u = Ψ * (k * tan(k * y) + y / (l^2))
  v = -Ψ * k * tan(k * x)

  ρ = FT(1)
  ρu = ρ * (U + ϵ * u)
  ρv = ρ * (ϵ * v)
  ρθ = ρ * sin(k * y)

  SVector(ρ, ρu, ρv, ρθ)
end

function run(A, FT, N, K; volume_form=WeakForm(), outputvtk=true)
  Nq = N + 1

  law = ShallowWaterLaw{FT, 2}()

  cell = LobattoCell{FT, A}(Nq, Nq)
  v1d = range(FT(-2π), stop=FT(2π), length=K+1)
  grid = brickgrid(cell, (v1d, v1d); periodic = (true, false))

  dg = DGSEM(; law, grid, volume_form,
               surface_numericalflux = RoeFlux())

  cfl = FT(1 // 8)
  dt = cfl * step(v1d) / N / sqrt(constants(law).grav)
  timeend = @isdefined(_testing) ? 10dt : FT(200)

  q = bickleyjet.(Ref(law), points(grid))

  if outputvtk
    vtkdir = joinpath("output", "shallow_water", "bickleyjet")
    mkpath(vtkdir)
    pvd = paraview_collection(joinpath(vtkdir, "timesteps"))
  end

  do_output = function(step, time, q)
    if outputvtk && step % ceil(Int, timeend / 100 / dt) == 0
      filename = "step$(lpad(step, 6, '0'))"
      vtkfile = vtk_grid(joinpath(vtkdir, filename), grid)
      P = Bennu.toequallyspaced(cell)
      ρθ = last(components(q))
      vtkfile["ρθ"] = vec(Array(P * ρθ))
      vtk_save(vtkfile)
      pvd[time] = vtkfile
    end
  end

  odesolver = LSRK54(dg, q, dt)

  outputvtk && do_output(0, FT(0), q)
  solve!(q, timeend, odesolver; after_step=do_output)
  outputvtk && vtk_save(pvd)
end

let
  A = Array
  FT = Float64
  N = 3

  K = 16
  run(A, FT, N, K)
end


text:
include("nonexist.jl")

foo(s)


text:
using DataFrames

module TextAnalysis
    using Languages
    using DataFrames
    using BinDeps
    using Compat

    importall Base
    import DataFrames.DataFrame

    export AbstractDocument, Document
    export FileDocument, StringDocument, TokenDocument, NGramDocument
    export GenericDocument
    export Corpus, DirectoryCorpus
    export stemmer_types, Stemmer
    export DocumentTermMatrix
    export text, tokens, ngrams
    export text!, tokens!, ngrams!
    export documents
    export language, name, author, timestamp
    export language!, name!, author!, timestamp!
    export ngram_complexity
    export lexicon, update_lexicon!, lexical_frequency, lexicon_size
    export inverse_index, update_inverse_index!, index_size
    export remove_corrupt_utf8
    export remove_corrupt_utf8!
    export remove_punctuation, remove_numbers, remove_case, remove_whitespace
    export remove_punctuation!, remove_numbers!, remove_case!, remove_whitespace!
    export remove_nonletters, remove_nonletters!
    export remove_words, remove_stop_words, remove_articles
    export remove_words!, remove_stop_words!, remove_articles!
    export remove_definite_articles, remove_indefinite_articles
    export remove_definite_articles!, remove_indefinite_articles!
    export remove_prepositions, remove_pronouns, stem, tag_pos
    export remove_prepositions!, remove_pronouns!, stem!, tag_pos!
    export remove_html_tags, remove_html_tags!
    export prepare!
    export frequent_terms, sparse_terms
    export remove_frequent_terms!, remove_sparse_terms!
    export dtv, each_dtv, dtm, tdm
    export TextHashFunction, index_hash, cardinality, hash_function, hash_function!
    export hash_dtv, each_hash_dtv, hash_dtm, hash_tdm
    export standardize!
    export tf_idf, lsa, lda
    export tf_idf!, lsa!, lda!
    export remove_patterns!, remove_patterns

    export strip_patterns, strip_corrupt_utf8, strip_case, stem_words, tag_part_of_speech, strip_whitespace, strip_punctuation
    export strip_numbers, strip_non_letters, strip_indefinite_articles, strip_definite_articles, strip_articles
    export strip_prepositions, strip_pronouns, strip_stopwords, strip_sparse_terms, strip_frequent_terms, strip_html_tags

    include("tokenizer.jl")
    include("ngramizer.jl")
    include("document.jl")
    include("hash.jl")
    include("corpus.jl")
    include("metadata.jl")
    include("preprocessing.jl")
    include("stemmer.jl")
    include("dtm.jl")
    include("tf_idf.jl")
    include("lsa.jl")
    include("lda.jl")
    include("show.jl")
end


text:
@testset "Symmetry" begin
    # chwck type instability
    f() = Symmetry(3,3)
    (@inferred f())::Symmetry{Tuple{3,3}}

    # constructors
    @test Symmetry(3,3) == Symmetry{Tuple{3,3}}()
    @test Symmetry((3,3)) == Symmetry{Tuple{3,3}}()
    @test Symmetry(Tuple{3,3}) == Symmetry{Tuple{3,3}}()

    # ncomponents
    @test Tensorial.ncomponents(Symmetry(3,3)) == 6
    @test Tensorial.ncomponents(Symmetry(2,2)) == 3

    # Dims
    @test Dims(Symmetry(3,3)) == (3,3)
    @test Dims(Symmetry(2,2)) == (2,2)

    # length
    @test length(Symmetry(3,3,3)) == 3
    @test length(Symmetry(2,2)) == 2

    # getindex
    S = Symmetry(3,3)
    @test (@inferred S[1])::Int == 3
    @test (@inferred S[2])::Int == 3
    @test_throws Exception S[3]

    # macro
    @test @Symmetry{3,3} == Symmetry{Tuple{3,3}}

    # errors
    @test_throws Exception Symmetry(Tuple{3,Int})
    @test_throws Exception Symmetry(3,2)
end


text:
# Autogenerated wrapper script for Deno_jll for x86_64-w64-mingw32
export deno

JLLWrappers.@generate_wrapper_header("Deno")
JLLWrappers.@declare_executable_product(deno)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        deno,
        "bin\\deno.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "Automata" begin

    using KnuthBendix

    abt = KnuthBendix.Alphabet(['a', 'e', 'b', 'p'])
    ta = KnuthBendix.Automaton(abt)
    σ = KnuthBendix.initialstate(ta)
    @test sprint(show, ta) isa String

    @test ta isa KnuthBendix.AbstractAutomaton
    @test ta isa KnuthBendix.Automaton
    @test ta isa KnuthBendix.Automaton{4, Word{UInt16}}
    @test σ isa KnuthBendix.AbstractState
    @test σ isa KnuthBendix.State
    @test σ isa KnuthBendix.State{4, Word{UInt16}}

    @test isone(KnuthBendix.name(σ))
    @test !KnuthBendix.isterminal(σ)
    @test isone(KnuthBendix.rightrule(σ))

    @test length(KnuthBendix.states(ta)) == 1
    @test length(KnuthBendix.inedges(σ)) == 0
    @test length(KnuthBendix.outedges(σ)) == 4
    @test length(KnuthBendix.stateslengths(ta)) == 1
    @test length(σ) == 0

    KnuthBendix.declarerightrule!(σ, one(KnuthBendix.name(σ)))
    @test KnuthBendix.isterminal(σ)
    @test isone(KnuthBendix.rightrule(σ))

    push!(ta, Word([1]))
    KnuthBendix.addedge!(ta, 1, 1, 2)
    @test length(KnuthBendix.inedges(KnuthBendix.states(ta)[2])) == 1
    @test KnuthBendix.walk(ta, Word([1])) == (1, KnuthBendix.states(ta)[2])
    @test_throws BoundsError KnuthBendix.addedge!(ta, 5, 1, 2)
    @test_throws BoundsError KnuthBendix.addedge!(ta, 1, 3, 2)
    @test_throws BoundsError KnuthBendix.addedge!(ta, 1, 1, 3)
    @test_throws BoundsError KnuthBendix.removeedge!(ta, 5, 1, 2)
    @test_throws BoundsError KnuthBendix.removeedge!(ta, 1, 3, 2)
    @test_throws BoundsError KnuthBendix.removeedge!(ta, 1, 1, 3)
    KnuthBendix.removeedge!(ta, 1, 1, 2)
    @test KnuthBendix.walk(ta, Word([1])) == (0, KnuthBendix.states(ta)[1])

    KnuthBendix.addedge!(ta, 1, 1, 2)
    deleteat!(ta, 2)
    @test length(KnuthBendix.states(ta)) == 1
    @test length(KnuthBendix.stateslengths(ta)) == 1
    @test KnuthBendix.isfailstate(KnuthBendix.outedges(σ)[1])

    A = KnuthBendix.Alphabet(['a', 'e', 'b', 'p'])
    KnuthBendix.set_inversion!(A, 'a', 'e')
    KnuthBendix.set_inversion!(A, 'b', 'p')

    a = Word([1,2])
    b = Word([2,1])
    c = Word([3,4])
    d = Word([4,3])
    ε = one(a)
    ba = Word([3,1])
    ab = Word([1,3])
    be = Word([3,2])
    eb = Word([2,3])
    pa = Word([4,1])
    ap = Word([1,4])
    pe = Word([4,2])
    ep = Word([2,4])

    lenlexord = KnuthBendix.LenLex(A)
    rsc = KnuthBendix.RewritingSystem([a=>ε, b=>ε, c=>ε, d=>ε, ba=>ab, be=>eb, pa=>ap, pe=>ep], lenlexord)
    ia = KnuthBendix.makeindexautomaton(rsc)

    testword = Word([1,1,1,1,1,2,2,2,3,4,2,2,3,3,3,4,4,4,4,3,4,3,4,1,2,1,1,1,1,1,1,1,2,1,3,4])
    @test KnuthBendix.rewrite_from_left(testword, rsc) == KnuthBendix.rewrite_from_left(testword, ia)

    w = Word([1,3,4,1,4,4,1,1,4,2,3,2,4,2,2,3,1,2,1])
    @test KnuthBendix.rewrite_from_left(w, rsc) == KnuthBendix.rewrite_from_left(w, ia)

    @test !isempty(ia)
    @test isempty(empty!(ia))

    rsd = KnuthBendix.RewritingSystem([a=>ε, b=>ε], lenlexord)
    iad = KnuthBendix.makeindexautomaton(rsd)
    @test KnuthBendix.rewrite_from_left(testword, rsd) == KnuthBendix.rewrite_from_left(testword, iad)
end


text:
module RegistryCLI
using ArgParse: ArgParse, ArgParseSettings, @add_arg_table!
using LibGit2: LibGit2, GitRepo, GitRemote, GitCommit, GitHash, GitObject
using LocalRegistry
using Mustache

# Create a new local Registry:
#
# - make a local clone of the public Registry (which should be initially empty)
# - prepare a commit in the local clone, initializing the Registry
#
# Everything made in the local clone of the Registry should be manually
# inspected before being pushed to the public Registry.

function arg_table_create!(s)
    @add_arg_table! s["create"] begin
        "URL"
        #==# help = "public URL of the Registry (should be an empty git repo)"
        #==# required = true

        "PATH"
        #==# help = "where a local clone of the Registry should be created"
        #==# required = true

        "--description"
        #==# help = "description of the Registry"
    end
end

function create(args)
    verbose = args[:verbose]
    url     = args[:create][:URL]
    path    = args[:create][:PATH] |> abspath
    desc    = args[:create][:description]

    verbose && @info "Creating and populating local clone for registry" url path
    LocalRegistry.create_registry(path, url; description=desc)
end

# Publish a new package version in the Registry.
#
# This function should be called from a local clone of the public Registry; a
# commit will be added to it to register the new package version.
#
# Everything made in the local clone of the Registry should be manually inspected
# before being pushed to the public Registry.

function arg_table_add!(s)
    @add_arg_table! s["add"] begin
        "URL"
        #==# help = "public URL of the package"
        #==# required = true

        "--ref"
        #==# help = "git reference identifying the version to publish"
    end
end

function add(args)
    verbose = args[:verbose]
    url     = args[:add][:URL]
    gitref  = args[:add][:ref]

    # Check that pwd() is a clone of a Julia Registry
    try
        repo = GitRepo(pwd())
        reg_url = LibGit2.get(GitRemote, repo, "origin") |> LibGit2.url
        @assert isfile("Registry.toml")
        verbose && @info "Registry" url=reg_url path=pwd()
    catch
        @error "Not in the git repository of a Registry"
        return 1
    end

    verbose && @info "Releasing new version of package" url gitref

    Base.Filesystem.mktempdir() do dir
        path = joinpath(dir, "package")

        verbose && @info "Cloning package repository"
        repo = LibGit2.clone(url, path)

        if gitref !== nothing
            verbose && @info "Checking out requested commit"
            commit = LibGit2.peel(GitCommit, GitObject(repo, gitref))
            LibGit2.checkout!(repo, string(GitHash(commit)))
        end

        verbose && @info "Registering new version in the registry"
        LocalRegistry.register(path, pwd())
    end

    return 0
end

# Main entry point

function parse_args(args)
    function exc_handler(settings::ArgParseSettings, err, err_code::Int = 1)
        @error err.text
        ArgParse.show_help(stderr, settings; exit_when_done=false)
        exit(err_code)
    end

    s = ArgParseSettings(
        prog = "jlreg",
        description = "Manage a local Julia Registry from the command line.",
        exit_after_help = false,
        exc_handler = exc_handler
    )

    @add_arg_table! s begin
        "create"
        #==# help = "Create and populate a local clone of a new, empty Registry"
        #==# action = :command

        "add"
        #==# help = "Publish a new version of a package in the Registry"
        #==# action = :command

        "--verbose", "-v"
        #==# help = "Print informative messages at each step"
        #==# action = :store_true
    end

    arg_table_create!(s)
    arg_table_add!(s)

    ArgParse.parse_args(args, s, as_symbols=true)
end

function main(clargs :: Vector{<:AbstractString})
    args = parse_args(clargs)

    # --help: nothing more to do
    args === nothing && return

    # command: call the relevant function
    getfield(@__MODULE__, args[:_COMMAND_])(args)
end

function main(clargs...)
    main(collect(clargs))
end

"""
    RegistryCLI.install(; kwargs)

Install the `jlreg` script for use in a terminal.

*Keyword arguments:*
- `julia`: path to the julia executable. Defaults to the currently running julia.
- `flags`: command-line flags for julia. Defaults to `--color=yes --startup-file=no --compile=min -O0`.
- `destdir`: directory where to install the executable script. Should be writable and available in PATH. Defaults to `~/.julia`.
- `command`: name of the executable script. Defaults to `jlreg`.
- `perms`: permissions of the executable script. Defaults to `0o755` (i.e. "-rwxr-xr-x").
- `force`: allow overwriting an existing file. Defaults to `false`.
"""
function install(;
                 julia::String   = first(Base.julia_cmd()),
                 flags::String   = "--color=yes --startup-file=no --compile=min -O0",
                 destdir::String = joinpath(first(DEPOT_PATH), "bin"),
                 command::String = "jlreg",
                 perms::Integer  = 0o755,
                 force::Bool     = false,
                 )
    if Sys.iswindows()
        @error "RegistryCLI does not (yet) support installing the command-line script on Windows systems."
    end

    destdir = abspath(expanduser(destdir))
    script  = joinpath(destdir, command)
    if ispath(script) && !force
        @error("Refusing to overwrite existing file. Use `RegistryCLI.install(force=true)` to force overwriting.",
               dest=script)
        return
    end

    mkpath(destdir)
    open(script, "w") do f
        tmpl = Mustache.load(joinpath(@__DIR__, "..", "bin", "jlreg-bash"))
        write(f, render(tmpl, Dict(
            :JULIA   => julia,
            :FLAGS   => flags,
            :COMMAND => command,
            :DESTDIR => destdir,
            :PROJECT => joinpath(@__DIR__, ".."),
        )))
    end
    chmod(script, perms)

    @info "Installed Registry CLI to `$script`"
end

end # module


text:
using BinaryBuilder

# Collection of sources required to build SundialsBuilder
sources = [
    "https://computation.llnl.gov/projects/sundials/download/sundials-3.1.1.tar.gz" =>
    "a24d643d31ed1f31a25b102a1e1759508ce84b1e4739425ad0e18106ab471a24",
    "http://faculty.cse.tamu.edu/davis/SuiteSparse/SuiteSparse-4.5.3.tar.gz" =>
    "6199a3a35fbce82b155fd2349cf81d2b7cddaf0dac218c08cb172f9bc143f37a",
    "./patches",
]

# Bash recipe for building across all platforms
script = raw"""
# SuiteSparse for KLU
cd $WORKSPACE/srcdir/SuiteSparse*/

# Patches for windows build system
patch -p0 < $WORKSPACE/srcdir/patches/SuiteSparse_windows.patch

for proj in SuiteSparse_config AMD COLAMD BTF KLU; do
    cd $WORKSPACE/srcdir/SuiteSparse/$proj
    make -j${nproc} library
    INSTALL=$WORKSPACE/destdir/ make install
done

echo "KLU Includes"
ls $WORKSPACE/destdir/include
echo "KLU Lib"
ls $WORKSPACE/destdir/lib

# Now the full Sundials build
cd $WORKSPACE/srcdir/sundials-*/
patch -p0 < $WORKSPACE/srcdir/patches/Sundials_windows.patch

# On 64-bit, we need to patch the BLAS libraries to use the Julia name-mangling scheme.
if [[ ${nbits} == 64 ]]; then
    patch -p0 < $WORKSPACE/srcdir/patches/Sundials_ilp64.patch
fi

mkdir build
cd build

CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=Release -DEXAMPLES_ENABLE_C=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DKLU_ENABLE=ON -DKLU_INCLUDE_DIR=\"$prefix/include/\" -DKLU_LIBRARY_DIR=\"$prefix/lib\""
CMAKE_FLAGS="${CMAKE_FLAGS} -DBLAS_ENABLE=ON"
CMAKE_FLAGS="${CMAKE_FLAGS} -DENABLE_LAPACK=ON -DLAPACK_LIBRARIES=$LIBBLAS"

if [[ ${nbits} == 32 ]]; then
    echo "***   32-bit BUILD   ***"
    LIBBLAS="$prefix/lib/libopenblas.so"
    cmake ${CMAKE_FLAGS} -DBLAS_LIBRARIES=\"$LIBBLAS\" -DSUNDIALS_INDEX_TYPE=int32_t ..
else
    echo "***   64-bit BUILD   ***"
    LIBBLAS="$prefix/lib/libopenblas64_.so"
    cmake ${CMAKE_FLAGS} -DBLAS_LIBRARIES=\"$LIBBLAS\" ..
fi

make -j${nproc}
make install

# On windows, move all `.dll` files to `bin`. We don't want to follow symlinks
# because non-administrative users cannot create symlinks on Windows, so we
# use `cp -L` followed by `rm` instead of just `mv`.
if [[ ${target} == *-mingw32 ]]; then
    cp -L $WORKSPACE/destdir/lib/*.dll $WORKSPACE/destdir/bin
    rm -f $WORKSPACE/destdir/lib/*.dll
fi
"""

# We attempt to build for all defined platforms
platforms = [
    BinaryProvider.Windows(:i686),
    BinaryProvider.Windows(:x86_64),
    BinaryProvider.MacOS(),
    BinaryProvider.Linux(:x86_64, :glibc),
    BinaryProvider.Linux(:i686, :glibc),
    BinaryProvider.Linux(:aarch64, :glibc),
    BinaryProvider.Linux(:armv7l, :glibc),
    BinaryProvider.Linux(:powerpc64le, :glibc),
]

products(prefix) = [
    LibraryProduct(prefix, "libbtf", :libbtf),
    LibraryProduct(prefix, "libsundials_sunlinsolspfgmr", :libsundials_sunlinsolspfgmr),
    LibraryProduct(prefix, "libsundials_ida", :libsundials_ida),
    LibraryProduct(prefix, "libsundials_cvode", :libsundials_cvode),
    LibraryProduct(prefix, "libsundials_cvodes", :libsundials_cvodes),
    LibraryProduct(prefix, "libcolamd", :libcolamd),
    LibraryProduct(prefix, "libsundials_sunmatrixdense", :libsundials_sunmatrixdense),
    LibraryProduct(prefix, "libsundials_sunlinsolspbcgs", :libsundials_sunlinsolspbcgs),
    LibraryProduct(prefix, "libsundials_idas", :libsundials_idas),
    LibraryProduct(prefix, "libsundials_nvecserial", :libsundials_nvecserial),
    LibraryProduct(prefix, "libsundials_sunlinsoldense", :libsundials_sunlinsoldense),
    LibraryProduct(prefix, "libsundials_sunlinsolspgmr", :libsundials_sunlinsolspgmr),
    LibraryProduct(prefix, "libsundials_sunlinsolpcg", :libsundials_sunlinsolpcg),
    LibraryProduct(prefix, "libsundials_sunlinsolsptfqmr", :libsundials_sunlinsolsptfqmr),
    LibraryProduct(prefix, "libsundials_sunlinsolklu", :libsundials_sunlinsolklu),
    LibraryProduct(prefix, "libsundials_sunmatrixsparse", :libsundials_sunmatrixsparse),
    LibraryProduct(prefix, "libsundials_sunlinsolband", :libsundials_sunlinsolband),
    LibraryProduct(prefix, "libsundials_sunmatrixband", :libsundials_sunmatrixband),
    LibraryProduct(prefix, "libsundials_kinsol", :libsundials_kinsol),
    LibraryProduct(prefix, "libsundials_arkode", :libsundials_arkode),
    LibraryProduct(prefix, "libklu", :libklu),
    LibraryProduct(prefix, "libsuitesparseconfig", :libsuitesparseconfig),
    LibraryProduct(prefix, "libamd", :libamd),
]

dependencies = [
    "https://github.com/staticfloat/OpenBLASBuilder/releases/download/v0.2.20-7/build.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "Sundials", sources, script, platforms, products, dependencies)

max_stars_repo_name: tshort/SundialsBuilder

max_issues_repo_name: tshort/SundialsBuilder

max_forks_repo_name: tshort/SundialsBuilder


text:
function pretraverse(f::Function, p::AbstractRNode, x)
    x = f(p, x)
    for c in children(p)
        x = pretraverse(f, c, x)
    end

    return x
end
pretraverse(f::Function, tree::AbstractRTree, x) = pretraverse(f, getroot(tree), x)

function posttraverse(f::Function, p::AbstractRNode, x)
    for c in children(p)
        x = posttraverse(f, c, x)
    end
    x = f(p, x)

    return x
end
posttraverse(f::Function, tree::AbstractRTree, x) = posttraverse(f, getroot(tree), x)

function _set_traversal_index(p, x::NamedTuple{(:V, :i)})
    i = x.i + 1
    setindex!(x.V, p, i)

    return (V = x.V, i = i)
end
function _grow_traversal(p, x::NamedTuple{(:V, :i)})
    i = x.i + 1
    push!(x.V, p)

    return (V = x.V, i = i)
end

"""
    preorder(tree)
    preorder(node)

Return a vector with the nodes of a tree or subtree in preoder.
"""
function preorder(tree::AbstractRTree)
    V = similar(tree.nodes)
    p = getroot(tree)
    V, i = pretraverse(_set_traversal_index, p, (V=V, i=0))

    return V
end
preorder(p::T) where T <: AbstractRNode =
    pretraverse(_grow_traversal, p, (V=T[], i=0))[:V]

"""
    postorder(tree)
    postorder(node)

Return a vector with the nodes of a tree or subtree in preoder.
"""
function postorder(tree::AbstractRTree)
    V = similar(tree.nodes)
    p = getroot(tree)
    V, i = posttraverse(_set_traversal_index, p, (V=V, i=0))

    return V
end
postorder(p::T) where T <: AbstractRNode =
    posttraverse(_grow_traversal, p, (V=T[], i=0))[:V]

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
using NeuralQuantum, Test

num_types = [Float32, Float64]

machines = Dict()

ma = (T, N) -> NDMComplex(T, N, 2, 3)
machines["NDMComplex"] = ma

N = 4

@testset "Test Properties $name" for name=keys(machines)
    for T=num_types
        net = machines[name](T,N)
        cnet = cached(net)

        @test NeuralQuantum.out_type(net) == Complex{real(T)}
        @test NeuralQuantum.is_analytic(net) == false
    end
end

@testset "Test Cached Value $name" for name=keys(machines)
    for T=num_types
        net = machines[name](T,N)
        cnet = cached(net)

        hilb = HomogeneousFock(N, 2)
        if net isa NeuralQuantum.MatrixNet
            hilb = SuperOpSpace(hilb)
        end

        v = state(T, hilb, net)
        # compute exact
        vals = []; cvals = [];
        for i=1:spacedimension(hilb)
            set!(v, hilb, i)
            push!(vals, net(v))
            push!(cvals, cnet(v))
        end
        @test vals ≈ cvals
    end
end