text:
export BrickletLoadCellV2WeightCallbackConfiguration
struct BrickletLoadCellV2WeightCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletLoadCellV2Configuration
struct BrickletLoadCellV2Configuration
    rate::Integer
    gain::Integer
end

export BrickletLoadCellV2SPITFPErrorCount
struct BrickletLoadCellV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletLoadCellV2Identity
struct BrickletLoadCellV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletLoadCellV2
"""
Measures weight with a load cell
"""
mutable struct BrickletLoadCellV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletLoadCellV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_load_cell_v2")
        deviceInternal = package.BrickletLoadCellV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_weight
"""
    $(SIGNATURES)

Returns the currently measured weight.

If you want to get the value periodically, it is recommended to use the
:cb:`Weight` callback. You can set the callback configuration
with :func:`Set Weight Callback Configuration`.
"""
function get_weight(device::BrickletLoadCellV2)
    return device.deviceInternal.get_weight()
end

export set_weight_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Weight` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Weight` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_weight_callback_configuration(device::BrickletLoadCellV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_weight_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_weight_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Weight Callback Configuration`.
"""
function get_weight_callback_configuration(device::BrickletLoadCellV2)
    return device.deviceInternal.get_weight_callback_configuration()
end

export set_moving_average
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the weight value.

Setting the length to 1 will turn the averaging off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletLoadCellV2, average)
    device.deviceInternal.set_moving_average(average)
end

export get_moving_average
"""
    $(SIGNATURES)

Returns the length moving average as set by :func:`Set Moving Average`.
"""
function get_moving_average(device::BrickletLoadCellV2)
    return device.deviceInternal.get_moving_average()
end

export set_info_led_config
"""
    $(SIGNATURES)

Configures the info LED to be either turned off, turned on, or blink in
heartbeat mode.
"""
function set_info_led_config(device::BrickletLoadCellV2, config)
    device.deviceInternal.set_info_led_config(config)
end

export get_info_led_config
"""
    $(SIGNATURES)

Returns the LED configuration as set by :func:`Set Info LED Config`
"""
function get_info_led_config(device::BrickletLoadCellV2)
    return device.deviceInternal.get_info_led_config()
end

export calibrate
"""
    $(SIGNATURES)

To calibrate your Load Cell Bricklet 2.0 you have to

* empty the scale and call this function with 0 and
* add a known weight to the scale and call this function with the weight.

The calibration is saved in the flash of the Bricklet and only
needs to be done once.

We recommend to use the Brick Viewer for calibration, you don't need
to call this function in your source code.
"""
function calibrate(device::BrickletLoadCellV2, weight)
    device.deviceInternal.calibrate(weight)
end

export tare
"""
    $(SIGNATURES)

Sets the currently measured weight as tare weight.
"""
function tare(device::BrickletLoadCellV2)
    device.deviceInternal.tare()
end

export set_configuration
"""
    $(SIGNATURES)

The measurement rate and gain are configurable.

The rate can be either 10Hz or 80Hz. A faster rate will produce more noise.
It is additionally possible to add a moving average
(see :func:`Set Moving Average`) to the measurements.

The gain can be 128x, 64x or 32x. It represents a measurement range of
±20mV, ±40mV and ±80mV respectively. The Load Cell Bricklet uses an
excitation voltage of 5V and most load cells use an output of 2mV/V. That
means the voltage range is ±15mV for most load cells (i.e. gain of 128x
is best). If you don't know what all of this means you should keep it at
128x, it will most likely be correct.
"""
function set_configuration(device::BrickletLoadCellV2, rate, gain)
    device.deviceInternal.set_configuration(rate, gain)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletLoadCellV2)
    return device.deviceInternal.get_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletLoadCellV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletLoadCellV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletLoadCellV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletLoadCellV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletLoadCellV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletLoadCellV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletLoadCellV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletLoadCellV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletLoadCellV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletLoadCellV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletLoadCellV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletLoadCellV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletLoadCellV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
# Full-dimensional local search quantization
using HDF5, Clustering

include("../utils.jl")
include("../initializations.jl")
include("../codebook_update.jl")
include("../encodings/encode_icm.jl")

function train_lsq{T <: AbstractFloat}(
  X::Matrix{T},         # d-by-n matrix of data points to train on.
  m::Integer,           # number of codebooks
  h::Integer,           # number of entries per codebook
  R::Matrix{T},         # init rotation
  B::Matrix{Int16},     # init codes
  C::Vector{Matrix{T}}, # init codebooks
  niter::Integer,       # number of optimization iterations
  ilsiter::Integer,     # number of ILS iterations to use during encoding
  icmiter::Integer,     # number of iterations in local search
  randord::Bool,        # whether to use random order
  npert::Integer,       # The number of codes to perturb
  V::Bool=false)        # whether to print progress

  # if V
  println("**********************************************************************************************");
  println("Doing local search with $m codebooks, $npert perturbations, $icmiter icm iterations and random order = $randord");
  println("**********************************************************************************************");
  # end

  d, n = size( X );

  # Update RX
  RX = R' * X;

  # Initialize C
  C = update_codebooks( RX, B, h, V, "lsqr" )

  # Apply the rotation to the codebooks
  for i = 1:m
    C[i] = R * C[i]
  end
  @printf("%3d %e \n", -2, qerror( X, B, C ));

  # Initialize B
  for i = 1:ilsiter
    B = encoding_icm( X, B, C, icmiter, randord, npert, V );
    @everywhere gc()
  end
  @printf("%3d %e \n", -1, qerror( X, B, C ));

  obj = zeros( Float32, niter );

  for iter = 1:niter

    obj[iter] = qerror( X, B, C  );
    @printf("%3d %e \n", iter, obj[iter]);

    # Update the codebooks
    C = update_codebooks( X, B, h, V, "lsqr" )

    # Update the codes with local search
    for i = 1:ilsiter
      B = encoding_icm( X, B, C, icmiter, randord, npert, V );
      @everywhere gc()
    end

  end

  # Get the codebook for norms
  CB = reconstruct(B, C);

  dbnorms = zeros(Float32, 1, n);
  for i = 1:n
     for j = 1:d
        dbnorms[i] += CB[j,i].^2;
     end
  end

  # Quantize the norms with plain-old k-means
  dbnormsq = kmeans(dbnorms, h);
  cbnorms  = dbnormsq.centers;

  # Add the dbnorms to the codes
  B_norms  = reshape( dbnormsq.assignments, 1, n )

  return C, B, cbnorms, B_norms, obj

end

ext: jl

max_issues_repo_name: jltmtz/local-search-quantization

max_forks_repo_name: jltmtz/local-search-quantization


text:
using KUparser, KUnet, HDF5, JLD, ArgParse, Compat, CUDArt
VERSION < v"0.4-" && eval(Expr(:using,:Dates))

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "datafiles"
        help = "JLD files for input: first one will be used for training"
        nargs = '+'
        required = true
        "--out"
        help = "JLD file for output trained network"
        default = ""
        "--nogpu"
        help = "Do not use gpu"
        action = :store_true
        "--parser"
        help = "Parsing algorithm to use: g(reedy)parser, b(eam)parser, or o(racle)parser for static training"
        default = "oparser"
        "--arctype"
        help = "Move set to use: ArcEager{R1,13}, ArcHybrid{R1,13}"
        default = "ArcEager13"
        "--feats"
        help = "Features to use from KUparser.Flist"
        default = "tacl13hybrid"
        "--hidden"
        help = "One or more hidden layer sizes"
        nargs = '+'
        arg_type = Int
        default = [20000]
        "--epochs"
        help = "Minimum number of epochs to train"
        arg_type = Int
        default = 100
        "--pepochs"
        help = "Number of epochs between parsing"
        arg_type = Int
        default = 1

        # This version does not work with multi cpu
        # "--ncpu"
        # help = "Number of workers for multithreaded parsing"
        # arg_type = Int
        # default = 12

        "--pbatch"
        help = "Minibatch size for parsing"
        arg_type = Int
        default = 2000
        "--tbatch"
        help = "Minibatch size for training"
        arg_type = Int
        default = 128
        "--nbeam"
        help = "Beam size for bparser"
        arg_type = Int
        default = 10
        "--adagrad"
        help = "If nonzero apply adagrad using arg as epsilon"
        arg_type = Float32
        nargs = '+'
        default = [1f-8]
        "--dropout"
        help = "Dropout probability"
        arg_type = Float32
        nargs = '+'
        default = [0.2f0, 0.7f0]
        "--l1reg"
        help = "L1 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--l2reg"
        help = "L2 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--learningRate"
        help = "Learning rate"
        arg_type = Float32
        nargs = '+'
        default = [2f-2]
        "--maxnorm"
        help = "If nonzero upper limit on weight matrix row norms"
        arg_type = Float32
        nargs = '+'
        "--momentum"
        help = "Momentum"
        arg_type = Float32
        nargs = '+'
        "--nesterov"
        help = "Apply nesterov's accelerated gradient"
        arg_type = Float32
        nargs = '+'
        "--seed"
        help = "Random number seed (-1 leaves the default seed)"
        arg_type = Int
        default = -1
    end
    args = parse_args(s)
    for (arg,val) in args
        print("$arg:$val ")
    end
    println("")
    return args
end

function main()
    args = parse_commandline()
    KUnet.gpu(!args["nogpu"])
    args["nogpu"] && blas_set_num_threads(20)
    args["seed"] >= 0 && (srand(args["seed"]); KUnet.gpuseed(args["seed"]))

    data = Corpus[]
    deprel = nothing
    ndeps = 0
    for f in args["datafiles"]
        @date d=load(f)
        push!(data, d["corpus"])
        if deprel == nothing
            deprel = d["deprel"]
            ndeps = length(deprel)
        else
            @assert deprel == d["deprel"]
        end
    end

    feats = eval(parse("Flist."*args["feats"]))
    pt = eval(parse(args["arctype"]))
    s1 = data[1][1]
    p1 = pt(wcnt(s1),length(deprel))
    xrows=KUparser.flen(p1, s1, feats)
    yrows=p1.nmove

    net=newnet(relu, [xrows; args["hidden"]; yrows]...)
    net[end].f=KUnet.logp
    for k in [fieldnames(UpdateParam); :dropout]
        haskey(args, string(k)) || continue
        v = args[string(k)]
        if isempty(v)
            continue
        elseif length(v)==1
            setparam!(net, k, v[1])
        else
            @assert length(v)==length(net) "$k should have 1 or $(length(net)) elements"
            for i=1:length(v)
                setparam!(net[i], k, v[i])
            end
        end
    end
    @show net

    pa = Array(Vector{pt}, length(data))
    x = y = nothing
    for i=1:length(data)
        # Initialize training set using oparse on first corpus
        if i==1
            @date (pa[i],x,y) = oparse(pt, data[i], ndeps, feats)
            @show map(size, (pa[i], x, y))
        else
            @date pa[i] = oparse(pt, data[i], ndeps)
        end
        @show evalparse(pa[i], data[i])
    end
    accuracy = Array(Float32, length(data))
    (bestscore,bestepoch,epoch)=(0,0,0)

    while true
        @show epoch += 1
        @date train(net, x, y; batch=args["tbatch"], loss=KUnet.logploss, shuffle=true)
        if epoch % args["pepochs"] == 0
            for i=1:length(data)
                @show i
                @date reset!(pa[i])
                if (args["parser"] == "oparser")
                    # We never change the training set with oparser, just report accuracy
                    @date gparse(pa[i], data[i], ndeps, feats, net, args["pbatch"])
                elseif (args["parser"] == "gparser")
                    # The first corpus gives us the new training set
                    if i==1
                        @date gparse(pa[i], data[i], ndeps, feats, net, args["pbatch"], x, y)
                    else
                        @date gparse(pa[i], data[i], ndeps, feats, net, args["pbatch"])
                    end
                elseif (args["parser"] == "bparser")
                    # The first corpus gives us the new training set
                    if i==1
                        @date bparse(pa[i], data[i], ndeps, feats, net, args["nbeam"], args["pbatch"], x, y)
                    else
                        @date bparse(pa[i], data[i], ndeps, feats, net, args["nbeam"], args["pbatch"])
                    end
                else
                    error("Unknown parser "*args["parser"])
                end
                @show e = evalparse(pa[i], data[i])
                accuracy[i] = e[1]  # e[1] is UAS including punct
            end
            println("DATA:\t$epoch\t"*join(accuracy, '\t')); flush(STDOUT)

            # We look at accuracy[2] (dev score) by default,
            # accuracy[1] (trn score) if there is no dev.

            score = (length(accuracy) > 1 ? accuracy[2] : accuracy[1])

            # Update best score and save best net

            if (score > bestscore)
                @show (bestscore,bestepoch)=(score,epoch)
                !isempty(args["out"]) && save(args["out"], net)
            end

            # Quit if we have the minimum number of epochs and have
            # not made progress in the last half.

            (epoch >= args["epochs"]) && (epoch >= 2*bestepoch) && break

        end
    end
end

main()


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

# struct _enum_TExecuteMode
#   GPU::Int32
#   CPU::Int32
# end
# const TExecuteMode = _enum_TExecuteMode(Int32(1), Int32(2))

struct _enum_TFileType
  DELIMITED::Int32
  POLYGON::Int32
  PARQUET::Int32
end
const TFileType = _enum_TFileType(Int32(0), Int32(1), Int32(2))

struct _enum_TPartitionDetail
  DEFAULT::Int32
  REPLICATED::Int32
  SHARDED::Int32
  OTHER::Int32
end
const TPartitionDetail = _enum_TPartitionDetail(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_TGeoFileLayerContents
  EMPTY::Int32
  GEO::Int32
  NON_GEO::Int32
  UNSUPPORTED_GEO::Int32
end
const TGeoFileLayerContents = _enum_TGeoFileLayerContents(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_TImportHeaderRow
  AUTODETECT::Int32
  NO_HEADER::Int32
  HAS_HEADER::Int32
end
const TImportHeaderRow = _enum_TImportHeaderRow(Int32(0), Int32(1), Int32(2))

struct _enum_TRole
  SERVER::Int32
  AGGREGATOR::Int32
  LEAF::Int32
  STRING_DICTIONARY::Int32
end
const TRole = _enum_TRole(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_TMergeType
  UNION::Int32
  REDUCE::Int32
end
const TMergeType = _enum_TMergeType(Int32(0), Int32(1))

struct _enum_TExpressionRangeType
  INVALID::Int32
  INTEGER::Int32
  FLOAT::Int32
  DOUBLE::Int32
end
const TExpressionRangeType = _enum_TExpressionRangeType(Int32(0), Int32(1), Int32(2), Int32(3))

struct _enum_TDBObjectType
  AbstractDBObjectType::Int32
  DatabaseDBObjectType::Int32
  TableDBObjectType::Int32
  DashboardDBObjectType::Int32
  ViewDBObjectType::Int32
end
const TDBObjectType = _enum_TDBObjectType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

#const TRowDescriptor = Vector{TColumnType}

#const TTableDescriptor = Dict{String,TColumnType}

const TSessionId = String

const TKrb5Token = String

#const TRenderPassMap = Dict{Int32,TRawRenderPassDataResult}
const TQueryId = Int64

# const TRenderAggDataMap = Dict{String,Dict{String,Dict{String,Dict{String,Vector{TRenderDatum}}}}}

# mutable struct TDatumVal <: Thrift.TMsg
#   int_val::Int64
#   real_val::Float64
#   str_val::String
#   arr_val::Vector{TDatum}
#   TDatumVal() = (o=new(); fillunset(o); o)
# end # mutable struct TDatumVal

mutable struct TDatumVal <: Thrift.TMsg
  int_val::Int64
  real_val::Float64
  str_val::String
  arr_val::Vector{Any}
  TDatumVal() = (o=new(); fillunset(o); o)
end # mutable struct TDatumVal

#https://github.com/tanmaykm/Thrift.jl/issues/37#issuecomment-428024388
function Thrift.meta(t::Type{TDatumVal})
    ThriftMeta(t, [
        ThriftMetaAttribs(1, :int_val, 10, true, Any[], ThriftMeta[])
        ThriftMetaAttribs(2, :real_val, 4, true, Any[], ThriftMeta[])
        ThriftMetaAttribs(3, :str_val, 11, true, Any[], ThriftMeta[])
        ThriftMetaAttribs(4, :arr_val, 15, true, Any[], ThriftMeta[meta(Core.eval(Main, Meta.parse("TDatum")))])
    ])
end

mutable struct TDatum <: Thrift.TMsg
  val::TDatumVal
  is_null::Bool
  TDatum() = (o=new(); fillunset(o); o)
end # mutable struct TDatum

mutable struct TStringValue <: Thrift.TMsg
  str_val::String
  is_null::Bool
  TStringValue() = (o=new(); fillunset(o); o)
end # mutable struct TStringValue

mutable struct TColumnType <: Thrift.TMsg
  col_name::String
  col_type::TTypeInfo
  is_reserved_keyword::Bool
  src_name::String
  is_system::Bool
  is_physical::Bool
  col_id::Int64
  TColumnType() = (o=new(); fillunset(o); o)
end # mutable struct TColumnType

const TRowDescriptor = Vector{TColumnType} #move these after TColumnType definition

const TTableDescriptor = Dict{String,TColumnType}

mutable struct TRow <: Thrift.TMsg
  cols::Vector{TDatum}
  TRow() = (o=new(); fillunset(o); o)
end # mutable struct TRow

# mutable struct TColumnData <: Thrift.TMsg
#   int_col::Vector{Int64}
#   real_col::Vector{Float64}
#   str_col::Vector{String}
#   arr_col::Vector{TColumn}
#   TColumnData() = (o=new(); fillunset(o); o)
# end # mutable struct TColumnData
#
# mutable struct TColumn <: Thrift.TMsg
#   data::TColumnData
#   nulls::Vector{Bool}
#   TColumn() = (o=new(); fillunset(o); o)
# end # mutable struct TColumn

mutable struct TColumnData <: Thrift.TMsg
  int_col::Vector{Int64}
  real_col::Vector{Float64}
  str_col::Vector{String}
  arr_col::Vector{Any}
  TColumnData() = (o=new(); fillunset(o); o)
end # mutable struct TColumnData

mutable struct TColumn <: Thrift.TMsg
  data::TColumnData
  nulls::Vector{Bool}
  TColumn() = (o=new(); fillunset(o); o)
end # mutable struct TColumn

function Thrift.meta(t::Type{TColumnData})
    ThriftMeta(t, [
        ThriftMetaAttribs(1, :int_col, 15, false, Any[], ThriftMeta[])
        ThriftMetaAttribs(2, :real_col,15, false, Any[], ThriftMeta[])
        ThriftMetaAttribs(3, :str_col, 15, false, Any[], ThriftMeta[])
        ThriftMetaAttribs(4, :arr_col, 15, false, Any[], ThriftMeta[meta(Core.eval(Main, Meta.parse("TColumn")))])
    ])
end

mutable struct TStringRow <: Thrift.TMsg
  cols::Vector{TStringValue}
  TStringRow() = (o=new(); fillunset(o); o)
end # mutable struct TStringRow

mutable struct TKrb5Session <: Thrift.TMsg
  sessionId::TSessionId
  krbToken::TKrb5Token
  TKrb5Session() = (o=new(); fillunset(o); o)
end # mutable struct TKrb5Session

mutable struct TStepResult <: Thrift.TMsg
  serialized_rows::TSerializedRows
  execution_finished::Bool
  merge_type::Int32
  sharded::Bool
  row_desc::TRowDescriptor
  node_id::Int32
  TStepResult() = (o=new(); fillunset(o); o)
end # mutable struct TStepResult

mutable struct TRowSet <: Thrift.TMsg
  row_desc::TRowDescriptor
  rows::Vector{TRow}
  columns::Vector{TColumn}
  is_columnar::Bool
  TRowSet() = (o=new(); fillunset(o); o)
end # mutable struct TRowSet

mutable struct TQueryResult <: Thrift.TMsg
  row_set::TRowSet
  execution_time_ms::Int64
  total_time_ms::Int64
  nonce::String
  TQueryResult() = (o=new(); fillunset(o); o)
end # mutable struct TQueryResult

mutable struct TDataFrame <: Thrift.TMsg
  sm_handle::Vector{UInt8}
  sm_size::Int64
  df_handle::Vector{UInt8}
  df_size::Int64
  TDataFrame() = (o=new(); fillunset(o); o)
end # mutable struct TDataFrame

mutable struct TDBInfo <: Thrift.TMsg
  db_name::String
  db_owner::String
  TDBInfo() = (o=new(); fillunset(o); o)
end # mutable struct TDBInfo

mutable struct TMapDException <: Exception
  error_msg::String
  TMapDException() = (o=new(); fillunset(o); o)
end # mutable struct TMapDException

mutable struct TCopyParams <: Thrift.TMsg
  delimiter::String
  null_str::String
  has_header::Int32
  quoted::Bool
  _quote::String
  escape::String
  line_delim::String
  array_delim::String
  array_begin::String
  array_end::String
  threads::Int32
  file_type::Int32
  s3_access_key::String
  s3_secret_key::String
  s3_region::String
  geo_coords_encoding::Int32
  geo_coords_comp_param::Int32
  geo_coords_type::Int32
  geo_coords_srid::Int32
  sanitize_column_names::Bool
  geo_layer_name::String
  s3_endpoint::String
  TCopyParams() = (o=new(); fillunset(o); o)
end # mutable struct TCopyParams
meta(t::Type{TCopyParams}) = meta(t, Symbol[], Int[], Dict{Symbol,Any}(:has_header => Int32(0), :file_type => Int32(0), :geo_coords_encoding => Int32(6), :geo_coords_comp_param => Int32(32), :geo_coords_type => Int32(18), :geo_coords_srid => Int32(4326), :sanitize_column_names => true))

mutable struct TCreateParams <: Thrift.TMsg
  is_replicated::Bool
  TCreateParams() = (o=new(); fillunset(o); o)
end # mutable struct TCreateParams

mutable struct TDetectResult <: Thrift.TMsg
  row_set::TRowSet
  copy_params::TCopyParams
  TDetectResult() = (o=new(); fillunset(o); o)
end # mutable struct TDetectResult

mutable struct TImportStatus <: Thrift.TMsg
  elapsed::Int64
  rows_completed::Int64
  rows_estimated::Int64
  rows_rejected::Int64
  TImportStatus() = (o=new(); fillunset(o); o)
end # mutable struct TImportStatus

mutable struct TFrontendView <: Thrift.TMsg
  view_name::String
  view_state::String
  image_hash::String
  update_time::String
  view_metadata::String
  TFrontendView() = (o=new(); fillunset(o); o)
end # mutable struct TFrontendView

mutable struct TDashboard <: Thrift.TMsg
  dashboard_name::String
  dashboard_state::String
  image_hash::String
  update_time::String
  dashboard_metadata::String
  dashboard_id::Int32
  dashboard_owner::String
  is_dash_shared::Bool
  TDashboard() = (o=new(); fillunset(o); o)
end # mutable struct TDashboard

mutable struct TServerStatus <: Thrift.TMsg
  read_only::Bool
  version::String
  rendering_enabled::Bool
  start_time::Int64
  edition::String
  host_name::String
  poly_rendering_enabled::Bool
  role::Int32
  TServerStatus() = (o=new(); fillunset(o); o)
end # mutable struct TServerStatus

mutable struct TPixel <: Thrift.TMsg
  x::Int64
  y::Int64
  TPixel() = (o=new(); fillunset(o); o)
end # mutable struct TPixel

mutable struct TPixelTableRowResult <: Thrift.TMsg
  pixel::TPixel
  vega_table_name::String
  table_id::Vector{Int64}
  row_id::Vector{Int64}
  row_set::TRowSet
  nonce::String
  TPixelTableRowResult() = (o=new(); fillunset(o); o)
end # mutable struct TPixelTableRowResult

mutable struct TRenderResult <: Thrift.TMsg
  image::Vector{UInt8}
  nonce::String
  execution_time_ms::Int64
  render_time_ms::Int64
  total_time_ms::Int64
  vega_metadata::String
  TRenderResult() = (o=new(); fillunset(o); o)
end # mutable struct TRenderResult

mutable struct TGpuSpecification <: Thrift.TMsg
  num_sm::Int32
  clock_frequency_kHz::Int64
  memory::Int64
  compute_capability_major::Int16
  compute_capability_minor::Int16
  TGpuSpecification() = (o=new(); fillunset(o); o)
end # mutable struct TGpuSpecification

mutable struct THardwareInfo <: Thrift.TMsg
  num_gpu_hw::Int16
  num_cpu_hw::Int16
  num_gpu_allocated::Int16
  start_gpu::Int16
  host_name::String
  gpu_info::Vector{TGpuSpecification}
  THardwareInfo() = (o=new(); fillunset(o); o)
end # mutable struct THardwareInfo

mutable struct TClusterHardwareInfo <: Thrift.TMsg
  hardware_info::Vector{THardwareInfo}
  TClusterHardwareInfo() = (o=new(); fillunset(o); o)
end # mutable struct TClusterHardwareInfo

mutable struct TMemoryData <: Thrift.TMsg
  slab::Int64
  start_page::Int32
  num_pages::Int64
  touch::Int32
  chunk_key::Vector{Int64}
  buffer_epoch::Int32
  is_free::Bool
  TMemoryData() = (o=new(); fillunset(o); o)
end # mutable struct TMemoryData

mutable struct TNodeMemoryInfo <: Thrift.TMsg
  host_name::String
  page_size::Int64
  max_num_pages::Int64
  num_pages_allocated::Int64
  is_allocation_capped::Bool
  node_memory_data::Vector{TMemoryData}
  TNodeMemoryInfo() = (o=new(); fillunset(o); o)
end # mutable struct TNodeMemoryInfo

mutable struct TTableMeta <: Thrift.TMsg
  table_name::String
  num_cols::Int64
  is_view::Bool
  is_replicated::Bool
  shard_count::Int64
  max_rows::Int64
  table_id::Int64
  max_table_id::Int64
  col_types::Vector{TTypeInfo}
  col_names::Vector{String}
  TTableMeta() = (o=new(); fillunset(o); o)
end # mutable struct TTableMeta
meta(t::Type{TTableMeta}) = meta(t, Symbol[], Int[1,2,4,5,6,7,8,9,10,11], Dict{Symbol,Any}())

mutable struct TTableDetails <: Thrift.TMsg
  row_desc::TRowDescriptor
  fragment_size::Int64
  page_size::Int64
  max_rows::Int64
  view_sql::String
  shard_count::Int64
  key_metainfo::String
  is_temporary::Bool
  partition_detail::Int32
  TTableDetails() = (o=new(); fillunset(o); o)
end # mutable struct TTableDetails

mutable struct TColumnRange <: Thrift.TMsg
  _type::Int32
  col_id::Int32
  table_id::Int32
  has_nulls::Bool
  int_min::Int64
  int_max::Int64
  bucket::Int64
  fp_min::Float64
  fp_max::Float64
  TColumnRange() = (o=new(); fillunset(o); o)
end # mutable struct TColumnRange

mutable struct TDictionaryGeneration <: Thrift.TMsg
  dict_id::Int32
  entry_count::Int64
  TDictionaryGeneration() = (o=new(); fillunset(o); o)
end # mutable struct TDictionaryGeneration

mutable struct TTableGeneration <: Thrift.TMsg
  table_id::Int32
  tuple_count::Int64
  start_rowid::Int64
  TTableGeneration() = (o=new(); fillunset(o); o)
end # mutable struct TTableGeneration

mutable struct TPendingQuery <: Thrift.TMsg
  id::TQueryId
  column_ranges::Vector{TColumnRange}
  dictionary_generations::Vector{TDictionaryGeneration}
  table_generations::Vector{TTableGeneration}
  TPendingQuery() = (o=new(); fillunset(o); o)
end # mutable struct TPendingQuery

mutable struct TVarLen <: Thrift.TMsg
  payload::Vector{UInt8}
  is_null::Bool
  TVarLen() = (o=new(); fillunset(o); o)
end # mutable struct TVarLen

mutable struct TDataBlockPtr <: Thrift.TMsg
  fixed_len_data::Vector{UInt8}
  var_len_data::Vector{TVarLen}
  TDataBlockPtr() = (o=new(); fillunset(o); o)
end # mutable struct TDataBlockPtr
meta(t::Type{TDataBlockPtr}) = meta(t, Symbol[:fixed_len_data,:var_len_data], Int[], Dict{Symbol,Any}())

mutable struct TInsertData <: Thrift.TMsg
  db_id::Int32
  table_id::Int32
  column_ids::Vector{Int32}
  data::Vector{TDataBlockPtr}
  num_rows::Int64
  TInsertData() = (o=new(); fillunset(o); o)
end # mutable struct TInsertData

mutable struct TPendingRenderQuery <: Thrift.TMsg
  id::TQueryId
  TPendingRenderQuery() = (o=new(); fillunset(o); o)
end # mutable struct TPendingRenderQuery

mutable struct TRenderParseResult <: Thrift.TMsg
  merge_type::Int32
  node_id::Int32
  execution_time_ms::Int64
  render_time_ms::Int64
  total_time_ms::Int64
  TRenderParseResult() = (o=new(); fillunset(o); o)
end # mutable struct TRenderParseResult

mutable struct TRawRenderPassDataResult <: Thrift.TMsg
  num_pixel_channels::Int32
  num_pixel_samples::Int32
  pixels::Vector{UInt8}
  row_ids_A::Vector{UInt8}
  row_ids_B::Vector{UInt8}
  table_ids::Vector{UInt8}
  accum_data::Vector{UInt8}
  TRawRenderPassDataResult() = (o=new(); fillunset(o); o)
end # mutable struct TRawRenderPassDataResult

const TRenderPassMap = Dict{Int32,TRawRenderPassDataResult} #move after TRawRenderPassDataResult definition

mutable struct TRawPixelData <: Thrift.TMsg
  width::Int32
  height::Int32
  render_pass_map::TRenderPassMap
  TRawPixelData() = (o=new(); fillunset(o); o)
end # mutable struct TRawPixelData

mutable struct TRenderDatum <: Thrift.TMsg
  _type::Int32
  cnt::Int32
  value::Vector{UInt8}
  TRenderDatum() = (o=new(); fillunset(o); o)
end # mutable struct TRenderDatum

const TRenderAggDataMap = Dict{String,Dict{String,Dict{String,Dict{String,Vector{TRenderDatum}}}}} #moved to after definition

mutable struct TRenderStepResult <: Thrift.TMsg
  merge_data::TRenderAggDataMap
  raw_pixel_data::TRawPixelData
  execution_time_ms::Int64
  render_time_ms::Int64
  total_time_ms::Int64
  TRenderStepResult() = (o=new(); fillunset(o); o)
end # mutable struct TRenderStepResult

mutable struct TDatabasePermissions <: Thrift.TMsg
  create_::Bool
  delete_::Bool
  view_sql_editor_::Bool
  access_::Bool
  TDatabasePermissions() = (o=new(); fillunset(o); o)
end # mutable struct TDatabasePermissions

mutable struct TTablePermissions <: Thrift.TMsg
  create_::Bool
  drop_::Bool
  select_::Bool
  insert_::Bool
  update_::Bool
  delete_::Bool
  truncate_::Bool
  alter_::Bool
  TTablePermissions() = (o=new(); fillunset(o); o)
end # mutable struct TTablePermissions

mutable struct TDashboardPermissions <: Thrift.TMsg
  create_::Bool
  delete_::Bool
  view_::Bool
  edit_::Bool
  TDashboardPermissions() = (o=new(); fillunset(o); o)
end # mutable struct TDashboardPermissions

mutable struct TViewPermissions <: Thrift.TMsg
  create_::Bool
  drop_::Bool
  select_::Bool
  insert_::Bool
  update_::Bool
  delete_::Bool
  TViewPermissions() = (o=new(); fillunset(o); o)
end # mutable struct TViewPermissions

mutable struct TDBObjectPermissions <: Thrift.TMsg
  database_permissions_::TDatabasePermissions
  table_permissions_::TTablePermissions
  dashboard_permissions_::TDashboardPermissions
  view_permissions_::TViewPermissions
  TDBObjectPermissions() = (o=new(); fillunset(o); o)
end # mutable struct TDBObjectPermissions
meta(t::Type{TDBObjectPermissions}) = meta(t, Symbol[:database_permissions_,:table_permissions_,:dashboard_permissions_,:view_permissions_], Int[], Dict{Symbol,Any}())

mutable struct TDBObject <: Thrift.TMsg
  objectName::String
  objectType::Int32
  privs::Vector{Bool}
  grantee::String
  privilegeObjectType::Int32
  TDBObject() = (o=new(); fillunset(o); o)
end # mutable struct TDBObject

mutable struct TDashboardGrantees <: Thrift.TMsg
  name::String
  is_user::Bool
  permissions::TDashboardPermissions
  TDashboardGrantees() = (o=new(); fillunset(o); o)
end # mutable struct TDashboardGrantees

mutable struct TLicenseInfo <: Thrift.TMsg
  claims::Vector{String}
  TLicenseInfo() = (o=new(); fillunset(o); o)
end # mutable struct TLicenseInfo

mutable struct TSessionInfo <: Thrift.TMsg
  user::String
  database::String
  start_time::Int64
  is_super::Bool
  TSessionInfo() = (o=new(); fillunset(o); o)
end # mutable struct TSessionInfo

mutable struct TGeoFileLayerInfo <: Thrift.TMsg
  name::String
  contents::Int32
  TGeoFileLayerInfo() = (o=new(); fillunset(o); o)
end # mutable struct TGeoFileLayerInfo

abstract type MapDClientBase end


text:
### Turing interface

import .Turing
using .Turing: TypedVarInfo, tonamedtuple, decondition, logprior, logjoint, VarName
using .Turing.DynamicPPL: evaluate!!

export TuringMuseProblem

struct TuringMuseProblem{A<:AD.AbstractBackend, M<:Turing.Model} <: AbstractMuseProblem

    autodiff :: A
    model :: M
    model_for_prior
    x
    observed_vars
    latent_vars
    hyper_vars

end

@doc doc"""

    TuringMuseProblem(
        model;
        observed_vars=(:x,), latent_vars=(:z,), hyper_vars=(:θ,),
        autodiff=MuseInference.ForwardDiffBackend()
    )

Wrap a Turing model to be ready to pass to [`muse`](@ref).

The names of variables within the model which are observed, latent, or
hyperparameters, can be customized via keyword arguments which accept
tuples of names. E.g.,

```julia
@model function demo()
    σ ~ Normal(0, 3)
    z ~ MvNormal(zeros(512), exp(σ/2))
    w ~ MvNormal(z, 1)
    x ~ MvNormal(w, 1)
    y ~ MvNormal(w, 2)
    (;σ,z,w,x,y)
end
truth = demo()()
model = demo() | (;truth.x, truth.y)
prob = TuringMuseProblem(
    model, observed_vars=(:x,:y), latent_vars=(:z,:w), hyper_vars=(:σ,)
)
```

!!! note

    When defining Turing models to be used with MuseInference, the new-style definition
    of Turing models is required, where the random variables do not appear as arguments
    to the function. This is because internally, MuseInference needs
    to [`condition`](https://turinglang.github.io/DynamicPPL.jl/stable/#AbstractPPL.condition-Tuple{Model})
    your model on various variables.

The `autodiff` parameter should be either
`MuseInference.ForwardDiffBackend()` or
`MuseInference.ZygoteBackend()`, specifying which library to use for
automatic differenation.

"""
function TuringMuseProblem(
    model;
    observed_vars = (:x,),
    latent_vars = (:z,),
    hyper_vars = (:θ,),
    autodiff = AD.ForwardDiffBackend()
)

    x = ComponentVector(select(model.context.values, observed_vars))
    model = decondition(model)
    vars = _extract_vars(evaluate!!(model)[2])
    model_for_prior = model | map(zero, select(vars, observed_vars)) | map(zero, select(vars, latent_vars))
    TuringMuseProblem(autodiff, model, model_for_prior, x, observed_vars, latent_vars, hyper_vars)

end

standardizeθ(prob::TuringMuseProblem, θ::NamedTuple) = ComponentVector(θ)
standardizeθ(prob::TuringMuseProblem, θ::Number) = length(prob.hyper_vars) == 1 ? ComponentVector(;θ) : error("Invalid θ type for this problem.")

function logPriorθ(prob::TuringMuseProblem, θ)
    logprior(prob.model_for_prior, θ)
end

function ∇θ_logLike(prob::TuringMuseProblem, x, z, θ)
    first(AD.gradient(prob.autodiff, θ -> logjoint(prob.model, (;NamedTupleView(x)..., NamedTupleView(z)..., NamedTupleView(θ)...)), θ))
end

function logLike_and_∇z_logLike(prob::TuringMuseProblem, x, z, θ)
    first.(AD.value_and_gradient(prob.autodiff, z -> logjoint(prob.model, (;NamedTupleView(x)..., NamedTupleView(z)..., NamedTupleView(θ)...)), z))
end

function sample_x_z(prob::TuringMuseProblem, rng::AbstractRNG, θ)
    model = (prob.model | NamedTupleView(θ))
    vars = _extract_vars(evaluate!!(model,rng)[2])
    (;x=ComponentVector(select(vars, prob.observed_vars)), z=ComponentVector(select(vars, prob.latent_vars)))
end

# helped to extract parameters from a sampled model. feels like there
# should be a less hacky way to do this...
function _extract_vars(vi)
    map(TypedVarInfo(vi).metadata) do m
        if m.vns[1] isa VarName{<:Any,Setfield.IdentityLens} && length(m.vals)==1
            m.vals[1]
        else
            m.vals
        end
    end
end

muse!(result::MuseResult, model::Turing.Model, args...; kwargs...) = muse!(result, TuringMuseProblem(model), args...; kwargs...)
get_J!(result::MuseResult, model::Turing.Model, args...; kwargs...) = get_J!(result, TuringMuseProblem(model), args...; kwargs...)
get_H!(result::MuseResult, model::Turing.Model, args...; kwargs...) = get_H!(result, TuringMuseProblem(model), args...; kwargs...)


text:
using Base

"""
    smld_sub = isolatesmld(smld::SMLD2D, subind::Int)

Isolate the `smld` localizations specified by `subind` index.

# Description
This method grabs the localization indexed by `subind` from `smld` and outputs
a new SMLD2D structure containing only those localization.
"""
function isolatesmld(smld::SMLD2D, subind::Int)
    # Loop through each field and add it to the output SMLD.  For scalar
    # fields, we don't have to apply the indexing of `subind`.
    fields = fieldnames(SMLMData.SMLD2D)
    nfields = Base.length(fields)
    indlength = Base.length(subind)
    smld_sub = deepcopy(smld)
    for ii = 1:nfields
        if isdefined(smld, fields[ii])
            currentfield = getfield(smld, fields[ii])
            if isa(currentfield, Vector) &&
            any(fields[ii] .== smld.datafields) &&
            (Base.length(currentfield) >= indlength)
                # If this field is a vector, we'll keep only the `subind` element.
                setfield!(smld_sub, fields[ii], [currentfield[subind]])
            end
        end
    end

    return smld_sub
end

"""
    smld_sub = isolatesmld(smld::SMLD2D, subind::Vector{Int})

Isolate the `smld` localizations specified by `subind` indices.

# Description
This method grabs the localization indexed by `subind` from `smld` and outputs
a new SMLD2D structure containing only those localization.
"""
function isolatesmld(smld::SMLD2D, subind::Vector{Int})
    # Loop through each field and add it to the output SMLD.  For scalar
    # fields, we don't have to apply the indexing of `subind`.
    fields = fieldnames(SMLMData.SMLD2D)
    nfields = Base.length(fields)
    indlength = Base.length(subind)
    smld_sub = deepcopy(smld)
    for ii = 1:nfields
        if isdefined(smld, fields[ii])
            currentfield = getfield(smld, fields[ii])
            if isa(currentfield, Vector) &&
            any(fields[ii] .== smld.datafields) &&
            (Base.length(currentfield) >= indlength)
                # If this field is a vector, we'll keep only the `subind` element.
                setfield!(smld_sub, fields[ii], currentfield[subind])
            end
        end
    end

    return smld_sub
end

"""
smld_sub = isolatesmld(smld::SMLD2D, keepbit::BitVector)

Isolate the `smld` localizations specified by `keepbit` BitVector.

# Description
This method grabs the localizations requested by the BitVector `keepbit`,
which should have the same length as the number of localizations in `smld`.
"""
function isolatesmld(smld::SMLD2D, keepbit::BitVector)
    # Loop through each field and add it to the output SMLD.  For scalar
    # fields, we don't have to apply the indexing of `keepbit`.
    fields = fieldnames(SMLMData.SMLD2D)
    nfields = Base.length(fields)
    bitlength = Base.length(keepbit)
    smld_sub = deepcopy(smld)
    for ii = 1:nfields
        if isdefined(smld, fields[ii])
            currentfield = getfield(smld, fields[ii])
            if isa(currentfield, Vector) &&
            any(fields[ii] .== smld.datafields) &&
            (Base.length(currentfield) >= bitlength)
                # If this field is a vector, we'll keep only the `keepbit` elements.
                setfield!(smld_sub, fields[ii], currentfield[keepbit])
            end
        end
    end

    return smld_sub
end

"""
    smld_sub = isolatesmld(smld::SMLD2D, subind::UnitRange)

Isolate the `smld` localizations specified by `subind` indices.

# Description
This method grabs the localizations indexed by `subind` from `smld` and outputs
a new SMLD2D structure containing only those localization.
"""
function isolatesmld(smld::SMLD2D, subind::UnitRange)
    return SMLMData.isolatesmld(smld, collect(subind))
end

"""
    smld_connected, keepbool = isolateconnected(smld::SMLMData.SMLD2D)

Prepare individual SMLD2D structures for each precluster in `smld`.

# Description
Connected localizations in `smld` (based on their connectID field) are
reorganized into distinct SMLD2D structures, with each structure corresponding
to a single connectID.

# Inputs
-`smld`: SMLD2D structure containing the localization data, with the field
         connectID populated with a meaningful set of connection indices.

# outputs
-`smld_connected`: Vector of SMLD2D structures with each structure containing
                   localizations that share the same `connectID`.
-`keepbool`: Vector of BitVector defining the `smld` indices sharing the same
             `connectID`.
"""
function isolateconnected(smld::SMLMData.SMLD2D)
    # Ensure that the provided `smld` has a meaningful connectID field,
    # returning if not.
    nloc = Base.length(smld)
    if Base.length(smld.connectID) != nloc
        @warn "Input `smld.connectID` is not valid."
        return [smld]
    end

    # Loop through clusters in `smld` and prepare the output.
    connectIDs = unique(smld.connectID)
    smld_connected = Vector{SMLMData.SMLD2D}(undef, Base.length(connectIDs))
    keepbool = Vector{BitVector}(undef, Base.length(connectIDs))
    for cc in connectIDs
        keepbool[cc] = smld.connectID .== cc
        smld_connected[cc] = SMLMData.isolatesmld(smld, keepbool[cc])
    end

    return smld_connected, keepbool
end


text:
module OP2019

import sbp

using sbp.Operator: read_operator, to_matrix, OperatorData
using sbp.Grid:grid_xp, grid_xm

default_path = string(@__DIR__, "/")

struct Description
        Dp::OperatorData
        Hp::OperatorData
        Pp::OperatorData
        Dm::OperatorData
        Hm::OperatorData
        Pm::OperatorData
end

function load_description(;path=""::String, order=4::Int64, verbose=false::Bool)
        if path == ""
                path = default_path
        end

        if order==4
                desc_Dp = sbp.Operator.read_operator(string(path, "Dp_42.txt"))
                desc_Hp = sbp.Operator.read_operator(string(path, "Hp_42.txt"))
                desc_Pp = sbp.Operator.read_operator(string(path, "Pp_42.txt"))
                desc_Dm = sbp.Operator.read_operator(string(path, "Dm_42.txt"))
                desc_Hm = sbp.Operator.read_operator(string(path, "Hm_42.txt"))
                desc_Pm = sbp.Operator.read_operator(string(path, "Pm_42.txt"))
        else
           error("No operators implemented for order = ", order)
        end

        return Description(desc_Dp, desc_Hp, desc_Pp, desc_Dm, desc_Hm, desc_Pm)
end

function build_operators(desc, m::Int64, is_sparse::Bool=true)

        xm, h = grid_xm(m+1)
        xp, h = grid_xp(m)

        Dp = to_matrix(desc.Dp, m, m + 1, is_sparse) / h
        Dm = to_matrix(desc.Dm, m + 1, m, is_sparse) / h

        Hp = to_matrix(desc.Hp, m, m, is_sparse) * h
        Hm = to_matrix(desc.Hm, m + 1, m + 1, is_sparse) * h

        Pp = to_matrix(desc.Pp, m, m + 1, is_sparse)
        Pm = to_matrix(desc.Pm, m + 1, m, is_sparse)

        return xp, xm, h, Dp, Dm, Hp, Hm, Pp, Pm

end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# All matrices are DxN, where N is the number of positions and D is the dimensionality

# Here, P is the probe (to be rotated) and Q is the refereence
# https://en.wikipedia.org/wiki/Kabsch_algorithm
# This has been generalized to support weighted points:
# https://igl.ethz.ch/projects/ARAP/svd_rot.pdf

# assuming P and Q are already centered at the origin
# returns the rotation for alignment
function kabsch_centered(P,Q,w)
    @assert size(P) == size(Q)
    W = diagm(w/sum(w)) # here, the weights are assumed to sum to 1
    H = P*W*Q'
    D = Matrix{Float64}(I,size(H,1), size(H,2))
    U,Σ,V = svd(H)
    D[end] = sign(det(V*U'))
    return LinearMap(V * D * U')
end

# transform DxN matrices
function (tform::Translation)(A::AbstractMatrix)
    return hcat([tform(A[:,i]) for i=1:size(A,2)]...)
end

# P and Q are not necessarily centered
# returns the transformation for alignment
function kabsch(P, Q, w::AbstractVector=ones(size(P,2)))
    @assert !any(w.<0) && sum(w)>0
    wn = w/sum(w) # weights should sum to 1 for computing the centroid
    centerP, centerQ = center_translation(P,wn), center_translation(Q,wn)
    R = kabsch_centered(centerP(P), centerQ(Q), wn)
    return inv(centerQ) ∘ R ∘ centerP
end

function kabsch(P,Q,matches::AbstractVector{<:Tuple{Int,Int}},wp=ones(size(P,2)),wq=ones(size(Q,2)))
    matchedP, matchedQ = matched_points(P,Q,matches)
    w = [wp[i]*wq[j] for (i,j) in matches]
    return kabsch(matchedP, matchedQ, w)
end

# centroid of positions in A, weighted by weights in w (assumed to sum to 1)
centroid(A, w=fill(1/size(A,2), size(A,2))) = A*w

# translation moving centroid to origin
center_translation(A, w=fill(1/size(A,2), size(A,2))) = Translation(-centroid(A,w))


text:
#=
Simple 2D ellipse with convex combination initial state.
=#
using Reachability, Plots

function compute(input_options::Pair{Symbol,<:Any}...)
    # =====================
    # Problem specification
    # =====================
    A = sparse([1, 1, 2, 2], [1, 2, 1, 2], [3., -9., 4., -3.])

    # initial set
    X0_1 = Singleton([1.0, 0.0])
    X0_2 = Singleton([1.5, 0.0])
    alpha = 0.8
    Z0 = Singleton(alpha * X0_1.element + (1.-alpha) * X0_2.element)

    # instantiate continuous LTI system
    S1 = ContinuousSystem(A, X0_1)
    S2 = ContinuousSystem(A, X0_2)
    S3 = ContinuousSystem(A, Z0)

    # ===============
    # Problem solving
    # ===============
    # define solver-specific options
    options = Options(
        :δ => 0.01,
        :approx_model => "nobloating",
        :plot_vars => [1, 2]
        )

    t1 = 1.2
    t2 = 1.0
    t3 = 1.09

    o1 = merge(options, Options(:T => t1))
    o2 = merge(options, Options(:T => t2))
    o3 = merge(options, Options(:T => t3))

    # solve
    result1 = solve(S1, o1)
    result2 = solve(S2, o2)
    result3 = solve(S3, o3)

    # ========
    # Plotting
    # ========
    plot(result1,color="red")
    plot!(result2,color="blue")
    plot!(result3,color="green")
    @eval(savefig(@filename_to_png))
end # function

compute()

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module UnstructuredGridsTests

using Test
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Io
using Gridap.TensorValues

# Unstructured grid from raw data

trian = GridMock()

node_coordinates = get_node_coordinates(trian)
cell_node_ids = get_cell_node_ids(trian)
reffes = get_reffes(trian)
cell_types = get_cell_type(trian)
facet_normal = fill(VectorValue(0,1),num_cells(trian)) #dummy value

grid = UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types,Oriented(),facet_normal)
test_grid(grid)
@test is_oriented(grid) == true
@test get_facet_normal(grid) === facet_normal

grid = UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types,Oriented())
test_grid(grid)
@test is_oriented(grid) == true

grid = UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types)
test_grid(grid)

cell_map1 = get_cell_map(grid)
cell_map2 = get_cell_map(grid)
@test cell_map1 === cell_map2

q1i = Point(0.5,0.5)
np1 = 4
q1 = fill(q1i,np1)

q2i = Point(0.25,0.25)
np2 = 3
q2 = fill(q2i,np2)

q = CompressedArray([q1,q2],get_cell_type(grid))

cell_map = get_cell_map(grid)
x = lazy_map(evaluate,cell_map,q)

x1i = Point(0.5, 0.5)
x2i = Point(1.25, 0.25)
x3i = Point(1.75, 0.5)
x4i = Point(0.5, 1.5)
x5i = Point(1.5, 1.5)
x1 = fill(x1i,np1)
x2 = fill(x2i,np2)
x3 = fill(x3i,np2)
x4 = fill(x4i,np1)
x5 = fill(x5i,np1)
x = [x1,x2,x3,x4,x5]

# UnstructuredGrid from Grid

grid = UnstructuredGrid(trian)
test_grid(grid)
@test grid === UnstructuredGrid(grid)

# UnstructuredGrid from LagrangianRefFE

quad8 = LagrangianRefFE(Float64,QUAD,2)
grid = UnstructuredGrid(quad8)
@test num_nodes(grid) == num_nodes(quad8)

# UnstructuredGrid from Polytope

grid = UnstructuredGrid(ReferenceFE{2},WEDGE)
@test num_cells(grid) == 5
@test num_cell_dims(grid) == 2
@test num_point_dims(grid) == 3

grid = UnstructuredGrid(ReferenceFE{3},WEDGE)
@test num_cells(grid) == 1
@test num_cell_dims(grid) == 3
@test num_point_dims(grid) == 3

#UnstructuredGrid from coordinate vector

grid = UnstructuredGrid(get_node_coordinates(quad8))
@test num_nodes(grid) == num_nodes(quad8)
@test num_cells(grid) == num_nodes(quad8)
@test num_cell_dims(grid) == 0
@test num_point_dims(grid) == 2

## get low dim grid
#
#grid1 = UnstructuredGrid(ReferenceFE{1},trian)
#@test num_cell_dims(grid1) == 1

grid =  UnstructuredGrid(GridMock())
dict = to_dict(grid)
g = from_dict(UnstructuredGrid,dict)
test_grid(g)

s = to_json(grid)
g = from_json(UnstructuredGrid,s)
test_grid(g)

d = mktempdir()
f = joinpath(d,"grid.jld2")

to_jld2_file(grid,f)
grid2 = from_jld2_file(typeof(grid),f)

@test grid.cell_node_ids == grid2.cell_node_ids
@test grid.cell_types == grid2.cell_types
@test grid.node_coordinates == grid2.node_coordinates
@test grid.reffes == grid2.reffes

rm(d,recursive=true)

end # module


text:
using PyCall, Conda

# try importing
# catch the installation
# 1) Julia built-in miniconda
# 2) global pip installer

#try
using Conda
#    Conda.add_channel("conda-forge")
#    Conda.add("meshio")
Conda.add("scipy")
#    pyimport("meshio")
#    pyimport("scipy")
#catch
cmd = `pip3 install meshio scipy --user`
run(cmd)
#    pyimport("meshio")
#    pyimport("scipy")
#end


text:
using Base.Rounding: setrounding_raw, to_fenv
using Printf

function compare_calc_raw(op, updown, calc, raw, args...)
    if isequal(calc, raw)
        true
    else
        @info("Erorr", op, updown)
        for (i, v) in enumerate(args...)
            @info(@sprintf("a%d = %0.18e, bit rep : %s", i, v, bitstring(v)))
        end

        @info(@sprintf("calc = %0.18e, bit rep : %s", calc, bitstring(calc)))
        @info(@sprintf("raw = %0.18e, bit rep : %s", raw, bitstring(raw)))
        false
    end
end

function rounding_check(op, base_op, arrays...)
    elt = eltype(first(arrays))
    setrounding_raw(elt, to_fenv(RoundNearest))

    @eval begin
        up_calc = broadcast($(Symbol(op, "_up")), $(arrays...))
        down_calc = broadcast($(Symbol(op, "_down")), $(arrays...))

        setrounding_raw($elt, to_fenv(RoundUp))
        up_raw = broadcast($base_op, $(arrays...))

        setrounding_raw($elt, to_fenv(RoundDown))
        down_raw = broadcast($base_op, $(arrays...))
    end

    # Compare
    for (calc, raw, args) in zip(up_calc, up_raw, zip(arrays...))
        @test compare_calc_raw(op, "up", calc, raw, args)
    end

    for (calc, raw, args) in zip(down_calc, down_raw, zip(arrays...))
        @test compare_calc_raw(op, "down", calc, raw, args)
    end

    setrounding_raw(elt, to_fenv(RoundNearest))
end

rounding_check_unary(a::AbstractVector) = rounding_check(:sqrt, :sqrt, a)
rounding_check_unary(a) = rounding_check_unary([a])

function rounding_check_binary(a::T, b::T) where {T<:AbstractVector}
    for (op, base_op) in zip((:add, :sub, :mul, :div), (:+, :-, :*, :/))
        rounding_check(op, base_op, a, b)
    end
end
rounding_check_binary(a, b) = rounding_check_binary([a], [b])

special_value_list(T::Type) = [
    zero(T), -zero(T),                                  # 0.0, -0.0
    one(T), -one(T),                                    # 1.0, -1.0
    nextfloat(zero(T)), prevfloat(zero(T)),             # N_min^s, -N_min^s
    prevfloat(floatmin(T)), nextfloat(-floatmin(T)),    # N_max^s, -N_max^s
    floatmin(T), -floatmin(T),                          # N_min^n, -N_min^n
    floatmax(T), -floatmax(T),                          # N_max^n, -N_max^n
    eps(T), -eps(T),                                    # machine epsilon
    typemax(T), typemin(T),                             # Inf, -Inf
    T(NaN)                                              # NaN
]

for T in (Float64, Float32)
    @testset "$(T), Special Cases" begin
        special_values = special_value_list(T)
        len = Base.length(special_values)
        a = repeat(special_values, len)
        b = sort(a)
        rounding_check_unary(filter(x->x ≥ zero(x), special_values)) # sqrt
        rounding_check_binary(a, b)
    end
end

for n in 3:6
    N = 10^n
    for T in (Float64, Float32)
        @testset "$(T), Random Sampling, 10^$(n)" begin
            rand_a = reinterpret.(T, rand(Base.uinttype(T), N))
            rand_b = reinterpret.(T, rand(Base.uinttype(T), N))
            rounding_check_unary(abs.(rand_a))
            rounding_check_unary(abs.(rand_b))
            rounding_check_binary(rand_a, rand_b)
            rounding_check_binary(rand_b, rand_a)
        end
    end
end


text:
# Used for when a lazy version should be constructed on materialize
abstract type AbstractQuasiArrayApplyStyle <: ApplyStyle end
struct QuasiArrayApplyStyle <: AbstractQuasiArrayApplyStyle end

##
# LazyQuasiArray
##

"""
A `LazyQuasiArray` is an abstract type for any array that is
stored and should be manipulated lazily.
"""
abstract type LazyQuasiArray{T,N} <: AbstractQuasiArray{T,N} end

const LazyQuasiVector{T} = LazyQuasiArray{T,1}
const LazyQuasiMatrix{T} = LazyQuasiArray{T,2}

struct LazyQuasiArrayApplyStyle <: AbstractQuasiArrayApplyStyle end

struct QuasiLazyLayout <: AbstractLazyLayout end

MemoryLayout(::Type{<:LazyQuasiArray}) = QuasiLazyLayout()
ndims(M::Applied{LazyQuasiArrayApplyStyle,typeof(*)}) = ndims(last(M.args))

combine_mul_styles(::QuasiLazyLayout) = LazyQuasiArrayApplyStyle()
result_mul_style(::LazyQuasiArrayApplyStyle, ::LazyQuasiArrayApplyStyle) = LazyQuasiArrayApplyStyle()
result_mul_style(::LazyQuasiArrayApplyStyle, ::MulAddStyle) = LazyQuasiArrayApplyStyle()
result_mul_style(::MulAddStyle, ::LazyQuasiArrayApplyStyle) = LazyQuasiArrayApplyStyle()
result_mul_style(::LazyQuasiArrayApplyStyle, ::LazyArrayApplyStyle) = LazyQuasiArrayApplyStyle()
result_mul_style(::LazyArrayApplyStyle, ::LazyQuasiArrayApplyStyle) = LazyQuasiArrayApplyStyle()
result_mul_style(::LazyQuasiArrayApplyStyle, _) = LazyQuasiArrayApplyStyle()
result_mul_style(_, ::LazyQuasiArrayApplyStyle) = LazyQuasiArrayApplyStyle()

###
# ApplyQuasiArray
###

"""
A `ApplyQuasiArray` is a lazy realization of a function that
creates a quasi-array.
"""
struct ApplyQuasiArray{T, N, F, Args<:Tuple} <: LazyQuasiArray{T,N}
    f::F
    args::Args
end

const ApplyQuasiVector{T, F, Args<:Tuple} = ApplyQuasiArray{T, 1, F, Args}
const ApplyQuasiMatrix{T, F, Args<:Tuple} = ApplyQuasiArray{T, 2, F, Args}

QuasiLazyArray(A::Applied) = ApplyQuasiArray(A)

ApplyQuasiArray{T,N,F,Args}(M::Applied) where {T,N,F,Args} = ApplyQuasiArray{T,N,F,Args}(M.f, M.args)
ApplyQuasiArray{T,N}(M::Applied{Style,F,Args}) where {T,N,Style,F,Args} = ApplyQuasiArray{T,N,F,Args}(instantiate(M))
ApplyQuasiArray{T}(M::Applied) where {T} = ApplyQuasiArray{T,ndims(M)}(M)
ApplyQuasiArray(M::Applied) = ApplyQuasiArray{eltype(M)}(M)
ApplyQuasiVector(M::Applied) = ApplyQuasiVector{eltype(M)}(M)
ApplyQuasiMatrix(M::Applied) = ApplyQuasiMatrix{eltype(M)}(M)

ApplyQuasiArray(f, factors...) = ApplyQuasiArray(applied(f, factors...))
ApplyQuasiArray{T}(f, factors...) where T = ApplyQuasiArray{T}(applied(f, factors...))
ApplyQuasiArray{T,N}(f, factors...) where {T,N} = ApplyQuasiArray{T,N}(applied(f, factors...))

ApplyQuasiVector(f, factors...) = ApplyQuasiVector(applied(f, factors...))
ApplyQuasiMatrix(f, factors...) = ApplyQuasiMatrix(applied(f, factors...))

@inline Applied(A::AbstractQuasiArray) = Applied(call(A), arguments(A)...)
@inline ApplyQuasiArray(A::AbstractQuasiArray) = ApplyQuasiArray(call(A), arguments(A)...)

axes(A::ApplyQuasiArray) = axes(Applied(A))
size(A::ApplyQuasiArray) = map(length, axes(A))
copy(A::ApplyQuasiArray) = copy(Applied(A))

IndexStyle(::ApplyQuasiArray{<:Any,1}) = IndexLinear()

@propagate_inbounds getindex(A::ApplyQuasiArray{T,N}, kj::Vararg{Number,N}) where {T,N} =
    Applied(A)[kj...]

MemoryLayout(M::Type{ApplyQuasiArray{T,N,F,Args}}) where {T,N,F,Args} =
    applylayout(F, tuple_type_memorylayouts(Args)...)

copy(A::Applied{LazyQuasiArrayApplyStyle}) = ApplyQuasiArray(A)
copy(A::Applied{<:AbstractQuasiArrayApplyStyle}) = QuasiArray(A)
QuasiArray(A::Applied) = QuasiArray(ApplyQuasiArray(A))

####
# broadcasting
###

struct LazyQuasiArrayStyle{N} <: AbstractQuasiArrayStyle{N} end
LazyQuasiArrayStyle(::Val{N}) where N = LazyQuasiArrayStyle{N}()
LazyQuasiArrayStyle{M}(::Val{N}) where {N,M} = LazyQuasiArrayStyle{N}()
quasisubbroadcaststyle(::LazyQuasiArrayStyle{N}, _) where N = LazyQuasiArrayStyle{N}()
subbroadcaststyle(::LazyQuasiArrayStyle{N}, _) where N = LazyArrayStyle{N}()
subbroadcaststyle(::AbstractQuasiArrayStyle{N}, _) where N = DefaultArrayStyle{N}()

struct BroadcastQuasiArray{T, N, F, Args} <: LazyQuasiArray{T, N}
    f::F
    args::Args
end

const BroadcastQuasiVector{T,F,Args} = BroadcastQuasiArray{T,1,F,Args}
const BroadcastQuasiMatrix{T,F,Args} = BroadcastQuasiArray{T,2,F,Args}

LazyQuasiArray(bc::Broadcasted) = BroadcastQuasiArray(bc)

BroadcastQuasiArray{T,N,F,Args}(bc::Broadcasted) where {T,N,F,Args} = BroadcastQuasiArray{T,N,F,Args}(bc.f,bc.args)
BroadcastQuasiArray{T,N}(bc::Broadcasted{Style,Axes,F,Args}) where {T,N,Style,Axes,F,Args} = BroadcastQuasiArray{T,N,F,Args}(bc.f,bc.args)
BroadcastQuasiArray{T}(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Tuple{Vararg{Any,N}},<:Any,<:Tuple}) where {T,N} =
    BroadcastQuasiArray{T,N}(bc)

_broadcast2broadcastarray(a, b...) = tuple(a, b...)
_broadcast2broadcastarray(a::Broadcasted{<:LazyQuasiArrayStyle}, b...) = tuple(BroadcastQuasiArray(a), b...)

_BroadcastQuasiArray(bc::Broadcasted) = BroadcastQuasiArray{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastQuasiArray(bc::Broadcasted{S}) where S =
    _BroadcastQuasiArray(instantiate(Broadcasted{S}(bc.f, _broadcast2broadcastarray(bc.args...))))
BroadcastQuasiArray(b::BroadcastQuasiArray) = b
BroadcastQuasiArray(f, A, As...) = BroadcastQuasiArray(broadcasted(f, A, As...))

Broadcasted(A::BroadcastQuasiArray) = instantiate(broadcasted(A.f, A.args...))

axes(A::BroadcastQuasiArray) = axes(Broadcasted(A))
size(A::BroadcastQuasiArray) = map(length, axes(A))

IndexStyle(::BroadcastQuasiArray{<:Any,1}) = IndexLinear()

@propagate_inbounds getindex(A::BroadcastQuasiArray, kj::Number...) = Broadcasted(A)[kj...]

@propagate_inbounds _broadcast_getindex_range(A::Union{Ref,AbstractQuasiArray{<:Any,0},Number}, I) = A[] # Scalar-likes can just ignore all indices
# Everything else falls back to dynamically dropping broadcasted indices based upon its axes
@propagate_inbounds _broadcast_getindex_range(A, I) = A[I]

getindex(B::BroadcastQuasiArray{<:Any,1}, kr::AbstractVector{<:Number}) =
    BroadcastArray(B.f, map(a -> _broadcast_getindex_range(a,kr), B.args)...)

copy(bc::Broadcasted{<:LazyQuasiArrayStyle}) = BroadcastQuasiArray(bc)

BroadcastStyle(::Type{<:LazyQuasiArray{<:Any,N}}) where N = LazyQuasiArrayStyle{N}()

MemoryLayout(M::Type{BroadcastQuasiArray{T,N,F,Args}}) where {T,N,F,Args} =
    broadcastlayout(F, tuple_type_memorylayouts(Args)...)

###
# sub of *
###
call(a::AbstractQuasiArray) = call(MemoryLayout(typeof(a)), a)
call(::ApplyLayout{typeof(*)}, V::SubQuasiArray) = *

arguments(a::AbstractQuasiArray) = arguments(MemoryLayout(typeof(a)), a)
arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2}) = _mat_mul_arguments(V)
arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,1}) = _vec_mul_arguments(V)


text:
# test/test_all_mirt.jl

export test_all_mirt

using Test: @testset, detect_ambiguities

srclist = (
"../data",
"algorithm",
"fbp",
"io",
"mri",
"nufft",
"plot",
"regularize",
"system",
"utility",
)

"""
`test_all_mirt()`

Run all MIRT tests.
To avoid plotting during tests, first do `default(show=false)`
"""
function test_all_mirt()

	basedir = dirname(pathof(MIRT))
#	@testset "all MIRT" begin
		for root in srclist
		#	@show root
			@testset "MIRT:$root" begin
				tmp = joinpath(basedir, root, "z-test.jl")
			#	printstyled(tmp * "\n", color=:blue, bold=true)
				include(tmp)
			end
		end
#	end

	@test length(detect_ambiguities(MIRT)) == 0
	true
end


text:
module Structure
using  ...PatternStructure.Trees
export MetaMethod, MetaMethodTable, MetaMethodError, unlabeled

struct MetaMethod
  label   :: Symbol
  matcher :: Function
  method  :: Function
  tree    :: PatternRoot
end

struct MetaMethodTable
  name    :: AbstractString
  labels  :: Dict{Symbol, MetaMethod}
  methods :: Vector{MetaMethod}

  MetaMethodTable(name) =
    new(name, Dict{Symbol, MetaMethod}(), MetaMethod[])
end

struct MetaMethodError <: Exception
  name :: AbstractString
  expr :: Any
end

const unlabeled = gensym("unlabeled")

end


text:
"""
```
```

The dictionary `diffop` map each variable to a number of derivative operators,
each of which should be applied to the `diffvar` dictionary.

I should also create another differentiate method that uses the same operators
for all the variables (so a higher level differentiate function).

You trigger this by adding a setting to your model object.

Use derivative order of the derivative operator to infer
the correct symbol name.
"""
function differentiate!(stategrid::StateGrid, diffvars::AbstractDict{Symbol, Vector{S}}, derivs::AbstractDict{Symbol, Vector{S}},
                        drift::Vector{S} = Vector{S}(undef, 0);
                        L₁ = nothing, L₂ = nothing, populate::Bool = false, uniform::Bool = false,
                        skipvar::Vector{Symbol} = Vector{Symbol}(undef, 0)) where {S <: Real}

    # Initialize finite difference operators with homogeneous Neumann boundary conditions
    if isnothing(L₁)
        drift = isempty(drift) ? vcat(ones(length(stategrid) - 1), -1.) : sign.(drift)

        # Default to a first-order forward finite difference w/backward finite difference at the end
        dx = uniform ? diff(values(stategrid.x)[1][1:2]) :
            vcat(values(stategrid.x[1])[1], diff(values(stategrid.x)[1]), diff(values(stategrid.x)[1][end - 1:end]))
        L₁ = if uniform
            UpwindDifference(1, 1, dx, length(stategrid), drift) *
                RobinBC((0., 1., 0.), (0., 1., 0.), dx)
        else
            UpwindDifference(1, 1, dx, length(stategrid), drift) *
                RobinBC((0., 1., 0.), (0., 1., 0.), dx)
        end
    else
        drift = isempty(drift) ? sign.(L₁.L.coefficients) : sign.(drift)
    end
    if isnothing(L₂)
        # Default to a second-order centered finite difference
        dx = uniform ? diff(values(stategrid.x)[1][1:2]) :
            vcat(stategrid.x[1], diff(values(stategrid.x)[1]), diff(stategrid.x[end - 1:end]))
        L₂ = CenteredDifference(2, 2, dx, length(stategrid)) * RobinBC((0., 1., 0.), (0., 1., 0.), dx)
    end

    # Differentiate all requested derivatives
    state_name = keys(stategrid.x)[1] # derivs maps to a vector ⇒ one-dimensional model
    for (k, v) in diffvars
        if k in skipvar
            continue
        end
        deriv1 = Symbol(:∂, k, :_∂, state_name)
        deriv2 = Symbol(:∂², k, :_∂, state_name, Symbol("²"))

        # Check if first derivative is in derivs or requested by populate
        if haskey(derivs, deriv1)
            derivs[deriv1] .= (drift * L₁) * v
        elseif populate
            derivs[deriv1]  = (drift * L₁) * v
        end

        # Check if second derivative is in derivs or requested by populate
        if haskey(derivs, deriv2)
            derivs[deriv2] .= L₂ * v
        elseif populate
            derivs[deriv2]  = L₂ * v
        end
    end
end

#=# Add another version of the above function with a dict that which maps variables to drifts that informs the DiffEqOperator
function differentiate(stategrid::StateGrid, diffvar::AbstractDict{Symbol, Vector{S}},
                       diffop::AbstractDict{Symbol, Vector{DiffEqOperators.AbstractDerivativeOperator}}) where {S <: Real}
    for (k, v) in diffop
        for state in keys(stategrid.x)
	        diffvar[Symbol(:∂, k :_∂, state)] = diffop[:A] * v
        end
    end
end
=#

"""
```
function differentiate(x::AbstractVector{S}, y::AbstractVector{S}) where {S <: Real}
```

first-order differentiates a variable using second-order central differences in the interior
and first-order forward/backward differences on the boundaries.
"""
function differentiate(x::AbstractVector{S}, y::AbstractVector{S}) where {S <: Real}

    # Set up
    N        = length(x)
    y_left   = y[1:(N - 2)]
    y_right  = y[3:N]
    y_middle = y[2:(N - 1)]
    x_left   = x[1:(N - 2)]
    x_right  = x[3:N]
    x_middle = x[2:(N - 1)]
    yx       = Vector{S}(undef, N)

    # Differentiate
    yx[2:(N - 1)] = (x_right - x_middle) ./ (x_right-x_left) .* (y_middle - y_left) ./ (x_middle - x_left) +
        (x_middle - x_left) ./ (x_right-x_left) .* (y_right - y_middle) ./ (x_right - x_middle)
    yx[1]   = (y[2] - y[1]) / (x[2] - x[1] );
    yx[end] = (y[end] - y[end - 1] ) / ( x[end] - x[end - 1])

    return yx
end


text:
#=
Defining various quadrature routines.

Based on the quadrature routines found in the CompEcon toolbox by
Miranda and Fackler.

@author: Spencer Lyon

@date: 2014-08-15

References
----------
Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
and Finance, MIT Press, 2002.
=#

const qnw_func_notes = """
##### Notes

If any of the parameters to this function are scalars while others are
`Vector`s of length `n`, the the scalar parameter is repeated `n` times.
"""

const qnw_returns = """
##### Returns

- `nodes::Array{Float64}` : An array of quadrature nodes
- `weights::Array{Float64}` : An array of corresponding quadrature weights
"""

const qnw_refs = """
##### References

Miranda, Mario J, and Paul L Fackler. Applied Computational Economics and
Finance, MIT Press, 2002.
"""

"""
Computes multivariate Guass-Legendre  quadrature nodes and weights.

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwlege

## 1d versions
function qnwlege(n::Int, a::Real, b::Real)
    maxit = 10000
    m = fix((n + 1) / 2.0)
    xm = 0.5 * (b+a)
    xl = 0.5 * (b-a)
    nodes = zeros(n)

    weights = copy(nodes)
    i = 1:m

    z = cos(pi * (i - 0.25) ./ (n + 0.5))

    # allocate memory for loop arrays
    p3 = similar(z)
    pp = similar(z)

    its = 0
    for its=1:maxit
        p1 = ones(z)
        p2 = zeros(z)
        for j=1:n
            p3 = p2
            p2 = p1
            p1 = ((2*j-1)*z.*p2-(j-1)*p3)./j
        end

        pp = n*(z.*p1-p2)./(z.*z-1)
        z1 = z
        z = z1 - p1./pp

        err = Base.maxabs(z - z1)
        if err < 1e-14
            break
        end
    end

    if its == maxit
        error("Maximum iterations in _qnwlege1")
    end

    nodes[i] = xm - xl * z
    nodes[n+1-i] = xm + xl * z

    weights[i] = 2*xl./((1-z.*z).*pp.*pp)
    weights[n+1-i] = weights[i]

    return nodes, weights
end

"""
Computes multivariate Guass-Checbychev quadrature nodes and weights.

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwcheb

function qnwcheb(n::Int, a::Real, b::Real)
    nodes = (b+a)/2 - (b-a)/2 .* cos(pi/n .* (0.5:(n-0.5)))
    weights = ((b-a)/n) .* (cos(pi/n .* (collect(1:n)-0.5)*collect(0:2:n-1)')
                            *vcat(1, -2./(collect(1:2:n-2).*collect(3:2:n))))
    return nodes, weights
end

"""
Computes nodes and weights for multivariate normal distribution

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `mu::Union{Real, Vector{Real}}` : Mean along each dimension
- `sig2::Union{Real, Vector{Real}, Matrix{Real}}(eye(length(n)))` : Covariance
structure

$(qnw_returns)

##### Notes

This function has many methods. I try to describe them here.

`n` or `mu` can be a vector or a scalar. If just one is a scalar the other is
repeated to match the length of the other. If both are scalars, then the number
of repeats is inferred from `sig2`.

`sig2` can be a matrix, vector or scalar. If it is a matrix, it is treated as
the covariance matrix. If it is a vector, it is considered the diagonal of a
diagonal covariance matrix. If it is a scalar it is repeated along the diagonal
as many times as necessary, where the number of repeats is determined by the
length of either n and/or mu (which ever is a vector).

If all 3 are scalars, then 1d nodes are computed. `mu` and `sig2` are treated as
the mean and variance of a 1d normal distribution

$(qnw_refs)
"""
:qnwnorm

function qnwnorm(n::Int)
    maxit = 100
    pim4 = 1 / pi^(0.25)
    m = floor(Int, (n + 1) / 2)
    nodes = zeros(n)
    weights = zeros(n)

    z = sqrt(2n+1) - 1.85575 * ((2n+1).^(-1/6))

    for i=1:m
        # Reasonable starting values for root finding
        if i == 1
            z = sqrt(2n+1) - 1.85575 * ((2n+1).^(-1/6))
        elseif i == 2
            z = z - 1.14 * (n.^0.426)./z
        elseif i == 3
            z = 1.86z + 0.86nodes[1]
        elseif i == 4
            z = 1.91z + 0.91nodes[2]
        else
            z = 2z + nodes[i-2]
        end

        # root finding iterations
        it = 0
        pp = 0.0  # initialize pp so it is available outside while
        while it < maxit
            it += 1
            p1 = pim4
            p2 = 0

            for j=1:n
                p3 = p2
                p2 = p1
                p1 = z .* sqrt(2/j) .*p2 - sqrt((j-1)/j).*p3
            end

            pp = sqrt(2n).*p2
            z1 = z
            z = z1 - p1./pp

            if abs(z - z1) < 1e-14
                break
            end
        end

        if it >= maxit
            error("Failed to converge in qnwnorm")
        end

        nodes[n+1-i] = z
        nodes[i] = -z
        weights[i] = 2 ./ (pp.*pp)
        weights[n+1-i] = weights[i]
    end

    weights ./= sqrt(pi)
    nodes *= sqrt(2)

    nodes = ndims(nodes) == 2 && size(nodes, 2) == 1 ? squeeze(nodes, 2): nodes

    return nodes, weights
end

"""
Computes multivariate Simpson quadrature nodes and weights.

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwsimp

function qnwsimp(n::Int, a::Real, b::Real)
    if n<=1
        error("In qnwsimp: n must be integer greater than one.")
    end

    if n % 2 ==0
        warn("In qnwsimp: n must be odd integer - increasing by 1.")
        n += 1
    end

    dx = (b - a) / (n - 1)
    nodes = collect(a:dx:b)
    weights = repeat([2.0, 4.0], outer=Int[(n + 1.0) / 2.0])
    weights = weights[1:n]
    weights[1] = 1
    weights[end] = 1
    weights = (dx / 3) * weights
    return nodes, weights
end

"""
Computes multivariate trapezoid quadrature nodes and weights.

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwtrap

function qnwtrap(n::Int, a::Real, b::Real)
    if n < 1
        error("n must be at least 1")
    end

    dx = (b - a) / (n - 1)
    nodes = collect(a:dx:b)
    weights = fill(dx, n)
    weights[[1, n]] .*= 0.5
    return nodes, weights
end

"""
Computes nodes and weights for beta distribution

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : First parameter of the beta distribution,
along each dimension
- `b::Union{Real, Vector{Real}}` : Second parameter of the beta distribution,
along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwbeta

function qnwbeta(n::Int, a::Real, b::Real)
    a -= 1
    b -= 1
    ab = a + b

    maxit = 25

    x = zeros(n)
    w = zeros(n)

    z::Float64 = 0.0

    for i=1:n
        if i == 1
            an = a / n
            bn = b / n
            r1 = (1 + a) * (2.78 / (4 + n * n) + 0.768an / n)
            r2 = 1 + 1.48 * an + 0.96bn + 0.452an*an + 0.83an*bn
            z = 1 - r1 / r2

        elseif i == 2
            r1 = (4.1 + a) / ((1 + a) * (1 + 0.156a))
            r2 = 1 + 0.06 * (n - 8) * (1 + 0.12a) / n
            r3 = 1 + 0.012b * (1 + 0.25 * abs(a)) / n
            z = z - (1 - z) * r1 * r2 * r3

        elseif i == 3
            r1 = (1.67 + 0.28a) / (1 + 0.37a)
            r2 = 1 + 0.22 * (n - 8) / n
            r3 = 1 + 8 * b / ((6.28 + b) * n * n)
            z = z - (x[1] - z) * r1 * r2 * r3

        elseif i == n - 1
            r1 = (1 + 0.235b) / (0.766 + 0.119b)
            r2 = 1 / (1 + 0.639 * (n - 4) / (1 + 0.71 * (n - 4)))
            r3 = 1 / (1 + 20a / ((7.5+ a ) * n * n))
            z = z + (z - x[n-3]) * r1 * r2 * r3

        elseif i == n
            r1 = (1 + 0.37b) / (1.67 + 0.28b)
            r2 = 1 / (1 + 0.22 * (n - 8) / n)
            r3 = 1 / (1 + 8 * a / ((6.28+ a ) * n * n))
            z = z + (z - x[n-2]) * r1 * r2 * r3

        else
            z = 3 * x[i-1] - 3 * x[i-2] + x[i-3]
        end

        its = 1
        temp = 0.0
        pp, p2 = 0.0, 0.0
        for its = 1:maxit
            temp = 2 + ab
            p1 = (a - b + temp * z) / 2
            p2 = 1
            for j=2:n
              p3 = p2
              p2 = p1
              temp = 2 * j + ab
              aa = 2 * j * (j + ab) * (temp - 2)
              bb = (temp - 1) * (a * a - b * b + temp * (temp - 2) * z)
              c = 2 * (j - 1 + a) * (j - 1 + b) * temp
              p1 = (bb * p2 - c * p3) / aa
            end
            pp = (n * (a - b - temp * z) * p1 +
                  2 * (n + a) * (n + b) * p2) / (temp * (1 - z * z))
            z1 = z
            z = z1 - p1 ./ pp
            if abs(z - z1) < 3e-14 break end
        end

        if its >= maxit
            error("Failure to converge in qnwbeta")
        end

        x[i] = z
        w[i] = temp / (pp * p2)
    end

    x = (1 - x) ./ 2
    w = w * exp(lgamma(a + n) +
                lgamma(b + n) -
                lgamma(n + 1) -
                lgamma(n + ab + 1))
    w = w / (2 * exp(lgamma(a + 1) +
                     lgamma(b + 1) -
                     lgamma(ab + 2)))

    return x, w
end

"""
Computes nodes and weights for beta distribution

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : First parameter of the gamma distribution,
along each dimension
- `b::Union{Real, Vector{Real}}` : Second parameter of the gamma distribution,
along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
:qnwgamma

function qnwgamma(n::Int, a::Real=1.0, b::Real=1.0)
    a < 0 && error("shape parameter must be positive")
    b < 0 && error("scale parameter must be positive")

    a -= 1
    maxit = 10
    fact = -exp(lgamma(a+n)-lgamma(n)-lgamma(a+1))
    nodes = zeros(n)
    weights = zeros(n)

    z = (1+a)*(3 + 0.92a)/(1 + 2.4n + 1.8a)

    for i=1:n
        # get starting values
        if i==1
            z = (1+a)*(3 + 0.92a)/(1 + 2.4n + 1.8a)
        elseif i==2
            z += (15 + 6.25a)./(1 + 0.9a + 2.5n)
        else
            j = i-2
            z += ((1+2.55j)./(1.9j) + 1.26j*a./(1+3.5j)) * (z-nodes[j])./(1+0.3a)
        end

        # rootfinding iterations
        it = 0
        pp = 0.0
        p2 = 0.0
        for it = 1:maxit
            p1 = 1.0
            p2 = 0.0
            for j=1:n
                p3 = p2
                p2 = p1
                p1 = ((2j - 1 + a - z) * p2 - (j - 1 + a) * p3) ./ j
            end
            pp = (n*p1-(n+a)*p2)./z
            z1 = z
            z = z1-p1 ./ pp
            if abs(z - z1) < 3e-14
                break
            end
        end

        if it >= maxit
            error("failure to converge.")
        end

        nodes[i] = z
        weights[i] = fact / (pp * n * p2)
    end

    return nodes .* b, weights
end

## Multidim versions
for f in [:qnwlege, :qnwcheb, :qnwsimp, :qnwtrap, :qnwbeta, :qnwgamma]
    @eval begin
        function ($f)(n::Vector{Int}, a::Real, b::Real)
            n_n = length(n)
            ($f)(n, fill(a, n_n), fill(b, n_n))
        end

        function ($f)(n::Int, a::Vector, b::Real)
            n_a = length(a)
            ($f)(fill(n, n_a), a, fill(b, n_a))
        end

        function ($f)(n::Int, a::Real, b::Vector)
            n_b = length(b)
            ($f)(fill(n, n_b), fill(a, n_b), b)
        end

        function ($f)(n::Vector{Int}, a::Vector, b::Real)
            n_n, n_a = length(n), length(a)
            if n_n != n_a
                msg = "Cannot construct nodes/weights. n and a have different"
                msg *= " lengths"
                error(msg)
            end
            ($f)(n, a, fill(b, n_a))
        end

        function ($f)(n::Vector{Int}, a::Real, b::Vector)
            n_n, n_b = length(n), length(b)
            if n_n != n_b
                msg = "Cannot construct nodes/weights. n and b have different"
                msg *= " lengths"
                error(msg)
            end
            ($f)(n, fill(a, n_b), b)
        end

        function ($f)(n::Real, a::Vector, b::Vector)
            n_a, n_b = length(a), length(b)
            if n_a != n_b
                msg = "Cannot construct nodes/weights. a and b have different"
                msg *= " lengths"
                error(msg)
            end
            ($f)(fill(n, n_a), a, b)
        end

        function ($f)(n::Vector{Int}, a::Vector, b::Vector)
            n_n, n_a, n_b = length(n), length(a), length(b)

            if !(n_n == n_a == n_b)
                error("n, a, and b must have same number of elements")
            end

            nodes = Vector{Float64}[]
            weights = Vector{Float64}[]

            for i=1:n_n
                _1d = $f(n[i], a[i], b[i])
                push!(nodes, _1d[1])
                push!(weights, _1d[2])
            end
            weights = ckron(weights[end:-1:1]...)
            nodes = gridmake(nodes...)
            return nodes, weights
        end
    end  # @eval
end

## Multidim version for qnworm
function qnwnorm(n::Vector{Int}, mu::Vector, sig2::Matrix=eye(length(n)))
    n_n, n_mu = length(n), length(mu)

    if !(n_n == n_mu)
        error("n and mu must have same number of elements")
    end

    nodes = Vector{Float64}[]
    weights = Vector{Float64}[]

    for i=1:n_n
        _1d = qnwnorm(n[i])
        push!(nodes, _1d[1])
        push!(weights, _1d[2])
    end

    weights = ckron(weights[end:-1:1]...)
    nodes = gridmake(nodes...)

    new_sig2 = chol(sig2)

    nodes = nodes * new_sig2 .+ mu'

    nodes = size(nodes, 2) == 1 ? squeeze(nodes, 2) : nodes

    return nodes, weights
end

# other types of args
qnwnorm(n::Vector{Int}, mu::Vector, sig2::Real) =
    qnwnorm(n, mu, diagm(fill(convert(Float64, sig2), length(n))))

qnwnorm(n::Vector{Int}, mu::Real, sig2::Matrix=eye(length(n))) =
    qnwnorm(n, fill(mu, length(n)), sig2)

qnwnorm(n::Vector{Int}, mu::Real, sig2::Real) =
    qnwnorm(n, fill(mu, length(n)), diagm(fill(convert(Float64, sig2), length(n))))

qnwnorm(n::Int, mu::Vector, sig2::Matrix=eye(length(mu))) =
    qnwnorm(fill(n, length(mu)), mu, sig2)

qnwnorm(n::Int, mu::Vector, sig2::Real) =
    qnwnorm(fill(n, length(mu)), mu, diagm(fill(convert(Float64, sig2), length(mu))))

qnwnorm(n::Int, mu::Real, sig2::Matrix=eye(length(mu))) =
    qnwnorm(fill(n, size(sig2, 1)), fill(mu, size(sig2, 1)), sig2)

qnwnorm(n::Int, mu::Real, sig2::Real) =
    qnwnorm([n], [mu], fill(convert(Float64, sig2), 1, 1))

qnwnorm(n::Vector{Int}, mu::Vector, sig2::Vector) =
    qnwnorm(n, mu, diagm(convert(Array{Float64}, sig2)))

qnwnorm(n::Vector{Int}, mu::Real, sig2::Vector) =
    qnwnorm(n, fill(mu, length(n)), diagm(convert(Array{Float64}, sig2)))

qnwnorm(n::Int, mu::Vector, sig2::Vector) =
    qnwnorm(fill(n, length(mu)), mu, diagm(convert(Array{Float64}, sig2)))

qnwnorm(n::Int, mu::Real, sig2::Vector) =
    qnwnorm(fill(n, length(sig2)), fill(mu, length(sig2)), diagm(convert(Array{Float64}, sig2)))

"""
Computes quadrature nodes and weights for multivariate uniform distribution

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)
"""
function qnwunif(n, a, b)
    nodes, weights = qnwlege(n, a, b)
    weights ./= prod(b - a)
    return nodes, weights
end

"""
Computes quadrature nodes and weights for multivariate uniform distribution

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `mu::Union{Real, Vector{Real}}` : Mean along each dimension
- `sig2::Union{Real, Vector{Real}, Matrix{Real}}(eye(length(n)))` : Covariance
structure

$(qnw_returns)

##### Notes

See also the documentation for `qnwnorm`

$(qnw_refs)
"""
function qnwlogn(n, mu, sig2)
    nodes, weights = qnwnorm(n, mu, sig2)
    return exp(nodes), weights
end

## qnwequi
const equidist_pp = sqrt(primes(7920))  # good for d <= 1000

"""
Generates equidistributed sequences with property that averages
value of integrable function evaluated over the sequence converges
to the integral as n goes to infinity.

##### Arguments

- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension
- `kind::AbstractString("N")`: One of the following:
    - N - Neiderreiter (default)
    - W - Weyl
    - H - Haber
    - R - pseudo Random

$(qnw_returns)

$(qnw_func_notes)

$(qnw_refs)

"""
:qnwequi

function qnwequi(n::Int, a::Vector, b::Vector, kind::AbstractString="N")
    # error checking
    n_a, n_b = length(a), length(b)
    if !(n_a == n_b)
        error("a and b must have same number of elements")
    end

    d = n_a
    i = reshape(1:n, n, 1)
    if kind == "N"
        j = 2.^((1:d)/(d+1))
        nodes = i*j'
        nodes -= fix(nodes)

    elseif kind == "W"
        j = equidist_pp[1:d]
        nodes = i*j'
        nodes -= fix(nodes)

    elseif kind == "H"
        j = equidist_pp[1:d]
        nodes = (i.*(i+1)./2)*j'
        nodes -= fix(nodes)

    elseif kind == "R"
        nodes = rand(n, d)

    else
        error("Unknown `kind` specified. Valid choices are N, W, H, R")
    end

    r = b - a
    nodes = a' .+ nodes .* r'  # use broadcasting here.
    weights = fill((prod(r) / n), n)

    nodes = size(nodes, 2) == 1 ? squeeze(nodes, 2) : nodes

    return nodes, weights
end

# Other argument types
qnwequi(n::Vector{Int}, a::Vector, b::Vector, kind::AbstractString="N") =
    qnwequi(prod(n), a, b, kind)

qnwequi(n::Vector{Int}, a::Real, b::Vector, kind::AbstractString="N") =
    qnwequi(prod(n), fill(a, length(b)), b, kind)

qnwequi(n::Vector{Int}, a::Vector, b::Real, kind::AbstractString="N") =
    qnwequi(prod(n), a, fill(b, length(a)), kind)

qnwequi(n::Vector{Int}, a::Real, b::Real, kind::AbstractString="N") =
    qnwequi(prod(n), fill(a, length(n)), fill(b, length(n)), kind)

qnwequi(n::Int, a::Real, b::Vector, kind::AbstractString="N") =
    qnwequi(n, fill(a, length(b)), b, kind)

qnwequi(n::Int, a::Vector, b::Real, kind::AbstractString="N") =
    qnwequi(n, a, fill(b, length(a)), kind)

qnwequi(n::Int, a::Real, b::Real, kind::AbstractString="N") =
    qnwequi(n, [a], [b], kind)

## Doing the quadrature
"""
Approximate the integral of `f`, given quadrature `nodes` and `weights`

##### Arguments

- `f::Function`: A callable function that is to be approximated over the domain
spanned by `nodes`.
- `nodes::Array`: Quadrature nodes
- `weights::Array`: Quadrature nodes
- `args...(Void)`: additional positional arguments to pass to `f`
- `;kwargs...(Void)`: additional keyword arguments to pass to `f`

##### Returns

- `out::Float64` : The scalar that approximates integral of `f` on the hypercube
formed by `[a, b]`

"""
:do_quad

function do_quad(f::Function, nodes::Array, weights::Vector, args...;
                 kwargs...)
    return dot(f(nodes, args...; kwargs...), weights)
end
do_quad(f::Function, nodes::Array, weights::Vector) = dot(f(nodes), weights)

"""
Integrate the d-dimensional function f on a rectangle with lower and upper bound
for dimension i defined by a[i] and b[i], respectively; using n[i] points.

##### Arguments

- `f::Function` The function to integrate over. This should be a function that
accepts as its first argument a matrix representing points along each dimension
(each dimension is a column). Other arguments that need to be passed to the
function are caught by `args...` and `kwargs...``
- `n::Union{Int, Vector{Int}}` : Number of desired nodes along each dimension
- `a::Union{Real, Vector{Real}}` : Lower endpoint along each dimension
- `b::Union{Real, Vector{Real}}` : Upper endpoint along each dimension
- `kind::AbstractString("lege")` Specifies which type of integration to perform. Valid
values are:
    - `"lege"` : Gauss-Legendre
    - `"cheb"` : Gauss-Chebyshev
    - `"trap"` : trapezoid rule
    - `"simp"` : Simpson rule
    - `"N"` : Neiderreiter equidistributed sequence
    - `"W"` : Weyl equidistributed sequence
    - `"H"` : Haber  equidistributed sequence
    - `"R"` : Monte Carlo
    - `args...(Void)`: additional positional arguments to pass to `f`
    - `;kwargs...(Void)`: additional keyword arguments to pass to `f`

##### Returns

- `out::Float64` : The scalar that approximates integral of `f` on the hypercube
formed by `[a, b]`

$(qnw_refs)

"""
:quadrect

function quadrect(f::Function, n, a, b, kind="lege", args...; kwargs...)
    if lowercase(kind)[1] == 'l'
        nodes, weights = qnwlege(n, a, b)
    elseif lowercase(kind)[1] == 'c'
        nodes, weights = qnwcheb(n, a, b)
    elseif lowercase(kind)[1] == 't'
        nodes, weights = qnwtrap(n, a, b)
    elseif lowercase(kind)[1] == 's'
        nodes, weights = qnwsimp(n, a, b)
    else
        nodes, weights = qnwequi(n, a, b, kind)
    end

    return do_quad(f, nodes, weights, args...; kwargs...)
end

function quadrect(f::Function, n, a, b, kind="lege")
    if lowercase(kind)[1] == 'l'
        nodes, weights = qnwlege(n, a, b)
    elseif lowercase(kind)[1] == 'c'
        nodes, weights = qnwcheb(n, a, b)
    elseif lowercase(kind)[1] == 't'
        nodes, weights = qnwtrap(n, a, b)
    elseif lowercase(kind)[1] == 's'
        nodes, weights = qnwsimp(n, a, b)
    else
        nodes, weights = qnwequi(n, a, b, kind)
    end

    return do_quad(f, nodes, weights)
end


text:
export Hilbert, fermion

# 2x2 Annhilation Operator
a = complex([0  1.0; 0  0])
Id = complex(I + zeros(2,2))

"""
A generic Hilbert space N-particle Hamiltonian. User defines
creation/annhilation operator strings applied to .mat
"""
mutable struct Hilbert
    N::Int
    dim::Int
    mat::Array{Complex}

    function Hilbert(N::Int)
        if N >= 20
            @warn "Full Matrix has more than 1 trillion elements."
        end

        # N-site full dimension
        dim = 2^N
        mat = zeros(Complex,2^N,2^N)

        # Construct
        new(N,dim,mat)
    end
end

"""
Implements fermionic matrix element to a Hamiltonian using a naive
Jordan-Wigner string. For example a creation operator on site 3 out of 6 sites
would produce the following matrix.

a'_3 = Z ⊗ Z ⊗ a' ⊗ 1 ⊗ 1 ⊗ 1

FIXME: very inefficient
"""
function fermion(Ham::Hilbert, s::String, value::AbstractFloat)

    # Get full dimension
    N = Ham.N

    # Temporary matrix storage
    working_mat = zeros(Complex,2^N,2^N)

    # Get sites for the term
    operators = parse_operator_string(s)

    isFirst = true

    # Loop over single creation/annhilation operators
    # in the fermion operator.
    for term in operators

        # Get site and operator type
        idx = term[1]
        create = term[2]

        # Input checks
        @assert idx >= 1 && idx <= N
        @assert create == 0 || create == 1

        # Populate with either creation or annihilation
        if create == 1
            op = a'
        else
            op = a
        end

        # Build Matrix for each operator, then multiply
        if idx == 1
            mat = op
            for i = idx+1:N
                mat = mat ⊗ Id
            end
        else
            mat = Z
            for i = 2:idx-1
                mat = mat ⊗ Z
            end

            mat = mat ⊗ op

            for i = idx+1:N
                mat = mat ⊗ Id
            end
        end

        # Check if first term in list
        if isFirst
            working_mat = mat
            isFirst = false
        else
            working_mat = working_mat * mat
        end

    end

    # Add Matrix element to Hamiltonian
    Ham.mat += value * working_mat
end


text:
#
# This script uses Clang.jl to generate a C binding from the HSA header files
#
using Clang.cindex
using Clang.wrap_c

script_path = dirname(Base.source_path())

include("$script_path/../src/discovery.jl") # pull in include path

hsa_hdrs = map(x->joinpath(hsa_include_path, x), [
    "hsa.h",
    "hsa_ext_finalize.h",
    "hsa_ext_image.h"
])

clang_includes = [
    hsa_include_path
]

if haskey(ENV, "INCLUDE_PATH")
    unshift!(clang_includes, ENV["INCLUDE_PATH"])
else
    warn("environment variable INCLUDE_PATH is not set, generator may not be able to find all headers included by hsa.h")
end

excluded_symbols = Set([
   # ignore noise
   "HSA_API",
   "HSA_IMPORT",
   # Custom wrappers for
   "hsa_system_get_info",
   "hsa_iterate_agents",
   "hsa_executable_iterate_symbols"
])

const wc = wrap_c.init(
    common_file = "hsa_datatypes.jl",
    output_file = "hsa_wrappers.jl",
    output_dir = "$script_path/../src/binding/generated"
)

wc.clang_includes = clang_includes

wc.header_library = x-> "libhsa"

wc.header_wrapped = (x,y)-> begin
    contains(y,"hsa")
end

# This function determines for each clang cursor encountered
# whether it should be included in the generated julia files
wc.cursor_wrapped = (name, cursor)-> begin
    if in(name, excluded_symbols)
        return false
    elseif endswith(name, "_s") || endswith(name, "_u")
        # Names ending in _s are the underlying struct types
        # those ending in _u underlying union types
        # typedefd to something ending in _t
        return false
    elseif cindex.name(cursor) == ""
        return false
    else
        return true
    end
end

wc.headers = hsa_hdrs

include("$script_path/../src/binding/custom/def.jl") # pull in type mappings

include("gen_argtypes.jl")
include("gen_atomics.jl")
include("gen_constants.jl")
include("gen_getters.jl")
include("gen_fieldgetters.jl")

wc.rewriter = function(obuf)
    map_argtypes(obuf)
    map_memspec_funcs(obuf)
    map_constants(obuf)
    gen_fieldgetters(obuf)
    gen_getters(obuf)

    return obuf
end

if !isinteractive()
    run(wc)
end


text:
module sparsePolyChaos

using PolyChaos
using LinearAlgebra: diag, cond

include("error_estimation.jl")
include("regression.jl")
include("sparse_pce.jl")
include("sparse_pce_mult.jl")

end # module

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
function grabscrab(anagram::String, dictionary::Array{String})::Array{String}
    result = String[]
    sorted_anagram = sort!([c for c in anagram])
    for word in dictionary
        if sorted_anagram == sort!([c for c in word])
            push!(result,word)
        end
    end
    result
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas