text:
using Halite3
using Logging

# Initialize the game.
log, constants, game = init_game()
# log is a function which takes a string and allows you to save messages for yourself.
# log is not available to use inside /src/*jl instead use @info to log from inside Halite3 module.
# This is required because the regular STDOUT (print statements) are reserved for the engine-bot communication.

# At this point "game" variable is populated with initial map data.
# This is a good place to do computationally expensive start-up pre-processing.
# As soon as you call "start_game" function below, the 2 second per turn timer will start.

start_game("MyJuliaBot")
# Now that your bot is initialized, save a message to yourself in the log file with some important information.
#   Here, you log here your id, which you can always fetch from the game object by using game.id.
log("Successfully created bot! My Player ID is $(game.id)")

while true
    # This loop handles each turn of the game. The game object changes every turn, and you refresh that state by
    #   running update_frame!(game).
    update_frame!(game)
    log("=============== TURN $(game.turn_number) ================")
    # You extract player metadata and the updated map metadata here for convenience.
    me = game.me
    game_map = game.game_map

    # A command queue holds all the commands you will run this turn. You build this list up and submit it at the
    #   end of the turn.
    command_queue = Vector()

    for ship in get_ships(me)
        # For each of your ships, move randomly if the ship is on a low halite location or the ship is full.
        #   Else, collect halite.
        if game_map[ship.position].halite_amount < constants.MAX_HALITE / 10 || is_full(ship, constants.MAX_HALITE)
            push!(command_queue, move(ship, random_direction()))
        else
            push!(command_queue, stay_still(ship))
        end
    end

    # If the game is in the first 200 turns and you have enough halite, spawn a ship.
    # Don't spawn a ship if you currently have a ship at port, though - the ships will collide.
    if game.turn_number <= 200 && me.halite_amount >= constants.SHIP_COST && !is_occupied(game_map[me.shipyard.position])
        push!(command_queue, spawn(me.shipyard))
    end

    # Send your moves back to the game environment, ending this turn.
    end_turn(command_queue)
end

max_stars_repo_name: johnnykwwang/Halite-III

max_issues_repo_name: johnnykwwang/Halite-III

max_forks_repo_name: johnnykwwang/Halite-III


text:
using Distributions
using JSON

p = 0.6
d = Geometric( p )

x = [ -0.5, -0.25, 0, 0.25, 0.5 ]

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("p", p),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/number.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

ext: jl

max_stars_repo_name: distributions-io/geometric-mgf

max_issues_repo_name: distributions-io/geometric-mgf

max_forks_repo_name: distributions-io/geometric-mgf


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

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

export BrickletLaserRangeFinderV2Configuration
struct BrickletLaserRangeFinderV2Configuration
    acquisition_count::Integer
    enable_quick_termination::Bool
    threshold_value::Integer
    measurement_frequency::Integer
end

export BrickletLaserRangeFinderV2MovingAverage
struct BrickletLaserRangeFinderV2MovingAverage
    distance_average_length::Integer
    velocity_average_length::Integer
end

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

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

export BrickletLaserRangeFinderV2
"""
Measures distance up to 40m with laser light
"""
mutable struct BrickletLaserRangeFinderV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_distance
"""
    $(SIGNATURES)

Returns the measured distance.

The laser has to be enabled, see :func:`Set Enable`.

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

export set_distance_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Distance` 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:`Distance` 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_distance_callback_configuration(device::BrickletLaserRangeFinderV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_distance_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_distance_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Distance Callback Configuration`.
"""
function get_distance_callback_configuration(device::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_distance_callback_configuration()
end

export get_velocity
"""
    $(SIGNATURES)

Returns the measured velocity. The value has a range of -12800 to 12700
and is given in 1/100 m/s.

The velocity measurement only produces stables results if a fixed
measurement rate (see :func:`Set Configuration`) is configured. Also the laser
has to be enabled, see :func:`Set Enable`.

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

export set_velocity_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Velocity` 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:`Velocity` 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_velocity_callback_configuration(device::BrickletLaserRangeFinderV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_velocity_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_velocity_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Velocity Callback Configuration`.
"""
function get_velocity_callback_configuration(device::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_velocity_callback_configuration()
end

export set_enable
"""
    $(SIGNATURES)

Enables the laser of the LIDAR if set to *true*.

We recommend that you wait 250ms after enabling the laser before
the first call of :func:`Get Distance` to ensure stable measurements.
"""
function set_enable(device::BrickletLaserRangeFinderV2, enable)
    device.deviceInternal.set_enable(enable)
end

export get_enable
"""
    $(SIGNATURES)

Returns the value as set by :func:`Set Enable`.
"""
function get_enable(device::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_enable()
end

export set_configuration
"""
    $(SIGNATURES)

The **Acquisition Count** defines the number of times the Laser Range Finder Bricklet
will integrate acquisitions to find a correlation record peak. With a higher count,
the Bricklet can measure longer distances. With a lower count, the rate increases. The
allowed values are 1-255.

If you set **Enable Quick Termination** to true, the distance measurement will be terminated
early if a high peak was already detected. This means that a higher measurement rate can be achieved
and long distances can be measured at the same time. However, the chance of false-positive
distance measurements increases.

Normally the distance is calculated with a detection algorithm that uses peak value,
signal strength and noise. You can however also define a fixed **Threshold Value**.
Set this to a low value if you want to measure the distance to something that has
very little reflection (e.g. glass) and set it to a high value if you want to measure
the distance to something with a very high reflection (e.g. mirror). Set this to 0 to
use the default algorithm. The other allowed values are 1-255.

Set the **Measurement Frequency** to force a fixed measurement rate. If set to 0,
the Laser Range Finder Bricklet will use the optimal frequency according to the other
configurations and the actual measured distance. Since the rate is not fixed in this case,
the velocity measurement is not stable. For a stable velocity measurement you should
set a fixed measurement frequency. The lower the frequency, the higher is the resolution
of the calculated velocity. The allowed values are 10Hz-500Hz (and 0 to turn the fixed
frequency off).

The default values for Acquisition Count, Enable Quick Termination, Threshold Value and
Measurement Frequency are 128, false, 0 and 0.
"""
function set_configuration(device::BrickletLaserRangeFinderV2, acquisition_count, enable_quick_termination, threshold_value, measurement_frequency)
    device.deviceInternal.set_configuration(acquisition_count, enable_quick_termination, threshold_value, measurement_frequency)
end

export get_configuration
"""
    $(SIGNATURES)

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

export set_moving_average
"""
    $(SIGNATURES)

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

Setting the length to 0 will turn the averaging completely off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletLaserRangeFinderV2, distance_average_length, velocity_average_length)
    device.deviceInternal.set_moving_average(distance_average_length, velocity_average_length)
end

export get_moving_average
"""
    $(SIGNATURES)

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

export set_offset_calibration
"""
    $(SIGNATURES)

The offset is added to the measured distance.
It is saved in non-volatile memory, you only have to set it once.

The Bricklet comes with a per-sensor factory-calibrated offset value,
you should not have to call this function.

If you want to re-calibrate the offset you first have to set it to 0.
Calculate the offset by measuring the distance to a known distance
and set it again.
"""
function set_offset_calibration(device::BrickletLaserRangeFinderV2, offset)
    device.deviceInternal.set_offset_calibration(offset)
end

export get_offset_calibration
"""
    $(SIGNATURES)

Returns the offset value as set by :func:`Set Offset Calibration`.
"""
function get_offset_calibration(device::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_offset_calibration()
end

export set_distance_led_config
"""
    $(SIGNATURES)

Configures the distance LED to be either turned off, turned on, blink in
heartbeat mode or show the distance (brighter = object is nearer).
"""
function set_distance_led_config(device::BrickletLaserRangeFinderV2, config)
    device.deviceInternal.set_distance_led_config(config)
end

export get_distance_led_config
"""
    $(SIGNATURES)

Returns the LED configuration as set by :func:`Set Distance LED Config`
"""
function get_distance_led_config(device::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_distance_led_config()
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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2, 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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2, 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::BrickletLaserRangeFinderV2, 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::BrickletLaserRangeFinderV2, 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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2, 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::BrickletLaserRangeFinderV2)
    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::BrickletLaserRangeFinderV2)
    return device.deviceInternal.get_identity()
end

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


text:
using PyCall

function __init__()
    py"""
    import os
    import shutil
    import re

    def main(model_path):
        ''' Convert fitparam/n/*.dat -> out/n/*.npy

        Parameters
        ----------
        model_path : str
            Path to your model written in Julia.

        Usage
        -----
        $ python dat2npy.py
        $ mv dat2npy/out/ path_to_biomass/

        '''
        try:
            import numpy as np
        except ImportError:
            print("numpy: Not installed")
        n_file = []
        fitparam_files = os.listdir(
            os.path.join(
                model_path,
                "fitparam",
            )
        )
        for file in fitparam_files:
            if re.match(r"\d", file):
                n_file.append(int(file))
        for nth_paramset in n_file:
            os.makedirs(
                os.path.join(
                    model_path,
                    "dat2npy",
                    "out",
                    f"{nth_paramset:d}",
                ),
                exist_ok=True,
            )
            nth_fitparam_files = os.listdir(
                os.path.join(
                    model_path,
                    "fitparam",
                    f"{nth_paramset:d}",
                )
            )
            for dat_file in nth_fitparam_files:
                if os.path.splitext(dat_file)[-1] == ".dat":
                    if "fit" in dat_file:
                        '''
                        - fit_param%d.dat -> fit_param%d.npy
                        - best_fitness.dat -> best_fitness.npy
                        '''
                        try:
                            data = np.loadtxt(
                                os.path.join(
                                    model_path,
                                    "fitparam",
                                    f"{nth_paramset:d}",
                                    f"{dat_file}",
                                ),
                                dtype="float",
                            )
                        except ValueError:
                            pass
                    else:
                        '''
                        - count_num.dat -> count_num.npy
                        - generation.dat -> generation.npy
                        '''
                        data = np.loadtxt(
                            os.path.join(
                                model_path,
                                "fitparam",
                                f"{nth_paramset:d}",
                                f"{dat_file}",
                            ),
                            dtype="int",
                        )
                np.save(
                    os.path.join(
                        model_path,
                        "dat2npy",
                        "out",
                        f"{nth_paramset:d}",
                        dat_file.replace(".dat", ".npy"),
                    ),
                    data,
                )
            if os.path.isfile(
                os.path.join(
                    model_path,
                    "logs",
                    f"{nth_paramset:d}.log",
                )
            ):
                shutil.copyfile(
                    os.path.join(
                        model_path,
                        "logs",
                        f"{nth_paramset:d}.log",
                    ),
                    os.path.join(
                        model_path,
                        "dat2npy",
                        "out",
                        f"{nth_paramset:d}",
                        "optimization.log",
                    ),
                )
    """
end

param2biomass(model_path::String) = py"main"(model_path)


text:
## Testing Evaluate Monomial
@testset "Test evaluate monomials" begin
    # Test single point
    x = collect(1.0:3.0)
    m = NarmaxLasso.Monomial([2, 2, 2])
    @test NarmaxLasso.evaluate(x, m) ≈ 36.0
    # Test vector of points
    x = [1.0 2.0 3.0; 1.0 2.0 3.0; 1.0 2.0 3.0]
    m = NarmaxLasso.Monomial([2, 2, 2])
    @test NarmaxLasso.evaluate(x, m) ≈ [36.0, 36.0, 36.0]
    # Test vector of monomials
    m = [NarmaxLasso.Monomial([2, 2, 2]), NarmaxLasso.Monomial([2, 2, 2])]
    @test NarmaxLasso.evaluate(x, m) ≈ 36*ones(3, 2)
end

## Test print monomials
@testset "Test print monomials list" begin
    @test string(NarmaxLasso.Monomial([1, 0])) == "x_1"
    @test string(NarmaxLasso.Monomial([1, 2])) == "x_1 * x_2^2"
    @test string(NarmaxLasso.Monomial([1, 0, 3, 0])) == "x_1 * x_3^3"
    @test string(NarmaxLasso.Monomial([1, 0, 3, 1])) == "x_1 * x_3^3 * x_4"
end

## Test generate monomials
@testset "Test generate monomials" begin
    # Test first order monomials
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 3, 1)
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 0, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 1, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 0, 1]))
    @test ~contains(isequal, l, NarmaxLasso.Monomial([0, 1, 1]))
    # Test second order monomials
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 3, 2)
    @test contains(isequal, l, NarmaxLasso.Monomial([2, 0, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 1, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 2, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 1, 1]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 0, 2]))
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 0, 1]))
    @test ~contains(isequal, l, NarmaxLasso.Monomial([0, 0, 1]))
    # Test maximum order 2
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 3, 1:2)
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 0, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 1, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 0, 1]))
    @test contains(isequal, l, NarmaxLasso.Monomial([2, 0, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 1, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 2, 0]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 1, 1]))
    @test contains(isequal, l, NarmaxLasso.Monomial([0, 0, 2]))
    @test contains(isequal, l, NarmaxLasso.Monomial([1, 0, 1]))
end

## Test split monomial list
@testset "Test detach terms" begin
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 3, 1:2)
    l1, l2 = NarmaxLasso.detach_terms(l, BitVector([false, true, true]))
    @test length(l2) == 2
    @test length(l1) == 7
    @test contains(isequal, l2, NarmaxLasso.Monomial([1, 0, 0]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([0, 1, 0]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([0, 0, 1]))
    @test contains(isequal, l2, NarmaxLasso.Monomial([2, 0, 0]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([1, 1, 0]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([0, 2, 0]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([0, 1, 1]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([0, 0, 2]))
    @test contains(isequal, l1, NarmaxLasso.Monomial([1, 0, 1]))
end


text:
"""
    CovNWFn(g0,m=0)

Calculates covariance matrix of sample average.

# Input
- `g0::Matrix`: Txq Matrix of q moment conditions
- `m:int`:     scalar, number of lags to use

# Output
- `S::Matrix`: qxq covariance matrix(average g0)

"""
function CovNWFn(g0,m=0)

    T = size(g0,1)                    #g0 is Txq
    m = min(m,T-1)                    #number of lags

    g = g0 .- mean(g0,dims=1)         #normalizing to zero means

    S = g'g                           #(qxT)*(Txq)
    for s = 1:m
        Λ_s = g[s+1:T,:]'g[1:T-s,:]   #same as Sum[g_t*g_{t-s}',t=s+1,T]
        S   = S  +  (1 - s/(m+1))*(Λ_s + Λ_s')
    end

    return S

end

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
module TestError

using Test
using Configurations

Base.@kwdef struct OptionA
    x::Int = 1
end

Configurations.is_option(::Type{OptionA}) = true
Configurations.is_option(::OptionA) = true

@testset "error" begin
    @test_throws ErrorException from_dict(OptionA, Dict{String, Any}("x" => 1))
    @test_throws ErrorException Configurations.type_alias(Float32)
    @test_throws ErrorException Configurations.type_alias(OptionA)
    @test_throws ErrorException Configurations.get_type_alias_map(Float32)
end

end


text: mkpath("/tmp/unusefuldir/helloworld.d/test123")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Cliffords

@test Clifford(1/sqrt(2) * [1 1; 1 -1]) == H
@test Clifford(eye(2)) == RI
@test Clifford(eye(4)) == kron(RI,RI)
@test Clifford([1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]) == CNOT
@test Clifford(diagm(0 => [1,1,1,-1])) == CZ
@test Clifford(exp(-im*pi/4*[1 0; 0 -1])) == S

function eq_upto_phase(A, B)
    idx = findfirst(x -> x != 0, A)
    rel_phase = (B[idx] == 0) ? 1.0 : (A[idx] / B[idx])
    return rel_phase * A ≈ B
end
@test eq_upto_phase(complex(RI), eye(2))
@test eq_upto_phase(complex(H), 1/sqrt(2) * [1 1; 1 -1])
@test eq_upto_phase(complex(CNOT), [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0])

# Cliffords * Paulis
@test H * Id == Id
@test H * X == Z
@test H * Z == X
@test H * Y == -Y

@test kron(RI, H) * II == II
@test kron(RI, H) * IX == IZ
@test kron(RI, H) * IZ == IX
@test kron(RI, H) * IY == -IY

@test CNOT * II == II
@test CNOT * IX == IX
@test CNOT * XI == XX
@test CNOT * YY == -XZ

# Cliffords \ Paulis
@test H \ Id == Id
@test H \ X == Z
@test H \ Z == X
@test H \ Y == -Y

@test kron(RI, H) \ II == II
@test kron(RI, H) \ IX == IZ
@test kron(RI, H) \ IZ == IX
@test kron(RI, H) \ IY == -IY

@test CNOT \ II == II
@test CNOT \ IX == IX
@test CNOT \ XI == XX
@test CNOT \ YY == -XZ

# Cliffords * Cliffords

@test H * S * H == Clifford(Dict(+X=>+X,+Z=>-Y),Dict(+X=>+X,+Z=>+Y))
@test H * S * H == Clifford(exp(-im*pi/4*[0 1; 1 0]))
@test RY * RX * RY == SelfInverseClifford(Dict(+X=>+X,+Z=>-Z))
@test kron(RI, H) * CNOT * kron(RI, H) == SelfInverseClifford(Dict(+IZ=>+IZ,+ZI=>+ZI,+XI=>+XZ,+IX=>+ZX))
@test kron(RI, H) * CNOT * kron(RI, H) == CZ
@test S * S == RZ

# 3 CNOTs is a SWAP
CNOT21 = Cliffords.expand(CNOT, [2,1], 2)
@test CNOT * CNOT21 * CNOT == SWAP

# iSWAP
iSWAP = CZ*kron(S,S)*SWAP
@test iSWAP == Clifford([1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1])

# Cliffords \ Cliffords
@test H \ H == H * H
@test CNOT \ CNOT == CNOT * CNOT
@test RY * RX * RY == RY \ (RX * RY)
@test S \ S == RI
@test S \ RZ == S

# Local Cliffords
@test localclifford(1) == RI

# Cliffords -> Pauli
@test Clifford(Id) == RI
@test Clifford(X)  == RX
@test Clifford(Y)  == RY
@test Clifford(Z)  == RZ
@test Clifford(XX) == kron(RX,RX)
@test Clifford(XI) == kron(RX,RI)

# Clifford and Matrix interaction
@test typeof(promote(H, eye(2))) == Tuple{Array{Float64,2}, Array{Float64,2}}
@test typeof(promote(H, eye(ComplexF64, 2))) == Tuple{Array{ComplexF64,2}, Array{ComplexF64,2}}


text:
#####
##### AlphaZero.jl
##### Jonathan Laurent, Carnegie Mellon University (2019-2020)
#####

module AlphaZero

# Submodules
export MCTS, MinMax, GameInterface, GI, Report, Benchmark
export Network, KnetLib, FluxLib, NetLib
# AlphaZero parameters
export Params, SelfPlayParams, LearningParams, ArenaParams
export MctsParams, MemAnalysisParams
export SamplesWeighingPolicy, CONSTANT_WEIGHT, LOG_WEIGHT, LINEAR_WEIGHT
export AbstractSchedule, ConstSchedule, PLSchedule, StepSchedule
# Players and games
export AbstractGame, AbstractPlayer, TwoPlayers, Trace
export think, select_move, reset_player!, player_temperature, apply_temperature
export play_game, interactive!, total_reward
export MctsPlayer, RandomPlayer, EpsilonGreedyPlayer, NetworkPlayer, Human
export ColorPolicy, ALTERNATE_COLORS, BASELINE_WHITE, CONTENDER_WHITE
# Networks
export AbstractNetwork, OptimiserSpec, Nesterov, CyclicNesterov, Adam
# Training environments
export Env, train!, get_experience
# User interface
export UserInterface
export Session, resume!, save, play_interactive_game
export Explorer, start_explorer

include("util.jl")
import .Util
using .Util: Option, apply_temperature

include("game.jl")
using .GameInterface
const GI = GameInterface

include("mcts.jl")
import .MCTS

include("networks/network.jl")
using .Network

include("batchifier.jl")
import .Batchifier

using Formatting
using Base: @kwdef
using DataStructures: CircularBuffer
using Distributions: Categorical, Dirichlet
using Statistics: mean

import Distributed

include("schedule.jl")
include("params.jl")
include("report.jl")
include("trace.jl")
include("memory.jl")
include("learning.jl")
include("play.jl")
include("training.jl")

include("minmax.jl")
import .MinMax

include("benchmark.jl")
using .Benchmark

include("networks/knet.jl")
# import .KnetLib

include("networks/flux.jl")
# import .FluxLib

# The default user interface is included here for convenience but it could be
# replaced or separated from the main AlphaZero.jl package (which would also
# enable dropping some dependencies such as Crayons or JSON3).

include("ui/ui.jl")
using .UserInterface

# Choose the default DL framework based on an environment variable
function __init__()
  @eval begin
    const DEFAULT_DL_FRAMEWORK = get(ENV, "ALPHAZERO_DEFAULT_DL_FRAMEWORK", "FLUX")
    const NetLib =
      if DEFAULT_DL_FRAMEWORK == "FLUX"
        @info "Using Flux."
        FluxLib
      elseif DEFAULT_DL_FRAMEWORK == "KNET"
        @info "Using Knet."
        KnetLib
      else
        error("Unknown DL framework: $(DEFAULT_DL_FRAMEWORK)")
      end
  end
end

end


text:
function sum_border(A::AbstractArray{Float64,3})
    s = 0.
    ndim, h, w = size(A)
    for i = 1:h, j = 1:w
        if i == 1 || j == 1 || i == h || j == w
            for d = 1:ndim
                s += abs(A[d, i, j])
            end
        end
    end
    s
end

@testset "uniform_field" begin
    @test_throws UndefVarError uniform_field
    @test typeof(Augmentor.uniform_field) <: Function

    A = @inferred Augmentor.uniform_field(5, 6)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test 2.5 > mapreduce(abs, +, A) > 0.

    A = @inferred Augmentor.uniform_field(5, 6, .5, false, true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 2.5

    A = Augmentor.uniform_field(5, 6, scale=.5)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 2.5

    A = @inferred Augmentor.uniform_field(5, 6, .5, true, true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) != 0.
    @test mapreduce(abs, +, A) > 2.5

    A = Augmentor.uniform_field(5, 6, scale=.5, border=true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) != 0.
    @test mapreduce(abs, +, A) > 2.5

    A = @inferred Augmentor.uniform_field(5, 6, .5, false, false)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 5.

    A = Augmentor.uniform_field(5, 6, scale=.5, normalize=false)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 5.

    # TODO: test normalization
end

@testset "gaussian_field" begin
    @test_throws UndefVarError gaussian_field
    @test typeof(Augmentor.gaussian_field) <: Function

    A = @inferred Augmentor.gaussian_field(5, 6)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test 2.2 > mapreduce(abs, +, A) > 0.

    A = @inferred Augmentor.gaussian_field(5, 6, .5, 2, 1, false, true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 2.2

    A = Augmentor.gaussian_field(5, 6, scale=.5)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.
    @test mapreduce(abs, +, A) > 2.2

    A = @inferred Augmentor.gaussian_field(5, 6, .5, 2, 1, true, true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) != 0.
    @test mapreduce(abs, +, A) > 2.2

    A = Augmentor.gaussian_field(5, 6, scale=.5, border=true)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) != 0.
    @test mapreduce(abs, +, A) > 2.2

    A = Augmentor.gaussian_field(5, 6, scale=.5, normalize=false)
    @test typeof(A) <: Array{Float64}
    @test size(A) === (2, 5, 6)
    @test sum_border(A) == 0.

    # TODO: test normalization
    # TODO: test iterations
end

@testset "DistortedView" begin
    @test_throws UndefVarError DistortedView
    @test Augmentor.DistortedView <: AbstractArray

    A1 = [0. 0 0; 0 -.5 0; 0 0 0]
    A2 = [0. 0 0; 0  .5 0; 0 0 0]
    A = zeros(2, 3, 3)
    A[1,:,:] = A1
    A[2,:,:] = A2
    dv = @inferred Augmentor.DistortedView(camera, A)
    @test parent(dv) === camera
    @test size(dv) == size(camera)
    @test eltype(dv) == eltype(camera)

    @test summary(dv) == "512×512 Augmentor.DistortedView(::Array{Gray{N0f8},2}, ::Array{Float64,3} as 3×3 vector field) with eltype Gray{Normed{UInt8,8}}"
    @test_reference "reference/distort_static.txt" dv

    camerao = OffsetArrays.no_offset_view(OffsetArray(camera, (-5,-10)))
    dv2 = @inferred Augmentor.DistortedView(camerao, A)
    @test size(dv2) == size(camera)
    @test eltype(dv2) == eltype(camera)
    @test summary(dv2) == "512×512 Augmentor.DistortedView(::Array{Gray{N0f8},2}, ::Array{Float64,3} as 3×3 vector field) with eltype Gray{Normed{UInt8,8}}"
    @test_reference "reference/distort_static.txt" dv2
    v = view(Augmentor.DistortedView(rand(10,10), A), 2:8, 3:10)
    @test summary(v) == "7×8 view(Augmentor.DistortedView(::Array{Float64,2}, ::Array{Float64,3} as 3×3 vector field), 2:8, 3:10) with eltype Float64"
end


text:
"""
	modelframe(formula::EconFormula, data::AbstractDataFrame;
		contrasts::Dict = Dict())

	This formula is akin to ModelFrame for Formula, but rather than being a
	struct it returns
		data::DataFrame
		varlist::Vector{String}
		response::Vector{T}
		exogenous::Matrix{T}
		endogenous::Matrix{T}
		instruments::Matrix{T}
		absorb::Vector{Vector{Vector{Int64}}}
		clusters::Vector{Vector{Vector{Int64}}}
"""
function modelframe(formula::EconFormula,
	data::AbstractDataFrame;
	contrasts::Dict = Dict(),
	θ::Vector{T} = Vector{Real}()) where T <: Real
	response = formula.exogenous.lhs
	exogenous = StatsModels.Terms(formula.exogenous)
	endogenous = StatsModels.Terms(formula.endogenous)
	instruments = StatsModels.Terms(formula.instruments)
	absorb = StatsModels.Terms(formula.absorb)
	clusters = StatsModels.Terms(formula.clusters)
	vars = Symbol.(reduce(union, getfield.((exogenous, endogenous, instruments, absorb, clusters), :eterms)))
	df = data[vars]
	dropmissing!(df)
	mf = ModelFrame(exogenous, df, contrasts = contrasts)
	y = Vector{Float64}(model_response(mf))
	varlist = coefnames(mf)
	MM = ModelMatrix(mf)
	assign = MM.assign
	X = MM.m
	if length(endogenous.eterms) > 1
		mf = ModelFrame(endogenous, df)
		append!(varlist, coefnames(mf)[2:end])
		mm = ModelMatrix(mf)
		if unique(mm.assign) ≠ mm.assign
			@assert false "Endogenous variables must not be categorical variables with more than two levels."
		end
		append!(assign, mm.assign .+ length(unique(assign)))
		z = mm.m[:, map(elem -> elem > 0, mm.assign)]
		mf = ModelFrame(instruments, df)
		Z = ModelMatrix(mf).m[:,2:end]
	else
		z = zeros(length(y),0)
		Z = zeros(length(y),0)
	end
	D = groups(formula.absorb, data)
	R = groups(formula.clusters, data)
	if isempty(D)
		D = Vector{Vector{Vector{Int64}}}()
	elseif !isempty(θ)
		@assert length(D) == 1 "Partial demeaning is implemented for only one dimension."
		@assert length(θ) == length(D[1]) "θ must be empty or same length as D[1]."
		y = partialwithin(y, D[1], θ)
		X = partialwithin(X, D[1], θ)
		if size(z, 2) > 0
			z = partialwithin(z, D[1], θ)
			Z = partialwithin(Z, D[1], θ)
		end
	else
		(m, singletons) = dropsingletons!(D)
		remapper = makeremapper(m, singletons)
		remapping!(D, remapper)
		validobs = setdiff(1:length(y), singletons)
		y = within(y[validobs], D)
		X = within(X[validobs,:], D)
		if size(z, 2) > 0
			z = within(z[validobs,:], D)
			Z = within(Z[validobs,:], D)
		end
	end
	if isempty(R)
		R = Vector{Vector{Vector{Int64}}}()
	end
	return df, varlist, assign, y, X, z, Z, D, R
end


text:
abstract type ClassificationLikelihood{T<:Real} <: Likelihood{T} end

include("logistic.jl")
include("bayesiansvm.jl")

export AugmentedLogisticLikelihood, LogisticLikelihood

""" Return the labels in a vector of vectors for multiple outputs"""
function treat_labels!(y::AbstractArray{T,N},likelihood::L) where {T,N,L<:ClassificationLikelihood}
    @assert T<:Real "For classification target(s) should be real valued (Bool,Integer or Float)"
    @assert N <= 2 "Target should be a matrix or a vector"
    labels = Int64.(unique(y))
    @assert (length(labels) <= 2) && ((sort(labels) == [0;1]) || (sort(labels) == [-1;1])) "Labels of y should be binary {-1,1} or {0,1}"
    if N == 1
        return [y],1,likelihood
    else
        return [y[:,i] for i in 1:size(y,2)],size(y,2),likelihood
    end
end


text:
struct Propagator{T,D} <: QuMatrix
    data::Matrix{Complex{T}}
    Δt::T
    dims::Dims{D}
end

function (U::Propagator)(ψ::Ket)
    dimsmatch(U,ψ)
    return Ket(data(U)*data(ψ),dims(ψ))
end

function (U::Propagator)(ψ::Ket, n::Integer)
    dimsmatch(U,ψ)
    return Ket(data(U)^n*data(ψ),dims(ψ))
end

function (U::Propagator)(ρ::Operator)
    dimsmatch(U,ρ)
    d = prod(dims(ρ))
    return Operator(reshape(data(U)*vec(data(ρ)),(d,d)),dims(ρ))
end

function (U::Propagator)(ρ::Operator, n::Integer)
    dimsmatch(U,ρ)
    d = prod(dims(ρ))
    return Operator(reshape(data(U)^n*vec(data(ρ)),(d,d)),dims(ρ))
end

# Convert a Propagator to an Operator
Operator(U::Propagator) = Operator(data(U),dims(U))

# Convert a Liouvillian to a Propagator
function SchrodingerProp(L::Liouvillian{T},tspan,alg=Vern8();kwargs...) where T
    d = prod(dims(L))
    prob = ODEProblem(LiouvillianODE(L),Matrix{Complex{T}}(I,d,d),tspan)
    sol  = __solve(prob,alg;save_start=false,saveat=[tspan[end]],abstol=1E-8,reltol=1E-6,kwargs...)
    return Propagator(sol.u[],(tspan[2]-tspan[1]),dims(L))
end

function LindbladProp(L::Liouvillian{T},tspan,alg=Tsit5();kwargs...) where T
    d = prod(dims(L))^2
    prob = ODEProblem(LiouvillianODE(L),Matrix{Complex{T}}(I,d,d),tspan)
    sol  = __solve(prob,alg;save_start=false,saveat=[tspan[end]],abstol=1E-7,reltol=1E-5,kwargs...)
    return Propagator(sol.u[],(tspan[2]-tspan[1]),dims(L))
end


text:
# test stacks and queues

using DataStructures
using Base.Test

# Stack

s = Stack(Int, 5)
n = 100

@test length(s) == 0
@test isempty(s)
@test_throws ErrorException top(s)
@test_throws ErrorException pop!(s)

for i = 1 : n
    push!(s, i)
    @test top(s) == i
    @test !isempty(s)
    @test length(s) == i
end

for i = 1 : n
    x = pop!(s)
    @test x == n - i + 1
    if i < n
        @test top(s) == n - i
    else
        @test_throws ErrorException top(s)
    end
    @test isempty(s) == (i == n)
    @test length(s) == n - i
end

# Queue

s = Queue(Int, 5)
n = 100

@test length(s) == 0
@test isempty(s)
@test_throws ErrorException front(s)
@test_throws ErrorException back(s)
@test_throws ErrorException dequeue!(s)

for i = 1 : n
    enqueue!(s, i)
    @test front(s) == 1
    @test back(s) == i
    @test !isempty(s)
    @test length(s) == i
end

for i = 1 : n
    x = dequeue!(s)
    @test x == i
    if i < n
        @test front(s) == i + 1
        @test back(s) == n
    else
        @test_throws ErrorException front(s)
        @test_throws ErrorException back(s)
    end
    @test isempty(s) == (i == n)
    @test length(s) == n - i
end


text:
@testset "Requests" begin

@testset "getentry for component" begin
    martinez = Data.martinez()

    @test getentry(martinez, "methane") isa Dict
    @test_throws CubicEoSDatabase.NotFoundError getentry(martinez, "philosophic_stone")
    @test getentry(martinez, "methane") == Dict(
        :name => "methane",
        :molecular_mass => 0.01604276,
        :critical_temperature => 191.0,
        :critical_pressure => 4600000.0,
        :acentric_factor => 0.01,
        :critical_compressibility => 0.288,
        :number_carbons => 1,
    )

    eos = Data.brusilovsky_comp()
    @test getentry(eos, "methane") isa Dict
    @test_throws CubicEoSDatabase.NotFoundError getentry(eos, "philosophic_stone")
    @test getentry(eos, "methane") == Dict(
        :name => "methane",
        :eos_critical_compressibility => 0.33294,
        :eos_critical_omega => 0.7563,
        :eos_psi => 0.37447
    )
    @test getentry(eos, "n-octane") == Dict(
        :name => "n-octane",
        :eos_critical_compressibility => "C5+",
        :eos_critical_omega => "C5+",
        :eos_psi => "C5+"
    )
end
@testset "getentry for binary" begin
    binary = Data.brusilovsky_mix()
    @test getentry(binary, "methane", "n-pentane") isa Dict
    @test_throws CubicEoSDatabase.NotFoundError getentry(binary, "methane", "philosophic_stone")

    @testset "Ordered requests" begin
        @test_throws CubicEoSDatabase.NotFoundError getentry(binary, "n-pentane", "methane", keeporder=true)
        @test getentry(binary, "methane", "n-pentane") == getentry(binary, "methane", "n-pentane", keeporder=true)
    end

    @test getentry(binary, "methane", "n-pentane") == Dict(
        :comp1 => "methane",
        :comp2 => "n-pentane",
        :constant => 0.001,
        :linear => 0.000604,
        :quadratic => 0
    )
end

@testset "getmatrix" begin
    db = Data.brusilovsky_mix()
    @test getmatrix(db, ("methane", "propane", "ethane")) isa Dict{Symbol,Matrix{Float64}}
    @test_throws CubicEoSDatabase.NotFoundError getmatrix(db, ("methane", "philosophic_stone"))

    response = getmatrix(db, ("methane", "propane", "ethane"))
    @test response[:constant] == [0 0.019 -0.015; 0.019 0 -0.015; -0.015 -0.015 0]
    @test response[:linear] == [0 0.000502 0.000123; 0.000502 0 0; 0.000123 0 0]
    let q=-0.41e-5
        @test response[:quadratic] == [0 0 q; 0 0 0; q 0 0]
    end
end

end  # @testset


text:
#########
# Tasks #
#########

mutable struct DelayedTask
    # Fields for use in processed task ready to be recorded
    used_modules::Vector
    code::String
    value_names::Vector{Tuple{ValueId,String}}
    effects::Dict{ValueId,String}
    pa_union::Vector{PartitionAnnotation} # Enumeration of applicable PAs
    # Fields for use in task yet to be processed in a call to `compute`
    partitioned_using_func::Union{Function,Nothing}
    partitioned_with_func::Union{Function,Nothing}
    mutation::Dict{Future,Future} # This gets converted to `effects`
end

DelayedTask() = DelayedTask([], "", [], Dict(), [PartitionAnnotation()], nothing, nothing, Dict())

function to_jl(task::DelayedTask)
    return Dict(
        "code" => task.code,
        "value_names" => task.value_names,
        "effects" => task.effects,
        "pa_union" => [to_jl(pa) for pa in task.pa_union],
    )
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
mutable struct PathSample
    path::Array
    g::Real # likelihood
    length::Real # path length
    N::Integer # number of samples
end

################################################################################
# The Main Functions
################################################################################

# The ultimate form of the function
function monte_carlo_path_sampling(origin::T, destination::T, adj_mtx::Array{T,2}, link_length_dict::Dict, N1=5000, N2=10000) where {T<:Integer}
    samples = length_distribution_method(origin, destination, adj_mtx, link_length_dict, N1, N2)

    # no_path_est, x_data, y_data = estimate_cumulative_count(samples, N2)
    # return no_path_est, x_data, y_data

    return samples
end

# - - - - - -

# Useful for testing the instances in
# Roberts, B., & Kroese, D. P. (2007). Estimating the Number of st Paths in a Graph. J. Graph Algorithms Appl., 11(1), 195-214.
# http://dx.doi.org/10.7155/jgaa.00142

# This function returns samples of paths.
function monte_carlo_path_sampling(origin::T, destination::T, adj_mtx::Array{T,2}, N1=5000, N2=10000) where {T<:Integer}
    link_length_dict = getLinkLengthDict(adj_mtx)
    return monte_carlo_path_sampling(origin, destination, adj_mtx, link_length_dict, N1, N2)
end

# This function returns the estimated number of paths, based on sampled paths.
function monte_carlo_path_number(origin::T, destination::T, adj_mtx::Array{T,2}, N1=5000, N2=10000) where {T<:Integer}
    link_length_dict = getLinkLengthDict(adj_mtx)
    samples = monte_carlo_path_sampling(origin, destination, adj_mtx, link_length_dict, N1, N2)
    return estimate_number(samples)
end

# - - - - - -

# Useful for experimenting many realistic road networks
function monte_carlo_path_sampling(origin::T, destination::T, start_node::Array{T,1},  end_node::Array{T,1}, link_length, N1=5000, N2=10000) where {T<:Integer}
    link_length_dict = getLinkLengthDict(start_node, end_node, link_length)
    adj_mtx = getAdjacency(start_node, end_node)

    return monte_carlo_path_sampling(origin, destination, adj_mtx, link_length_dict, N1, N2)
end

################################################################################
################################################################################

function estimate_number(samples::Array{PathSample,1})
    no_path_est = 0
    for s in samples
        no_path_est += 1 / s.g
    end
    no_path_est = no_path_est / samples[1].N
    return no_path_est
end

function estimate_cumulative_count(samples::Array{PathSample,1}, option=:uniform)
    no_path_est = 0
    path_lengths = Array{Float64,1}(undef, 0)
    likelihoods = Array{Float64,1}(undef, 0)
    for s in samples
        no_path_est += 1 / s.g
        push!(path_lengths, s.length)
        push!(likelihoods, s.g)
    end
    no_path_est = no_path_est / samples[1].N

    N_data = 100
    # option == :uniform
    x_data = collect(range(minimum(path_lengths), stop=maximum(path_lengths), length=N_data))
    if option == :unique
        x_data = sort(unique(path_lengths))
    elseif option == :first_quarter
        x_q1 = range(minimum(path_lengths), stop=0.25*maximum(path_lengths), length=N_data/2)
        x_q234 = range(0.25*maximum(path_lengths), stop=maximum(path_lengths), length=N_data/2)
        x_data = append!(collect(x_q1), collect(x_q234))
    end

    y_data = similar(x_data)

    for i=1:length(x_data)
        idx = findall(x-> x<=x_data[i], path_lengths)
        y_data[i] = sum(1 ./ likelihoods[idx])
    end
    y_data =  y_data / samples[1].N

    return x_data, y_data
end

################################################################################
################################################################################

# Algorithm 1 (Naive Path Generation) of Roberts and Kroese (2007)
function naive_path_generation(origin, destination, adj_mtx, N1)

    naive_samples = PathSample[]
    for k=1:N1
        adj = copy(adj_mtx)

        # 1
        x = [origin]
        g = 1
        current = origin

        #2
        adj[:,origin] .= 0

        while current != destination
            #3
            V = Int[]
            for i=1:size(adj,1) # number of nodes
                if adj[current,i] == 1
                    push!(V, i)
                end
            end
            if length(V)==0
                break
            end

            #4
            next = rand(V)
            x = [x; next]

            #5
            current = next
            adj[:,next] .= 0
            g = g / length(V)
        end

        if x[end] == destination
            this_sample = PathSample(x, g, 0.0, N1)
            push!(naive_samples, this_sample)
        end
    end

    return naive_samples
end

# Computing \hat{l}_k as in Equation (4) of Roberts and Kroese (2007)
function length_distribution_vector(naive_samples::Array{PathSample}, adj::Array{Int,2}, no_node::Int, destination::Int)
    # Length-Distribution
    # l_hat = numerator / denominator

    # adj = getAdjacency(start_node, end_node)
    # no_node = getNoNode(start_node, end_node)

    l_hat = Array{Float64}(undef, no_node-1)
    for k = 1:length(l_hat)
        numerator = 0.0
        denominator = 0.0

        for s in naive_samples
            if length(s.path)-1 == k
                numerator += 1 / s.g
            end

            if length(s.path)-1 >= k
                if adj[ s.path[k], destination ]==1
                    denominator += 1 / s.g
                end
            end
        end

        l_hat[k] = numerator / denominator
    end

    return l_hat
end

# Algorithm 2 (Length-Distribution Method) of Roberts and Kroese (2007)
function length_distribution_method(origin, destination, adj_mtx, link_length_dict, N1, N2)

    no_node = size(adj_mtx,1)

    naive_samples = naive_path_generation(origin, destination, adj_mtx, N1)
    l_hat = length_distribution_vector(naive_samples, adj_mtx, no_node, destination)

    # Algorithm 2
    # Naive Path Generation
    better_samples = PathSample[]
    for k=1:N2
        adj = copy(adj_mtx)

        #2
        x = [origin]
        g = 1
        current = origin

        #3
        adj[:,origin] .= 0

        next = []
        while current != destination
            t = length(x)

            #4
            if adj[current, destination] == 1
            #4. A(c,n)=1

                if sum(adj, dims=2)[current] == 1
                #4. destination is the only available vertex adjacent to current
                    next = destination
                else
                #4. If there are other vertices adjacent to current
                    if rand() <= l_hat[t]
                    #   choose the next vertext to be destination with probability l_hat[t]
                        g = g * l_hat[t]
                        next = destination
                    else
                        g = g * (1 - l_hat[t])
                        adj[current, destination] = 0
                    end
                end

            end

            if next == destination
                x = [x; next]
                break

            else
            #5
                V = Int[]
                for i=1:size(adj,1) # number of nodes
                    if adj[current,i] == 1
                        push!(V, i)
                    end
                end
                if length(V)==0
                    break
                end

                #6
                next = rand(V)
                x = [x; next]

                #7
                current = next
                adj[:,next] .= 0
                g = g / length(V)

            end

        end

        if x[end] == destination
            this_sample = PathSample(x, g, getPathLength(x, link_length_dict), N2)
            push!(better_samples, this_sample)
        end
    end

    return better_samples
end


text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type RequestHeaderProto
    methodName::AbstractString
    declaringClassProtocolName::AbstractString
    clientProtocolVersion::UInt64
    RequestHeaderProto() = (o=new(); fillunset(o); o)
end #type RequestHeaderProto
const __req_RequestHeaderProto = Symbol[:methodName,:declaringClassProtocolName,:clientProtocolVersion]
meta(t::Type{RequestHeaderProto}) = meta(t, __req_RequestHeaderProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RequestHeaderProto) = ProtoBuf.protohash(v)
isequal(v1::RequestHeaderProto, v2::RequestHeaderProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RequestHeaderProto, v2::RequestHeaderProto) = ProtoBuf.protoeq(v1, v2)

export RequestHeaderProto


text:
# Author: Gabrio Rizzuti, grizzuti33@gatech.edu
# Date: September 2020

using InvertibleNetworks, LinearAlgebra, Test, Statistics

# Input
nx = 28
ny = 28*2
nc = 4
n_hidden = 64
batchsize = 5

# Initialization network
logdet = false
X1 = rand(Float32, nx, ny, nc, batchsize)
N1 = ActNorm(nc; logdet=logdet); Y1 = N1.forward(X1)
X2 = rand(Float32, nx, ny, nc, batchsize)
N2 = ActNorm(nc; logdet=logdet); Y2 = N2.forward(X2)
ΔX = X1-X2
Δθ = get_params(N1)-get_params(N2)
ΔY = Y1-Y2

# Jacobian initialization X×θ↦Y
N = ActNorm(nc; logdet=logdet); N.forward(rand(Float32, nx, ny, nc, batchsize))
X = rand(Float32, nx, ny, nc, batchsize)
J = Jacobian(N, X; io_mode="X×θ↦Y")

# Evaluation
JΔ = J*(ΔX, Δθ)

# Adjoint
Jadj = adjoint(J)
ΔX_, Δθ_ = Jadj*ΔY

# Adjoint Test
a = dot(ΔX, ΔX_)+dot(Δθ, Δθ_)
b = dot(ΔY, JΔ)
@test isapprox(a, b; rtol=1f-3)

# Jacobian initialization θ↦Y
N = ActNorm(nc; logdet=logdet); N.forward(rand(Float32, nx, ny, nc, batchsize))
X = rand(Float32, nx, ny, nc, batchsize)
J = Jacobian(N, X; io_mode="θ↦Y")

# Evaluation
JΔ = J*Δθ

# Adjoint
Jadj = adjoint(J)
Δθ_ = Jadj*ΔY

# Adjoint Test
a = dot(Δθ, Δθ_)
b = dot(ΔY, JΔ)
@test isapprox(a, b; rtol=1f-3)


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f1c9d31f-c565-4e96-8801-22352b014397
import Pkg

# ╔═╡ f6daaa32-fe10-4a20-b30e-ba31b63c6b90
Pkg.add("Distributions")

# ╔═╡ 338200b2-c10c-4e7d-8a24-bfd32e6757e5
Pkg.add("DataFrames")

# ╔═╡ 525ba450-c635-4bb9-a2f3-aeb12d5c30d3
Pkg.add("CSV")

# ╔═╡ 64adb0f4-b695-4be1-9b02-c42562951b39
Pkg.add("Plots")

# ╔═╡ 55405eea-1c12-4060-81f5-cbcc082a992c
Pkg.add("IterTools")

# ╔═╡ 7e392b1a-31a8-4364-ac93-742221501f1c
Pkg.add("PlutoUI")

# ╔═╡ b708f1d2-e8dd-40ee-92fd-cb4ae02a6a9a
Pkg.add("StatsBase")

# ╔═╡ 34a840a1-7050-4719-a81e-c4545551ec7d
using Random, Distributions, StatsBase

# ╔═╡ 83920316-034e-4ca1-94fa-5aea4d4d117b
using LinearAlgebra

# ╔═╡ 8287cc85-37d3-4b53-ab65-0aef5d996a72
using Plots

# ╔═╡ 1196f179-17ff-42d9-a582-1b8498b7c372
using IterTools

# ╔═╡ 2d8d5528-d64d-4026-b5a6-93b66e17e609
using PlutoUI

# ╔═╡ c3929637-0224-4f95-a8ce-b1a159084014
using DataFrames

# ╔═╡ 2f28edec-8461-4eda-a95f-05f4db5b3a81
using CSV

# ╔═╡ 0a435f6c-9f0a-11eb-114f-fb8c733dafbb
md"# Brute Force Simulation"

# ╔═╡ 5c181bf0-f92e-40ac-9193-62c48927f60f
md"""
## Simulation Configuration
"""

# ╔═╡ 0fe833e8-a359-41e2-ba4a-b8313a222182
md"""
### Classical field parameters

$(@bind bigN Slider(1:1000, show_value = true))

Number of emitting atoms

$(@bind temp Slider(293:3000, show_value = true))

Temperature (K)

$(@bind tres NumberField(10:500))

Time resolution (ps)

$(@bind tmax NumberField(1.0:0.1:10.0))

Maximum time (ns)

"""

# ╔═╡ c0e353c3-53e7-4419-b8cf-878f6fbc761d
md"""
### Photon count parameters

$(@bind nbar NumberField(1.0:1000.0;default=10))

Expected counts in time tmax

$(@bind ntot NumberField(1:1000000;default=1000))

Total photon counts overall (approximate)
"""

# ╔═╡ 402b4a2a-d22f-4ee7-ae3d-ca5fe7fda110
md"""
Number of classical simulations: $(convert(Integer,ceil(ntot/nbar)))
"""

# ╔═╡ 89b17676-6496-4d37-929d-3f626e329cdd
md"""
### Plotting options

**Available Plots:**
- Single instance of classical calculations $(@bind classicalPlots CheckBox())
- Average of all classical calculations $(@bind classicalAvgPlots CheckBox())
- Total counts from all time windows $(@bind countsPlot CheckBox())
- Autocorrelation of photon count time series  $(@bind corrPlot CheckBox())

#### Make Plots [$(@bind makePlots CheckBox())]

"""

# ╔═╡ 30f9f516-568b-4e73-8aa2-4f5bb48f9ab4
md"""
## Plots
"""

# ╔═╡ 57da1ea1-5a5d-49d1-88f2-2624a14f8dff
if makePlots && classicalPlots
	md"""
	### Single instance of classical calculations:
	"""
end

# ╔═╡ 1e48204e-23d8-4581-a160-7ba84b83b8d8
if makePlots && classicalAvgPlots
	md"""
	### Average of all classical calculations:
	"""
end

# ╔═╡ 0f36d9b7-6ba6-402d-9601-53b16e952110
if makePlots && countsPlot
	md"""
	### Total counts:
	"""
end

# ╔═╡ c54073d5-b58a-46ed-b086-d6347c1dbcda
if makePlots && corrPlot
	md"""
	### Autocorrelation of photon counts

	Note that the τ = 0 bin matches the number of simulations. This only happens because we are using the autocorrelation for a single beam. This will not be true when I add a beam splitter.
	"""
end

# ╔═╡ 13e70af4-146d-43e7-9fd6-ec8f901185d4
md"""
## Calculations
"""

# ╔═╡ 807e056a-3c89-46f9-8dc2-5dde04023ce2
md"""
### Classical calculations
"""

# ╔═╡ f2e2b9fd-9d21-410c-89d8-2286b163bb2a
# these parameters are needed for all simulations
begin
	# Balmer-α lines specified here
	ωM = 2*π*[456811.0, 456812.0]
	# magnitude of each line
	mag = convert(Vector{ComplexF64},ones(length(ωM)))
	# calculate line differences
	ΔM = ωM .- ωM[1]
	# generate times in tres ps intervals up to 2*tmax
	times = collect(0:tres*1e-3:2*tmax);
	# limit the window to tmax to avoid correlation cutoff
	window = convert(Integer,floor(length(times)/2));
	# τ is just the times up to our window
	τ = times[1:window];
end;

# ╔═╡ 7ac53f6b-15ee-440d-acae-005dc3acce51
md"""
f0 (GHz) = $(ωM[1]/(2*π))

Δf (GHz) = $(round((ωM[2]-ωM[1])/(2*π),digits=2))

σ (GHz) = $(round(sqrt(9.178e-14*temp)*ωM[1]/(2*π),digits=2))
"""

# ╔═╡ a3cd4171-5bc5-4988-a5ef-428dd6fa31ab
if makePlots
	md"""
	Size: $(@bind windowSize Slider(1:window,default=window))
	"""
end

# ╔═╡ e2fab25a-ccde-44e8-a3bf-e402fc32d2ba
if makePlots
	md"""
	Position: $(@bind windowPos Slider(0:(window-windowSize)))
	"""
end

# ╔═╡ 8b607526-47b0-4a8a-aa16-2aee351b590f
if makePlots
	md"""
	**Window:**

	Window size (ns) = $(tres*windowSize/1000)

	Window position (ns) = $(windowPos*tres/1000)
	"""
end

# ╔═╡ 21e16aae-1b50-4f31-936a-17525352afbe
if makePlots
	# our calculated g2τ
	g2τCalc = 1
	Emag2 = real.(mag .* conj(mag))
	Emag4 = Emag2 .* Emag2
	sumEmag2 = sum(Emag2)
	sumEmag4 = sum(Emag4)
	term2 = sumEmag4/(bigN*sumEmag2^2)
	g2τCalc -= term2

	term3 = sum(Emag2 .* exp.(-im*ΔM .* transpose(τ)),dims=1)/sumEmag2
	term3 = real.(term3 .* conj(term3))

	kbOverMhC2 = 9.178e-14;
	σ = sqrt(kbOverMhC2*temp)*ωM[1]
	stauAvg = transpose(bigN .+ bigN*(bigN-1)*exp.(-σ^2*τ .^2))
	term3 = term3 .* stauAvg/bigN^2

	g2τCalc =  g2τCalc .+ term3
end;

# ╔═╡ 173ca217-29e6-4c5f-8995-0bd5b62c32dd
md"""
- Our calculated version

$g^{(2)}(\tau) =1-\frac{1}{N}\frac{\sum_{m=1}^M|\mathcal{E}|_m^4}{\left(\sum_{m=1}^M|\mathcal{E}_m|^2\right)^2}+\left|\frac{\sum_{m=1}^M|\mathcal{E}_m|^2e^{-i\Delta_m\tau}}{\sum_{m=1}^M|\mathcal{E}|_m^2}\right|^2\frac{\langle S(\tau)\rangle_\omega}{N^2}$
"""

# ╔═╡ 689bc5a6-2469-41eb-99ce-c300dd73874e
md"""
### Calculate photon counts
"""

# ╔═╡ f2a04b1f-f61c-45d4-a67a-1609f53d452e
md"""
Photon counts are calculated by treating the intensity in each time bin as the average photon count rate, then sampling from a poisson distribution with that average count rate.
"""

# ╔═╡ 730e47f5-1c6a-4fd3-a83f-f877b14cde06
md"""
### Calculate autocorrelation of photon counts

Note that I only look at whether two bins both have counts or not when calculating the autocorrelation. I **do not** look at *how many* counts there are in each bin.
"""

# ╔═╡ 282a7d64-8d67-4a10-9aa9-d64be7c5dd54
md"""
## Functions
"""

# ╔═╡ d2ec427e-50f8-47d2-9a44-cb0d5c0dcb9b
"""
	function countDeltaTimes(τ::Vector,γCounts::Vector)

Returns an array of τ values for which the γCounts autocorrelation is non-zero.
"""
function countTimes(times::Vector,γCounts::Vector)
	out = Vector{Real}(undef,0)
	for (i,counts) in enumerate(γCounts)
		if counts != 0
			countTimes = times[i]*ones(counts)
			out = vcat(out,countTimes)
		end
	end
	return out
end

# ╔═╡ 26e08020-6775-431f-b33a-0efaad94189a
"""
	eFieldParams64(mag::Array{S},ΔM::Array{T},ωN::Array{T},ϕ::Array{S} ) where {T<:Real, S<:Complex}

Static parameters for the electric field
"""
struct eFieldParams
	mag::Vector
	ΔM::Vector
	ωN::Vector
	ϕ::Matrix

	function eFieldParams(mag::Vector,ΔM::Vector,ωN::Vector,ϕ::Matrix )
		@assert length(mag) == length(ΔM) "Number of magnitudes must match number of emission lines"
		@assert (size(ϕ)[1] == length(ΔM) && size(ϕ)[2] == length(ωN)) "Must have a unique phase for each n and m"
		new(
			convert(Vector{Complex},mag),
			convert(Vector{Real},ΔM),
			convert(Vector{Real},ωN),
			convert(Matrix{Real},ϕ)
		)

	end
end

# ╔═╡ 8175ac47-b261-4942-a70d-f127f22aee74
"""
	function electricField(t::Real,params::eFieldParams)

Returns the electric field value at time t
"""
function electricField(t::Real,params::eFieldParams)
	# generate frequencies
	ωNM = transpose(params.ωN) .+ params.ΔM
	# add the phase
	exponentnm = -im*(t*ωNM+params.ϕ)
	# put them in the exponent
	enm = exp.(exponentnm)
	# multiply by the field magnitude
	fieldnm = params.mag .* enm
	# add it all together
	return sum(ivec(fieldnm))
end

# ╔═╡ d8c63dd0-d6dd-4234-8515-caf5803af683
"""
	function ωnDoppler(ω0::Real,N::Integer,temp::Real,seed::Integer = -1)

Generates N doppler shifted frequencies around frequency ω0 for a source at temperature temp. Seed optional for reproducible results.
"""
function ωnDoppler(ω0::Real,N::Integer,temp::Real,seed::Integer = -1)
	rng = MersenneTwister()
	if seed != -1
		rng = MersenneTwister(seed)
	end
	kbOverMhC2 = 9.178e-14;
	σ = sqrt(kbOverMhC2*temp)*ω0
	d = Normal(ω0,σ)
	return rand(rng,d,N)
end

# ╔═╡ b237ca1e-1719-420c-820d-67b8e18c7a57
"""
	function correlate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

Calculates correlation between vectors u and v with given offset. Specify averaging window to limit range of correlation. If the window extends beyond the end of one vector, it treats out-of-bounds indices as zero.
"""
function correlate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

	@assert offset <= length(u) "Offset out of bounds"
	@assert window <= length(u) && window <= length(v) "Window must be smaller than input vector lengths"

	if window == -1
		window = length(u)
	end

	v1 = view(u,1:window)
	v2 = view(v,1+offset:min(window+offset,length(v)))
	if window+offset > length(v)
		v2 = vcat(v2,zeros(window+offset-length(v)))
	end

	return dot(v1,v2)/window
end

# ╔═╡ 3620f708-3eac-4cd8-90db-72ac383b5917
"""
	function autocorrelate(u::Vector{T},offset::Integer, window::Integer = -1) where {T<:Number}

Calculates correlation of vector u with itself.
"""
function autocorrelate(u::Vector{T},offset::Integer, window::Integer = -1) where {T<:Number}
	correlate(u,u,offset,window)
end

# ╔═╡ 18f7bc41-751f-4283-8a90-af3ad217b58b
# if ONLY the classical plots are desired, then this calculates a single instance of the classical fields and correlations
if makePlots && (classicalPlots && !(classicalAvgPlots || countsPlot || corrPlot))
	# generate N doppler broadened frequencies
	ω1Doppler = ωnDoppler(ωM[1],bigN,temp)
	# generate N*M random phases
	ϕ1nm = 2*π*rand(Float64,(length(ωM),bigN))
	# construct field parameter object
	testParams1 = eFieldParams(mag,ΔM,ω1Doppler,ϕ1nm)

	# calculate electric field vs time
	e1fieldt = map(t->electricField(t,testParams1),times)

	# calculate g1τ
	g1τ1Norm = correlate(e1fieldt,conj.(e1fieldt),0,window);
	g1τ1 = abs.(map(i->correlate(e1fieldt,conj.(e1fieldt),i,window),collect(0:window-1))/g1τ1Norm)

	# calculate the intensity vs time
	intensity1t = real.( e1fieldt .* conj(e1fieldt))

	# calculate g2τ
	g2τ1Norm = mean(intensity1t)^2
	g2τ1 = map(i->autocorrelate(intensity1t,i,window),collect(0:window-1))/g2τ1Norm
end;

# ╔═╡ 39e8d2c6-c8de-4e0a-b531-60dbaadd6aa1
# the classical average plots and photon-based plots require multiple instances of the classical field calculations
if makePlots && (classicalAvgPlots || countsPlot || corrPlot)
	# calculate number of trials from the total desired photon counts and the average photon count per trial
	nTrials = convert(Integer,ceil(ntot/nbar))
	# make an array of the average photon counts per trial for array broadcasting
	nPerTrial = bigN*ones(Integer,nTrials)
	# calculate nTrials instances of doppler broadened frequencies
	ωDoppler = ωnDoppler.(ωM[1],nPerTrial,temp);
	# calculate nTrials instances of random phases
	ϕnm = map(n->2*π*rand(Float64,(length(ωM),n)),nPerTrial);
	# generate nTrials instances of field parameters
	testParams = map((x,y)->eFieldParams(mag,ΔM,x,y),ωDoppler,ϕnm);

	# calculate nTrials instances of the time dependent electric field
	efieldt = map(x->map(t->electricField(t,x),times),testParams);

	# calculate nTrials instances of g1τ
	g1τNorm = map(eft->correlate(eft,conj.(eft),0,window),efieldt);
	g1τ = map((eft,normG1τ)->abs.(map(i->correlate(eft,conj.(eft),i,window),collect(0:window-1))/normG1τ),efieldt,g1τNorm);

	# calculate nTrials instances of the time dependent intensity
	intensityt = map(eft->real.( eft .* conj(eft)),efieldt);

	# calculate nTrials instances of g2τ
	g2τNorm = map(intens->mean(intens)^2,intensityt)
	g2τ = map((intens,normG2τ)->map(i->autocorrelate(intens,i,window),collect(0:window-1))/normG2τ,intensityt,g2τNorm)

	# pick out one instance of everything for plotting
	e1fieldtM = efieldt[1]
	g1τ1M = g1τ[1]
	intensity1tM = intensityt[1]
	g2τ1M = g2τ[1]
end;

# ╔═╡ 9e013576-7682-4d6c-8caa-e3f2c1034f73
if makePlots && classicalPlots
	if !(classicalAvgPlots || countsPlot || corrPlot)
		e1fieldPlot = plot(times[1+windowPos:windowSize+windowPos],real.(e1fieldt)[1+windowPos:windowSize+windowPos],label=false,title="E(t)")
		xlabel!(e1fieldPlot,"t (ns)")

		eCorr1Plot = plot(τ[1+windowPos:windowSize+windowPos],g1τ1[1+windowPos:windowSize+windowPos],label=false,title="g1(τ)")
		xlabel!(eCorr1Plot,"τ (ns)")

		intensity1tPlot = plot(times[1+windowPos:windowSize+windowPos],intensity1t[1+windowPos:windowSize+windowPos],label=false,color=2,title="I(t)")
		xlabel!(intensity1tPlot,"t (ns)")

		iCorr1Plot = plot(τ[1+windowPos:windowSize+windowPos],g2τ1[1+windowPos:windowSize+windowPos],label="data",color=2,title="g2(τ)")
		plot!(iCorr1Plot,τ[1+windowPos:windowSize+windowPos],transpose(g2τCalc)[1+windowPos:windowSize+windowPos],label="calc")
		xlabel!(iCorr1Plot,"τ (ns)")

		plot(e1fieldPlot,intensity1tPlot,eCorr1Plot,iCorr1Plot,layout=(2,2),size=(800,800))
	else
		e1fieldPlot = plot(times[1+windowPos:windowSize+windowPos],real.(e1fieldtM)[1+windowPos:windowSize+windowPos],label=false,title="E(t)")
		xlabel!(e1fieldPlot,"t (ns)")

		eCorr1Plot = plot(τ[1+windowPos:windowSize+windowPos],g1τ1M[1+windowPos:windowSize+windowPos],label=false,title="g1(τ)")
		xlabel!(eCorr1Plot,"τ (ns)")

		intensity1tPlot = plot(times[1+windowPos:windowSize+windowPos],intensity1tM[1+windowPos:windowSize+windowPos],label=false,color=2,title="I(t)")
		xlabel!(intensity1tPlot,"t (ns)")

		iCorr1Plot = plot(τ[1+windowPos:windowSize+windowPos],g2τ1M[1+windowPos:windowSize+windowPos],label="data",color=2,title="g2(τ)")
		plot!(iCorr1Plot,τ[1+windowPos:windowSize+windowPos],transpose(g2τCalc)[1+windowPos:windowSize+windowPos],label="calc",color=3)
		xlabel!(iCorr1Plot,"τ (ns)")

		plot(e1fieldPlot,intensity1tPlot,eCorr1Plot,iCorr1Plot,layout=(2,2),size=(800,800))
	end
end

# ╔═╡ b5342349-5b28-4918-bedb-effa46df10c3
if makePlots && classicalPlots
	if !(classicalAvgPlots || countsPlot || corrPlot)
		classicalFieldDf = DataFrame(
			time              = times,
			realElectricField = real.(e1fieldt),
			imagElectricField = imag.(e1fieldt),
			intensity         = intensity1t
			)
		classicalCorrDf = DataFrame(
			tau               = τ,
			g1tau             = g1τ1,
			g2tau             = g2τ1
			)
	else
		classicalFieldDf = DataFrame(
			time              = times,
			realElectricField = real.(e1fieldtM),
			imagElectricField = imag.(e1fieldtM),
			intensity         = intensity1tM
			)
		classicalCorrDf = DataFrame(
			tau               = τ,
			g1tau             = g1τ1M,
			g2tau             = g2τ1M
			)
	end
	classicalFieldString = string(collect(CSV.RowWriter(classicalFieldDf))...)
	classicalCorrString  = string(collect(CSV.RowWriter(classicalCorrDf))...)
	md"""
	E-field/Intensity:
	$(DownloadButton(classicalFieldString,"classical_field_single.csv"))

	Correlations:
	$(DownloadButton(classicalCorrString,"classical_correlations_single.csv"))
	"""
end

# ╔═╡ a3e82734-c1d3-4e9b-8603-ce6bc1cb5018
"""
	function singleDeltaTimes(τ::Vector,γCounts::Vector)

Returns an array of τ values for which the γCounts autocorrelation is non-zero.
"""
function singleDeltaTimes(τ::Vector,γCounts::Vector)
	return τ[map(i->autocorrelate(γCounts,i,length(τ)) > 0 ? true : false,collect(0:length(τ)-1)) ]
end

# ╔═╡ 48fd09eb-f5a0-4bc4-bdbf-3648366978a1
"""
	function γIntensity(intensity::Vector,nbar::Real)

Calculates the photon count rate in each bin of an intensity histogram
"""
function γIntensity(intensity::Vector,nbar::Real)
	nintensity = intensity/sum(intensity)
	return nbar*nintensity
end

# ╔═╡ 31b92f9f-6be0-48cb-9a7e-bfc5edb79bd2
"""
	function poissonCount(nbar::Real)

Returns Poisson distributed counts for average count rate nbar
"""
function poissonCount(nbar::Real)
	p = exp(-nbar)
	s = p
	r = rand()
	count = 0
	while r > s
		count += 1
		p *= nbar/count
		s += p
	end
	return count
end

# ╔═╡ 045fbba2-5ae7-4a41-9e8e-343a5ab41c7b
if makePlots && (countsPlot || corrPlot)
	γCounts = map(intens->poissonCount.(γIntensity(intens,nbar*2)),intensityt)
	γcountTimes = map(γCt->countTimes(times,γCt),γCounts)
	flatγCountTimes = vcat(γcountTimes...)
end;

# ╔═╡ 92d93e89-1d36-4bf0-b3f8-d4370e08857b
if makePlots && countsPlot
	countHist =  fit(Histogram,flatγCountTimes, vcat(τ,tres/10)[1+windowPos:windowSize+windowPos],closed=:left)
	countPlot = plot(countHist,label=false)
	title!(countPlot,"Total counts displayed = $(sum(countHist.weights))")
	xlabel!(countPlot,"t (ns)")
end

# ╔═╡ 9561be7c-ab38-493a-8b33-cd7fb7aec495
if makePlots && countsPlot
	countEdgesDf   = DataFrame(count_bin_edges = countHist.edges[1])
	countWeightsDf = DataFrame(count_bin_weights = countHist.weights)
	countDataDf    = DataFrame(time_ns = flatγCountTimes)

	countEdgesString   = string(collect(CSV.RowWriter(countEdgesDf))...)
	countWeightsString = string(collect(CSV.RowWriter(countWeightsDf))...)
	countDataString    = string(collect(CSV.RowWriter(countDataDf))...)

	md"""
	Bin edges:
	$(DownloadButton(countEdgesString,"photon-count_histogram_bin-edges.csv"))

	Bin weights:
	$(DownloadButton(countWeightsString,"photon-count_histogram_bin-weights.csv"))

	Arrival times:
	$(DownloadButton(countWeightsString,"photon-arrival-time_data.csv"))
	"""
end

# ╔═╡ a62a3b1c-4980-492c-95d5-20681ca343ae
if makePlots && corrPlot
	correlationTimes = map(γCt->singleDeltaTimes(τ,γCt),γCounts)
	flatCorrelationTimes = vcat(correlationTimes...)
end;

# ╔═╡ 25b20ebd-3fc1-4326-9848-1adf5afb768d
if makePlots && corrPlot
	corrHist = fit(Histogram, flatCorrelationTimes,vcat(τ,tres/10)[1+windowPos:windowSize+windowPos],closed=:left)
	plot(corrHist,label=false)

end

# ╔═╡ ac647fa5-3294-4328-9b9f-fae54ef54174
if makePlots && corrPlot
	corrEdgesDf   = DataFrame(corr_bin_edges = corrHist.edges[1])
	corrWeightsDf = DataFrame(corr_bin_weights = corrHist.weights)
	corrDataDf    = DataFrame(corr_times = flatCorrelationTimes)

	corrEdgesString = string(collect(CSV.RowWriter(corrEdgesDf))...)
	corrWeightsString = string(collect(CSV.RowWriter(corrWeightsDf))...)
	corrDataString = string(collect(CSV.RowWriter(corrDataDf))...)

	md"""
	Bin edges:
	$(DownloadButton(corrEdgesString,"photon-correlation_histogram_bin-edges.csv"))

	Bin weights:
	$(DownloadButton(corrWeightsString,"photon-correlation_histogram_bin-weights.csv"))

	Correlation times:
	$(DownloadButton(corrWeightsString,"photon-correlation-tau_data.csv"))
	"""

end

# ╔═╡ 3edfcf46-9f6d-4ff0-8179-f9a7ddb1897f
"""
	function beCount(nbar::Real)

Returns Bose-Einstein distributed counts for average count rate nbar
"""
function beCount(nbar::Real)
	p = 1/(nbar+1)
	fnbar = p*nbar
	f = p*nbar
	s = p
	r = rand()
	count = 0
	while r>s
		count += 1
		p *= f
		s += p
	end
	return count
end

# ╔═╡ 2328359b-8f27-4a18-8d63-0538dbdfcdff
function vectorAvg(someVector::Vector)
	return +(someVector...)/length(someVector)
end

# ╔═╡ 099f95dc-1dbc-48f6-ab3f-0920e86f0c0a
if makePlots && classicalAvgPlots
	efieldtAvg = vectorAvg(efieldt)
	g1τAvg = vectorAvg(g1τ)
	intensitytAvg = vectorAvg(intensityt)
	g2τAvg = vectorAvg(g2τ)
end;

# ╔═╡ 0e50e055-6c57-40d9-96ca-d8c795213999
if makePlots && classicalAvgPlots
	efieldtAvgPlot = plot(times[1+windowPos:windowSize+windowPos],real.(efieldtAvg)[1+windowPos:windowSize+windowPos],label=false,title="Average E(t): N = $(nTrials)")
	xlabel!(efieldtAvgPlot,"t (ns)")

	eCorrAvgPlot = plot(τ[1+windowPos:windowSize+windowPos],g1τAvg[1+windowPos:windowSize+windowPos],label=false,title="Average g1(τ): N = $(nTrials)")
	xlabel!(eCorrAvgPlot,"τ (ns)")

	intensitytAvgPlot = plot(times[1+windowPos:windowSize+windowPos],intensitytAvg[1+windowPos:windowSize+windowPos],label=false,color=2,title="Average I(t): N = $(nTrials)")
	xlabel!(intensitytAvgPlot,"t (ns)")

	iCorrAvgPlot = plot(τ[1+windowPos:windowSize+windowPos],g2τAvg[1+windowPos:windowSize+windowPos],label="data",color=2,title="Average g2(τ): N = $(nTrials)")
	plot!(iCorrAvgPlot,τ[1+windowPos:windowSize+windowPos],transpose(g2τCalc)[1+windowPos:windowSize+windowPos],label="calc",color=3)
	xlabel!(iCorrAvgPlot,"τ (ns)")

	plot(efieldtAvgPlot,intensitytAvgPlot,eCorrAvgPlot,iCorrAvgPlot,layout=(2,2),size=(800,800))

end

# ╔═╡ c023ca90-8bc8-4183-b0b9-9b1a8bfb4b59
if makePlots && classicalAvgPlots
	classicalFieldAvgDf = DataFrame(
		time                 = times,
		realElectricFieldAvg = real.(efieldtAvg),
		imagElectricFieldAvg = imag.(efieldtAvg),
		intensityAvg         = intensitytAvg
		)
	classicalCorrAvgDf = DataFrame(
		tau      = τ,
		g1tauAvg = g1τAvg,
		g2tauAvg = g2τAvg
		)

	classicalFieldAvgString = string(collect(CSV.RowWriter(classicalFieldAvgDf))...)
	classicalCorrAvgString  = string(collect(CSV.RowWriter(classicalCorrAvgDf))...)
	md"""
	E-field/intensity:
	$(DownloadButton(classicalFieldString,"classical_field_avg.csv"))

	Correlations:
	$(DownloadButton(classicalCorrString,"classical_correlations_avg.csv"))
	"""
end

# ╔═╡ 57b0a316-782e-4f85-8a77-8743abc601f1
md"""
## Load Prerequisites
"""

# ╔═╡ Cell order:
# ╟─0a435f6c-9f0a-11eb-114f-fb8c733dafbb
# ╟─5c181bf0-f92e-40ac-9193-62c48927f60f
# ╟─0fe833e8-a359-41e2-ba4a-b8313a222182
# ╟─7ac53f6b-15ee-440d-acae-005dc3acce51
# ╟─c0e353c3-53e7-4419-b8cf-878f6fbc761d
# ╟─402b4a2a-d22f-4ee7-ae3d-ca5fe7fda110
# ╟─89b17676-6496-4d37-929d-3f626e329cdd
# ╟─30f9f516-568b-4e73-8aa2-4f5bb48f9ab4
# ╟─8b607526-47b0-4a8a-aa16-2aee351b590f
# ╟─a3cd4171-5bc5-4988-a5ef-428dd6fa31ab
# ╟─e2fab25a-ccde-44e8-a3bf-e402fc32d2ba
# ╟─57da1ea1-5a5d-49d1-88f2-2624a14f8dff
# ╟─9e013576-7682-4d6c-8caa-e3f2c1034f73
# ╟─b5342349-5b28-4918-bedb-effa46df10c3
# ╟─1e48204e-23d8-4581-a160-7ba84b83b8d8
# ╟─0e50e055-6c57-40d9-96ca-d8c795213999
# ╟─c023ca90-8bc8-4183-b0b9-9b1a8bfb4b59
# ╟─0f36d9b7-6ba6-402d-9601-53b16e952110
# ╠═92d93e89-1d36-4bf0-b3f8-d4370e08857b
# ╠═9561be7c-ab38-493a-8b33-cd7fb7aec495
# ╠═c54073d5-b58a-46ed-b086-d6347c1dbcda
# ╠═25b20ebd-3fc1-4326-9848-1adf5afb768d
# ╟─ac647fa5-3294-4328-9b9f-fae54ef54174
# ╟─13e70af4-146d-43e7-9fd6-ec8f901185d4
# ╟─807e056a-3c89-46f9-8dc2-5dde04023ce2
# ╠═f2e2b9fd-9d21-410c-89d8-2286b163bb2a
# ╠═21e16aae-1b50-4f31-936a-17525352afbe
# ╟─173ca217-29e6-4c5f-8995-0bd5b62c32dd
# ╠═18f7bc41-751f-4283-8a90-af3ad217b58b
# ╠═39e8d2c6-c8de-4e0a-b531-60dbaadd6aa1
# ╠═099f95dc-1dbc-48f6-ab3f-0920e86f0c0a
# ╟─689bc5a6-2469-41eb-99ce-c300dd73874e
# ╟─f2a04b1f-f61c-45d4-a67a-1609f53d452e
# ╠═045fbba2-5ae7-4a41-9e8e-343a5ab41c7b
# ╟─730e47f5-1c6a-4fd3-a83f-f877b14cde06
# ╠═a62a3b1c-4980-492c-95d5-20681ca343ae
# ╟─282a7d64-8d67-4a10-9aa9-d64be7c5dd54
# ╠═a3e82734-c1d3-4e9b-8603-ce6bc1cb5018
# ╠═d2ec427e-50f8-47d2-9a44-cb0d5c0dcb9b
# ╠═26e08020-6775-431f-b33a-0efaad94189a
# ╠═8175ac47-b261-4942-a70d-f127f22aee74
# ╠═d8c63dd0-d6dd-4234-8515-caf5803af683
# ╠═b237ca1e-1719-420c-820d-67b8e18c7a57
# ╠═3620f708-3eac-4cd8-90db-72ac383b5917
# ╠═48fd09eb-f5a0-4bc4-bdbf-3648366978a1
# ╠═31b92f9f-6be0-48cb-9a7e-bfc5edb79bd2
# ╠═3edfcf46-9f6d-4ff0-8179-f9a7ddb1897f
# ╠═2328359b-8f27-4a18-8d63-0538dbdfcdff
# ╟─57b0a316-782e-4f85-8a77-8743abc601f1
# ╠═f1c9d31f-c565-4e96-8801-22352b014397
# ╠═f6daaa32-fe10-4a20-b30e-ba31b63c6b90
# ╠═338200b2-c10c-4e7d-8a24-bfd32e6757e5
# ╠═525ba450-c635-4bb9-a2f3-aeb12d5c30d3
# ╠═64adb0f4-b695-4be1-9b02-c42562951b39
# ╠═55405eea-1c12-4060-81f5-cbcc082a992c
# ╠═7e392b1a-31a8-4364-ac93-742221501f1c
# ╠═b708f1d2-e8dd-40ee-92fd-cb4ae02a6a9a
# ╠═34a840a1-7050-4719-a81e-c4545551ec7d
# ╠═83920316-034e-4ca1-94fa-5aea4d4d117b
# ╠═8287cc85-37d3-4b53-ab65-0aef5d996a72
# ╠═1196f179-17ff-42d9-a582-1b8498b7c372
# ╠═2d8d5528-d64d-4026-b5a6-93b66e17e609
# ╠═c3929637-0224-4f95-a8ce-b1a159084014
# ╠═2f28edec-8461-4eda-a95f-05f4db5b3a81

max_stars_repo_name: greenspaceexplorer/SpeckleSim

max_issues_repo_name: greenspaceexplorer/SpeckleSim

max_forks_repo_name: greenspaceexplorer/SpeckleSim


text:
@with_kw struct ParticleSystem{T}
    m::T = 87 * 1.661e-27  # in kg
    coupling::T = -32e9 * 6.626e-34 * 1e-12 * 1e12  # in kg um^6 / us^2
    k_b::T = 1.38e-23
    range::Int64 = 6
    temperature::T = 0.0

    c_m::T = coupling / m  # in um^6/us^2
end

function sum_forces!(force::AbstractVector{T}, pos_1::AbstractVector{T}, pos_2::AbstractVector{T},
    particle_system::ParticleSystem{T}) where T
    range = particle_system.range
    force .= force .- range * particle_system.c_m .* (pos_1 .- pos_2) ./ norm(pos_1 .- pos_2)^(range + 2)
end

function sum_forces!(force::AbstractVector{T}, dist::AbstractVector{T},
    particle_system::ParticleSystem{T}) where T
    range = particle_system.range
    force .= force .- range * particle_system.c_m .* dist ./ norm(dist)^(range + 2)
end

const dist_particle_system = zeros(3)
const force_particle_system = zeros(3)

function hamiltonian_explosian!(dv::Array{T, 2}, v::Array{T, 2}, u::Array{T, 2},
    system::ParticleSystem{T}, t::T) where T

    for i in 1:size(u, 2)
        force_particle_system .= 0
        for j in 1:size(u, 2)
            if i == j
                continue
            end
            pos_1 = @view u[:, i]
            pos_2 = @view u[:, j]
            dist_particle_system .= pos_1 .- pos_2
            sum_forces!(force_particle_system, dist_particle_system, system)
        end
        dv[:, i] .= force_particle_system
    end
end

function get_velocity(N::Int64; system=ParticleSystem())
    sigma = sqrt(system.k_b * system.temperature / system.m)
    gaussian_distribution = Normal(0.0, sigma)
    multivariate = product_distribution([gaussian_distribution, gaussian_distribution, gaussian_distribution])
    return rand(multivariate, N)
end

function solve_particle_system(
    positions::Matrix{T}, times::AbstractVector{T};
    system=ParticleSystem()
    ) where T
    N = size(positions)[2]
    velocities = get_velocity(N, system=system)
    prob = SecondOrderODEProblem(
        hamiltonian_explosian!, velocities, positions, (times[1], times[end]), system
        )
    return solve(prob, DPRKN6())
end

function get_J(sol::DifferentialEquations.ODESolution, t)
    pos = sol(t).x[2]
    J = get_J(pos, PowerLaw(6, 60e3))
    return @SMatrix[J]
end

function get_J_itp(
    positions, times;
    system::ParticleSystem=ParticleSystem()
    )

    sol = solve_particle_system(positions, times, system=system)
    itp = interpolate([get_J(sol, t) for t in times], BSpline(Linear()))
    f_itp = Interpolations.scale(itp, times)
    return (t) -> f_itp(t).data[1]
end

max_stars_repo_name: titusfranz/DTWA

max_issues_repo_name: titusfranz/DTWA

max_forks_repo_name: titusfranz/DTWA


text:
import SemialgebraicSets
import PolyJuMP

@testset "Non-symmetric matrix SOS constraint" begin
    @polyvar x
    model = SOSModel()
    err = ErrorException("In `@constraint(model, [1 x; -x 0] in SOSMatrixCone())`: The polynomial matrix constrained to be SOS must be symmetric.")
    @test_throws err @constraint(model, [1 x; -x 0] in SOSMatrixCone())
end

@testset "Printing" begin
    @polyvar x
    model = SOSModel()
    @variable(model, a)
    cref = @constraint(model, a * x^2 >= 1)
    @test sprint(show, MIME"text/plain"(), cref) == "(a)x² + (-1) is SOS"
    @test sprint(show, MIME"text/latex"(), cref) == "\$ (a)x^{2} + (-1) \\text{ is SOS} \$"
    sdref = @constraint(model, a * x^2 in SDSOSCone())
    @test sprint(show, MIME"text/plain"(), sdref) == "(a)x² is SDSOS"
    @test sprint(show, MIME"text/latex"(), sdref) == "\$ (a)x^{2} \\text{ is SDSOS} \$"
    dref = @constraint(model, a * x^2 in DSOSCone())
    @test sprint(show, MIME"text/plain"(), dref) == "(a)x² is DSOS"
    @test sprint(show, MIME"text/latex"(), dref) == "\$ (a)x^{2} \\text{ is DSOS} \$"
end


text:
using Distances, LinearAlgebra
using Test
using KernelFunctions
using PDMats
using Kronecker
dims = [10,5]

A = rand(dims...)
B = rand(dims...)
C = rand(8,9)
K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
Kdiag = [zeros(dims[1]),zeros(dims[2])]
s = rand()
k = SqExponentialKernel()
kt = transform(SqExponentialKernel(),s)
@testset "Kernel Matrix Operations" begin
    @testset "Inplace Kernel Matrix" begin
        for obsdim in [1,2]
            @test kernelmatrix!(K[obsdim],k,A,B,obsdim=obsdim) == kernelmatrix(k,A,B,obsdim=obsdim)
            @test kernelmatrix!(K[obsdim],k,A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
            @test kerneldiagmatrix!(Kdiag[obsdim],k,A,obsdim=obsdim) == kerneldiagmatrix(k,A,obsdim=obsdim)
            @test_throws DimensionMismatch kernelmatrix!(K[obsdim],k,A,C,obsdim=obsdim)
            @test_throws DimensionMismatch kernelmatrix!(K[obsdim],k,C,obsdim=obsdim)
            @test_throws DimensionMismatch kerneldiagmatrix!(Kdiag[obsdim],k,C,obsdim=obsdim)
        end
    end
    @testset "Kernel matrix" begin
        for obsdim in [1,2]
            @test kernelmatrix(k,A,B,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,B,dims=obsdim))
            @test kernelmatrix(k,A,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,dims=obsdim))
            @test kerneldiagmatrix(k,A,obsdim=obsdim) == diag(kernelmatrix(k,A,obsdim=obsdim))
            @test k(A,B,obsdim=obsdim) == kernelmatrix(k,A,B,obsdim=obsdim)
            @test k(A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
            @test KernelFunctions._kernel(k,1.0,2.0) == KernelFunctions._kernel(k,[1.0],[2.0])
            @test_throws DimensionMismatch kernelmatrix(k,A,C,obsdim=obsdim)
        end
    end
    @testset "Transformed Kernel Matrix Operations" begin
        @testset "Inplace Kernel Matrix" begin
            for obsdim in [1,2]
                @test kernelmatrix!(K[obsdim],kt,A,B,obsdim=obsdim) == kernelmatrix(k,s*A,s*B,obsdim=obsdim)
                @test kernelmatrix!(K[obsdim],kt,A,obsdim=obsdim) == kernelmatrix(k,s*A,obsdim=obsdim)
                @test kerneldiagmatrix!(Kdiag[obsdim],kt,A,obsdim=obsdim) == kerneldiagmatrix(k,s*A,obsdim=obsdim)
            end
        end
        @testset "Kernel matrix" begin
            for obsdim in [1,2]
                @test kernelmatrix(kt,A,B,obsdim=obsdim) == kernelmatrix(k,s*A,s*B,obsdim=obsdim)
                @test kernelmatrix(kt,A,obsdim=obsdim) == kernelmatrix(k,s*A,obsdim=obsdim)
                @test kerneldiagmatrix(kt,A,obsdim=obsdim) == kerneldiagmatrix(k,s*A,obsdim=obsdim)
            end
        end
    end
    @testset "KernelSum" begin
        k1 = SqExponentialKernel()
        k2 = LinearKernel()
        ks = k1 + k2
        w1 = 0.4; w2 = 1.2;
        ks2 = KernelSum([k1,k2],weights=[w1,w2])
        @test all(kernelmatrix(ks,A) .== kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(ks+k1,A) .≈ 2*kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(k1+ks,A) .≈ 2*kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(ks,A,B) .== kernelmatrix(k1,A,B) + kernelmatrix(k2,A,B))
        @test all(kerneldiagmatrix(ks,A) .== kerneldiagmatrix(k1,A) + kerneldiagmatrix(k2,A))
        @test all(kernelmatrix(ks2,A) .== w1*kernelmatrix(k1,A) + w2*kernelmatrix(k2,A))
    end
    @testset "KernelProduct" begin
        k1 = SqExponentialKernel()
        k2 = LinearKernel()
        k3 = RationalQuadraticKernel()
        kp = k1 * k2
        kp2 = k1 * k3
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp*k1,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
        @test all(kernelmatrix(k1*kp,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp,A,B) .≈ kernelmatrix(k1,A,B) .* kernelmatrix(k2,A,B))
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kerneldiagmatrix(kp,A) .== kerneldiagmatrix(k1,A) .* kerneldiagmatrix(k2,A))
    end
    @testset "PDMat" begin
        for obsdim in [1,2]
            @test all(Matrix(kernelpdmat(k,A,obsdim=obsdim)) .≈ Matrix(PDMat(kernelmatrix(k,A,obsdim=obsdim))))
            # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
        end
    end
    @testset "Kronecker" begin
        x = range(0,1,length=10)
        X = vcat(collect.(Iterators.product(x,x))'...)
        @test all(collect(kernelkronmat(k,collect(x),2)).≈kernelmatrix(k,X,obsdim=1))
        @test all(collect(kernelkronmat(k,[x,x])).≈kernelmatrix(k,X,obsdim=1))
        @test_throws AssertionError kernelkronmat(LinearKernel(),collect(x),2)
    end
end


text:
using Distributions
using ArrayViews
using Base.LinAlg.BLAS

import Base.length

abstract RBM

typealias Mat{T} AbstractArray{T, 2}
typealias Vec{T} AbstractArray{T, 1}

const UNIT_CLASSES = [:bernoulli, :gaussian]

type Units
    bias::Vector{Float64}
    class::Symbol

    function Units(num::Int, class::Symbol=:bernoulli)
        #@assert class in UNIT_CLASSES
        new(zeros(num), class)
    end

    function Units(params::Dict)
        @assert haskey(params, "num")
        if haskey(params, "class")
            return new(zeros(params["num"]), symbol(params["class"]))
        else
            return new(zeros(params["num"]), :bernoulli)
        end
    end

    function Base.show(io::IO, units::Units)
        num = length(units)
        print(io, "Units(class=$(units.class), num=$num)")
    end
end

length(units::Units) = length(units.bias)
insert!(units::Units, index::Int, item::Any) = insert!(units.bias, index, item)
deleteat!(units::Units, index::Int) = deleteat!(units.bias, index)

type BaseRBM <: RBM
    W::Matrix{Float64}
    visible::Units
    hidden::Units
    dW_prev::Matrix{Float64}
    persistent_chain::Matrix{Float64}
    momentum::Float64
    directory

    function BaseRBM(visible::Units, hidden::Units; sigma=0.001, momentum=0.9, directory="./")
        n_hid = length(hidden)
        n_vis = length(visible)

        new(
            rand(Normal(0, sigma), (n_hid, n_vis)),
            visible, hidden,
            zeros(n_hid, n_vis),
            Array(Float64, 0, 0),
            momentum,
            directory
        )
    end

    function BaseRBM(params::Dict)
        @assert haskey(params, "visible")
        @assert haskey(params, "hidden")

        visible = Units(params["visible"])
        hidden = Units(params["hidden"])
        sigma = 0.001
        momentum = 0.9
        directory="./"

        sigma = haskey(params, "sigma") ? params["sigma"] : sigma
        momentum = haskey(params, "momentum") ? params["momentum"] : momentum
        directory = haskey(params, "directory") ? params["directory"] : directory

        n_vis = length(visible)
        n_hid = length(hidden)

        new(
            rand(Normal(0, sigma), (n_hid, n_vis)),
            visible, hidden,
            zeros(n_hid, n_vis),
            Array(Float64, 0, 0),
            momentum,
            directory
        )
    end

    function Base.show(io::IO, rbm::BaseRBM)
        n_vis = length(rbm.visible)
        n_hid = length(rbm.hidden)
        print(io, "BaseRBM($n_vis, $n_hid, momentum=$(rbm.momentum))")
    end
end

function logistic(x)
    return 1 ./ (1 + exp(-x))
end

function mean_hiddens(rbm::RBM, vis::Mat{Float64})
    p = rbm.W * vis .+ rbm.hidden.bias
    return logistic(p)
end

function sample_hiddens(rbm::RBM, vis::Mat{Float64})
    p = mean_hiddens(rbm, vis)
    return float(rand(size(p)) .< p)
end

function sample_visibles(rbm::RBM, hid::Mat{Float64})
    if rbm.visible.class == :bernoulli
        p = rbm.W' * hid .+ rbm.visible.bias
        p = logistic(p)
        return float(rand(size(p)) .< p)
    elseif rbm.visible.class == :gaussian
        mu = logistic(rbm.W' * hid .+ rbm.visible.bias)
        sigma2 = 0.01                   # using fixed standard diviation
        samples = zeros(size(mu))
        for j=1:size(mu, 2), i=1:size(mu, 1)
            samples[i, j] = rand(Normal(mu[i, j], sigma2))
        end
        return samples
    end
end

function gibbs(rbm::RBM, vis::Mat{Float64}; n_times=1)
    v_pos = vis
    h_pos = sample_hiddens(rbm, v_pos)
    v_neg = sample_visibles(rbm, h_pos)
    h_neg = sample_hiddens(rbm, v_neg)
    for i=1:n_times-1
        v_neg = sample_visibles(rbm, h_neg)
        h_neg = sample_hiddens(rbm, v_neg)
    end
    return v_pos, h_pos, v_neg, h_neg
end

function free_energy(rbm::RBM, vis::Mat{Float64})
    vb = sum(vis .* rbm.visible.bias, 1)
    Wx_b_log = sum(log(1 + exp(rbm.W * vis .+ rbm.hidden.bias)), 1)
    return - vb - Wx_b_log
end

# This function is very memory intensive. Force subsampling of the input
# space with vectors should allow people to trade of itereations of memory consumption
function score_samples(rbm::RBM, X::Mat{Float64}, scoring_ratio)
    # Randomly select the scoring sample view
    sample_size = int(size(X, 2) * scoring_ratio)
    start = rand(1:size(X,2)-sample_size)
    vis = view(X, :, start:start+sample_size)

    # Still have to do a copy
    n_feat, n_samples = size(vis)
    vis_corrupted = copy(vis)
    idxs = rand(1:n_feat, n_samples)
    for (i, j) in zip(idxs, 1:n_samples)
        vis_corrupted[i, j] = 1 - vis_corrupted[i, j]
    end

    # These are expensive, but necessary for now.
    fe = free_energy(rbm, vis)
    fe_corrupted = free_energy(rbm, vis_corrupted)

    return n_feat * log(logistic(fe_corrupted - fe))
end

function update_weights!(rbm::RBM, h_pos, v_pos, h_neg, v_neg, lr, buf)
    dW = buf
    # dW = (h_pos * v_pos') - (h_neg * v_neg')
    gemm!('N', 'T', 1.0, h_neg, v_neg, 0.0, dW)
    gemm!('N', 'T', 1.0, h_pos, v_pos, -1.0, dW)
    # rbm.W += lr * dW
    axpy!(lr, dW, rbm.W)
    # rbm.W += rbm.momentum * rbm.dW_prev
    axpy!(lr * rbm.momentum, rbm.dW_prev, rbm.W)
    # save current dW
    copy!(rbm.dW_prev, dW)
end

function contdiv(rbm::RBM, vis::Mat{Float64}, n_gibbs::Int)
    v_pos, h_pos, v_neg, h_neg = gibbs(rbm, vis, n_times=n_gibbs)
    return v_pos, h_pos, v_neg, h_neg
end

function persistent_contdiv(rbm::RBM, vis::Mat{Float64}, n_gibbs::Int)
    if size(rbm.persistent_chain) != size(vis)
        # persistent_chain not initialized or batch size changed, re-initialize
        rbm.persistent_chain = vis
    end
    # take positive samples from real data
    v_pos, h_pos, _, _ = gibbs(rbm, vis)
    # take negative samples from "fantasy particles"
    rbm.persistent_chain, _, v_neg, h_neg = gibbs(rbm, vis, n_times=n_gibbs)
    return v_pos, h_pos, v_neg, h_neg
end

function fit_batch!(rbm::RBM, vis::Mat{Float64};
                    persistent=true, buf=None, lr=0.1, n_gibbs=1)
    buf = buf == None ? zeros(size(rbm.W)) : buf
    # v_pos, h_pos, v_neg, h_neg = gibbs(rbm, vis, n_times=n_gibbs)
    sampler = persistent ? persistent_contdiv : contdiv
    v_pos, h_pos, v_neg, h_neg = sampler(rbm, vis, n_gibbs)
    lr = lr / size(v_pos, 1)
    update_weights!(rbm, h_pos, v_pos, h_neg, v_neg, lr, buf)
    rbm.hidden.bias += vec(lr * (sum(h_pos, 2) - sum(h_neg, 2)))
    rbm.visible.bias += vec(lr * (sum(v_pos, 2) - sum(v_neg, 2)))
    return rbm
end

function fit(rbm::RBM, X::Mat{Float64};
             persistent=true, lr=0.1, n_iter=10, batch_size=100, n_gibbs=1, scoring_ratio=0.1)
    @assert minimum(X) >= 0 && maximum(X) <= 1
    n_samples = size(X, 2)
    n_batches = int(ceil(n_samples / batch_size))
    w_buf = zeros(size(rbm.W))
    batch = Array(Float64, size(X, 1), batch_size)  # pre-allocate batch
    for itr=1:n_iter
        tic()
        for i=1:n_batches
            # println("fitting $(i)th batch")
            batch = X[:, ((i-1)*batch_size + 1):min(i*batch_size, size(X, 2))]
            batch = full(batch)
            fit_batch!(rbm, batch, persistent=persistent,
                       buf=w_buf, n_gibbs=n_gibbs)
        end
        toc()
        pseudo_likelihood = mean(score_samples(rbm, X, scoring_ratio))
        info("Iteration #$itr, pseudo-likelihood = $pseudo_likelihood")
    end
    return rbm
end

function transform(rbm::RBM, X::Mat{Float64})
    return mean_hiddens(rbm, X)
end

function generate(rbm::RBM, vis::Vec{Float64}; n_gibbs=1)
    return gibbs(rbm, reshape(vis, length(vis), 1); n_times=n_gibbs)[3]
end

function generate(rbm::RBM, X::Mat{Float64}; n_gibbs=1)
    return gibbs(rbm, X; n_times=n_gibbs)[3]
end

function features(rbm::RBM; transpose=true)
    return if transpose rbm.W' else rbm.W end
end
# synonym for backward compatibility
components(rbm::RBM; transpose=true) = features(rbm, transpose)


text:
using VLKeggSDK
using DataFrames

# specify a list of genes (notice this has the KEGG organism code prepended) -
gene_array = [
    "eco:b2388"
]

# DNA sequence for each gene in the list of genes -> piped to a DataFrame
# DNA sequence is a BioSequences.LongDNASeq type -
gene_table = get_sequence_for_gene(gene_array) |> check |> DataFrame

# Protein sequence for each gene in the list of genes -> piped to a DataFrame
# Protein sequence is a BioSequences.LongAminoAcidSeq  -
protein_table = get_sequence_for_protein(gene_array) |> check |> DataFrame