text:
export BrickIMUV2Acceleration
struct BrickIMUV2Acceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2MagneticField
struct BrickIMUV2MagneticField
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2AngularVelocity
struct BrickIMUV2AngularVelocity
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2Orientation
struct BrickIMUV2Orientation
    heading::Integer
    roll::Integer
    pitch::Integer
end

export BrickIMUV2LinearAcceleration
struct BrickIMUV2LinearAcceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2GravityVector
struct BrickIMUV2GravityVector
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2Quaternion
struct BrickIMUV2Quaternion
    w::Integer
    x::Integer
    y::Integer
    z::Integer
end

export BrickIMUV2AllData
struct BrickIMUV2AllData
    acceleration::Vector{Integer}
    magnetic_field::Vector{Integer}
    angular_velocity::Vector{Integer}
    euler_angle::Vector{Integer}
    quaternion::Vector{Integer}
    linear_acceleration::Vector{Integer}
    gravity_vector::Vector{Integer}
    temperature::Integer
    calibration_status::Integer
end

export BrickIMUV2SensorConfiguration
struct BrickIMUV2SensorConfiguration
    magnetometer_rate::Integer
    gyroscope_range::Integer
    gyroscope_bandwidth::Integer
    accelerometer_range::Integer
    accelerometer_bandwidth::Integer
end

export BrickIMUV2SPITFPBaudrateConfig
struct BrickIMUV2SPITFPBaudrateConfig
    enable_dynamic_baudrate::Bool
    minimum_dynamic_baudrate::Integer
end

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

export BrickIMUV2Protocol1BrickletName
struct BrickIMUV2Protocol1BrickletName
    protocol_version::Integer
    firmware_version::Vector{Integer}
    name::String
end

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

export BrickIMUV2
"""
Full fledged AHRS with 9 degrees of freedom
"""
mutable struct BrickIMUV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the calibrated acceleration from the accelerometer for the
x, y and z axis. The acceleration is in the range configured with
:func:`Set Sensor Configuration`.

If you want to get the acceleration periodically, it is recommended
to use the :cb:`Acceleration` callback and set the period with
:func:`Set Acceleration Period`.
"""
function get_acceleration(device::BrickIMUV2)
    return device.deviceInternal.get_acceleration()
end

export get_magnetic_field
"""
    $(SIGNATURES)

Returns the calibrated magnetic field from the magnetometer for the
x, y and z axis.

If you want to get the magnetic field periodically, it is recommended
to use the :cb:`Magnetic Field` callback and set the period with
:func:`Set Magnetic Field Period`.
"""
function get_magnetic_field(device::BrickIMUV2)
    return device.deviceInternal.get_magnetic_field()
end

export get_angular_velocity
"""
    $(SIGNATURES)

Returns the calibrated angular velocity from the gyroscope for the
x, y and z axis. The angular velocity is in the range configured with
:func:`Set Sensor Configuration`.

If you want to get the angular velocity periodically, it is recommended
to use the :cb:`Angular Velocity` acallback nd set the period with
:func:`Set Angular Velocity Period`.
"""
function get_angular_velocity(device::BrickIMUV2)
    return device.deviceInternal.get_angular_velocity()
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the IMU Brick.
The temperature is measured in the core of the BNO055 IC, it is not the
ambient temperature
"""
function get_temperature(device::BrickIMUV2)
    return device.deviceInternal.get_temperature()
end

export get_orientation
"""
    $(SIGNATURES)

Returns the current orientation (heading, roll, pitch) of the IMU Brick as
independent Euler angles. Note that Euler angles always
experience a `gimbal lock <https://en.wikipedia.org/wiki/Gimbal_lock>`__.
We recommend that you use quaternions instead, if you need the absolute
orientation.

If you want to get the orientation periodically, it is recommended
to use the :cb:`Orientation` callback and set the period with
:func:`Set Orientation Period`.
"""
function get_orientation(device::BrickIMUV2)
    return device.deviceInternal.get_orientation()
end

export get_linear_acceleration
"""
    $(SIGNATURES)

Returns the linear acceleration of the IMU Brick for the
x, y and z axis. The acceleration is in the range configured with
:func:`Set Sensor Configuration`.

The linear acceleration is the acceleration in each of the three
axis of the IMU Brick with the influences of gravity removed.

It is also possible to get the gravity vector with the influence of linear
acceleration removed, see :func:`Get Gravity Vector`.

If you want to get the linear acceleration periodically, it is recommended
to use the :cb:`Linear Acceleration` callback and set the period with
:func:`Set Linear Acceleration Period`.
"""
function get_linear_acceleration(device::BrickIMUV2)
    return device.deviceInternal.get_linear_acceleration()
end

export get_gravity_vector
"""
    $(SIGNATURES)

Returns the current gravity vector of the IMU Brick for the
x, y and z axis.

The gravity vector is the acceleration that occurs due to gravity.
Influences of additional linear acceleration are removed.

It is also possible to get the linear acceleration with the influence
of gravity removed, see :func:`Get Linear Acceleration`.

If you want to get the gravity vector periodically, it is recommended
to use the :cb:`Gravity Vector` callback and set the period with
:func:`Set Gravity Vector Period`.
"""
function get_gravity_vector(device::BrickIMUV2)
    return device.deviceInternal.get_gravity_vector()
end

export get_quaternion
"""
    $(SIGNATURES)

Returns the current orientation (w, x, y, z) of the IMU Brick as
`quaternions <https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation>`__.

You have to divide the return values by 16383 (14 bit) to get
the usual range of -1.0 to +1.0 for quaternions.

If you want to get the quaternions periodically, it is recommended
to use the :cb:`Quaternion` callback and set the period with
:func:`Set Quaternion Period`.
"""
function get_quaternion(device::BrickIMUV2)
    return device.deviceInternal.get_quaternion()
end

export get_all_data
"""
    $(SIGNATURES)

Return all of the available data of the IMU Brick.

* acceleration (see :func:`Get Acceleration`)
* magnetic field (see :func:`Get Magnetic Field`)
* angular velocity (see :func:`Get Angular Velocity`)
* Euler angles (see :func:`Get Orientation`)
* quaternion (see :func:`Get Quaternion`)
* linear acceleration (see :func:`Get Linear Acceleration`)
* gravity vector (see :func:`Get Gravity Vector`)
* temperature (see :func:`Get Temperature`)
* calibration status (see below)

The calibration status consists of four pairs of two bits. Each pair
of bits represents the status of the current calibration.

* bit 0-1: Magnetometer
* bit 2-3: Accelerometer
* bit 4-5: Gyroscope
* bit 6-7: System

A value of 0 means for "not calibrated" and a value of 3 means
"fully calibrated". In your program you should always be able to
ignore the calibration status, it is used by the calibration
window of the Brick Viewer and it can be ignored after the first
calibration. See the documentation in the calibration window for
more information regarding the calibration of the IMU Brick.

If you want to get the data periodically, it is recommended
to use the :cb:`All Data` callback and set the period with
:func:`Set All Data Period`.
"""
function get_all_data(device::BrickIMUV2)
    return device.deviceInternal.get_all_data()
end

export leds_on
"""
    $(SIGNATURES)

Turns the orientation and direction LEDs of the IMU Brick on.
"""
function leds_on(device::BrickIMUV2)
    device.deviceInternal.leds_on()
end

export leds_off
"""
    $(SIGNATURES)

Turns the orientation and direction LEDs of the IMU Brick off.
"""
function leds_off(device::BrickIMUV2)
    device.deviceInternal.leds_off()
end

export are_leds_on
"""
    $(SIGNATURES)

Returns *true* if the orientation and direction LEDs of the IMU Brick
are on, *false* otherwise.
"""
function are_leds_on(device::BrickIMUV2)
    return device.deviceInternal.are_leds_on()
end

export save_calibration
"""
    $(SIGNATURES)

A call of this function saves the current calibration to be used
as a starting point for the next restart of continuous calibration
of the IMU Brick.

A return value of *true* means that the calibration could be used and
*false* means that it could not be used (this happens if the calibration
status is not "fully calibrated").

This function is used by the calibration window of the Brick Viewer, you
should not need to call it in your program.
"""
function save_calibration(device::BrickIMUV2)
    return device.deviceInternal.save_calibration()
end

export set_acceleration_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Acceleration` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_acceleration_period(device::BrickIMUV2, period)
    device.deviceInternal.set_acceleration_period(period)
end

export get_acceleration_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Acceleration Period`.
"""
function get_acceleration_period(device::BrickIMUV2)
    return device.deviceInternal.get_acceleration_period()
end

export set_magnetic_field_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Magnetic Field` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_magnetic_field_period(device::BrickIMUV2, period)
    device.deviceInternal.set_magnetic_field_period(period)
end

export get_magnetic_field_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Magnetic Field Period`.
"""
function get_magnetic_field_period(device::BrickIMUV2)
    return device.deviceInternal.get_magnetic_field_period()
end

export set_angular_velocity_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Angular Velocity` callback is
triggered periodically. A value of 0 turns the callback off.
"""
function set_angular_velocity_period(device::BrickIMUV2, period)
    device.deviceInternal.set_angular_velocity_period(period)
end

export get_angular_velocity_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Angular Velocity Period`.
"""
function get_angular_velocity_period(device::BrickIMUV2)
    return device.deviceInternal.get_angular_velocity_period()
end

export set_temperature_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_temperature_period(device::BrickIMUV2, period)
    device.deviceInternal.set_temperature_period(period)
end

export get_temperature_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Temperature Period`.
"""
function get_temperature_period(device::BrickIMUV2)
    return device.deviceInternal.get_temperature_period()
end

export set_orientation_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Orientation` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_orientation_period(device::BrickIMUV2, period)
    device.deviceInternal.set_orientation_period(period)
end

export get_orientation_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Orientation Period`.
"""
function get_orientation_period(device::BrickIMUV2)
    return device.deviceInternal.get_orientation_period()
end

export set_linear_acceleration_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Linear Acceleration` callback is
triggered periodically. A value of 0 turns the callback off.
"""
function set_linear_acceleration_period(device::BrickIMUV2, period)
    device.deviceInternal.set_linear_acceleration_period(period)
end

export get_linear_acceleration_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Linear Acceleration Period`.
"""
function get_linear_acceleration_period(device::BrickIMUV2)
    return device.deviceInternal.get_linear_acceleration_period()
end

export set_gravity_vector_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Gravity Vector` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_gravity_vector_period(device::BrickIMUV2, period)
    device.deviceInternal.set_gravity_vector_period(period)
end

export get_gravity_vector_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Gravity Vector Period`.
"""
function get_gravity_vector_period(device::BrickIMUV2)
    return device.deviceInternal.get_gravity_vector_period()
end

export set_quaternion_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Quaternion` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_quaternion_period(device::BrickIMUV2, period)
    device.deviceInternal.set_quaternion_period(period)
end

export get_quaternion_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Quaternion Period`.
"""
function get_quaternion_period(device::BrickIMUV2)
    return device.deviceInternal.get_quaternion_period()
end

export set_all_data_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`All Data` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_all_data_period(device::BrickIMUV2, period)
    device.deviceInternal.set_all_data_period(period)
end

export get_all_data_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set All Data Period`.
"""
function get_all_data_period(device::BrickIMUV2)
    return device.deviceInternal.get_all_data_period()
end

export set_sensor_configuration
"""
    $(SIGNATURES)

Sets the available sensor configuration for the Magnetometer, Gyroscope and
Accelerometer. The Accelerometer Range is user selectable in all fusion modes,
all other configurations are auto-controlled in fusion mode.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function set_sensor_configuration(device::BrickIMUV2, magnetometer_rate, gyroscope_range, gyroscope_bandwidth, accelerometer_range, accelerometer_bandwidth)
    device.deviceInternal.set_sensor_configuration(magnetometer_rate, gyroscope_range, gyroscope_bandwidth, accelerometer_range, accelerometer_bandwidth)
end

export get_sensor_configuration
"""
    $(SIGNATURES)

Returns the sensor configuration as set by :func:`Set Sensor Configuration`.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function get_sensor_configuration(device::BrickIMUV2)
    return device.deviceInternal.get_sensor_configuration()
end

export set_sensor_fusion_mode
"""
    $(SIGNATURES)

If the fusion mode is turned off, the functions :func:`Get Acceleration`,
:func:`Get Magnetic Field` and :func:`Get Angular Velocity` return uncalibrated
and uncompensated sensor data. All other sensor data getters return no data.

Since firmware version 2.0.6 you can also use a fusion mode without magnetometer.
In this mode the calculated orientation is relative (with magnetometer it is
absolute with respect to the earth). However, the calculation can't be influenced
by spurious magnetic fields.

Since firmware version 2.0.13 you can also use a fusion mode without fast
magnetometer calibration. This mode is the same as the normal fusion mode,
but the fast magnetometer calibration is turned off. So to find the orientation
the first time will likely take longer, but small magnetic influences might
not affect the automatic calibration as much.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function set_sensor_fusion_mode(device::BrickIMUV2, mode)
    device.deviceInternal.set_sensor_fusion_mode(mode)
end

export get_sensor_fusion_mode
"""
    $(SIGNATURES)

Returns the sensor fusion mode as set by :func:`Set Sensor Fusion Mode`.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function get_sensor_fusion_mode(device::BrickIMUV2)
    return device.deviceInternal.get_sensor_fusion_mode()
end

export set_spitfp_baudrate_config
"""
    $(SIGNATURES)

The SPITF protocol can be used with a dynamic baudrate. If the dynamic baudrate is
enabled, the Brick will try to adapt the baudrate for the communication
between Bricks and Bricklets according to the amount of data that is transferred.

The baudrate will be increased exponentially if lots of data is sent/received and
decreased linearly if little data is sent/received.

This lowers the baudrate in applications where little data is transferred (e.g.
a weather station) and increases the robustness. If there is lots of data to transfer
(e.g. Thermal Imaging Bricklet) it automatically increases the baudrate as needed.

In cases where some data has to transferred as fast as possible every few seconds
(e.g. RS485 Bricklet with a high baudrate but small payload) you may want to turn
the dynamic baudrate off to get the highest possible performance.

The maximum value of the baudrate can be set per port with the function
:func:`Set SPITFP Baudrate`. If the dynamic baudrate is disabled, the baudrate
as set by :func:`Set SPITFP Baudrate` will be used statically.

.. versionadded:: 2.0.10\$nbsp;(Firmware)
"""
function set_spitfp_baudrate_config(device::BrickIMUV2, enable_dynamic_baudrate, minimum_dynamic_baudrate)
    device.deviceInternal.set_spitfp_baudrate_config(enable_dynamic_baudrate, minimum_dynamic_baudrate)
end

export get_spitfp_baudrate_config
"""
    $(SIGNATURES)

Returns the baudrate config, see :func:`Set SPITFP Baudrate Config`.

.. versionadded:: 2.0.10\$nbsp;(Firmware)
"""
function get_spitfp_baudrate_config(device::BrickIMUV2)
    return device.deviceInternal.get_spitfp_baudrate_config()
end

export get_send_timeout_count
"""
    $(SIGNATURES)

Returns the timeout count for the different communication methods.

The methods 0-2 are available for all Bricks, 3-7 only for Master Bricks.

This function is mostly used for debugging during development, in normal operation
the counters should nearly always stay at 0.

.. versionadded:: 2.0.7\$nbsp;(Firmware)
"""
function get_send_timeout_count(device::BrickIMUV2, communication_method)
    return device.deviceInternal.get_send_timeout_count(communication_method)
end

export set_spitfp_baudrate
"""
    $(SIGNATURES)

Sets the baudrate for a specific Bricklet port.

If you want to increase the throughput of Bricklets you can increase
the baudrate. If you get a high error count because of high
interference (see :func:`Get SPITFP Error Count`) you can decrease the
baudrate.

If the dynamic baudrate feature is enabled, the baudrate set by this
function corresponds to the maximum baudrate (see :func:`Set SPITFP Baudrate Config`).

Regulatory testing is done with the default baudrate. If CE compatibility
or similar is necessary in your applications we recommend to not change
the baudrate.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function set_spitfp_baudrate(device::BrickIMUV2, bricklet_port, baudrate)
    device.deviceInternal.set_spitfp_baudrate(bricklet_port, baudrate)
end

export get_spitfp_baudrate
"""
    $(SIGNATURES)

Returns the baudrate for a given Bricklet port, see :func:`Set SPITFP Baudrate`.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function get_spitfp_baudrate(device::BrickIMUV2, bricklet_port)
    return device.deviceInternal.get_spitfp_baudrate(bricklet_port)
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 Brick side. All
Bricklets have a similar function that returns the errors on the Bricklet side.

.. versionadded:: 2.0.5\$nbsp;(Firmware)
"""
function get_spitfp_error_count(device::BrickIMUV2, bricklet_port)
    return device.deviceInternal.get_spitfp_error_count(bricklet_port)
end

export enable_status_led
"""
    $(SIGNATURES)

Enables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.
"""
function enable_status_led(device::BrickIMUV2)
    device.deviceInternal.enable_status_led()
end

export disable_status_led
"""
    $(SIGNATURES)

Disables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.
"""
function disable_status_led(device::BrickIMUV2)
    device.deviceInternal.disable_status_led()
end

export is_status_led_enabled
"""
    $(SIGNATURES)

Returns *true* if the status LED is enabled, *false* otherwise.
"""
function is_status_led_enabled(device::BrickIMUV2)
    return device.deviceInternal.is_status_led_enabled()
end

export get_protocol1_bricklet_name
"""
    $(SIGNATURES)

Returns the firmware and protocol version and the name of the Bricklet for a
given port.

This functions sole purpose is to allow automatic flashing of v1.x.y Bricklet
plugins.
"""
function get_protocol1_bricklet_name(device::BrickIMUV2, port)
    return device.deviceInternal.get_protocol1_bricklet_name(port)
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 an
accuracy of ±15%. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickIMUV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Brick. Calling this function
on a Brick inside of a stack will reset the whole stack.

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

export write_bricklet_plugin
"""
    $(SIGNATURES)

Writes 32 bytes of firmware to the bricklet attached at the given port.
The bytes are written to the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_bricklet_plugin(device::BrickIMUV2, port, offset, chunk)
    device.deviceInternal.write_bricklet_plugin(port, offset, chunk)
end

export read_bricklet_plugin
"""
    $(SIGNATURES)

Reads 32 bytes of firmware from the bricklet attached at the given port.
The bytes are read starting at the position offset * 32.

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

export get_identity
"""
    $(SIGNATURES)

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

The position is the position in the stack from '0' (bottom) to '8' (top).

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

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


text:
# BBox

export BBox
"""
Bounding box for a region of integer points. Lower and upper bounds
are inclusive.
"""
struct BBox{D}
    min::SVector{D,Int}
    max::SVector{D,Int}
    BBox{D}() where {D} = new{D}(zero(SVector{D,Int}), -one(SVector{D,Int}))
    function BBox{D}(min::SVector{D,Int}, max::SVector{D,Int}) where {D}
        return new{D}(min, max)
    end
end

BBox(min::SVector{D,Int}, max::SVector{D,Int}) where {D} = BBox{D}(min, max)

Base.min(bbox::BBox) = bbox.min
Base.max(bbox::BBox) = bbox.max

Base.ndims(bbox::BBox{D}) where {D} = D

Base.eltype(bbox::BBox) = Int

Base.isempty(bbox::BBox) = any(max(bbox) > min(bbox))

function Base.size(bbox::BBox{D}) where {D}
    isempty(bbox.max) && return zero(SVector{D,Int})
    return max(bbox) - min(bbox) .+ 1
end

Base.length(bbox::BBox) = prod(size(bbox))

function Base.strides(bbox::BBox{D}) where {D}
    res = zero(SVector{D,Int})
    str = 1
    for d in 1:D
        res = setindex!(res, str, d)
        str *= size(bbox, d)
    end
    return res
end

export linear_index
function linear_index(bbox::BBox{D}, pos::SVector{D,Int}) where {D}
    return sum((pos - min(bbox)) .* strides(bbox))
end

export shift
function shift(bbox::BBox{D}, offset::SVector{D,Int}) where {D}
    isempty(bbox) && return BBox{D}()
    return BBox{D}(min(bbox) + offset, max(bbox) + offset)
end

export extend
function extend(bbox::BBox{D}, lo::SVector{D,Int}, hi::SVector{D,Int}) where {D}
    isempty(bbox) && return BBox{D}()
    return BBox{D}(min(bbox) - lo, max(bbox) + hi)
end

function Base.in(pos::SVector{D,Int}, bbox::BBox{D}) where {D}
    isempty(bbox) && return false
    return all(min(bbox) .≤ pos .≤ max(bbox))
end

function Base.intersect(bbox1::BBox{D}, bbox2::BBox{D}) where {D}
    (isempty(bbox1) || isempty(bbox2)) && return BBox{D}()
    return BBox{D}(max(min(bbox1), min(bbox2)), min(max(bbox1), max(bbox2)))
end

export intersects
function intersects(bbox1::BBox{D}, bbox2::BBox{D}) where {D}
    return !isempty(intersect(bbox1, bbox2))
end

ext: jl


text:
########
# TESTS
########
# To be completed

#################################
# PATHS
# change "path_main" if necessary
path_main = "/home/julien/Final-Project-Julien-Pascal"
path_table = string(path_main,"/tables/")  #the path tables
path_surpluses = string(path_main,"/surpluses/")#path to precalculated objects, to speed up the process

file_random_numbers =  jldopen(string(path_surpluses,"random_numbers.jld"), "r")
random_numbers = read(file_random_numbers, "random_numbers")

include(string(path_main,"/src/objective_function.jl"))
include(string(path_main,"/src/objective_function_2.jl"))

number_parameters = 6 #number of parameters
# Read the estimated paramaters
if isfile(string(path_table, "estimated_parameters.csv"))
	println("loading estimated parameters")
	coeff = readtable(string(path_table, "estimated_parameters.csv"))
else #if does not exist, use the parameters from the original paper
	parameters = [0.77, 0.023, 0.99, 2.00, 5.56, 0.73]
end

z0 = coeff[:Value][1];
sigma = coeff[:Value][2];
lambda0 = coeff[:Value][3];
eta =  coeff[:Value][4];
mu =  coeff[:Value][5];
x_lower_bar= coeff[:Value][6] #lower bound of the ability level

parameters = [z0, sigma, lambda0, eta, mu, x_lower_bar] #parameters to be estimated

b0 = [1.0; 0.023;  0.059;  0.210; 0.761; 0.085; 0.046] #moments to match
number_moments = 7

Weighting_matrix = eye(number_moments,number_moments)./(b0.*b0)

# 1. Make sure the objective functions used by NLopt and DE coincide:
context("Test the objective functions") do

	facts("Test objective = objective 2") do

	result_objective = objective_function(parameters,parameters, Weighting_matrix, random_numbers, b0)
	result_objective2 = objective_function_2(parameters, Weighting_matrix, random_numbers, b0)
	difference = result_objective - result_objective2

		#the difference should be close to zero
		@fact difference --> roughly(0; atol=0.001)

	end
end

# Look if the estimated parameters perform well compared the original ones:
context("Test estimated parameters vs origianl paper's parameters") do

	facts("estimated parameters vs original paper's parameters") do

	parameters_original = [0.77, 0.023, 0.99, 2.00, 5.56, 0.73]

	result_original = objective_function(parameters_original ,parameters_original , Weighting_matrix, random_numbers, b0)

	result_estimated = objective_function_2(parameters, Weighting_matrix, random_numbers, b0)

	difference = result_estimated - result_original

		# if difference is negative, then the estimated parameters are better wrt to the criteria used:
		@fact difference --> less_than_or_equal(0)

	end
end

# 2. With the estimated coefficient, make sure the system behave properly
# e.g. the unemployment rate cannot be negative

# 3. Make sure the wages are finite

# 4. check the accuracy of the minimum of the objective function
# Use a local mininmization algorithm (example: COBYLA) with the estimated parameters
# as an initial condition.
# A. If the local minimization find a better point, then the global
# maximum is not one. Use this new local minimum instead.
# B. Else, keep the global max


text:
using LinearAlgebra
using DynamicPolynomials
using SwitchOnSafety
using Combinatorics
using SparseArrays
using JuMP, Ipopt, MosekTools,NLopt
using SpecialFunctions
using ControlSystems

include("../src/RandomTrajectories.jl")
include("../src/AlgebraicLift.jl")
include("../src/ScenarioOpti.jl")
include("../src/ProbabilisticCertificates.jl")
include("../src/WhiteBoxAnalysis.jl")

dim = 3; numMode = 4; dimIn = 3

numScen_budget = 30000

c = 1375
R12 = 1.5
Ro12 = 3
Ro3 = 2.7
tau = 3*60
A = []
for R13 in [0.8 1.2]
    for R23 in [0.8 1.2]
        Ai = [1-tau/c*(1/R12+1/R13+1/Ro12) tau/c/R12 tau/c/R13;tau/c/R12 1-tau/c*(1/R12+1/R23+1/Ro12) tau/c/R23;tau/c/R13 tau/c/R23 1-tau/c*(1/R13+1/R23+1/Ro3)]
        push!(A,Ai)
    end
end

B = tau/c*Matrix(I,dimIn,dimIn)

Q = Matrix(I,dim,dim)
R = 0.02*Matrix(I,dimIn,dimIn)

(Kw,Pw) = white_box_LQR(A,B,Q,R)
println("K white-box: $Kw")

(state0_budget,state_budget) = generate_trajectories(1,A,numScen_budget)

K0 = zeros(dimIn,dim)

#=
K,xi,P,ratio = probabilistc_LQR_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,Q=Q,R=R,kappaU=1e2,K0=K0,beta=0.01,tol=1e-4)

state = [15;10;12]
X = zeros(dim,10)
for i in 1:10
    global state
    X[:,i] = state
    b = rand((1:numMode),1)
    u = K*state
    state = A[b[1]]*state+B*u
end

X = X.+24

using Plots
gr(size = (400, 300))
fn = plot(0:3:9*3, Any[X[1,:],X[2,:],X[3,:]],label = ["Zone 1" "Zone 2" "Zone 3"], line = [:solid :solid :solid], lw = 2)
xlabel!("Time (minute)")
ylabel!("Temperature (°C) ")

=#
xisq = []
ratiosq = []

for N in 1000:1000:numScen_budget
    K,xi,P,ratio = probabilistc_LQR_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,Q=Q,R=R,kappaU=1e2,K0=K0,beta=0.01,tol=1e-4)
    push!(xisq,xi)
    push!(ratiosq,ratio)
    println("xi sequence: $xisq")
    println("ratio sequence: $ratiosq")
    println("K: $K")
    println("K white-box: $Kw")
    println("P: $P")
    println("P white-box: $Pw")
end

#=
using Plots
using LaTeXStrings
gr(size = (400, 300))
fn = plot(500:500:12000, Any[feas,fill(1.0,24)],label = "", line = [:solid :dashdot], lw = 2)
xlabel!(L"N")
ylabel!(L"\bar{\xi}(\omega_N)")
savefig(fn,"fn.png")
=#

#feas = xisq./ratiosq
#[1.1462996266804206, 1.0304548488231406, 1.1839711666921258, 0.9796333488331512, 0.9668010693346112, 1.1983116266741682, 1.0747101169748587, 0.953405955179486, 0.9531081705886429, 0.9502858772790718]
#K: [-0.028447388063758478 -0.02397013534333861 -0.024868877374478904; -0.022049588263423403 -0.029909164375603235 -0.025289972953130273; -0.02557500118128832 -0.02521458873586027 -0.025649666991379057]

#xisq =[0.9506955125447197, 0.9500814227738398, 0.9496688890425853, 0.9506770820573243, 0.9491450126708667, 0.9491309513543958, 0.9512678391859847, 0.9511263256916642, 0.9509734949035534, 0.9496753424296864, 0.9495745388568111, 0.9517758890790594, 0.9517405648377967, 0.9516439742658995, 0.9498230924269193, 0.9500617518717405, 0.9499383280065629, 0.9498282336931118, 0.9516072138900373, 0.9518346791078194, 0.9518833690400422, 0.9518458883283956, 0.95181760735376, 0.9520342911681299]
#ratiosq = [0.32929113429561285, 0.6359092714236014, 0.7458849122101749, 0.8080460152288503, 0.8388407832548765, 0.8651945928045293, 0.8807624981871738, 0.895048424191853, 0.9057136189419326, 0.9146100154266903, 0.9218573770210922, 0.9277461652693918, 0.9328689844736177, 0.9374298867002537, 0.9413125864511712, 0.9503108804475847, 0.9529825663527286, 0.9554076570766022, 0.9575460410445409, 0.9595593129229043, 0.9614069403433664, 0.963069357269896, 0.9646071007278422, 0.965914139694424]

# K = [-3.379476201830191 -0.5565284859098218 -0.6653948967484442; -0.5610849053980036 -3.374365413285184 -0.665561226645545; -0.664814150076262 -0.6724159301492334 -3.2436777438358977]
# P = [1.4041357466418602 0.11383342685787007 0.13332235372007795; 0.11383342685787007 1.4041274464054574 0.13332286522268716; 0.13332235372007795 0.13332286522268716 1.3787232916679992]

#P white-box: [1.384392543034221 0.10849589561899299 0.1270429070953619; 0.10849589561899299 1.3843925454395583 0.1270429068809821; 0.1270429070953619 0.1270429068809821 1.3602060587775264]
#K white-box: [-3.173628196889124 -0.4840291283071396 -0.5938068939821004; -0.48402965221411914 -3.1736302770962785 -0.5938067500136359; -0.5881971655090553 -0.588195839338544 -3.031968422173231]

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
_has_aarch64_sve() = (Base.libllvm_version ≥ v"11") && (Base.BinaryPlatforms.CPUID.test_cpu_feature(Base.BinaryPlatforms.CPUID.JL_AArch64_sve))

if Int === Int64
    @noinline vscale() = ccall("llvm.vscale.i64", llvmcall, Int64, ())
else
    @noinline vscale() = ccall("llvm.vscale.i32", llvmcall, Int32, ())
end

# TODO: find actually support SVE
# _dynamic_register_size() = _has_aarch64_sve() ? 16vscale() : 16
_dynamic_register_size() = 16

function _set_sve_vector_width!(bytes = _dynamic_register_size())
    @eval begin
        register_size() = StaticInt{$bytes}()
        simd_integer_register_size() = StaticInt{$bytes}()
    end
    nothing
end

if _has_aarch64_sve()# && !(Bool(has_feature(Val(:aarch64_sve))))
    has_feature(::Val{:aarch64_sve_cpuid}) = True()
    _set_sve_vector_width!()
else
    # has_feature(::Val{:aarch64_svejl}) = False()
    register_size() = StaticInt{16}()
    simd_integer_register_size() = StaticInt{16}()
end

function reset_extra_features!()
    drs = _dynamic_register_size()
    register_size() ≠ drs && _set_sve_vector_width!(drs)
    hassve = _has_aarch64_sve()
    if hassve ≠ has_feature(Val(:aarch64_sve_cpuid))
        @eval has_feature(::Val{:aarch64_sve_cpuid}) = $(Expr(:call, hassve ? :True : :False))
    end
end

fma_fast() = True()
register_count() = StaticInt{32}()
has_opmask_registers() = has_feature(Val(:aarch64_sve_cpuid))

fast_int64_to_double() = True()

fast_half() = False()


text:
#=
CMUTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/CMU.jl")
matrix =  GBMatrix(sparse([[0, 1, 1, 0, 0] [1, 0, 1, 1, 1] [1, 1, 0, 0, 1] [0, 1, 0, 0, 1] [0, 1, 1, 1, 0]]))
@test cmu(matrix, 5) == 3


text:
abstract type MTTranslationModel <: TranslationModel end

struct MTTranslationParam <: EoSParam
    acentricfactor::SingleParam{Float64}
end

@newmodelsimple MTTranslation MTTranslationModel MTTranslationParam

export MTTranslation
function MTTranslation(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = MTTranslationParam(acentricfactor)
    model = MTTranslation(packagedparams, verbose=verbose)
    return model
end

function translation(model::CubicModel,V,T,z,translation_model::MTTranslation)
    Tc = model.params.Tc.values
    Pc = model.params.Pc.values
    ω  = translation_model.params.acentricfactor.values

    Zc = @. 0.289-0.0701*ω-0.0207*ω^2
    β  = @. -10.2447-28.6312*ω

    t0 = @. R̄*Tc/Pc*(-0.014471+0.067498*ω-0.084852*ω^2+0.067298*ω^3-0.017366*ω^4)
    tc = @. R̄*Tc/Pc*(0.3074-Zc)

    Tr = @. T/Tc
    return @. t0+(tc-t0)*exp(β*abs(1-Tr))
end


text:
using Test, DocumentationGenerator
using Pkg

const julia = first(Base.julia_cmd())

@test length(DocumentationGenerator.installable_on_version()) > 1500

@testset "Running code with a timeout" begin
    let
        tempfile = tempname()
        str = """
        sleep(4)
        sleep(4)
        write("$tempfile", "hi")
        sleep(3)
        """
        proc, _ = DocumentationGenerator.run_with_timeout(`$julia -e $str`, timeout=7)
        wait(proc)
        @test !success(proc)
        @test !isfile(tempfile)
    end

    let
        tempfile = tempname()
        str = """
        for i = 1:10
        sleep(1)
        println(i)
        end
        write("$tempfile", "hi")
        """
        proc, _ = DocumentationGenerator.run_with_timeout(`$julia -e $str`, timeout=6)
        wait(proc)
        @test success(proc)
        @test String(read(tempfile)) == "hi"
    end

    let
        tempfile = tempname()
        logfile = tempname()
        str = """
        for i in 1:10
        println(i)
        sleep(1)
        end
        write("$tempfile", "hi")
        """
        proc, task = DocumentationGenerator.run_with_timeout(`$julia -e $str`, timeout=3, log=logfile)
        wait(proc)
        @test success(proc)
        @test String(read(tempfile)) == "hi"
        logstr = ""
        for i in 1:10
            logstr *= "$i"
        end
        wait(task)
        @test isfile(logfile)
        @test replace(String(read(logfile)), '\n' => "") == logstr
    end
end

@testset "Documentation Generation" begin
    packages = [
        # without docs
        (
            name = "FixedPointNumbers",
            uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93",
            url = "https://github.com/JuliaMath/FixedPointNumbers.jl.git",
            versions = [v"0.5.3"],
            installs = [true],
            doctype = ["default"],
        ),
        (
            name = "ReactionNetworkImporters",
            url = "https://github.com/isaacsas/ReactionNetworkImporters.jl.git",
            uuid = "b4db0fb7-de2a-5028-82bf-5021f5cfa881",
            versions = [v"0.1.4"],
            installs = [true],
            doctype = ["default"],
        ),
        # with docs
        (
            name = "Example",
            url = "https://github.com/JuliaLang/Example.jl.git",
            uuid = "7876af07-990d-54b4-ab0e-23690620f79a",
            versions = [v"0.5.1", v"0.5.2"],
            installs = [true, false],
            doctype = ["default", "real"],
        ),
        (
            name = "DynamicHMC",
            url = "https://github.com/tpapp/DynamicHMC.jl.git",
            uuid = "bbc10e6e-7c05-544b-b16e-64fede858acb",
            versions = [v"2.0.1"],
            installs = [true],
            doctype = ["real"],
        ),
        (
            name = "jlpkg",
            url = "https://github.com/fredrikekre/jlpkg",
            uuid = "c4c688b2-6cc8-11e9-1c12-6d20b663313d",
            versions = [v"1.0.2"],
            installs = [true],
            doctype = ["default"],
        ),
        # with fancy docs
        # not installable (wrong julia version)
        (
            name = "Flux",
            url = "https://github.com/FluxML/Flux.jl.git",
            uuid = "587475ba-b771-5e3f-ad9e-33799f191a9c",
            versions = [v"0.2.2", v"0.7.3", v"0.8.3"],
            installs = [false, true, true],
            doctype = [nothing, "real", "real"],
        ),
        # with hosted docs
        (
            name = "Juno",
            url = "https://github.com/JunoLab/Juno.jl.git",
            uuid = "e5e0dc1b-0480-54bc-9374-aad01c23163d",
            versions = [v"0.7.0"],
            installs = [true],
            doctype = ["real"]
        ),
    ]

    basepath = @__DIR__
    rm(joinpath(basepath, "logs"), force = true, recursive = true)
    rm(joinpath(basepath, "build"), force = true, recursive = true)

    DocumentationGenerator.build_documentations(
        packages, basepath = basepath, filter_versions = identity
    )

    build = joinpath(basepath, "build")
    @testset "build folder" begin
        for pkg in packages
            pkgbuild = joinpath(build, DocumentationGenerator.get_docs_dir(pkg.name, pkg.uuid))
            @test isdir(pkgbuild)
            @testset "$(pkg.name): $(version)" for (i, version) in enumerate(pkg.versions)
                @test isfile(basepath, "logs", string(pkg.name, "-", pkg.uuid, " ", version, ".log"))

                versiondir = joinpath(pkgbuild, string(version))
                @test isdir(versiondir)
                toml_path = joinpath(versiondir, "meta.toml")
                @test isfile(toml_path)
                toml = Pkg.TOML.parsefile(toml_path)
                pkginstalls = get(toml, "installs", false)
                @test pkginstalls == pkg.installs[i]
                if pkginstalls
                    doctype = get(toml, "doctype", nothing)
                    @test doctype == pkg.doctype[i]
                    @test isfile(joinpath(versiondir, "index.html"))
                    if doctype == "default"
                        @test isdir(joinpath(versiondir, "autodocs"))
                    end
                end
            end
        end
    end
end

@testset "Ruby Gems" begin
    @test isfile(DocumentationGenerator.find_ruby_gem("licensee"))
    @test isfile(DocumentationGenerator.find_ruby_gem("commonmarker"))
end

@testset "Documentation Registry" begin
    mktempdir() do dir
        reg = DocumentationGenerator.download_registry(dir)
        @test isfile(reg)
        toml = Pkg.TOML.parsefile(reg)
        @test length(keys(toml)) > 0
        @test haskey(toml, "e5e0dc1b-0480-54bc-9374-aad01c23163d")
        junosettings = toml["e5e0dc1b-0480-54bc-9374-aad01c23163d"]
        @test junosettings["method"] == "hosted"
    end
end


text:
module CG
  import Base: +, -, *, getindex, inv;

  export Vec4, Vec4f, Mat4, Mat4f, Object, homo, xproduct, unitize;

  # Basic types
  type Vec4{T<:Number}
    e1::T
    e2::T
    e3::T
    e4::T
  end

  type Mat4{T<:Vec4}
    v1::T
    v2::T
    v3::T
    v4::T
  end

  getindex(v::Vec4, index) =
  begin
    if (index == 1)
      return v.e1;
    elseif (index == 2)
      return v.e2;
    elseif (index == 3)
      return v.e3;
    elseif (index == 4)
      return v.e4;
    elseif (index == :)
      return v;
    end
  end

  getindex(m::Mat4, row, col) =
  begin
    if ((row == :) && (col == :))
      return m;
    end

    r;
    if (row == 1)
      r = m.v1;
    elseif (row == 2)
      r = m.v2;
    elseif (row == 3)
      r = m.v3;
    elseif (row == 4)
      r = m.v4;
    elseif (row == :)
      return Vec4(m.v1[col], m.v2[col], m.v3[col], m.v4[col]);
    end

    return r[col];
  end

  #typedefs
  Vec4f = Vec4{Float32};
  Mat4f = Mat4{Vec4f};

  ##operators
  # 1) vector addition
  + (v1::Vec4, v2::Vec4) =
  begin
    return Vec4f(
      v1.e1 + v2.e1,
      v1.e2 + v2.e2,
      v1.e3 + v2.e3,
      v1.e4 + v2.e4
    );
  end

  -(v1::Vec4) =
  begin
    return Vec4(-v1.e1, -v1.e2, -v1.e3, -v1.e4);
  end

  -(v1::Vec4, v2::Vec4) =
  begin
    return v1 + - v2;
  end

  *(a::Number, v::Vec4) =
  begin
    return Vec4f(a * v.e1, a * v.e2, a * v.e3, a * v.e4);
  end

  # dot product
  *(v1::Vec4, v2::Vec4) =
  begin
    return v1.e1 * v2.e1 + v1.e2 * v2.e2 + v1.e3 * v2.e3;
  end

  homo = function(v::Vec4)
    return 1.0/v[4] * v;
  end

  unitize = function(v::Vec4)
    return 1.0/sqrt(v * v) * v;
  end

  xproduct = function(v1::Vec4, v2::Vec4)
    v1n = homo(v1);
    v2n = homo(v2);
    return Vec4(v1n[2] * v2n[3] - v1n[3] * v2n[2],
                v1n[3] * v2n[1] - v1n[1] * v2n[3],
                v1n[1] * v2n[2] - v1n[2] * v2n[1],
                v1n[4]);
  end

  # 2) matrix operations
  *(m::Mat4, v::Vec4) =
  begin
    return Vec4f(
      m.v1 * v,
      m.v2 * v,
      m.v3 * v,
      m.v4 * v
    );
  end

  *(m1::Mat4f, m2::Mat4f) =
  begin
    return Mat4f(
      Vec4f(m1[1, :] * m2[:, 1], m1[1, :] * m2[:, 2], m1[1, :] * m2[:, 3], m1[1, :] * m2[:, 4]),
      Vec4f(m1[2, :] * m2[:, 1], m1[2, :] * m2[:, 2], m1[2, :] * m2[:, 3], m1[2, :] * m2[:, 4]),
      Vec4f(m1[3, :] * m2[:, 1], m1[3, :] * m2[:, 2], m1[3, :] * m2[:, 3], m1[3, :] * m2[:, 4]),
      Vec4f(m1[4, :] * m2[:, 1], m1[4, :] * m2[:, 2], m1[4, :] * m2[:, 3], m1[4, :] * m2[:, 4])
    )
  end

  inv(m::Mat4f) =
  begin
    mat = zeros(4, 4);
    for r = 1:4
      for c = 1:4
        mat[r, c] = m[r, c];
      end
    end

    matinv = inv(mat);

    return Mat4f(
        Vec4f(matinv[1, 1], matinv[1, 2], matinv[1, 3], matinv[1, 4]),
        Vec4f(matinv[2, 1], matinv[2, 2], matinv[2, 3], matinv[2, 4]),
        Vec4f(matinv[3, 1], matinv[3, 2], matinv[3, 3], matinv[3, 4]),
        Vec4f(matinv[4, 1], matinv[4, 2], matinv[4, 3], matinv[4, 4])
      )
  end

  ## Object
  type Object
    vertices::Vector{Vec4f}
    Object(x::Vector{Vec4f}) = new(x)
    Object(x::Vec4f...) = new(collect(Vec4f, x))
  end

  ## Transformation
  export Transformation;
  type Transformation
    M::Mat4f

    Transformation(v1::Vec4f, v2::Vec4f, v3::Vec4f, v4::Vec4f) = new(Mat4f(v1, v2, v3, v4));
    Transformation(m::Mat4f) = new(m);
  end

  *(T::Transformation, v::Vec4f) =
  begin
    return T.M * v;
  end

  *(T::Transformation, o::Object) =
  begin
    return Object(map(x -> T * x, o.vertices))
  end

  # 1) Translation
  export translation;
  function translation(t::Vec4f)
    return Transformation(
      Vec4f(1, 0, 0, t.e1),
      Vec4f(0, 1, 0, t.e2),
      Vec4f(0, 0, 1, t.e3),
      Vec4f(0, 0, 0, 1)
    )
  end

  # 2) scaling
  export scaling;
  function scaling(t::Vec4f)
    return Transformation(
      Vec4f(t.e1, 0, 0, 0),
      Vec4f(0, t.e2, 0, 0),
      Vec4f(0, 0, t.e3, 0),
      Vec4f(0, 0, 0, 1)
    )
  end

  # 3) rotation
  export rotx;
  function rotx(a)
    return Transformation(
      Vec4f(1, 0, 0, 0),
      Vec4f(0, cos(a), -sin(a), 0),
      Vec4f(0, sin(a), cos(a), 0),
      Vec4f(0, 0, 0, 1)
    )
  end

  export roty;
  function roty(a)
    return Transformation(
      Vec4f(cos(a), 0, sin(a), 0),
      Vec4f(0, 1, 0, 0),
      Vec4f(-sin(a), 0, cos(a), 0),
      Vec4f(0, 0, 0, 1)
    )
  end

  export rotz;
  function rotz(a)
    return Transformation(
      Vec4f(cos(a), -sin(a), 0, 0),
      Vec4f(sin(a), cos(a), 0, 0),
      Vec4f(0, 0, 1, 0),
      Vec4f(0, 0, 0, 1)
    )
  end

  export euler;
  function euler(x, y, z)
    return rotz(z) * roty(y) * rotx(x);
  end

  ## Transformation chaining
  *(t1::Transformation, t2::Transformation) =
  begin
    return Transformation(t1.M * t2.M)
  end

  inv(t::Transformation) =
  begin
    return Transformation(inv(t.M));
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
function getBestMove(M,player,score,moves,movePossible,makeMove!,gameOver)
    # uses minimax algorithm to find best move

    # initialise algorithm values
    bestVal = -1000
    bestMove = -1
    bestDepth = 1000
    moveVal = 0

    for i in moves # search over all candidate moves
        if movePossible(M,i)
            N = deepcopy(M)
            N = makeMove!(N,i,player) # play the move
            moveVal, depth = minimax(N,0,false,3-player,moves,score,movePossible,makeMove!,gameOver) # evaluate it

            print("$i: $moveVal $depth\n")
            # if the move is better than the current best, keep it
            # if it's equal, sometimes keep it
            if (moveVal > bestVal) || (moveVal==bestVal && depth<=bestDepth && rand()<0.75)
                # or if it's equivalent, sometimes keep it
                bestVal = moveVal
                bestMove = i
                bestDepth = depth
            end
        end
    end
    print("Choosing $bestMove\n")
    return bestMove # we are left with one of the optimal moves
end

function minimax(M,depth,isMax,player,moves,score,movePossible,makeMove!,gameOver)
    # minimax step
    depthOut = depth
    if gameOver(M)
        return isMax ? score(M,player) : -score(M,player), depthOut
    else
        best = isMax ? -1000 : 1000
        for i in moves # serach over candidate moves
            if movePossible(M,i)
                N = deepcopy(M)
                N = makeMove!(N,i,player) # play the move

                # pass on to the minimiser
                value, depthOut = minimax(N,depth+1,!isMax,3-player,moves,score,movePossible,makeMove!,gameOver)

                # did this move do better?
                best = isMax ? max(best, value) : min(best, value)
            end
        end
        return best, depthOut
    end
end

function getBestMoveAB(M,player,score,moves,movePossible,makeMove!,gameOver,maxDepth)
    # uses minimax algorithm to find best move

    # initialise algorithm values
    alpha = -1000
    beta = 1000
    bestVal = -1000
    bestMove = [-1]
    bestDepth = [1000]
    moveVal = 0

    for i in moves # search over all candidate moves
        if movePossible(M,i)
            N = deepcopy(M)
            N = makeMove!(N,i,player) # play the move
            moveVal, depth = minimaxAB(alpha,beta,N,0,false,3-player,
                moves,score,movePossible,makeMove!,gameOver,maxDepth) # evaluate it

            print("$i: $moveVal $depth ($bestVal)\n")
            # if the move is better than the current best, keep it
            # if it's equal, sometimes keep it

            if moveVal > bestVal
                bestVal = moveVal
                bestMove = [i]
                bestDepth = [depth]
            elseif moveVal==bestVal
                push!(bestMove,i)
                push!(bestDepth,depth)
            end
        end
    end
    if bestVal<0
        goodIdxs = bestDepth.==maximum(bestDepth)
    else
        goodIdxs = bestDepth.==minimum(bestDepth)
    end
    bestMoveOptions = bestMove[goodIdxs]
    bestMove = length(bestMoveOptions)>1 ? rand(bestMoveOptions) : bestMoveOptions[1]
    print("From $bestMoveOptions, choosing $bestMove\n")
    return bestMove # we are left with one of the optimal moves
end

function minimaxAB(alpha,beta,M,depth,isMax,player,moves,score,
        movePossible,makeMove!,gameOver,maxDepth)
    # minimax step
    depthOut = depth
    if depth>maxDepth
        return score(M,player), depth
    elseif gameOver(M)
        return isMax ? score(M,player) : -score(M,player), depth
    else
        bestValue = isMax ? -1000 : 1000
        for i in moves # serach over candidate moves
            if movePossible(M,i)
                N = deepcopy(M)
                N = makeMove!(N,i,player) # play the move
                valueTemp, depthOut = minimaxAB(alpha,beta,N,depth+1,!isMax,
                    3-player,moves,score,movePossible,makeMove!,gameOver,maxDepth)

                if isMax
                    bestValue = max(bestValue,valueTemp)
                    alpha = max(alpha,bestValue)
                else
                    bestValue = min(bestValue,valueTemp)
                    beta = min(beta,bestValue)
                end
                if alpha>=beta
                    return bestValue, depthOut
                end
            end
        end

        return bestValue, depthOut
    end
end


text:
# return recordings of particle velocity by injecting windowed source cube
function MultiStepForward(nz::Int64, nx::Int64, irz::Array{Int64,1}, irx::Array{Int64,1},
                          wsc::WSrcCube, fidMtx::FidMtx; tmax=1.0)
    dt = wsc.dt; ext=wsc.ext; iflag=wsc.iflag;
    nt = round(Int64, tmax/dt)+1
    shot = InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddWsc2spt!(spt1, wsc)
    Spt2ShotV!(shot, spt1)
    tmp = zeros(length(spt1.Tzzz))
    tmp1= zeros(tmp);
    for it = 2 : wsc.nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        AddWsc2spt!(spt2, wsc)
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shot, spt1)
    end
    for it = wsc.nt+1 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shot, spt1)
    end
    return shot
end

# return records (vx, vz), inject single source
function MultiStepForward(irz::Array{Int64,1}, irx::Array{Int64,1}, src::Source, fidMtx::FidMtx; tmax=0.5)
    nz  =  src.nz;  nx = src.nx    ;
    ext = src.ext;  iflag=src.iflag; dt = src.dt;
    stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
    nt  = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    shotv = InitShotV(src.isz, src.isx, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    Spt2ShotV!(shotv, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddSource!(spt2, src)
        end
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shotv, spt1)
    end
    return shotv
end

# return records (vx, vz), inject multiple sources
function MultiStepForward(irz::Array{Int64,1}, irx::Array{Int64,1}, srcs::Array{Source,1}, fidMtx::FidMtx; tmax=0.5)
    nz  = srcs[1].nz ; nx = srcs[1].nx    ;
    ext = srcs[1].ext; iflag=srcs[1].iflag; dt = srcs[1].dt;
    (stl, stu) = SourcesTimeRange(srcs)
    nt  = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddMultiSources!(spt1, srcs)
    shotv = InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    Spt2ShotV!(shotv, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddMultiSources!(spt2, srcs)
        end
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shotv, spt1)
    end
    return shotv
end

# write wave field(default, 5 components) or snapshot(10 components), inject single source
function MultiStepForward(path::String, src::Source, fidMtx::FidMtx; tmax=0.5, otype="wfd")
    nz  =  src.nz;  nx = src.nx    ;
    ext = src.ext;  iflag=src.iflag; dt  = src.dt;
    stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
    nt   = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    if otype == "wfd"
       fid = WriteWfd(path, spt1)
    elseif otype == "spt"
       fid = WriteSnapShot(path, spt1)
    end
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddSource!(spt2, src)
        end
        CopySnapShot!(spt1, spt2)
        if otype == "wfd"
           WriteWfd(fid, spt1)
        elseif otype == "spt"
           WriteSnapShot(fid, spt1)
        end
    end
    close(fid)
    return nothing
end

# save the boundary of source side wave field
function MultiStepForward(path::String, src::Source, fidMtx::FidMtx, tmax::Float64)
    nz  =  src.nz;  nx = src.nx    ;
    ext = src.ext;  iflag=src.iflag; dt  = src.dt;
    stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
    nt   = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    bv = zeros((4*nz+(nx-4)*4)*5);
    fid = WriteWfdBoundary(path, spt1, bv)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddSource!(spt2, src)
        end
        CopySnapShot!(spt1, spt2)
        WriteWfdBoundary(fid, spt1, bv)
    end
    close(fid)
    return nothing
end

# write full wavefield, inject multiple souces
function MultiStepForward(path::String, srcs::Array{Source,1}, fidMtx::FidMtx; tmax=0.5, otype="wfd")
    nz  = srcs[1].nz ; nx   = srcs[1].nx;
    ext = srcs[1].ext; iflag= srcs[1].iflag; dt = srcs[1].dt;
    (stl, stu) = SourcesTimeRange(srcs)
    nt   = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddMultiSources!(spt1, srcs)
    if otype == "wfd"
       fid = WriteWfd(path, spt1)
    elseif otype == "spt"
       fid = WriteSnapShot(path, spt1)
    end
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddMultiSources!(spt2, srcs)
        end
        CopySnapShot!(spt1, spt2)
        if otype == "wfd"
           WriteWfd(fid, spt1)
        elseif otype == "spt"
           WriteSnapShot(fid, spt1)
        end
    end
    close(fid)
    return nothing
end

# write (vxx, vxz, vzx, vzz); inject single source
function MultiStepForward(path::String, src::Source, fidMtx::FidMtx, dz::Float64, dx::Float64; tmax=0.5)
    nz  =  src.nz;  nx = src.nx    ;
    ext = src.ext;  iflag=src.iflag; dt  = src.dt;
    stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
    nt  = round(Int64, tmax/dt)+1
    spt1= InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    vxx=zeros(Nz,Nx); vxz=zeros(Nz,Nx);
    vzx=zeros(Nz,Nx); vzz=zeros(Nz,Nx);
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1);
    fid = WritePv(path, vxx, vxz, vzx, vzz, nz, nx, ext, iflag, dt);
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddSource!(spt2, src)
        end
        CopySnapShot!(spt1, spt2)
        partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1); WritePv(fid, vxx, vxz, vzx, vzz);
    end
    close(fid)
    return nothing
end

# write (vxx, vxz, vzx, vzz); inject multiple souces
function MultiStepForward(path::String, srcs::Array{Source,1}, fidMtx::FidMtx, dz::Float64, dx::Float64; tmax=0.5)
    nz  = srcs[1].nz ; nx   = srcs[1].nx   ;
    ext = srcs[1].ext; iflag= srcs[1].iflag; dt = srcs[1].dt;
    (stl, stu) = SourcesTimeRange(srcs)
    nt  = round(Int64, tmax/dt)+1
    spt1= InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddMultiSources!(spt1, srcs)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    vxx=zeros(Nz,Nx); vxz=zeros(Nz,Nx);
    vzx=zeros(Nz,Nx); vzz=zeros(Nz,Nx);
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1);
    fid = WritePv(path, vxx, vxz, vzx, vzz, nz, nx, ext, iflag, dt);
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if stl <= (it-1)*dt <= stu
           AddMultiSources!(spt2, srcs)
        end
        CopySnapShot!(spt1, spt2)
        partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1); WritePv(fid, vxx, vxz, vzx, vzz);
    end
    close(fid)
    return nothing
end
# function MultiStepForward(path::String, src::Source, fidMtx::FidMtx, dz::Float64, dx::Float64; tmax=0.5)
#     nz  =  src.nz;  nx = src.nx    ;
#     ext = src.ext;  iflag=src.iflag; dt  = src.dt;
#     stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
#     nt  = round(Int64, tmax/dt)+1
#     spt1= InitSnapShot(nz, nx, ext, iflag, dt, 1)
#     spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
#     AddSource!(spt1, src)
#     vxx=zeros(nz,nx); vxz=zeros(nz,nx);
#     vzx=zeros(nz,nx); vzz=zeros(nz,nx);
#     tmp = zeros(length(spt1.Vxx))
#     tmp1= zeros(tmp)
#     partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1);
#     fid = WritePv(path, vxx, vxz, vzx, vzz, nz, nx, ext, iflag, dt);
#     for it = 2 : nt
#         OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
#         if stl <= (it-1)*dt <= stu
#            AddSource!(spt2, src)
#         end
#         CopySnapShot!(spt1, spt2)
#         partialV!(vxx, vxz, vzx, vzz, spt1, dz, dx, tmp, tmp1); WritePv(fid, vxx, vxz, vzx, vzz);
#     end
#     close(fid)
#     return nothing
# end

# output shotv(Vx, Vz), Born approximation
function MultiStepForward(irz::Array{Int64,1}, irx::Array{Int64,1}, path::String, dm::Array{Float64,1}, du::Array{Float64,1}, fidMtx::FidMtx; tmax=0.5)
    (nz, nx, ext, iflag, dt, ns) = InfoPv(path);
    nt  = round(Int64, tmax/dt)+1
    spt1= InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSourceBorn!(spt1, path, dm, du)
    shotv= InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    Spt2ShotV!(shotv, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if spt2.it <= ns
           AddSourceBorn!(spt2, path, dm, du)
        end
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shotv, spt1)
    end
    return shotv
end

# write wavefield to hard drive, do born approximation
function MultiStepForward(path_wfd, path::String, dm::Array{Float64,2}, du::Array{Float64,2}, fidMtx::FidMtx; tmax=0.5)
    (nz, nx, ext, iflag, dt, ns) = InfoPv(path);
    nt  = round(Int64, tmax/dt)+1
    spt1= InitSnapShot(nz, nx, ext, iflag, dt, 1)
    spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
    AddSourceBorn!(spt1, path, dm, du)
    fid = WriteWfd(path_wfd, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if spt2.it <= ns
           AddSourceBorn!(spt2, path, dm, du)
        end
        CopySnapShot!(spt1, spt2)
        WriteWfd(fid, spt1)
    end
    close(fid)
    return nothing
end

# output shotv(Vx, Vz), input snapshot
function MultiStepForward(irz::Array{Int64,1}, irx::Array{Int64,1}, path::String, fidMtx::FidMtx; tmax=0.5)
    (nz, nx, ext, iflag, dt, ns) = InfoSnapShots(path);
    nt  = round(Int64, tmax/dt)+1
    spt1= ReadSnapShot(path, 1)
    spt2= InitSnapShot(nz, nx, ext, iflag, dt, 2)
    shotv= InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    Spt2ShotV!(shotv, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = 2 : nt
        OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
        if spt2.it <= ns
           spt = ReadSnapShot(path, it)
           Add2SnapShots!(spt2, spt)
        end
        CopySnapShot!(spt1, spt2)
        Spt2ShotV!(shotv, spt1)
    end
    return shotv
end

# return shot(vx, vz, Txx, Tzz, Txz), inject multiple souces
# function MultiStepForward(irz::Array{Int64,1}, irx::Array{Int64,1}, srcs::Array{Source,1}, fidMtx::FidMtx; tmax=0.5)
#     nz  = srcs[1].nz ; nx = srcs[1].nx    ;
#     ext = srcs[1].ext; iflag=srcs[1].iflag; dt = srcs[1].dt;
#     (stl, stu) = SourcesTimeRange(srcs)
#     nt  = round(Int64, tmax/dt)+1
#     spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
#     spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
#     AddMultiSources!(spt1, srcs)
#     shot = InitShot(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
#     Spt2Shot!(shot, spt1)
    # tmp = zeros(length(spt1.Vxx))
    # tmp1= zeros(tmp)
#     for it = 2 : nt
#         OneStepForward!(spt2, spt1, fidMtx, tmp, tmp1)
#         if stl <= (it-1)*dt <= stu
#            AddMultiSources!(spt2, srcs)
#         end
#         CopySnapShot!(spt1, spt2)
#         Spt2Shot!(shot, spt1)
#     end
#     return shot
# end

max_stars_repo_name: WneleiGao/ElasticWave


text:
export distibution

export GaussianProduct

struct GaussianSum

    a  :: Float64
    b  :: Float64
    xa :: Float64
    xb :: Float64
    ra :: Float64
    rb :: Float64

end

"""

    distribution( mesh1, mesh2, gauss)

    ```math
    f(x,υ) = a \exp ( \frac{x-xa}{ra}^2 ) +
             b \exp ( \frac{υ-υb}{rb}^2 )
    ```

"""
function distribution( mesh1  :: UniformMesh,
                       mesh2  :: UniformMesh,
		       gauss  :: GaussianSum)

    x = mesh1.points
    v = mesh2.points
    a , b  = gauss.a,  gauss.b
    xa, xb = gauss.xa, gauss.xb
    ra, rb = gauss.ra, gauss.rb

    @. a * exp(-((x-xa)/ra)^2) * transpose(exp(-((v-xb)/rv)^2))

end

export GaussianProduct

struct GaussianProduct

    a  :: Float64
    ra :: Float64
    rb :: Float64

end

"""

    distribution( mesh1, mesh2, gauss)

    ```math
    f(x,υ) = a \exp (\frac{x}{ra}^2) ⋅ \exp (\frac{υ}{rb}^2)
    ```

"""
function distribution( mesh1  :: UniformMesh,
                       mesh2  :: UniformMesh,
		       gauss  :: GaussianProduct)

    x = mesh1.points
    v = mesh2.points
    rx, rv = gauss.rx, gauss.rv

    @. a * exp(-(x/rx)^2) * transpose(exp(-(v/rv)^2))

end


text:
using Random
Random.seed!(172)
using LinearAlgebra: Matrix
using CompressingSolvers
using SparseArrays
using LinearAlgebra
using MKLSparse

# Setting up the test domains
ρ = 4.0
# ρ = Inf
q = 4
# A, coarse_domains, scales, basis_functions, multicolor_ordering, fine_domains, tree_function = CompressingSolvers.FD_Laplacian_subdivision_2d(q, ρ)
α(x) = rand()
A, coarse_domains, scales, basis_functions, multicolor_ordering, fine_domains, tree_function = CompressingSolvers.FD_periodic_Laplacian_subdivision_2d(q, ρ, α)

measurement_matrix = CompressingSolvers.form_measurement_matrix(multicolor_ordering)
@show size(measurement_matrix)
@time measurement_results = CompressingSolvers.measure(cholesky(A), measurement_matrix)

@time L = CompressingSolvers.reconstruct(multicolor_ordering, CompressingSolvers.center.(fine_domains), measurement_matrix, measurement_results, tree_function)

# @show opnorm(inv(Matrix(A)) - L * L') / opnorm(inv(Matrix(A)))
@show CompressingSolvers.compute_relative_error(L, cholesky(A), 200, 200)
# 𝐅 = CompressingSolvers.SupernodalFactorization(multicolor_ordering, domain_supernodes);
#
# 𝐌 = CompressingSolvers.create_measurement_matrix(multicolor_ordering, 𝐅.row_supernodes);
#
# # 𝐎 = CompressingSolvers.measure(inv(Matrix(A)), 𝐌, 𝐅.row_supernodes);
# @time 𝐎 = CompressingSolvers.measure(cholesky(A), 𝐌, 𝐅.row_supernodes);
#
# 𝐇 = hcat(SparseMatrixCSC.(vcat(𝐅.column_supernodes...))...);
#
# @time CompressingSolvers.reconstruct!(𝐅, 𝐎, 𝐌, multicolor_ordering);
# L = SparseMatrixCSC(𝐅)
#
# # @show norm(L * L' - inv(Matrix(A))) / norm(inv(Matrix(A)))
# @show length(multicolor_ordering)
# @show length(reduce(vcat, multicolor_ordering))
# @show CompressingSolvers.compute_relative_error(L, cholesky(A))


text:
#   Author: V. Vitale
#   Feb 2022

using TensorOperations
using KrylovKit
using LinearAlgebra
using Random
using RandomMatrices
using JLD2, FileIO
using NPZ
abstract type AbstractTN end

"""
    length(::MPS/MPO)
The number of sites of an MPS/MPO.
"""

Base.:length(m::AbstractTN) = m.N

data(m::AbstractTN) = m.data

Base.:size(m::AbstractTN) = m.N

dims(m::AbstractTN) = Dict([key=>Base.size(m.data[key]) for key in keys(data(m))])

function writeTN(filename::String,m::AbstractTN)
    save(filename*".jld2", "data", m)
end

function readTN(filename::String)
    m=load(filename*".jld2")["data"]
    return m
end

function copy(A::AbstractTN)
    N=length(A)
    if typeof(A)==MPS
        temp=MPS()
    elseif typeof(A)==MPO
        temp=MPO()
    end
    temp.N=N
    for i in 1:N
        temp.data[i]=Base.copy(A.data[i])
    end
    return temp
end

++(A::AbstractArray, B::AbstractArray) = cat(A, B,dims=(1,2))
const ⨁=++

function direct_sum(A::AbstractTN, B::AbstractTN)
    N=length(A)
    temp=copy(A)
    for i in 1:L
        temp.data[i]=++(A.data[i],B.data[i])
    end
    return temp
end

function direct_sub(A::AbstractTN, B::AbstractTN)
    N=length(A)
    temp=copy(A)
    println(typeof(temp))
    for i in 1:L
        temp.data[i]=++(A.data[i],-B.data[i])
    end
    return temp
end

Base.:+(A::AbstractTN, B::AbstractTN) = direct_sum(A, B)
Base.:-(A::AbstractTN, B::AbstractTN) = direct_sub(A, B)

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
using Plots
using Reproduce
using JLD2
using FileIO

function plot_single_reward(results, log_interval, logger_key)
    previous_GKSwstype = get(ENV, "GKSwstype", "")
    ENV["GKSwstype"] = "100"

    save_file = "plotting/chunlok/generated_plots/single_reward.svg"

    p = plot()

    error = results[logger_key]

    for i in 1:4
        # println(size(error[i, :]))
        plot!(p, 1:log_interval:size(error)[2] * log_interval, error[i, :])
    end

    # print(error[3,10])
    savefig(save_file)
    println("Single reward plot saved to $(save_file)")

    ENV["GKSwstype"] = previous_GKSwstype
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
import Pkg;
Pkg.activate(".")
# %%
path = pwd() * "/functions"
push!(LOAD_PATH,path)
# %%
using PyCall
sym = pyimport("sympy")
# %%
p, ξ, δ, ω, μᵖ, μᶠ, c₀ᶜ , c₀⁰ , n₀ᶠ, Nᶜ , N⁰ , Nᶠ, α₁, α₂, nᶜ, ρ₀ᶜ, ρ₁ᶜ, ρ₂ᶜ, n⁰, ρ₀⁰, ρ₁⁰, ρ₂⁰, nᶠ, ρ₁ᶠ, ρ₂ᶠ, E₀ᶜ, Eᶜ, E₀⁰, E⁰, α₀, β₀, Wᶜ, W⁰, n₀ᶜ, n₀⁰, ν₀ᶜ, ν₁ᶜ, ν₂ᶜ, ν₀⁰, ν₁⁰, ν₂⁰, W₀ᶜ, W₀⁰, W∞ᶜ, W∞⁰ = sym.symbols("p, ξ, δ, ω , μᵖ , μᶠ, c₀ᶜ , c₀⁰ , n₀ᶠ, Nᶜ , N⁰ , Nᶠ, α₁, α₂, nᶜ, ρ₀ᶜ, ρ₁ᶜ, ρ₂ᶜ, n⁰, ρ₀⁰, ρ₁⁰, ρ₂⁰, nᶠ, ρ₁ᶠ, ρ₂ᶠ, E₀ᶜ, Eᶜ, E₀⁰, E⁰, α₀, β₀, Wᶜ, W⁰, n₀ᶜ, n₀⁰, ν₀ᶜ, ν₁ᶜ, ν₂ᶜ, ν₀⁰, ν₁⁰, ν₂⁰, W₀ᶜ, W₀⁰, W∞ᶜ, W∞⁰")
# %%
# ====================================================== #
# ================= Dual to dual ======================= #
# ====================================================== #
system = [(p+ξ)*Nᶜ + (ω+ξ)*N⁰ - c₀ᶜ, μᵖ*Nᶜ + (p + μᵖ - ω)*N⁰ + μᵖ*Nᶠ - ( μᵖ/p + c₀⁰ ), μᶠ*Nᶜ + μᶠ*N⁰ + (p+δ+μᶠ)*Nᶠ - (μᶠ/p + n₀ᶠ) ]
sol = sym.linsolve(system, [Nᶜ, N⁰, Nᶠ])
# %%
sol = sym.factor(sol, p)
# %%
Nᶜ = sol.args[1][1]
N⁰ = sol.args[1][2]
Nᶠ = sol.args[1][3]
# %%
top_Nᶜ, bot_Nᶜ = sym.fraction(Nᶜ)
top_Nᶜ =  sym.Poly(top_Nᶜ,p)
for c in top_Nᶜ.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# %%
top_N⁰, bot_N⁰ = sym.fraction(N⁰)
top_N⁰ =  sym.Poly(top_N⁰,p)
for c in top_N⁰.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# %%
top_Nᶠ, bot_Nᶠ = sym.fraction(Nᶠ)
top_Nᶠ =  sym.Poly(top_Nᶠ,p)
for c in top_Nᶠ.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# Laplace transforms seem correct !
# %%
# ======================= Nᶜ ========================= #
bot_Nᶜ = sym.Poly(bot_Nᶜ, p)
roots = sym.roots(bot_Nᶜ)
# %%
expr = sym.cancel(top_Nᶜ*p/(p*(p-ω)*(p-α₁)*(p-α₂)))
expr = expr.subs(p , 0)
# nᶜ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_Nᶜ/(p*(p-ω)*(p-α₁)*(p-α₂)) - nᶜ/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀ᶜ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀ᶜ/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p-α₁)/bot)
expr = sym.factor(expr.subs(p , α₁),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ])
# ρ₁ᶜ : Checks out !
# %%
expr = sym.cancel(top/bot - ρ₁ᶜ/(p-α₁))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(α₁)
expr = sym.cancel(top*(p-α₂)/bot)
expr = sym.factor(expr.subs(p , α₂),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ, ρ₁ᶜ])
# %%
# ======================= N⁰ ========================= #
bot_N⁰ = sym.Poly(bot_N⁰, p)
roots = sym.roots(bot_N⁰)
# %%
expr = sym.cancel(top_N⁰*p/(p*(p-ω)*(p-α₁)*(p-α₂)))
expr = expr.subs(p , 0)
# n⁰ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_N⁰/(p*(p-ω)*(p-α₁)*(p-α₂)) - n⁰/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀⁰ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀⁰/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p-α₁)/bot)
expr = sym.factor(expr.subs(p , α₁),[c₀ᶜ, c₀⁰, n₀ᶠ, n⁰, ρ₀⁰])
# ρ₁⁰ : Checks out !
# %%
expr = sym.cancel(top/bot - ρ₁⁰/(p-α₁))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(α₁)
expr = sym.cancel(top*(p-α₂)/bot)
expr = sym.factor(expr.subs(p , α₂),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ, ρ₁ᶜ])
# ρ₂⁰ : Checks out !
# %%
# ======================= Nᶠ ========================= #
bot_Nᶠ = sym.Poly(bot_Nᶠ, p)
roots = sym.roots(bot_Nᶠ)
# %%
expr = sym.cancel(top_Nᶠ*p/(p*(p-α₁)*(p-α₂)))
expr = expr.subs(p , 0)
# nᶠ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_Nᶠ/(p*(p-α₁)*(p-α₂)) - nᶠ/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-α₁)/bot)
expr = sym.factor(expr.subs(p, α₁),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₁ᶠ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₁ᶠ/(p-α₁))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(α₁)
expr = sym.cancel(top*(p-α₂)/bot)
expr = sym.factor(expr.subs(p , α₂),[c₀ᶜ, c₀⁰, n₀ᶠ, n⁰, ρ₀⁰])
# ρ₂ᶠ : Checks out!
# %%
# ====================================================== #
# ================= Dual to classic ==================== #
# ====================================================== #
system = [(p+ξ)*Nᶜ + (ω+ξ)*N⁰ - c₀ᶜ, μᵖ*Nᶜ + (p + μᵖ - ω)*N⁰ + μᵖ*Nᶠ - ( μᵖ/p + c₀⁰ ), (p+δ)*Nᶠ - n₀ᶠ]
sol = sym.linsolve(system, [Nᶜ, N⁰, Nᶠ])
# %%
sol = sym.factor(sol, p)
# %%
Nᶜ = sol.args[1][1]
N⁰ = sol.args[1][2]
Nᶠ = sol.args[1][3]
# %%
top_Nᶜ, bot_Nᶜ = sym.fraction(Nᶜ)
top_Nᶜ =  sym.Poly(top_Nᶜ,p)
for c in top_Nᶜ.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# %%
top_N⁰, bot_N⁰ = sym.fraction(N⁰)
top_N⁰ =  sym.Poly(top_N⁰,p)
for c in top_N⁰.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# %%
top_Nᶠ, bot_Nᶠ = sym.fraction(Nᶠ)
top_Nᶠ =  sym.Poly(top_Nᶠ,p)
for c in top_Nᶠ.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# Laplace transforms seem correct !
# %%
# ======================= Nᶜ ========================= #
bot_Nᶜ = sym.Poly(bot_Nᶜ, p)
roots = sym.roots(bot_Nᶜ)
# %%
expr = sym.cancel(top_Nᶜ*p/(p*(p-ω)*(p + μᵖ + ξ)*(p + δ)))
expr = expr.subs(p , 0)
# nᶜ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_Nᶜ/(p*(p-ω)*(p + μᵖ + ξ)*(p + δ)) - nᶜ/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀ᶜ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀ᶜ/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p + μᵖ + ξ)/bot)
expr = sym.factor(expr.subs(p , -(μᵖ + ξ)),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ])
# ρ₁ᶜ : Checks out !
# %%
expr = sym.cancel(top/bot - ρ₁ᶜ/(p + μᵖ + ξ))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(-(μᵖ + ξ))
expr = sym.cancel(top*(p + δ)/bot)
expr = sym.factor(expr.subs(p , -δ),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ, ρ₁ᶜ])
# %%
# ======================= N⁰ ========================= #
bot_N⁰ = sym.Poly(bot_N⁰, p)
roots = sym.roots(bot_N⁰)
# %%
expr = sym.cancel(top_N⁰*p/(p*(p-ω)*(p + μᵖ + ξ)*(p + δ)))
expr = expr.subs(p , 0)
# n⁰ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_N⁰/(p*(p-ω)*(p + μᵖ + ξ)*(p + δ)) - n⁰/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀⁰ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀⁰/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p + μᵖ + ξ)/bot)
expr = sym.factor(expr.subs(p , -(μᵖ + ξ)),[c₀ᶜ, c₀⁰, n₀ᶠ, n⁰, ρ₀⁰])
# ρ₁⁰ : Checks out !
# %%
expr = sym.cancel(top/bot - ρ₁⁰/(p + μᵖ + ξ))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(-(μᵖ + ξ))
expr = sym.cancel(top*(p+δ)/bot)
expr = sym.factor(expr.subs(p , -δ),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ, ρ₁ᶜ])
# ρ₂⁰ : Checks out !
# %%
# ====================================================== #
# ================= Classic to classic ================= #
# ====================================================== #
system = [(p+ξ)*Nᶜ + (ω+ξ)*N⁰ - c₀ᶜ, μᵖ*Nᶜ + (p + μᵖ - ω)*N⁰ - ( μᵖ/p + c₀⁰ )]
sol = sym.linsolve(system, [Nᶜ, N⁰])
# %%
sol = sym.factor(sol, p)
# %%
Nᶜ = sol.args[1][1]
N⁰ = sol.args[1][2]
# %%
top_Nᶜ, bot_Nᶜ = sym.fraction(Nᶜ)
top_Nᶜ =  sym.Poly(top_Nᶜ,p)
for c in top_Nᶜ.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# %%
top_N⁰, bot_N⁰ = sym.fraction(N⁰)
top_N⁰ =  sym.Poly(top_N⁰,p)
for c in top_N⁰.coeffs()
    println(sym.factor(c, [c₀ᶜ, c₀⁰, n₀ᶠ]))
end
# Laplace transforms seem correct !
# %%
# ======================= Nᶜ ========================= #
bot_Nᶜ = sym.Poly(bot_Nᶜ, p)
roots = sym.roots(bot_Nᶜ)
# %%
expr = sym.cancel(top_Nᶜ*p/(p*(p-ω)*(p + μᵖ + ξ)))
expr = expr.subs(p , 0)
# nᶜ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_Nᶜ/(p*(p-ω)*(p + μᵖ + ξ)) - nᶜ/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀ᶜ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀ᶜ/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p + μᵖ + ξ)/bot)
expr = sym.factor(expr.subs(p , -(μᵖ + ξ)),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ])
# ρ₁ᶜ : Checks out !
# %%
# ======================= N⁰ ========================= #
bot_N⁰ = sym.Poly(bot_N⁰, p)
roots = sym.roots(bot_N⁰)
# %%
expr = sym.cancel(top_N⁰*p/(p*(p-ω)*(p + μᵖ + ξ)))
expr = expr.subs(p , 0)
# n⁰ : Checks out!
# %%
expr = sym.factor(sym.cancel(top_N⁰/(p*(p-ω)*(p + μᵖ + ξ)) - n⁰/p),p)
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(0)
expr = sym.cancel(top*(p-ω)/bot)
expr = sym.factor(expr.subs(p , ω),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ])
# ρ₀⁰ : Checks out!
# %%
expr = sym.cancel(top/bot - ρ₀⁰/(p-ω))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(ω)
expr = sym.cancel(top*(p + μᵖ + ξ)/bot)
expr = sym.factor(expr.subs(p , -(μᵖ + ξ)),[c₀ᶜ, c₀⁰, n₀ᶠ, n⁰, ρ₀⁰])
# ρ₁⁰ : Checks out !
# %%
expr = sym.cancel(top/bot - ρ₁⁰/(p + μᵖ + ξ))
top,bot = sym.fraction(expr)
top = sym.Poly(top, p)
top = top - top(-(μᵖ + ξ))
expr = sym.cancel(top*(p+δ)/bot)
expr = sym.factor(expr.subs(p , -δ),[c₀ᶜ, c₀⁰, n₀ᶠ, nᶜ, ρ₀ᶜ, ρ₁ᶜ])
# ρ₂⁰ : Checks out !
# %%
# ====================================================== #
# ======================= Wages ======================== #
# ====================================================== #
system = [(p-ω)*Wᶜ - W₀ᶜ + (ω+ξ)*Eᶜ*(nᶜ/p + ρ₀ᶜ/(p-ω) + ρ₁ᶜ/(p-α₁) + ρ₂ᶜ/(p-α₂) + n⁰/p + ρ₀⁰/(p-ω) + ρ₁⁰/(p-α₁) + ρ₂⁰/(p-α₂)) + α₀*E₀ᶜ*n₀ᶜ, (p-ω)*W⁰ - W₀⁰ - μᵖ*E⁰*(1/p - (nᶜ/p + ρ₀ᶜ/(p-ω) + ρ₁ᶜ/(p-α₁) + ρ₂ᶜ/(p-α₂)) - (n⁰/p + ρ₀⁰/(p-ω) + ρ₁⁰/(p-α₁) + ρ₂⁰/(p-α₂)) - (nᶠ/p + ρ₁ᶠ/(p-α₁) + ρ₂ᶠ/(p-α₂))) + β₀*E₀⁰*n₀⁰]
sol = sym.linsolve(system, [Wᶜ, W⁰])
# %%
Wᶜ = sol.args[1][1]
W⁰ = sol.args[1][2]
# %%
top_Wᶜ, bot_Wᶜ = sym.fraction(Wᶜ)
top_Wᶜ =  sym.Poly(top_Wᶜ,p)
# %%
top_W⁰, bot_W⁰ = sym.fraction(W⁰)
top_W⁰ =  sym.Poly(top_W⁰,p)
# Laplace transforms seem correct !
# %%
# ======================= Wᶜ ========================= #
bot_Wᶜ = sym.Poly(bot_Wᶜ, p)
roots = sym.roots(bot_Wᶜ)
# %%
expr = sym.cancel(top_Wᶜ/((p-ω)^2*(p-α₁)*(p-α₂)))
W∞ᶜ = sym.cancel(expr.subs(p , 0))
# %%
expr = sym.cancel(top_Wᶜ/(p*(p-α₁)*(p-α₂)))
ν₀ᶜ = sym.cancel(expr.subs(p , ω))
# %%
expr = sym.cancel(top_Wᶜ/(p*(p-ω)*(p-α₁)*(p-α₂)) - ν₀ᶜ/(p-ω))
expr = expr.subs(p, ω)
sym.apart(expr, ω)
# %%
expr = sym.cancel(top_Wᶜ/(p*(p-ω)^2*(p-α₂)))
ν₁ᶜ = sym.cancel(expr.subs(p , α₁))
# ν₁ᶜ : Checks out !
# %%
expr = sym.cancel(top_Wᶜ/(p*(p-ω)^2*(p-α₁)))
ν₂ᶜ = sym.cancel(expr.subs(p , α₂))
# ν₂ᶜ : Checks out !
# %%
# ======================= W⁰ ========================= #
bot_W⁰ = sym.Poly(bot_W⁰, p)
roots = sym.roots(bot_W⁰)
# %%
expr = sym.cancel(top_W⁰/((p-ω)^2*(p-α₁)*(p-α₂)))
sym.cancel(expr.subs(p , 0))
# %%
expr = sym.cancel(top_W⁰/(p*(p-α₁)*(p-α₂)))
ν₀⁰ = sym.cancel(expr.subs(p , ω))
# %%
expr = sym.cancel(top_W⁰/(p*(p-ω)*(p-α₁)*(p-α₂)) - ν₀⁰/(p-ω))
expr = expr.subs(p, ω)
sym.apart(expr, ω)
# %%
expr = sym.cancel(top_W⁰/(p*(p-ω)^2*(p-α₂)))
ν₁⁰ = sym.cancel(expr.subs(p , α₁))
# ν₁ᶜ : Checks out !
# %%
expr = sym.cancel(top_W⁰/(p*(p-ω)^2*(p-α₁)))
ν₂⁰ = sym.cancel(expr.subs(p , α₂))
#


text:
"""
$(SIGNATURES)

Convert a judoc html string into a html string (i.e. replace `{{ ... }}` blocks).
"""
function convert_html(hs::AS, allvars::PageVars; isoptim::Bool=false)::String
    # Tokenize
    tokens = find_tokens(hs, HTML_TOKENS, HTML_1C_TOKENS)

    # Find hblocks ({{ ... }})
    hblocks, tokens = find_all_ocblocks(tokens, HTML_OCB)
    filter!(hb -> hb.name != :COMMENT, hblocks)

    # Find qblocks (qualify the hblocks)
    qblocks = qualify_html_hblocks(hblocks)
    # Find overall conditional blocks (if ... elseif ... else ...  end)
    cblocks, qblocks = find_html_cblocks(qblocks)
    # Find conditional def blocks (isdef / isnotdef)
    cdblocks, qblocks = find_html_cdblocks(qblocks)
    # Find conditional page blocks (ispage / isnotpage)
    cpblocks, qblocks = find_html_cpblocks(qblocks)

    # Get the list of blocks to process
    hblocks = merge_blocks(qblocks, cblocks, cdblocks, cpblocks)

    # construct the final html
    htmls = IOBuffer()
    head = 1
    for (i, hb) ∈ enumerate(hblocks)
        fromhb = from(hb)
        (head < fromhb) && write(htmls, subs(hs, head, prevind(hs, fromhb)))
        write(htmls, convert_html_block(hb, allvars))
        head = nextind(hs, to(hb))
    end
    strlen = lastindex(hs)
    (head < strlen) && write(htmls, subs(hs, head, strlen))

    fhs = String(take!(htmls))

    # See issue #204, basically not all markdown links are processed  as
    # per common mark with the JuliaMarkdown, so this is a patch that kind
    # of does
    if haskey(allvars, "reflinks") && allvars["reflinks"].first
        fhs = find_and_fix_md_links(fhs)
    end

    # if it ends with </p>\n but doesn't start with <p>, chop it off
    # this may happen if the first element parsed is an ocblock (not text)
    δ = ifelse(endswith(fhs, "</p>\n") && !startswith(fhs, "<p>"), 5, 0)

    isempty(fhs) && return ""

    if !isempty(GLOBAL_PAGE_VARS["prepath"].first) && isoptim
        fhs = fix_links(fhs)
    end

    return String(chop(fhs, tail=δ))
end


text:
function norm2_metric(x)
	m, n = size(x)
	D = zeros(eltype(x), n, n)
	@inbounds for i in 1:n
		for j in i+1:n
			D[i,j] = D[j,i] = norm(view(x, :, i) .- view(x, :, j))
		end
	end
	D
end

function norm2_metric_back(x, D, sym_dD)
	dD = 0.5*(sym_dD + sym_dD')  # Cholesky decomposition only use lower triangle, reconstruct the true dD
	m, n = size(x)
	T = eltype(x)
	dx = zeros(T, m, n)
	@inbounds for i in 1:m
		for j in 1:n
			for k in 1:n
				dx[i,j] += 2.0*(dD[j,k]/(D[j,k]+eps(T)))*(x[i,j]-x[i,k])
			end
		end
	end
	dx
end

function norm2_metric(x, xo)
	N1 = size(x, 2)
	N2 = size(xo, 2)
	D = zeros(eltype(x), N1, N2)
	@inbounds for i in 1:N1
		for j in 1:N2
			D[i,j] = norm(view(x, :, i) .- view(xo, :, j))
		end
	end
	D
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
using Multilane
using MCTS
using RobustMCTS
using POMCP
using GenerativeModels
using POMDPToolbox
using JLD
using POMDPs
using DataFrames
using Plots

N=500

# filename = "initials_Jul_7_13_02.jld"
filename = "initials_Jul_27_23_05.jld"
initials = load(filename)
problems = initials["problems"]
initial_states = initials["initial_states"]

representative_problem = Dict(take(problems,1))

# Calculate single points
point_solvers = Dict{String, Solver}(
    "random" => RandomSolver(),
    "heuristic" => SimpleSolver()
)

point_results = evaluate(representative_problem, initial_states, point_solvers, parallel=true, N=N)

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

rsolver = RobustMCTSSolver(
    depth=20,
    c=50.0,
    n_iterations=2500,
    k_state=4.0,
    alpha_state=1/8,
    k_nature=4.0,
    alpha_nature=1/8,
    rollout_solver=SimpleSolver(),
    rollout_nature=StochasticBehaviorNoCrashMDP(collect(values(problems))[1]))

pomcps = POMCPDPWSolver(
    eps=0.001,
    max_depth=20,
    c=50.0,
    tree_queries=2500,
    k_observation=4.0,
    alpha_observation=1/8,
    rollout_solver=SimpleSolver()
)

curve_solvers = Dict{String, Solver}(
    "upper_bound" => dpws,
    "robust" => RobustMLSolver(rsolver),
    "single_behavior" =>
        SingleBehaviorSolver(dpws,
                 IDMMOBILBehavior("normal", 30.0, 10.0, 1)),
    "pomcp" => MLPOMDPSolver(pomcps, BehaviorRootUpdaterStub(0.1))
)

curve_results = evaluate(problems, initial_states, curve_solvers, parallel=true, N=N)

results = merge_results!(curve_results, point_results)

println(results["stats"])

filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results["histories"] = nothing
save(filename, results)
println("results saved to $filename")

stats = results["stats"]
mean_performance = by(stats, :solver_key) do df
    by(df, :lambda) do df
        DataFrame(steps_in_lane=mean(df[:steps_in_lane]),
                  nb_brakes=mean(df[:nb_brakes]),
                  time_to_lane=mean(df[:time_to_lane]),
                  steps=mean(df[:steps]),
                  )
    end
end

unicodeplots()
for g in groupby(mean_performance, :solver_key)
    if size(g,1) > 1
        # plot!(g, :steps_in_lane, :nb_brakes, group=:solver_key)
        plot!(g, :time_to_lane, :nb_brakes, group=:solver_key)
    else
        # scatter!(g, :steps_in_lane, :nb_brakes, group=:solver_key)
        scatter!(g, :time_to_lane, :nb_brakes, group=:solver_key)
    end
end
gui()


text:
#! /usr/bin/julia

# Rosetta Code, Vogel's approximation method

immutable TProblem{T<:Integer,U<:String}
    sd::Array{Array{T,1},1}
    toc::Array{T,2}
    labels::Array{Array{U,1},1}
    tsort::Array{Array{T,2}, 1}
end

function TProblem{T<:Integer,U<:String}(s::Array{T,1},
                                        d::Array{T,1},
                                        toc::Array{T,2},
                                        slab::Array{U,1},
                                        dlab::Array{U,1})
    scnt = length(s)
    dcnt = length(d)
    size(toc) = (scnt,dcnt) || error("Supply, Demand, TOC Size Mismatch")
    length(slab) == scnt || error("Supply Label Size Labels")
    length(dlab) == dcnt || error("Demand Label Size Labels")
    0 <= minimum(s) || error("Negative Supply Value")
    0 <= minimum(d) || error("Negative Demand Value")
    sd = Array{T,1}[]
    push!(sd, s)
    push!(sd, d)
    labels = Array{U,1}[]
    push!(labels, slab)
    push!(labels, dlab)
    tsort = Array{T,2}[]
    push!(tsort, mapslices(sortperm, toc, 2))
    push!(tsort, mapslices(sortperm, toc, 1))
    TProblem(sd, toc, labels, tsort)
end
isbalanced(tp::TProblem) = sum(tp.sd[1]) == sum(tp.sd[2])

type Resource{T<:Integer}
    dim::T
    i::T
    quant::T
    l::T
    m::T
    p::T
    q::T
end
function Resource{T<:Integer}(dim::T, i::T, quant::T)
    zed = zero(T)
    Resource(dim, i, quant, zed, zed, zed, zed)
end

isavailable(r::Resource) = 0 < r.quant
Base.isless(a::Resource, b::Resource) = a.p < b.p || (a.p == b.p && b.q < a.q)

function penalize!{T<:Integer,U<:String}(sd::Array{Array{Resource{T},1},1},
                                         tp::TProblem{T,U})
    avail = BitArray{1}[]
    for dim in 2:-1:1
        push!(avail, bitpack(map(isavailable, sd[dim])))
    end
    for dim in 1:2, r in sd[dim]
        if r.quant == 0
            r.l = r.m = r.p = r.q = 0
            continue
        end
        r.l == 0 || !avail[dim][r.l] || !avail[dim][r.m] || continue
        rsort = filter(x->avail[dim][x], vec(slicedim(tp.tsort[dim],dim,r.i)))
        rcost = vec(slicedim(tp.toc, dim, r.i))[rsort]
        if length(rsort) == 1
            r.l = r.m = rsort[1]
            r.p = r.q = rcost[1]
        else
            r.l, r.m = rsort[1:2]
            r.p = rcost[2] - rcost[1]
            r.q = rcost[1]
        end
    end
    nothing
end

function vogel{T<:Integer,U<:String}(tp::TProblem{T,U})
    sdcnt = collect(size(tp.toc))
    sol = spzeros(T, sdcnt[1], sdcnt[2])
    sd = Array{Resource{T},1}[]
    for dim in 1:2
        push!(sd, [Resource(dim, i, tp.sd[dim][i]) for i in 1:sdcnt[dim]])
    end
    while any(map(isavailable, sd[1])) && any(map(isavailable, sd[2]))
        penalize!(sd, tp)
        a = maximum([sd[1], sd[2]])
        b = sd[rem1(a.dim+1,2)][a.l]
        if a.dim == 2 # swap to make a supply and b demand
            a, b = b, a
        end
        expend = min(a.quant, b.quant)
        sol[a.i, b.i] = expend
        a.quant -= expend
        b.quant -= expend
    end
    return sol
end

sup = [50, 60, 50, 50]
slab = ["W", "X", "Y", "Z"]
dem = [30, 20, 70, 30, 60]
dlab = ["A", "B", "C", "D", "E"]
c = [16 16 13 22 17;
     14 14 13 19 15;
     19 19 20 23 50;
     50 12 50 15 11]

tp = TProblem(sup, dem, c, slab, dlab)
sol = vogel(tp)
cost = sum(tp.toc .* sol)

println("The solution is:")
print("        ")
for s in tp.labels[2]
    print(@sprintf "%4s" s)
end
println()
for i in 1:size(tp.toc)[1]
    print(@sprintf "    %4s" tp.labels[1][i])
    for j in 1:size(tp.toc)[2]
        print(@sprintf "%4d" sol[i,j])
    end
println()
end
println("The total cost is:  ", cost)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
Return either the statically known Val() or runtime length()
"""
@inline _size(a) = Val(length(a))
@inline _size(a::TupleVector{n}) where n = Val(n)

# Return static array from a set of arrays
@inline _first_static(a1::TupleVector, as...) = a1
@inline _first_static(a1, as...) = _first_static(as...)
@inline _first_static() = throw(ArgumentError("No TupleVector found in argument list"))

"""
Returns the common Val of the inputs (or else throws a DimensionMismatch)
"""
@inline function same_size(as...)
    s = Val(length(_first_static(as...)))
    _sizes_match(s, as...) || _throw_size_mismatch(as...)
    s
end
@inline _sizes_match(s::Val, a1, as...) = ((s == _size(a1)) ? _sizes_match(s, as...) : false)
@inline _sizes_match(s::Val) = true
@noinline function _throw_size_mismatch(as...)
    throw(DimensionMismatch("Sizes $(map(_size, as)) of input arrays do not match"))
end


text:
# Sequence Record
# ===============
#
# Sequence record with name and attached metadata.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
`SeqRecord{S,T}` is a type holding a named sequence of type `S`, along with
some arbitrary metadata of type `T`.
"""
type SeqRecord{S<:Sequence,T}
    name::StringField
    seq::S
    metadata::T
end

function SeqRecord(name::AbstractString, seq::Sequence, metadata=nothing)
    return SeqRecord(StringField(name), seq, metadata)
end

function (::Type{SeqRecord{S,T}}){S,T}()
    return SeqRecord{S,T}(StringField(), S(), T())
end

"""
    seqname(record)

Return the sequence name of `record`.
"""
function seqname(rec::SeqRecord)
    return rec.name
end

"""
    sequence(record)

Return the sequence of `record`.
"""
function sequence(rec::SeqRecord)
    return rec.seq
end

"""
    metadata(record)

Return the metadata of `record`.
"""
function metadata(rec::SeqRecord)
    return rec.metadata
end

function seqtype{S,T}(::Type{SeqRecord{S,T}})
    return S
end

function Base.:(==){T<:SeqRecord}(a::T, b::T)
    return a.name == b.name && a.seq == b.seq && a.metadata == b.metadata
end

function Base.copy{T <: SeqRecord}(seqrec::T)
    if seqrec.metadata == nothing
        # no copy method for the Void type
        return T(copy(seqrec.name), copy(seqrec.seq), nothing)
    else
        return T(copy(seqrec.name), copy(seqrec.seq), copy(seqrec.metadata))
    end
end

function Base.show(io::IO, rec::SeqRecord)
    println(io, summary(rec), ':')
    println(io, "  name: ", seqname(rec))
    println(io, "  sequence: ", string_compact(sequence(rec)))
      print(io, "  metadata: ", metadata(rec))
end


text:
include("Quad/Genus.jl")
include("Quad/GenusRep.jl")


text:
include("basic_setup.jl")

t_s = 100
t_x = 1
κ = mean_s / t_s
λ = 1 / t_s
ρ = 1 / t_x
μ = 1 / t_x
mean_x = mean_s * ρ / μ
gen = Trajectories.configuration_generator(sn, rn, [κ, λ], [ρ, μ], mean_s, mean_x)
system, initial = Trajectories.generate_configuration(gen, duration=100.0)

chain = Trajectories.chain(system, 1.0)
samples, acceptance = Trajectories.generate_mcmc_samples(initial, chain, 2^8, 2^18)

using Plots
using Statistics
using StatsBase
using LaTeXStrings

hist_acc = fit(Histogram, chain.accepted_list, nbins=30)
hist_rej = fit(Histogram, chain.rejected_list, hist_acc.edges[1])

acc_rate = hist_acc.weights ./ (hist_rej.weights .+ hist_acc.weights)

gr()
p = plot(
    hist_acc.edges[1],
    acc_rate,
    seriestype=:barbins,
    ylim=(0,1),
    legend=false,
    xlabel="duration of regrown segment",
    ylabel="ratio of accepted MCMC moves",
    size=(270,270),
)

savefig(p, projectdir("plots", "ratio_acceptance_rates.png"))
savefig(p, projectdir("plots", "ratio_acceptance_rates.tex"))


text:
#=
Author: Daniel Carrera (dcarrera@gmail.com)

Compute the accretion rate onto a small (up to Neptune size) planet
embedded in a protoplanetary disk.
=#

"""
Compute the gas accretion rate onto a planet, up to Neptune size, embedded
in a protoplanetary disk. This function implements Equation (B36) derived in
Carrera et al. (2018) which itself is adapted from Ginzburg et al. (2016).

**NOTE**: Equation (B36) of Carrera et al. is in units of M_earth/Myr but
this function returns values in untis of M_earth/year.

**NOTE**: The accretion rate for a planet with zero atmosphere diverges. You
need to initialize the planet's atmosphere to a non-zero value. A reasonable
choice is the gas mass contained inside the planet's Bondi radius. Furthermore,
the accretion rate is initially very high and small timesteps are required to
resolve the accretion correctly.

Example:

	#
	# Gas accretion rate (units: M_earth / year) for a small planet embedded
	# in a protoplanetary disk.
	#
	# Inputs:
	#
	M_core = 8.0  # Planet's core mass in Earth masses.
	M_atm = 0.02  # Planet's current H2 mass in Earth masses.
	T_disk = 500  # Local disk temperature in Kelvin.

	M_atm_dot = accretion_rate(M_core, M_atm, T_disk) # M_earth / year.

Citation:

	Carrera et al. (2018) and Ginzburg et al. (2016)

	http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1804.05069
	http://adsabs.harvard.edu/abs/2016ApJ...825...29G

### Initializing the planet's atmosphere:

	ρ = ... # Local disk gas density in g/cm^3
	T = ... # Local disk gas temperature in Kelvin

	M_core = ... # Core mass in M_earth

	Rb = 2.861e7 * M_core / T # Bondi radius in cm

	M_atm = ρ * Rb^3 * 1.6744e-28 # Initial atmosphere mass in M_earth
"""
function accretion_rate(M_core::Number, M_atm::Number, T::Number)
	return (4.5e-10 / M_atm) * M_core^3.625 / sqrt(T/1000)
end


text:
@testset "Finite Diff Dynamics" begin
## Test with ForwardDiff
RD.diffmethod(::RobotZoo.Cartpole) = RD.ForwardAD()
prob,opts = Problems.Cartpole()
Z0 = copy(prob.Z)
solver = Altro.iLQRSolver(prob, opts, show_summary=false)
solve!(solver)
@test solver.cache.fx1 == zeros(size(prob)[1])
@test (@allocated solve!(solver)) == 0
J = cost(solver)

solver = ALTROSolver(prob, opts, show_summary=false)
initial_trajectory!(solver, Z0)
solve!(solver)
initial_trajectory!(solver, Z0)
allocs = @allocated solve!(solver)

## Test with finite differencing
RD.diffmethod(::RobotZoo.Cartpole) = RD.FiniteDifference()
prob,opts = Problems.Cartpole()
solver = Altro.iLQRSolver(prob, opts, show_summary=false)
solve!(solver)
@test solver.cache.fx1 != zeros(size(prob)[1])
@test (@allocated solve!(solver)) == 0

# make sure finite diff doesn't add any allocs for the entire solve
solver = ALTROSolver(prob, opts)
initial_trajectory!(solver, Z0)
solve!(solver)
initial_trajectory!(solver, Z0)
@test (@allocated solve!(solver)) == allocs
end


text:
export measure, dirac
export variables, monomials, moments

# If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined
struct Measure{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasure{T}
    a::Vector{T}
    x::MVT

    function Measure{T, MT, MVT}(a::Vector{T}, x::MVT) where {T, MT, MVT}
        @assert length(a) == length(x)
        new(a, x)
    end
end
Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike} = Measure{T, monomialtype(TT), monovectype(x)}(monovec(a, x)...)

"""
    measure(a, X::AbstractVector{<:AbstractMonomial})

Creates a measure with moments `moment(a[i], X[i])` for each `i`.
"""
measure(a, X) = Measure(a, X)

"""
    variables(μ::AbstractMeasureLike)

Returns the variables of `μ` in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial.
"""
MP.variables(μ::Measure) = variables(μ.x)

"""
    monomials(μ::AbstractMeasureLike)

Returns an iterator over the monomials of `μ` sorted in the decreasing order.
"""
MP.monomials(μ::Measure) = μ.x

"""
    moments(μ::AbstractMeasureLike)

Returns an iterator over the moments of `μ` sorted in decreasing order of monomial.
"""
moments(μ::Measure) = map((α, x) -> moment(α, x), μ.a, μ.x)

Base.:(*)(α, μ::Measure) = measure(α * μ.a, μ.x)
Base.:(*)(μ::Measure, α) = measure(μ.a * α, μ.x)
Base.:(-)(μ::Measure) = measure(-μ.a, μ.x)
function Base.:(+)(μ::Measure, ν::Measure)
    @assert μ.x == ν.x
    measure(μ.a + ν.a, μ.x)
end

"""
    dirac(X::AbstractVector{<:AbstractMoment}, s::AbstractSubstitution...)

Creates the dirac measure by evaluating the moments of `X` using `s`.

## Examples

Calling `dirac([x*y, x*y^2], x=>3, y=>2)` should the measure with moment `x*y` of value `6` and moment `x*y^2` of value `12`.
"""
function dirac(x::AbstractVector{MT}, s::MP.AbstractSubstitution...) where {MT <: AbstractMonomial}
    Measure([m(s...) for m in x], x)
end


text:
//test 2 in julia
function a()
	x = 0
    y = 1
	while <= x 5 do
        print(x)
		x = + x 1
	end
end

max_stars_repo_name: tylerholmes/Interpreter

max_issues_repo_name: tylerholmes/Interpreter

max_forks_repo_name: tylerholmes/Interpreter