text:
@testset "Linear" begin
    xmax = 10
    g1(x) = sin((x-3)*2pi/(xmax-1)-1)
    f(x) = g1(x)
    A1 = Float64[f(x) for x in 1:xmax]
    fr(x) = (x^2) // 40 + 2

    ymax = 10
    g2(y) = cos(y/6)
    f(x,y) = g1(x)*g2(y)
    A2 = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]

    for (constructor, copier) in ((interpolate, identity), (interpolate!, copy))
        isinplace = constructor == interpolate!
        itp1 = @inferred(constructor(copier(A1), BSpline(Linear())))
        itp2 = @inferred(constructor(copier(A2), BSpline(Linear())))

        @test parent(itp1) === itp1.coefs
        @test Interpolations.lbounds(itp1) == (1,)
        @test Interpolations.ubounds(itp1) == (xmax,)

        for (itp, A) in ((itp1, A1), (itp2, A2))
            check_axes(itp, A, isinplace)
            check_inbounds_values(itp, A)
            check_oob(itp)
            can_eval_near_boundaries(itp)
            I = first(eachindex(itp))
            @test itp(I) == itp(Tuple(I)...)
        end

        # Just interpolation
        for x in 1:.2:xmax
            @test f(x) ≈ itp1(x) atol=abs(0.1 * f(x))
        end

        # 2D
        for x in 2.1:.2:xmax-1, y in 1.9:.2:ymax-.9
            @test ≈(f(x,y),itp2(x,y),atol=abs(0.25 * f(x,y)))
        end

        # Rational element types
        A1R = Rational{Int}[fr(x) for x in 1:10]
        itp1r = @inferred(constructor(copier(A1R), BSpline(Linear())))
        @test @inferred(size(itp1r)) == size(A1R)
        @test itp1r(23 // 10) ≈ fr(23 // 10) atol=abs(0.1 * fr(23 // 10))
        @test typeof(itp1r(23//10)) == Rational{Int}
        @test eltype(itp1r) == Rational{Int}
    end

    # Issue #183
    x = rand(3,3,3)
    itp = interpolate(x, BSpline(Linear()))
    @test itp(1.5, CartesianIndex((2, 3))) === itp(1.5, 2, 3)
    @test itp(CartesianIndex((1, 2)), 1.5) === itp(1, 2, 1.5)

    # Issue #289
    for (x, it) in ((ones(1, 3), (NoInterp(), BSpline(Linear()))),
                    (ones(3, 1), (BSpline(Linear()), NoInterp())))
        @test_throws ArgumentError interpolate(x, BSpline(Linear()))
        itp = interpolate(x, it)
        etp = extrapolate(itp, Flat())
        @test sum(isnan.(etp)) == 0
    end

    @testset "Linear periodic" begin
        # Constructors
        @test Linear() === Linear(Throw(OnGrid())) === Linear{Throw{OnGrid}}()
        @test Linear() isa Linear{Throw{OnGrid}}
        @test Linear(Periodic()) === Linear(Periodic(OnCell())) == Linear{Periodic{OnCell}}()
        @test Linear(Periodic()) isa Linear{Periodic{OnCell}}

        xmax = 10
        f2(x) = sin((x-3)*2pi/xmax - 1) # Periodicity is xmax, not xmax-1
        A1 = Float64[f2(x) for x in 1:xmax]

        for (constructor, copier) in ((interpolate, identity), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp = @inferred(constructor(copier(A1), BSpline(Linear())))
            itp_periodic = @inferred(constructor(copier(A1), BSpline(Linear(Periodic()))))

            @test all(Interpolations.lbounds(itp_periodic) .≈ (0.5,))
            @test all(Interpolations.ubounds(itp_periodic) .≈ (10.5,))

            for x in 1:.2:xmax
                @test itp(x) ≈ itp_periodic(x)
            end

            # Interpolation range for periodic should be [0.5, xmax + 0.5]
            for x in 0.5:.1:0.9
                @test f2(x) ≈ itp_periodic(x) atol=abs(0.1 * f2(x))
                @test_throws BoundsError itp(x)
            end
            for x in xmax+0.1:.1:xmax+0.5
                @test f2(x) ≈ itp_periodic(x) atol=abs(0.1 * f2(x))
                @test_throws BoundsError itp(x)
            end
        end
    end
end


text:
using Images, FixedPointNumbers, Unitful, Test, Random
import ImagineFormat

μm = u"μm"
img = load("test.imagine")
@test eltype(img) == N2f14
@test ndims(img) == 4
@test size(img) == (5,7,3,4)
@test timedim(img) == 4
@test axisnames(img) == (:x, :l, :z, :time)
@test pixelspacing(img)[1:2] == (0.71μm, 0.71μm)
@test pixelspacing(img)[3] ≈ 100μm
@test img.properties["imagineheader"]["bidirectional"] == false
imgb = load("test_bidi.imagine")
@test imgb.properties["imagineheader"]["bidirectional"] == true

bn = joinpath(tempdir(), randstring())
ifn = string(bn, ".imagine")
cfn = string(bn, ".cam")
A = rand(Float32,2,3,4,5)
open(cfn, "w") do io
    write(io, A)
end
ImagineFormat.save_header(ifn, "test.imagine", A)
img2 = load(ifn)
@test eltype(img2) == Float32
@test data(img2) == A

using ImagineFormat
io = IOBuffer()
imagine2nrrd(io, img["imagineheader"])
str = String(take!(io))
@test str == "NRRD0001\ntype: uint16\ndimension: 4\nsizes: 5 7 3 4\nkinds: space space space time\nencoding: raw\nendian: little\n"

# Optional fields
MHz = u"MHz"
μs = u"μs"
@test isnan(img["imagineheader"]["readout rate"]) && isnan(img2["imagineheader"]["readout rate"]) #marked as NA in test.imagine
@test img["imagineheader"]["vertical shift speed"] == img2["imagineheader"]["vertical shift speed"] == 1.9176μs

h = ImagineFormat.parse_header("test_noshift.imagine")
h["byte order"] = ENDIAN_BOM == 0x04030201 ? "l" : "b"
h["stimulus scan hi"] = [3, 5]
ImagineFormat.save_header(ifn, h; misc=("stimulus scan hi",))
h2 = ImagineFormat.parse_header(ifn)
@test isnan(h["readout rate"]) && isnan(h2["readout rate"])
@test isnan(h["vertical shift speed"]) && isnan(h2["vertical shift speed"])
@test h2["stimulus scan hi"] == [3, 5]
rm(ifn)
img2 = 0
GC.gc(); GC.gc(); GC.gc()
rm(cfn)

#test BidiImageArrays
A = ones(2,2,4,5)
zsize = size(A,3)
for t = 1:size(A,4)
    for z = 1:zsize
        if isodd(t)
            A[:,:,z,t] .= z
        else
            A[:,:,z,t] .= zsize-z+1
        end
    end
end

B = BidiImageArray(A)
for z = 1:zsize
    @test all(B[:,:,z,:].==z) #getindex
    B[:,:,z,:] .= zsize-z+1 #setindex!
    @test all(B[:,:,z,:].==zsize-z+1)
end
A_copy = copy(A)
B[1,1,1,2] = -1.0
@test B[1,1,1,2] == -1.0
@test A[1,1,4,2] == -1.0
@test all(A[:,:,1:4,[1;3;4;5]] .== A_copy[:,:,1:4,[1;3;4;5]]) #make sure we only changed one index

img = load("no_z.imagine")
@test eltype(img) == N2f14
@test ndims(img) == 3
@test size(img) == (6,5,2)
@test timedim(img) == 3
@test axisnames(img) == (:x, :l, :time)
@test pixelspacing(img) == (-1μm, -1μm)

# Fix deprecations in loading imagine-cam file incompatibility
@info "Two warnings are expected"
img = load("test_nocam.imagine")
@test isempty(img)

img = load("test_wrong_nstack.imagine")
@test nimages(img) == 4

ext: jl


text:
using Serialization

"""
There are two newick trees and two julia objects containing these respective trees in
the current data format.

Newick:
tree1.nwk
tree2.nwk

Julia objects
tree1.jls
tree2.jls

You can read in the julia objects as follows:
"""
tree = deserialize("tree1.jls")

 # please supply the correct name of the function
parsed_tree = newick_parser("tree1.nwk")

"""
You can use the deserialized Julia objects to compare them against your trees.

There is a function called 'RF' which you can call on two trees. If the function
returns 0 the trees have an identical topology.
"""
topological_distance = RF(tree, parsed_tree)

"""
If you properly numbered all of the nodes of the tree you parsed from the newick
file, the function 'get_branchlength_vector' returns a vector of all branch lengths
in the tree. You can compare your branch lengths to the original branch lenghts.
Be aware, that the numbering of the nodes might be different and thus the order
of the two vectors. If you abstract from the ordering difference it mights still be
a reasonable indicator.
"""
original_branch = get_branchlength_vector(tree)
parsed_branch = get_branchlength_vector(parsed_tree)


text:
using NeuroIO, MAT, BinDeps, Base.Test

# Unfortunately, these test files come from the SDK provided by Plexon Inc. and
# are assumed not to be redistributable
testdir = joinpath(dirname(@__FILE__), "Matlab Offline Files SDK", "mexPlex", "tests")
if !isdir(testdir)
    parentdir = dirname(@__FILE__)
    bundlezip = joinpath(parentdir, "Plexon Offline SDKs.zip")
    run(download_cmd("http://www.plexon.com/sites/default/files/downloads/OmniPlex%20and%20MAP%20Offline%20SDK%20Bundle_0.zip",
                     bundlezip))
    cd(parentdir) do
        run(`unzip $bundlezip`)
    end
    cd(parentdir) do
        run(`unzip $(joinpath(parentdir, "Plexon Offline SDKs", "Matlab Offline Files SDK.zip"))`)
    end
end

to_vec(x::Vector) = x
to_vec(x::Array) = vec(x)
to_vec(x) = [x]

m = matread(joinpath(testdir, "mexPlexData1.dat"))
for plx in m["data"]["plxs"]
    p = readplx(open(joinpath(testdir, plx["FileName"])))
    @test p.header.Version == plx["Version"]
    @test p.header.ADFrequency == plx["Freq"]
    @test p.header.Comment == plx["Comment"]
    @test p.header.Trodalness == plx["Trodalness"]
    @test p.header.NumPointsWave == plx["NPW"]
    @test p.header.NumPointsPreThr == plx["PreThresh"]
    @test p.header.SpikeMaxMagnitudeMV == plx["SpikePeakV"]
    @test p.header.BitsPerSpikeSample == plx["SpikeADResBits"]
    @test p.header.SlowMaxMagnitudeMV == plx["SlowPeakV"]
    @test p.header.BitsPerSlowSample == plx["SlowADResBits"]
    @test p.header.LastTimestamp/p.header.ADFrequency == plx["Duration"]

    # Spikes
    n = plx["ts"]["n"]
    ts = plx["ts"]["ts"]
    wf = plx["wf"]["wf"]
    wf_v = plx["wf_v"]["wf"]
    for iunit = 1:30, ich = -1:130
        myts = ts[iunit+2, ich+2]
        if myts == -1
            @test !haskey(p.spike_channels, ich) ||
                  sum(p.spike_channels[ich].unit_numbers .== iunit) .== 0
        else
            ch = p.spike_channels[ich]
            inunit = find(ch.unit_numbers .== iunit)
            @test vec(myts) == NeuroIO.times(ch)[inunit]
            readwf = NeuroIO.data(ch)[:, inunit]
            @test wf[iunit+2, ich+2]' == readwf
            @test_approx_eq wf_v[iunit+2, ich+2]' readwf/3276.8
        end
    end

    # Continuous channels
    freq = plx["ad"]["freq"]
    ts = plx["ad"]["ts"]
    val = plx["ad"]["val"]
    val_v = plx["ad_v"]["val"]
    for ich = -1:400
        myts = ts[ich+2]
        if myts == -1
            @test !haskey(p.continuous_channels, ich+1) ||
                  isempty(NeuroIO.data(p.continuous_channels[ich+1]))
        else
            ch = p.continuous_channels[ich+1]
            mydata = vec(val[ich+2])
            @test first(myts) == first(ch.times)
            @test samplerate(ch) == freq[ich+2]
            @test NeuroIO.data(ch) == vec(val[ich+2])
            @test NeuroIO.data(ch)*(ch.voltage_multiplier*1000) == vec(val_v[ich+2])
        end
    end

    # Events
    tsevs = plx["tsevs"]
    for ich = -1:310
        myts = tsevs[ich+2]
        if myts == -1
            @test !haskey(p.event_channels, ich) ||
                  isempty(p.event_channels[ich].times)
        else
            ch = p.event_channels[ich]
            @test_approx_eq NeuroIO.times(ch) to_vec(myts)
            if ich == 257
                v = to_vec(plx["strobed"]["v"])
                @test v == NeuroIO.data(ch)
            end
        end
    end
end


text:
export BrickletXMC1400BreakoutADCValuesCallbackConfiguration
struct BrickletXMC1400BreakoutADCValuesCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

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

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

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

export BrickletXMC1400Breakout
"""
Breakout for Infineon XMC1400 microcontroller
"""
mutable struct BrickletXMC1400Breakout <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export set_gpio_config
"""
    $(SIGNATURES)

Example for a setter function. The values are the values that can be given to
the XMC_GPIO_Init function. See communication.c in the firmware.
"""
function set_gpio_config(device::BrickletXMC1400Breakout, port, pin, mode, input_hysteresis, output_level)
    device.deviceInternal.set_gpio_config(port, pin, mode, input_hysteresis, output_level)
end

export get_gpio_input
"""
    $(SIGNATURES)

Example for a getter function. Returns the result of a
XMC_GPIO_GetInput call for the given port/pin.
"""
function get_gpio_input(device::BrickletXMC1400Breakout, port, pin)
    return device.deviceInternal.get_gpio_input(port, pin)
end

export set_adc_channel_config
"""
    $(SIGNATURES)

Enables a ADC channel for the ADC driver example (adc.c/adc.h).

There are 8 ADC channels and they correspond to the following pins:

* Channel 0: P2_6
* Channel 1: P2_8
* Channel 2: P2_9
* Channel 3: P2_10
* Channel 4: P2_11
* Channel 5: P2_0
* Channel 6: P2_1
* Channel 7: P2_2
"""
function set_adc_channel_config(device::BrickletXMC1400Breakout, channel, enable)
    device.deviceInternal.set_adc_channel_config(channel, enable)
end

export get_adc_channel_config
"""
    $(SIGNATURES)

Returns the config for the given channel as set by :func:`Set ADC Channel Config`.
"""
function get_adc_channel_config(device::BrickletXMC1400Breakout, channel)
    return device.deviceInternal.get_adc_channel_config(channel)
end

export get_adc_channel_value
"""
    $(SIGNATURES)

Returns the 12-bit value of the given ADC channel of the ADC driver example.
"""
function get_adc_channel_value(device::BrickletXMC1400Breakout, channel)
    return device.deviceInternal.get_adc_channel_value(channel)
end

export get_adc_values
"""
    $(SIGNATURES)

Returns the values for all 8 ADC channels of the adc driver example.

This example function also has a corresponding callback.

See :func:`Set ADC Values Callback Configuration` and :cb:`ADC Values` callback.
"""
function get_adc_values(device::BrickletXMC1400Breakout)
    return device.deviceInternal.get_adc_values()
end

export set_adc_values_callback_configuration
"""
    $(SIGNATURES)

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

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

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_adc_values_callback_configuration(device::BrickletXMC1400Breakout, period, value_has_to_change)
    device.deviceInternal.set_adc_values_callback_configuration(period, value_has_to_change)
end

export get_adc_values_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set ADC Values Callback Configuration`.
"""
function get_adc_values_callback_configuration(device::BrickletXMC1400Breakout)
    return device.deviceInternal.get_adc_values_callback_configuration()
end

export get_count
"""
    $(SIGNATURES)

Returns the value of the example count (see example.c).

This example function uses the "add_callback_value_function"-helper in the
generator. The getter as well as the callback and callback configuration
functions are auto-generated for the API as well as the firmware.

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

export set_count_callback_configuration
"""
    $(SIGNATURES)

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

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

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

It is furthermore possible to constrain the callback with thresholds.

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

The following options are possible:

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

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

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

export get_count_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Count Callback Configuration`.
"""
function get_count_callback_configuration(device::BrickletXMC1400Breakout)
    return device.deviceInternal.get_count_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

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

The errors are divided into

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

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

export set_bootloader_mode
"""
    $(SIGNATURES)

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

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

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

export get_bootloader_mode
"""
    $(SIGNATURES)

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

export set_write_firmware_pointer
"""
    $(SIGNATURES)

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

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

export write_firmware
"""
    $(SIGNATURES)

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

You can only write firmware in bootloader mode.

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

export set_status_led_config
"""
    $(SIGNATURES)

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

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

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

export get_status_led_config
"""
    $(SIGNATURES)

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

export get_chip_temperature
"""
    $(SIGNATURES)

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

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

export reset
"""
    $(SIGNATURES)

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

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

export write_uid
"""
    $(SIGNATURES)

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

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

export read_uid
"""
    $(SIGNATURES)

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

export get_identity
"""
    $(SIGNATURES)

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

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

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

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


text:
##############################################################################
#
# REFERENCES: "Statistical Distributions"
#
##############################################################################

immutable Erlang <: ContinuousUnivariateDistribution
    shape::Int
    scale::Float64
    nested_gamma::Gamma

    function Erlang(shape::Real, scale::Real)
        isinteger(shape) || error("Erlang shape parameter must be an integer")
        new(int(shape), float64(scale), Gamma(shape, scale))
    end

    Erlang(scale::Real) = Erlang(scale, 1.0)
    Erlang() = Erlang(1, 1.0)
end

@continuous_distr_support Erlang 0.0 Inf

cdf(d::Erlang, x::Real) = cdf(d.nested_gamma, x)

entropy(d::Erlang) = entropy(d.nested_gamma)

kurtosis(d::Erlang) = kurtosis(d.nested_gamma)

mean(d::Erlang) = d.shape * d.scale

median(d::Erlang) = median(d.nested_gamma)

mgf(d::Erlang, t::Real) = mgf(d.nested_gamma, t)
cf(d::Erlang, t::Real) = cf(d.nested_gamma, t)

mode(d::Erlang) = mode(d.nested_gamma)
modes(d::Erlang) = modes(d.nested_gamma)

function pdf(d::Erlang, x::Real)
    b, c = d.scale, d.shape
    ((x / b)^(c - 1.0) * exp(-x / b)) / (b * gamma(c))
end

quantile(d::Erlang, p::Real) = quantile(d.nested_gamma, p)

function rand(d::Erlang)
    b, c = d.scale, d.shape
    z = 1.0
    for i in 1:c
        z *= rand()
    end
    -b * log(z)
end

skewness(d::Erlang) = skewness(d.nested_gamma)

var(d::Erlang) = d.scale^2 * d.shape

show(io::IO, d::Erlang) = show(io, d, (:shape, :scale))


text:
module Networks

# using ..AmplTxt

include("./generator.jl")
include("./bus.jl")
include("./branch.jl")
include("./network.jl")

function get_id(obj)
    return obj.id
end

export
    # struct
    Generator, Bus, Branch, Network,
    PTDFDict, PTDFValues,
    # enum
    GeneratorType,
    # functions
    ## bus
    add_new_bus!, add_new_buses!, add_bus!, add_buses!,
    get_bus, safeget_bus, get_buses,
    get_generators,
    ## branch
    add_new_branch!, add_new_branches!, add_branch!, add_branches!,
    get_branch, safeget_branch, get_branches,
    get_limit,
    ## load and generator
    add_new_generator_to_bus!,
    get_generator, safeget_generator, get_generators,
    get_bus_id, get_type, get_p_min, get_p_max, get_start_cost, get_prop_cost, get_dmo, get_dp,
    ## ptdf
    add_ptdf_elt!,
    get_ptdf_component, safeget_ptdf_component,
    get_ptdf, safeget_ptdf,
    ## utils
    get_id,
    safeget_generator_or_bus, get_generator_or_bus,
    get_info

end #module Networks

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using ClangCompiler
using Test

@testset "Lookup | Decl" begin
    src = joinpath(@__DIR__, "code", "main.cpp")
    args = get_compiler_args()
    cc = IncrementalIRGenerator(src, args)
    decl_lookup = DeclFinder(cc.instance)
    @test decl_lookup(cc, "vector", "std::vector")
    @test_nowarn decl = get_decl(decl_lookup)
    dispose(decl_lookup)
    dispose(cc)
end


text:
function sub_quote(W, @nospecialize(T), fast::Bool)
    vtyp = vtype(W, T)
    instrs = "%res = fneg $(fast_flags(fast)) $vtyp %0\nret $vtyp %res"
    quote
        $(Expr(:meta, :inline))
        Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T}}, data(v)))
    end
end

@generated vsub(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, T, false)
@generated vsub_fast(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, T, true)

@inline vsub(v::Vec{<:Any,<:NativeTypes}) = vsub(zero(v), v)
@inline vsub_fast(v::Vec{<:Any,<:NativeTypes}) = vsub_fast(zero(v), v)

@inline vinv(v::AbstractSIMD{W,<:FloatingTypes}) where {W} = vfdiv(one(v), v)
@inline vinv(v::AbstractSIMD{W,<:IntegerTypesHW}) where {W} = inv(float(v))

@inline vabs(v::AbstractSIMD{W,<:Unsigned}) where {W} = v
@inline vabs(v::AbstractSIMD{W,<:Signed}) where {W} = ifelse(v > 0, v, -v)

@inline vround(v::AbstractSIMD{W,<:Integer}) where {W} = v
@inline vround(v::AbstractSIMD{W,<:Integer}, ::RoundingMode) where {W} = v


text:
SingleView{T} = SubArray{T,0,Array{T,1},Tuple{Int64},true} where T

abstract type AbstractTuneable <: Number end

struct Tuneable{T<:Number} <: AbstractTuneable
    val::T
    Tuneable(val::T) where T = new{T}(val)
end

struct InitializedTuneable{T<:Number} <: AbstractTuneable
    view::SingleView{T}
end

initialize!(vect, untuneable) = untuneable
function initialize!(vect::Vector{T}, tuneable::Tuneable) where T
    push!(vect, tuneable.val)
    return InitializedTuneable{T}(@view vect[end])
end
function initialize!(vect::Vector, array::AbstractArray{T}) where T
    new_array = T[]
    for elem in array
        push!(new_array, initialize!(vect, elem))
    end
    return new_array
end
function initialize!(vect::Vector, tup::NamedTuple)
    array = Pair[]
    for (key, val) in zip(keys(tup), values(tup))
        push!(array, key => initialize!(vect, val))
    end
    return (;array...)
end

struct TuneableModel{M, A}
    model::M
    array::Vector{<:A}
end
function TuneableModel{A}(model) where {A}
    arr = A[]
    for pname in propertynames(model)
        prop = getproperty(model, pname)
        initialize!(arr, prop)
    end
end

max_stars_repo_name: jonniedie/Tuneables

max_issues_repo_name: jonniedie/Tuneables

max_forks_repo_name: jonniedie/Tuneables


text:
using Formatting
using MKL
using Base.LinAlg.BLAS
using Gadfly
using DataFrames
using Colors

include("../utility.jl")

"""
In this benchmark (b1), we provide a comparison between functions from
Julia, OpenBLAS and MKL that runs on a n-dim vector.
The comparison is also made between Float32 / Float64.
"""
# BENCH PARAMETERS
path = Pkg.dir("MKL") * "/benchmark/final/"
timing = .5
α = 0.05

N = [10_000] # ne pas prendre de valeur inférieure
T = [Float32,Float64]

function ibench(f, args, timing=0.001 ; α=0.05)

    # TRIGGER JIT : first compilation will be ignored
    res = @timed f(args...)

    # eval nrun to match timing
    t = @elapsed f(args...)
    nrun = max(2,convert(Int64,div(timing,t))) # at least one run

    # allocate results arrays
    res_cpu = Array(Float64,nrun+1)
    res_gc = Array(Float64,nrun+1)
    res_alloc = Array(Int64,nrun+1)

    # write JIT results
    res_cpu[1] = res[2]
    res_gc[1] = res[4]
    res_alloc[1] = res[3]

    # write PROFILING results : call f n times
    gc()
    for i=1:nrun
        # gcdisable ? gc_enable(false) : ""
        res = @timed f(args...)
        # gcdisable ? gc_enable(true) : ""
        res_cpu[i+1] = res[2]
        res_gc[i+1] = res[4]
        res_alloc[i+1] = res[3]
    end

    return res_cpu, res_gc, res_alloc
end

function write_csv(fname,df)
    head = transpose([string(name) for name in names(df)])
    data = convert(Array,df)
    data = vcat(head,data)
    writedlm(fname,data,",")
end

# GENERATE RAW DATA
function bench_b1()

    df = DataFrame( Generic = ASCIIString[], Implementation = ASCIIString[], Type = ASCIIString[], Method = Any[],
                    N = Int64[],
                    CPU_mean = Float64[], CPU_conf = Float64[],
                    GC_mean = Float64[], GC_conf = Float64[],
                    ALLOC_mean = Float64[], ALLOC_conf = Float64[])

    for n in N

        println("#################################################################")
        println("N = $n")
        println("#################################################################")

        for t in T

            λ = convert(t,pi) # keep it to 1.0 otherwise scal is going inf
            a = rand(t,n)
            b = rand(t,n)
            y = zeros(t,n)
            z = zeros(t,n)

        # for (t,λ,a,b,y,z) in ((Float32,λ32,a32,b32,y32,z32),(Float64,λ64,a64,b64,y64,z64))

            tup = (

            ("{\$ X.Y \$}",         "BLAS"  ,    t,  BLAS.dot   ,   (n,a,1,b,1) ),
            ("{\$ X.Y \$}",         "MKL"   ,    t,  mkl_dot    ,   (n,a,1,b,1) ),
            ("{\$ X.Y \$}",         "Julia" ,    t,  vj_dot     ,   (a,b)       ),

            ("{\$ \\lambda X \$}",    "BLAS"  ,    t,  BLAS.scal! ,   (n,convert(t,0.9999),a,1)   ),
            ("{\$ \\lambda X \$}",    "MKL"   ,    t,  mkl_scal!  ,   (n,convert(t,0.9999),a,1)   ),
            ("{\$ \\lambda X \$}",    "Julia" ,    t,  *          ,   (λ,a)       ),

            ("{\$ X+Y \$}",           "BLAS"  ,    t,  BLAS.axpy! ,   (1.0,a,b)   ),
            ("{\$ X+Y \$}",           "MKL"   ,    t,  mkl_add!   ,   (n,a,b,y)   ),
            ("{\$ X+Y \$}",           "Julia" ,    t,  +          ,   (a,b)       ),

            ("{\$ X-Y \$}",           "BLAS"  ,    t,  BLAS.axpy! ,   (-1.0,a,b)  ),
            ("{\$ X-Y \$}",           "MKL"   ,    t,  mkl_sub!   ,   (n,a,b,y)   ),
            ("{\$ X-Y \$}",           "Julia" ,    t,  -          ,   (a,b)       ),
            #
            ("{\$ \\lambda X+Y \$}",  "BLAS"  ,    t,  BLAS.axpy! ,   (λ,a,b)     ),
            ("{\$ \\lambda X+Y \$}",  "MKL"   ,    t,  mkl_axpy!  ,   (n,λ,a,1,b,1)),
            ("{\$ \\lambda X+Y \$}",  "Julia" ,    t,  vj_axpy    ,   (λ,a,b)     ),

            ("{\$ X*Y \$}",           "BLAS"  ,    t,  nanf       ,   ()          ),
            ("{\$ X*Y \$}",           "MKL"   ,    t,  mkl_mul!   ,   (n,a,b,y)   ),
            ("{\$ X*Y \$}",           "Julia" ,    t,  .*         ,   (a,b)       ),

            ("{\$ X^2 \$}",           "BLAS"  ,    t,  nanf       ,   ()                  ),
            ("{\$ X^2 \$}",           "MKL"   ,    t,  mkl_powx!  ,   (n,a,convert(t,2),y)),
            ("{\$ X^2 \$}",           "Julia" ,    t,  .^         ,   (a,convert(t,2))    ),

            ("{\$ X/Y \$}",           "BLAS"  ,    t,  nanf       ,   ()          ),
            ("{\$ X/Y \$}",           "MKL"   ,    t,  mkl_div!   ,   (n,a,b,y)   ),
            ("{\$ X/Y \$}",           "Julia" ,    t,  ./         ,   (a,b)       ),

            ("{\$ 1/X \$}",           "BLAS"  ,    t,  nanf       ,   ()              ),
            ("{\$ 1/X \$}",           "MKL"   ,    t,  mkl_inv!   ,   (n,a,y)         ),
            ("{\$ 1/X \$}",           "Julia" ,    t,  ./         ,   (convert(t,1),b)),

            ("{\$ |X| \$}",           "BLAS"  ,    t,  nanf       ,   ()         ),
            ("{\$ |X| \$}",           "MKL"   ,    t,  mkl_abs!   ,   (n,a,y)     ),
            ("{\$ |X| \$}",           "Julia" ,    t,  abs        ,   (a,)        ),

            ("{\$ \\sqrt{X} \$}",     "BLAS"  ,    t,  nanf        ,   ()         ),
            ("{\$ \\sqrt{X} \$}",     "MKL"   ,    t,  mkl_sqrt!   ,   (n,a,y)    ),
            ("{\$ \\sqrt{X} \$}",     "Julia" ,    t,  sqrt        ,   (a,)       ),

            ("{\$ 1/\\sqrt{X} \$}",   "BLAS"  ,    t,  nanf        ,   ()          ),
            ("{\$ 1/\\sqrt{X} \$}",   "MKL"   ,    t,  mkl_invsqrt!,   (n,a,y)     ),
            ("{\$ 1/\\sqrt{X} \$}",   "Julia" ,    t,  vj_invsqrt  ,   (a,)        ),
            )

            for (gm, imp, ty, f, args) in tup
                if length(args) == 0 # void function
                    push!(df, [gm, imp, string(ty), f, n, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
                else
                    res_cpu, res_gc, res_alloc = ibench(f, args,timing ; α=α)
                    stat_cpu = bstat(res_cpu[2:end].*(1e9/n), α) # in micro second per op * (1e9/n)
                    stat_gc = bstat(res_gc[2:end].*(1e9/n), α)
                    stat_alloc = bstat(res_alloc[2:end], α)
                    push!(df, [gm, imp, string(ty), f, n, stat_cpu[1], stat_cpu[5], stat_gc[1], stat_gc[5], stat_alloc[1], stat_alloc[5]])
                end
            end
        end
    end
    # writetable(path * "arithmetic_bench_data.csv", df)
    df
end

function plot_gadfly()

    df = DataFrames.readtable(path*"arithmetic_bench_data.csv")

    N = union(df[:N])
    T = union(df[:Type])
    imp = union(df[:Implementation])
    gen = union(df[:Generic])

    for n in N
        for t in T

            println(n)
            println(t)

            filter_set = (df[:Type] .== string(t)) & (df[:N] .== n)
            sdf = df[filter_set , :]
            println(sdf)

            set_default_plot_size(20cm, 20cm)
            p = Gadfly.plot(
                        sdf, xgroup="Generic", x="Implementation", y="CPU_mean",
                        color="Implementation",
                        Scale.color_discrete_manual(Tdarkgray,Tblue,Tlightgray),
                        Scale.y_continuous(minvalue=0, maxvalue=15),
                        Geom.subplot_grid(Geom.bar, Guide.xticks(orientation=:vertical, label=false)),
                        Guide.xlabel(""),
                        Guide.ylabel(""),
                        Guide.title("CPU time on n-dim vectors (n=1e"*string(convert(Int32,log10(n)))),
                        Theme(bar_spacing=0.1cm)
                        )

            draw(PNG(path*"b1_$t"*"_$n.png", 30cm, 20cm), p)
            # draw(SVG(path*"output/b1_$t"*"_$n.svg", 25cm, 12cm), p)
            p
        end
    end

    for n in N

        println(n)

        filter_set = (df[:N] .== n)
        sdf = df[filter_set , :]
        println(sdf)

        p = Gadfly.plot(
                    layer(
                    sdf, xgroup="Generic", ygroup="Type", x="Implementation", y="CPU_mean",
                    color="Implementation",
                    Geom.subplot_grid(Geom.bar, Guide.xticks(orientation=:vertical, label=false))
                    ),
                    Scale.color_discrete_manual(Tdarkgray,Tblue,Tlightgray),
                    Scale.y_continuous(minvalue=0, maxvalue=15),
                    Guide.xlabel(""),
                    Guide.ylabel(""),
                    Guide.title("CPU time on n-dim vectors (n=1e"*string(convert(Int32,log10(n)))*")"),
                    Theme(bar_spacing=0.1cm)
                    )

        # draw(PNG(path*"output/b1_$n.png", 25cm, 12cm), p)
        # draw(SVG(path*"b1_$n.svg", 29cm, 20cm), p)
        draw(PDF(path*"b1_$n.pdf", 30cm, 20cm), p)
        p
    end
end

df = bench_b1()
plot_gadfly()

write_csv(path*"arithmetic_bench_data.csv",df)
df
# .csv pour histogramme temps cpu absolu & interval de confiance
function histo_cpu_conf(df::DataFrame, n::Integer, T::DataType)
    df_type = df[(df[:Type] .== string(T)) & (df[:N] .== n),:]
    gens = union(df_type[:Generic])
    imps = union(df_type[:Implementation])
    res = DataFrame()
    res[:Generic] = gens
    for imp in imps
        df_imp = df_type[df_type[:Implementation] .== imp,:]
        res[symbol(string("CPU_mean_",imp))] = df_imp[:CPU_mean]
        res[symbol(string("CPU_conf_",imp))] = df_imp[:CPU_conf]
    end
    res
end

res = histo_cpu_conf(df, 10000, Float64)
write_csv(path * "arithmetic_cpu_conf_Float64.csv", res)
res = histo_cpu_conf(df, 10000, Float32)
write_csv(path * "arithmetic_cpu_conf_Float32.csv", res)

# .csv pour histogramme temps cpu relatif à MKL sans interval de confiance
function histo_cpu_rel(df::DataFrame, n::Integer, T::DataType)
    res = histo_cpu_conf(df, n, T)[[1,2,4,6]]
    res[2] = res[2] ./ res[3]
    res[3] = res[3] ./ res[3]
    res[4] = res[4] ./ res[3]
    res
end
res = histo_cpu_rel(df, 10000, Float64)
write_csv(path * "arithmetic_cpu_rel_Float64.csv", res)
res = histo_cpu_rel(df, 10000, Float32)
write_csv(path * "arithmetic_cpu_rel_Float32.csv", res)

# .csv pour histogramme temps cpu relatif Float64/Float32
function histo_cpu_type(df::DataFrame, n::Integer)
    df_32 = histo_cpu_conf(df, n, Float32)[[1,2,4,6]]
    df_64 = histo_cpu_conf(df, n, Float64)[[1,2,4,6]]
    df_64[2] = df_64[2] ./ df_32[2]
    df_64[3] = df_64[3] ./ df_32[3]
    df_64[4] = df_64[4] ./ df_32[4]
    # replace NaN values by zero
    for i in 1:nrow(df_64)
        for j in 2:ncol(df_64)
            if isnan(df_64[i,j])
                df_64[i,j]=0.0
            end
        end
    end
    df_64
end
res = histo_cpu_type(df, 10000)
write_csv(path * "arithmetic_cpu_type.csv", res)

# println(df_blas)


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

using Markdown
using InteractiveUtils

# ╔═╡ 8c1de468-b339-11eb-2c9a-fb5d7065bf78
html"""

Section 3.7

 Resistors and stencils

Lecture Video

"""

# ╔═╡ 89b4bb1b-0d49-4cf2-9013-3d320711577f
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="OffsetArrays", version="1")
    ])
    using PlutoUI, OffsetArrays
end

# ╔═╡ 7e939280-ccb3-4d64-8799-82630fbb7811
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 8b3da08b-d900-4736-955a-d25f8c7b70a9
md"""
Julia
* `CartesianIndex`
* `OffsetArray` from `OffsetArrays`
"""

# ╔═╡ 1945c286-c9a3-41f0-b3dc-e3a7c10334ab
md"""
#  Resistors, Equilibrium, and Poisson's equation. (17:48 minute video)
"""

# ╔═╡ 87ff0226-a249-419d-bd86-62331be6a538
md"""
This video is really about the interaction of the discrete and the continuous.
We set up a grid of resistor problem, and we show that solving this problem is the discretized version of solving a partial differential equation
known as Poisson's equation.

There are some references to fall 2020 which date this lecture, including
a reference to the Biden/Trump election, and a reference to John Urschel's lecture
which some of you might want to check out:

[John Urschel's video](https://youtu.be/rRCGNvMdLEY)
"""

# ╔═╡ a4d5fe96-5fed-4c26-b3ad-7637246cbb87
html"""

"""

# ╔═╡ ea8801ad-275e-4cb8-b383-e2d80fb920ec
md"""
# Stencils (first 2.5 minutes or so)

Don't miss the fun video effects!
"""

# ╔═╡ c9c62724-4c0e-4858-9419-828289ec1f45
html"""

"""

# ╔═╡ d0c01582-8499-4458-894c-cb23cf31094f
md"""
In the remainder of this notebook, we will show how to set up ghost cells in Julia so as to be able to apply a stencil.
"""

# ╔═╡ b1148b31-fb3b-434d-984f-dfc439c0e5c7
md"""
# Cartesian Indices
"""

# ╔═╡ fff37e11-d618-4271-b48c-12b7d1ef8472
md"""
Grab a 6x7 matrix that we'd like to apply a stencil to.
"""

# ╔═╡ afa4a77e-28fc-11eb-1ab0-bbba1b653e46
data = rand(1:9,6,7)

# ╔═╡ 17b7f5e6-83d5-43e0-8b56-93e67e3af747
md"""
A "Cartesian index" lets you access an element with one index variable:
"""

# ╔═╡ e970a572-12f6-4b72-8a60-5cd70e1b7260
i = CartesianIndex(2,3)

# ╔═╡ 2a0c666f-f3fa-49ab-8ffc-45323e7cba15
data[i]

# ╔═╡ af8c6b1e-5f72-4547-bd9f-1a1fc7264dc8
md"""
Obtain all the Cartesian Indices of your data matrix.
"""

# ╔═╡ 0019e726-28fd-11eb-0e86-31ec28b3c1a9
I = CartesianIndices(data)

# ╔═╡ e4e51762-7010-4afe-9599-3746acbe9143
md"""
# Offset Arrays
"""

# ╔═╡ 4d03b75f-e43a-484b-8020-43244b7700d5
md"""
An offset array lets you index in ways other than the standard 1:m, 1:n
"""

# ╔═╡ 5fb6e7b6-2901-11eb-0e94-aba290fd0bae
A = OffsetArray(zeros(Int,8,9), 0:7 ,0:8)

# ╔═╡ 87c260a2-2901-11eb-1060-b1e4b6b5b02b
for i ∈ I
	A[i] = data[i]  # copy data
end

# ╔═╡ a5c7693a-2901-11eb-1083-0da8138a73c2
A

# ╔═╡ f2903f3f-9697-4cac-af87-b2cfee362638
A[1,1]

# ╔═╡ 4fb21151-fc95-40e2-b2b7-7d0a05c5a60a
A[0,0]

# ╔═╡ 423f22c0-336a-4640-bbd2-2649e6021de6
md"""
# Neighborhood: a 3x3 window built from Cartesian Indices
"""

# ╔═╡ b6fde83c-2901-11eb-0e3b-4b3766579cc8
neighborhood = CartesianIndices((-1:1, -1:1))

# ╔═╡ 0eacc41c-89f7-4c11-b727-1769a6e7f5d5
md"""
Grab all the neighborhoods of `A`.
"""

# ╔═╡ babe3c24-2901-11eb-2d30-51256eb97e11
[ A[i.+neighborhood] for i ∈ I]

# ╔═╡ 7bd3671d-e59d-4d04-a60c-4524b2057972
md"""
# Stencil
"""

# ╔═╡ e6bd9dea-2901-11eb-1100-ad10705f41cc
stencil =  [ 0  -1   0
            -1   4  -1
            0  -1   0]

# ╔═╡ fe4f6df0-2901-11eb-1945-27e3f041ed1f
[  sum(A[i.+neighborhood].*stencil) for i ∈ I]

# ╔═╡ 48374720-6c79-4c2b-8b81-86565cbf19a2
md"""
Notice the result is the same size as the original data, and the stencil
"worked" on the edges.
"""

# ╔═╡ a7615570-0826-4ef1-80b2-da21c0c640b6
md"""
# Other boundary conditions.
We just saw 0 boundary conditions, what about periodic or zero derivative?
"""

# ╔═╡ 77c06ce6-2902-11eb-30a7-51f210dbd723
begin
 B = copy(A)

 B[0,:] = B[7,:]  ## periodic
 B[7,:] = B[1,:]
 B[:,0] = B[:,7]
 B[:,8] = B[:,1]

 #B[0,:] = B[1,:]  ## zero derivative
 #B[7,:] = B[7:]
 #B[:,0] = B[:,1]
 #B[:,8] = B[:,7]

 B
end

# ╔═╡ 4f342744-2902-11eb-1401-55e770d9d751

for i∈I
	B[i] = sum(A[i.+neighborhood].*stencil)
end

# ╔═╡ 6223e374-2902-11eb-3bb2-4d2d0d352801
B

# ╔═╡ Cell order:
# ╟─8c1de468-b339-11eb-2c9a-fb5d7065bf78
# ╠═89b4bb1b-0d49-4cf2-9013-3d320711577f
# ╠═7e939280-ccb3-4d64-8799-82630fbb7811
# ╟─8b3da08b-d900-4736-955a-d25f8c7b70a9
# ╟─1945c286-c9a3-41f0-b3dc-e3a7c10334ab
# ╟─87ff0226-a249-419d-bd86-62331be6a538
# ╟─a4d5fe96-5fed-4c26-b3ad-7637246cbb87
# ╟─ea8801ad-275e-4cb8-b383-e2d80fb920ec
# ╟─c9c62724-4c0e-4858-9419-828289ec1f45
# ╟─d0c01582-8499-4458-894c-cb23cf31094f
# ╟─b1148b31-fb3b-434d-984f-dfc439c0e5c7
# ╟─fff37e11-d618-4271-b48c-12b7d1ef8472
# ╠═afa4a77e-28fc-11eb-1ab0-bbba1b653e46
# ╟─17b7f5e6-83d5-43e0-8b56-93e67e3af747
# ╠═e970a572-12f6-4b72-8a60-5cd70e1b7260
# ╠═2a0c666f-f3fa-49ab-8ffc-45323e7cba15
# ╟─af8c6b1e-5f72-4547-bd9f-1a1fc7264dc8
# ╠═0019e726-28fd-11eb-0e86-31ec28b3c1a9
# ╟─e4e51762-7010-4afe-9599-3746acbe9143
# ╟─4d03b75f-e43a-484b-8020-43244b7700d5
# ╠═5fb6e7b6-2901-11eb-0e94-aba290fd0bae
# ╠═87c260a2-2901-11eb-1060-b1e4b6b5b02b
# ╠═a5c7693a-2901-11eb-1083-0da8138a73c2
# ╠═f2903f3f-9697-4cac-af87-b2cfee362638
# ╠═4fb21151-fc95-40e2-b2b7-7d0a05c5a60a
# ╟─423f22c0-336a-4640-bbd2-2649e6021de6
# ╠═b6fde83c-2901-11eb-0e3b-4b3766579cc8
# ╟─0eacc41c-89f7-4c11-b727-1769a6e7f5d5
# ╠═babe3c24-2901-11eb-2d30-51256eb97e11
# ╟─7bd3671d-e59d-4d04-a60c-4524b2057972
# ╠═e6bd9dea-2901-11eb-1100-ad10705f41cc
# ╠═fe4f6df0-2901-11eb-1945-27e3f041ed1f
# ╟─48374720-6c79-4c2b-8b81-86565cbf19a2
# ╟─a7615570-0826-4ef1-80b2-da21c0c640b6
# ╠═77c06ce6-2902-11eb-30a7-51f210dbd723
# ╠═4f342744-2902-11eb-1401-55e770d9d751
# ╠═6223e374-2902-11eb-3bb2-4d2d0d352801


text:
using PQEcolaPoint

# build a binary tree with three levels -
model = build(CRRLatticeModel; number_of_levels=84, σ=0.35, T=(28.0 / 365), μ=0.00174)

# compute -
a = price(model, 84; weights=true)


text:
"""
    NexNumber(nexSerial, name; pid=pid, cid=cid)

A Nextion Number UI control.

On Nextion side, it have the following properties:
- bco
- pco
- font
- xcen
- ycen
- val
- length
- format
- isbr
"""
struct NexNumber <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    numericalvalued::INumericalValued
    fontstyleable::IFontStyleable
    colourable::IColourable
    touchable::ITouchable

    function NexNumber(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), INumericalValued(nid, Int32, nothing), IFontStyleable(nid), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexNumber, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexNumber, property::Symbol)
    _getcommonproperty(obj, property)
end


text:
export Diagm

"C[:,:,k]=diagm(A[:,k])"
function Diagm(A::AFArray{Float32,2})
  C=zeros(AFArray{Float32},size(A,1),size(A,1),size(A,2))
  for k=1:1:size(A,2)
    @inbounds C[:,:,k]=diagm(A[:,k])
  end
  C
end


text:
function EulerIntegrator(arg0::jdouble)
    return EulerIntegrator((jdouble,), arg0)
end

function get_a(obj::EulerIntegrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::EulerIntegrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::EulerIntegrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end


text:
## -*-Julia-*-
## Test suite for Julia's WAV module
import WAV
using Test

# These float array comparison functions are from dists.jl
function absdiff(current::AbstractArray{T}, target::AbstractArray{T}) where T <: Real
    @test all(size(current) == size(target))
    maximum(abs.(current - target))
end

function reldiff(current::T, target::T) where T <: Real
    abs.((current - target)/(bool(target) ? target : 1))
end

function reldiff(current::AbstractArray{T}, target::AbstractArray{T}) where T <: Real
    @test all(size(current) == size(target))
    maximum([reldiff(current[i], target[i]) for i in 1:numel(target)])
end

## example from README, modified to use an IO buffer
@testset "1" begin
    x = [0:7999;]
    y = sin.(2 * pi * x / 8000)
    io = IOBuffer()
    WAV.wavwrite(y, io, Fs=8000)
    seek(io, 0)
    y, Fs = WAV.wavread(io)
    y = cos.(2 * pi * x / 8000)
    WAV.wavappend(y, io)
    seek(io, 0)
    y, Fs = WAV.wavread(io)
end

@testset "2" begin
    x = [0:7999;]
    y = sin.(2 * pi * x / 8000)
    WAV.wavwrite(y, "example.wav", Fs=8000)
    y, Fs = WAV.wavread("example.wav")
    y = cos.(2 * pi * x / 8000)
    WAV.wavappend(y, "example.wav")
    y, Fs = WAV.wavread("example.wav")
    @test length(y) == (2 * length(x))
    rm("example.wav")
end

## default arguments, GitHub Issue #10
@testset "3" begin
    tmp=rand(Float32,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; nbits=32)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @test typeof(y) == Array{Float32, 2}
    @test fs == 8000.0
    @test nbits == 32
    fmt = WAV.getformat(extra)
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
end

## malformed subchunk header, GitHub Issue #18
@testset "4" begin
    # Create a malformed WAV file
    samples = rand(Float32, (10, 1))
    io = IOBuffer()

    compression = WAV.get_default_compression(samples)
    nbits = WAV.get_default_precision(samples, compression)

    nchannels = size(samples, 2)
    sample_rate = 8000
    nbits = ceil(Integer, nbits / 8) * 8
    block_align = nbits / 8 * nchannels
    bps = sample_rate * block_align
    data_length::UInt32 = size(samples, 1) * block_align

    fmt = WAV.WAVFormat(compression,
                        nchannels,
                        sample_rate,
                        bps,
                        block_align,
                        nbits,
                        WAV.WAVFormatExtension())

    WAV.write_header(io, UInt32(data_length + 37)) # 37 instead of 36 is the broken part
    WAV.write_format(io, fmt)

    # write the data subchunk header
    WAV.write(io, b"data")
    WAV.write_le(io, data_length) # UInt32
    WAV.write_data(io, fmt, samples)

    seek(io, 0)
    y, fs, nbits, opt = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 32
    fmt = WAV.getformat(opt)
    @test fmt.compression_code == compression
    @test fmt.nchannels == nchannels
    @test fmt.sample_rate == sample_rate
    @test fmt.bytes_per_second == bps
    @test fmt.block_align == block_align
    @test WAV.isformat(fmt, compression)
    @test WAV.bits_per_sample(fmt) == nbits
    @test samples == y
end

@testset "5" begin
    tmp=rand(Float64,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; nbits=64)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @test typeof(y) == Array{Float64, 2}
    @test fs == 8000.0
    @test nbits == 64
    fmt = WAV.getformat(extra)
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_EXTENSIBLE)
end

@testset "6" begin
    tmp=rand(Float32,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; compression=WAV.WAVE_FORMAT_IEEE_FLOAT)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @test typeof(y) == Array{Float32, 2}
    @test fs == 8000.0
    @test nbits == 32
    @test WAV.isformat(WAV.getformat(extra), WAV.WAVE_FORMAT_IEEE_FLOAT)
end

function testread(io, ::Type{T}, sz) where T <: Real
    a = Array{T}(undef, sz)
    read!(io, a)
    return a
end

@testset "7"  begin
## Test wavread and wavwrite
## Generate some wav files for writing and reading
for fs = (8000,11025,22050,44100,48000,96000,192000), nbits = (1,7,8,9,12,16,20,24,32,64), nsamples = [0; 10 .^ (1:4)], nchans = 1:4
    ## Test wav files
    ## The tolerance is based on the number of bits used to encode the file in wavwrite
    tol = 2.0 / (2.0^(nbits - 1))

    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @test maximum(in_data) <= 1.0
        @test minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=WAV.WAVE_FORMAT_PCM)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @test testread(io, UInt8, 4) == b"RIFF"
    @test WAV.read_le(io, UInt32) == file_size - 8
    @test testread(io, UInt8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @test sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @test length(out_data) == nsamples * nchans
    @test size(out_data, 1) == nsamples
    @test size(out_data, 2) == nchans
    @test typeof(out_data) == Array{Float64, 2}
    @test out_fs == fs
    @test out_nbits == nbits
    fmt = WAV.getformat(out_extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
    @test fmt.nchannels == nchans
    if nsamples > 0
        @test absdiff(out_data, in_data) < tol
    end
    if nchans > 2
        @test WAV.isformat(fmt, WAV.WAVE_FORMAT_EXTENSIBLE)
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, trunc(Int, nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @test length(out_data) == subsamples * nchans
        @test size(out_data, 1) == subsamples
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == nbits
        fmt = WAV.getformat(out_extra)
        @test WAV.bits_per_sample(fmt) == nbits
        @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
        @test fmt.nchannels == nchans
        @test absdiff(out_data, in_data[1:subsamples, :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, trunc(Int, nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @test length(out_data) == length(sr) * nchans
        @test size(out_data, 1) == length(sr)
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == nbits
        fmt = WAV.getformat(out_extra)
        @test WAV.bits_per_sample(fmt) == nbits
        @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
        @test fmt.nchannels == nchans
        @test absdiff(out_data, in_data[sr, :]) < tol
    end
end
end

@testset "8" begin

## Test wavappend
x0=rand(2,3)
x1=rand(3,3)
x2=rand(4,3)
io = IOBuffer()
WAV.wavwrite(x0, io)
WAV.wavappend(x1, io)
WAV.wavappend(x2, io)
seek(io, 0)
x, fs, nbits, extra = WAV.wavread(io)
@test x == [x0; x1; x2]

## Test correct reading of 44-byte headers when using `wavappend()`
number_bits = 16
number_channels = 2
number_samples = 16000
sampling_frequency = 16000
x0 = rand(number_samples, number_channels)
x1 = rand(2*number_samples, number_channels)
io = IOBuffer()
### Write first chunk to wav-file
WAV.wavwrite(
    x0,
    io,
    Fs=sampling_frequency,
    nbits=number_bits,
    compression=WAV.WAVE_FORMAT_PCM)
seek(io, 4)
chunk_size_old = WAV.read_le(io, UInt32)
seek(io, 40)
data_length_old = WAV.read_le(io, UInt32)
### Append second chunk to wav-file
WAV.wavappend(
    x1,
    io)
seek(io, 4)
chunk_size_new = WAV.read_le(io, UInt32)
seek(io, 40)
data_length_new = WAV.read_le(io, UInt32)
### Compare data lengths.
data_length_old_in_samples =
    round(Int32, data_length_old/(number_channels*number_bits)*8)
data_length_new_in_samples =
    round(Int32, data_length_new/(number_channels*number_bits)*8)
@test data_length_new_in_samples == 3*data_length_old_in_samples
@test (chunk_size_new-36) == 3*(chunk_size_old-36)

## Test native encoding of 8 bits
for nchans = (1,2,4)
    in_data_8 = reshape(typemin(UInt8):typemax(UInt8), (trunc(Int, 256 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_8, io)

    seek(io, 0)
    out_data_8, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 8
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
    @test fmt.nchannels == nchans
    @test in_data_8 == out_data_8
end

## Test native encoding of 16 bits
for nchans = (1,2,4)
    in_data_16 = reshape(typemin(Int16):typemax(Int16), (trunc(Int, 65536 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_16, io)

    seek(io, 0)
    out_data_16, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 16
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
    @test fmt.nchannels == nchans
    @test in_data_16 == out_data_16
end

## Test native encoding of 24 bits
for nchans = (1,2,4)
    in_data_24 = convert(Array{Int32}, reshape(-63:64, trunc(Int, 128 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_24, io)

    seek(io, 0)
    out_data_24, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 24
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_PCM)
    @test fmt.nchannels == nchans
    @test in_data_24 == out_data_24
end

## Test encoding 32 bit values
for nchans = (1,2,4)
    in_data_single = convert(Array{Float32}, reshape(range(-1.0, stop=1.0, length=128), trunc(Int, 128 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 32
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test fmt.nchannels == nchans
    @test in_data_single == out_data_single
end

## Test encoding 32 bit values outside the [-1, 1] range
for nchans = (1,2,4)
    nsamps = trunc(Int, 128 / nchans)
    in_data_single = convert(Array{Float32}, reshape(-63:64, nsamps, nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 32
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test fmt.nchannels == nchans
    @test [in_data_single[i, j] for i = 1:nsamps, j = 1:nchans] == out_data_single
end

## Test encoding 64 bit values
for nchans = (1,2,4)
    in_data_single = convert(Array{Float64}, reshape(range(-1.0, stop=1.0, length=128), trunc(Int, 128 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 64
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test fmt.nchannels == nchans
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test in_data_single == out_data_single
end

## Test encoding 64 bit values outside the [-1, 1] range
for nchans = (1,2,4)
    nsamps = trunc(Int, 128 / nchans)
    in_data_single = convert(Array{Float64}, reshape(-63:64, nsamps, nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @test fs == 8000
    @test nbits == 64
    fmt = WAV.getformat(extra)
    @test WAV.bits_per_sample(fmt) == nbits
    @test WAV.isformat(fmt, WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test fmt.nchannels == nchans
    @test [in_data_single[i, j] for i = 1:nsamps, j = 1:nchans] == out_data_single
end

### Test A-Law and Mu-Law
for nbits = (8, 16), nsamples = [0; 10 .^ (1:4)], nchans = 1:4, fmt=(WAV.WAVE_FORMAT_ALAW, WAV.WAVE_FORMAT_MULAW)
    fs = 8000.0
    tol = 2.0 / (2.0^6)
    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @test maximum(in_data) <= 1.0
        @test minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=fmt)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @test testread(io, UInt8, 4) == b"RIFF"
    @test WAV.read_le(io, UInt32) == file_size - 8
    @test testread(io, UInt8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @test sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @test length(out_data) == nsamples * nchans
    @test size(out_data, 1) == nsamples
    @test size(out_data, 2) == nchans
    @test typeof(out_data) == Array{Float64, 2}
    @test out_fs == fs
    @test out_nbits == 8
    out_fmt = WAV.getformat(out_extra)
    @test WAV.bits_per_sample(out_fmt) == 8
    @test out_fmt.nchannels == nchans
    @test WAV.isformat(out_fmt, fmt)
    if nchans > 2
        @test WAV.isformat(out_fmt, WAV.WAVE_FORMAT_EXTENSIBLE)
    else
        @test !WAV.isformat(out_fmt, WAV.WAVE_FORMAT_EXTENSIBLE)
    end
    if nsamples > 0
        @test absdiff(out_data, in_data) < tol
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, trunc(Int, nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @test length(out_data) == subsamples * nchans
        @test size(out_data, 1) == subsamples
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == 8
        out_fmt = WAV.getformat(out_extra)
        @test WAV.bits_per_sample(out_fmt) == 8
        @test out_fmt.nchannels == nchans
        @test WAV.isformat(out_fmt, fmt)
        @test absdiff(out_data, in_data[1:subsamples, :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, trunc(Int, nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @test length(out_data) == length(sr) * nchans
        @test size(out_data, 1) == length(sr)
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == 8
        out_fmt = WAV.getformat(out_extra)
        @test WAV.bits_per_sample(out_fmt) == 8
        @test out_fmt.nchannels == nchans
        @test WAV.isformat(out_fmt, fmt)
        @test absdiff(out_data, in_data[sr, :]) < tol
    end
end

### Test float formatting
for nbits = (32, 64), nsamples = [0; 10 .^ (1:4)], nchans = 1:2, fmt=(WAV.WAVE_FORMAT_IEEE_FLOAT)
    fs = 8000.0
    tol = 1e-6
    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @test maximum(in_data) <= 1.0
        @test minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=fmt)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @test testread(io, UInt8, 4) == b"RIFF"
    @test WAV.read_le(io, UInt32) == file_size - 8
    @test testread(io, UInt8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @test sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @test length(out_data) == nsamples * nchans
    @test size(out_data, 1) == nsamples
    @test size(out_data, 2) == nchans
    @test typeof(out_data) == Array{Float64, 2}
    @test out_fs == fs
    @test out_nbits == nbits
    @test WAV.bits_per_sample(WAV.getformat(out_extra)) == nbits
    @test WAV.isformat(WAV.getformat(out_extra), WAV.WAVE_FORMAT_IEEE_FLOAT)
    @test WAV.getformat(out_extra).nchannels == nchans
    if nsamples > 0
        @test absdiff(out_data, in_data) < tol
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, trunc(Int, nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @test length(out_data) == subsamples * nchans
        @test size(out_data, 1) == subsamples
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == nbits
        @test WAV.bits_per_sample(WAV.getformat(out_extra)) == nbits
        @test WAV.getformat(out_extra).nchannels == nchans
        @test WAV.isformat(WAV.getformat(out_extra), WAV.WAVE_FORMAT_IEEE_FLOAT)
        @test absdiff(out_data, in_data[1:subsamples, :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, trunc(Int, nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @test length(out_data) == length(sr) * nchans
        @test size(out_data, 1) == length(sr)
        @test size(out_data, 2) == nchans
        @test typeof(out_data) == Array{Float64, 2}
        @test out_fs == fs
        @test out_nbits == nbits
        @test WAV.bits_per_sample(WAV.getformat(out_extra)) == nbits
        @test WAV.isformat(WAV.getformat(out_extra), WAV.WAVE_FORMAT_IEEE_FLOAT)
        @test WAV.getformat(out_extra).nchannels == nchans
        @test absdiff(out_data, in_data[sr, :]) < tol
    end
end
end

### Read unknown chunks, GitHub Issue #50
@testset "9" begin
    fs = 8000.0
    in_data = rand(1024, 2)
    io = IOBuffer()
    in_chunks = [WAV.WAVChunk(:test, [0x1, 0x2, 0x3])]
    WAV.wavwrite(in_data, io, Fs=fs, chunks=in_chunks)

    seek(io, 0)
    data, fs, nbits, ext = WAV.wavread(io)

    @test findfirst(c -> c.id == :test, ext) > 0
    test_chunk = ext[findfirst(c -> c.id == :test, ext)]
    @test test_chunk.id == in_chunks[1].id
    @test test_chunk.data == in_chunks[1].data
    @test length(data) == length(in_data)
    @test isapprox(data, in_data; atol=1.0e-6)
end

### Multiple LIST chunks, GitHub Issue #55
@testset "10" begin
    fs = 8000.0
    in_data = rand(1024, 2)
    io = IOBuffer()
    in_chunks = [WAV.WAVChunk(:LIST, [0x1, 0x2, 0x3]), WAV.WAVChunk(:LIST, [0x4, 0x5, 0x6])]
    WAV.wavwrite(in_data, io, Fs=fs, chunks=in_chunks)

    seek(io, 0)
    data, fs, nbits, ext = WAV.wavread(io)

    @test length(findall(c -> c.id == :LIST, ext)) == length(in_chunks)
    list_chunks = ext[findall(c -> c.id == :LIST, ext)]
    for (c, i) in zip(list_chunks, in_chunks)
        @test c.id == i.id
        @test c.data == i.data
    end
    @test length(data) == length(in_data)
    @test isapprox(data, in_data; atol=1.0e-6)
end

# Test WAVMarker and WavChunk
@testset "11" begin
    io = IOBuffer()
    fs = 16000
    samples = rand(fs)
    markers = Dict{UInt32, WAV.WAVMarker}()
    markers[1] = WAV.WAVMarker("Foo", 42, 1337)
    markers[2] = WAV.WAVMarker("Bar", 1337, 42)

    marker_chunks = WAV.wav_cue_write(markers)

    title = "Never Gonna Give You Up"
    artist = "Rick Astley"
    tags = Dict{Symbol, String}()
    tags[:INAM] = title
    tags[:IART] = artist
    tag_chunks = WAV.wav_info_write(tags)

    out_chunks = [tag_chunks; marker_chunks]

    WAV.wavwrite(samples, io, Fs=fs, nbits=16, compression=WAV.WAVE_FORMAT_PCM, chunks=out_chunks)
    seek(io, 0)

    x, fs, bits, in_chunks = WAV.wavread(io)
    in_markers = WAV.wav_cue_read(in_chunks)

    @test length(in_markers) == length(markers)
    for k in keys(in_markers)
        @test in_markers[k].label == markers[k].label
        @test in_markers[k].start_time == markers[k].start_time
        @test in_markers[k].duration == markers[k].duration
    end

    in_info = WAV.wav_info_read(in_chunks)
    @test in_info[:INAM] == title
    @test in_info[:IART] == artist
end

### WAVArray
struct TestHtmlDisplay <: AbstractDisplay
    io::IOBuffer
end
function display(d::TestHtmlDisplay, mime::MIME"text/html", x)
    print(d.io, repr(mime, x))
end

@testset "12" begin
    wa = WAV.WAVArray(8000, sin.(1:256 * 8000.0 / 1024));
    myio = IOBuffer()
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    wa = WAV.WAVArray(8000, sin.(1:256 * 8000.0 / 1024), "WAVArray caption test");
    seek(myio, 0)
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    x = [0:7999;]
    y = sin.(2 * pi * x / 8000)
    io = IOBuffer()
    WAV.wavwrite(y, io, Fs=8000)
    seek(io,0)
    wa = WAV.WAVArray(io)
    myio = IOBuffer()
    seek(myio, 0)
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    seek(io,0)
    wa = WAV.WAVArray(io, "WAVArray caption test")
    seek(myio, 0)
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    WAV.wavwrite(y, "example.wav", Fs=8000)
    wa = WAV.WAVArray("example.wav")
    seek(myio, 0)
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    wa = WAV.WAVArray("example.wav", "WAVArray caption test")
    seek(myio, 0)
    display(TestHtmlDisplay(myio), MIME"text/html"(), wa)
    @test occursin(r"audio controls", String(take!(copy(myio))))

    rm("example.wav")
end

### playback
if !haskey(ENV, "CI")
    @testset "13" begin
        fs = 44100.0
        t = 0.0:fs-1.0;
        left_data  = sin.(2pi * 500.0 * t / fs) * 1e-1;
        right_data = sin.(2pi * 800.0 * t / fs) * 1e-1;
        # 500 Hz mono
        @test WAV.wavplay(left_data, fs) === nothing;
        # 500 Hz left, 800 Hz right stereo
        @test WAV.wavplay([left_data right_data], fs) === nothing;
    end
end


text:
include("OpCodes2.jl")
using .OpCodes2

data = parse.(Int, split(readline("Dec2_data.txt"), ','))
data[2:3] .= (12, 2)
res = compute(data)[1]

# Question 2

function fun2(input)
    for i in 0:99, j in 0:99
        p = copy(input)
        p[2:3] .= (i, j)
        p = compute(p)
        p[1] == 19690720 && return(100i + j)
    end
    error("no input gave the right result")
end
fun2(data)


text:
# ============================================================
# Initialize Simulation
# ============================================================

"""
Initialize solver

"""
function initialize(configfilename::AbstractString, structure = "sol"::AbstractString)
    println("==============================================================")
    println("SKS.jl: A Software Package for Stochastic Kinetic Simulation")
    println("==============================================================")
    println("")
    println("reading configurations from $configfilename")
    println("")
    println("initializeing solver:")

    allowed = ["uqMethod", "nr", "nRec", "opType", "parameter1", "parameter2"]
    D = read_dict(configfilename, allowed)

    nr = D["nr"]
    nRec = D["nRec"]
    opType = D["opType"]
    parameter1 = D["parameter1"]
    parameter2 = D["parameter2"]
    uqMethod = D["uqMethod"]

    ks = SolverSet(configfilename)
    uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)

    if structure == "ctr"
        ctr, face = init_fvm(ks, uq)
        return ks, ctr, face, uq, 0.0
    elseif structure == "sol"
        sol, flux = init_sol(ks, uq)
        return ks, sol, flux, uq, 0.0
    else
        throw(ArgumentError("no data structure available for $structure"))
    end
end

"""
Initialize collocation solution structs

"""
function init_collo_sol(KS::AbstractSolverSet, uq::AbstractUQ)
    if KS.set.space[1:2] == "1d"
        w0 = zeros(axes(KS.ib.fw(KS.ps.x0, KS.ib.p), 1), uq.op.quad.Nquad)
        prim0 = uq_conserve_prim(w0, KS.gas.γ, uq)
        w = [deepcopy(w0) for i in axes(KS.ps.x, 1)]
        prim = deepcopy(w)
        facefw = [zeros(axes(w0)) for i = 1:KS.ps.nx+1]

        for i in axes(w, 1)
            for j in axes(w0, 2)
                w0[:, j] .= KS.ib.fw(KS.ps.x[i], KS.ib.p)
            end

            w[i] .= w0
            prim[i] .= uq_conserve_prim(w0, KS.gas.γ, uq)
        end

        if KS.set.space[3:4] == "1f"
            if KS.set.space[5:6] == "1v"
                f0 = uq_maxwellian(KS.vs.u, prim0, uq)
            elseif KS.set.space[5:6] == "3v"
                f0 = uq_maxwellian(KS.vs.u, KS.vs.v, KS.vs.w, prim0, uq)
            end

            f = [deepcopy(f0) for i in axes(KS.ps.x, 1)]
            for i in axes(w, 1)
                f[i] .= uq_maxwellian(KS.vs.u, prim[i], uq)
            end
            faceff = [zeros(axes(f0)) for i = 1:KS.ps.nx+1]

            sol = Solution1D(w, prim, f)
            flux = Flux1D(facefw, faceff)

            return sol, flux
        elseif KS.set.space[3:4] == "2f"
            if KS.set.space[5:6] == "1v"
                h0 = uq_maxwellian(KS.vs.u, prim0, uq)
                b0 = deepcopy(h0)
                for j in axes(b0, 2)
                    b0[:, j] .= h0[:, j] .* KS.gas.K ./ (2.0 * prim0[end, j])
                end
            elseif KS.set.space[5:6] == "2v"
                h0 = uq_maxwellian(KS.vs.u, KS.vs.v, prim0, uq)
                b0 = deepcopy(h0)
                for j in axes(b0, 3)
                    b0[:, :, j] .= h0[:, :, j] .* KS.gas.K ./ (2.0 * prim0[end, j])
                end
            end

            h = [deepcopy(h0) for i in axes(KS.ps.x, 1)]
            b = [deepcopy(b0) for i in axes(KS.ps.x, 1)]
            if KS.set.space[5:6] == "1v"
                for i in axes(w, 1)
                    h[i] .= uq_maxwellian(KS.vs.u, prim[i], uq)
                    b[i] .= uq_energy_distribution(h[i], prim[i], KS.gas.K, uq)
                end
            elseif KS.set.space[5:6] == "2v"
                for i in axes(w, 1)
                    h[i] .= uq_maxwellian(KS.vs.u, KS.vs.v, prim[i], uq)
                    b[i] .= uq_energy_distribution(h[i], prim[i], KS.gas.K, uq)
                end
            end

            facefh = [zeros(axes(h0)) for i = 1:KS.pSpace.nx+1]
            facefb = [zeros(axes(b0)) for i = 1:KS.pSpace.nx+1]

            sol = Solution1D(w, prim, h, b)
            flux = Flux1D(facefw, facefh, facefb)

            return sol, flux
        end
    end
end

"""
Initialize solution structs

"""
function init_sol(KS::AbstractSolverSet, uq::AbstractUQ)
    if uq.method == "galerkin"
        # upstream
        primL = zeros(axes(ks.ib.bc(ks.ps.x0, ks.ib.p), 1), uq.nr + 1)
        primL[:, 1] .= ks.ib.bc(ks.ps.x0, ks.ib.p)

        # downstream
        primR = zeros(axes(ks.ib.bc(ks.ps.x1, ks.ib.p), 1), uq.nr + 1)
        primR[:, 1] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
    elseif uq.method == "collocation"
        # upstream
        primL = zeros(axes(ks.ib.bc(ks.ps.x0, ks.ib.p), 1), uq.op.quad.Nquad)
        for j in axes(primL, 2)
            primL[:, j] .= ks.ib.bc(ks.ps.x0, ks.ib.p)
        end

        # downstream
        primR = zeros(axes(ks.ib.bc(ks.ps.x1, ks.ib.p), 1), uq.op.quad.Nquad)
        for j in axes(primR, 2)
            primR[:, j] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
        end
    end

    wL = uq_prim_conserve(primL, KS.gas.γ, uq)
    wR = uq_prim_conserve(primR, KS.gas.γ, uq)
    facew = [zeros(axes(wL)) for i = 1:KS.pSpace.nx+1]
    facefw = [zeros(axes(wL)) for i = 1:KS.pSpace.nx+1]

    if KS.set.space[1:2] == "1d"

        w = [deepcopy(wL) for i in axes(KS.pSpace.x, 1)]
        prim = [deepcopy(primL) for i in axes(KS.pSpace.x, 1)]
        for i in axes(w, 1)
            if i > KS.pSpace.nx ÷ 2
                w[i] .= deepcopy(wR)
                prim[i] .= deepcopy(primR)
            end
        end

        if KS.set.space[3:4] == "1f"

            if KS.set.space[5:6] == "1v"
                fL = uq_maxwellian(KS.vSpace.u, primL, uq)
                fR = uq_maxwellian(KS.vSpace.u, primR, uq)
            elseif KS.set.space[5:6] == "3v"
                fL = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, KS.vSpace.w, primL, uq)
                fR = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, KS.vSpace.w, primR, uq)
            end

            f = [deepcopy(fL) for i in axes(KS.pSpace.x, 1)]
            for i in axes(w, 1)
                if i > KS.pSpace.nx ÷ 2
                    f[i] .= deepcopy(fR)
                end
            end
            faceff = [zeros(axes(fL)) for i = 1:KS.pSpace.nx+1]

            sol = Solution1D1F(w, prim, f)
            flux = Flux1D1F(facew, facefw, faceff)

            return sol, flux

        elseif KS.set.space[3:4] == "2f"

            if KS.set.space[5:6] == "1v"
                hL = uq_maxwellian(KS.vSpace.u, primL, uq)
                hR = uq_maxwellian(KS.vSpace.u, primR, uq)
                bL = deepcopy(hL)
                bR = deepcopy(hR)
                for j in axes(bL, 2)
                    bL[:, j] .= hL[:, j] .* KS.gas.K ./ (2.0 * primL[end, j])
                    bR[:, j] .= hR[:, j] .* KS.gas.K ./ (2.0 * primR[end, j])
                end
            elseif KS.set.space[5:6] == "2v"
                hL = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primL, uq)
                hR = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primR, uq)
                bL = deepcopy(hL)
                bR = deepcopy(hR)
                for j in axes(bL, 3)
                    bL[:, :, j] .= hL[:, :, j] .* KS.gas.K ./ (2.0 * primL[end, j])
                    bR[:, :, j] .= hR[:, :, j] .* KS.gas.K ./ (2.0 * primR[end, j])
                end
            end

            h = [deepcopy(hL) for i in axes(KS.pSpace.x, 1)]
            b = [deepcopy(bL) for i in axes(KS.pSpace.x, 1)]
            for i in axes(w, 1)
                if i > KS.pSpace.nx ÷ 2
                    h[i] .= deepcopy(hR)
                    b[i] .= deepcopy(bR)
                end
            end
            facefh = [zeros(axes(hL)) for i = 1:KS.pSpace.nx+1]
            facefb = [zeros(axes(bL)) for i = 1:KS.pSpace.nx+1]

            sol = Solution1D2F(w, prim, h, b)
            flux = Flux1D2F(facew, facefw, facefh, facefb)

            return sol, flux

        end

    elseif KS.set.space[1:2] == "2d"

        #--- cell ---#
        w = [deepcopy(wL) for i in axes(KS.pSpace.x, 1), j in axes(KS.pSpace.x, 2)]
        prim = [deepcopy(primL) for i in axes(KS.pSpace.x, 1), j in axes(KS.pSpace.x, 2)]
        for j in axes(w, 2), i in axes(w, 1)
            if j > KS.pSpace.ny ÷ 2
                w[i, j] .= deepcopy(wR)
                prim[i, j] .= deepcopy(primR)
            end
        end

        #--- interface ---#
        n1 = [[1.0, 0.0] for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
        n2 = [[0.0, 1.0] for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]
        facew1 = [zeros(axes(wL)) for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
        facefw1 = [zeros(axes(wL)) for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
        facew2 = [zeros(axes(wL)) for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]
        facefw2 = [zeros(axes(wL)) for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]

        if KS.set.space[3:4] == "1f"

            #--- cell ---#
            fL = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primL, uq)
            fR = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primR, uq)

            f = [deepcopy(fL) for i in axes(KS.pSpace.x, 1), j in axes(KS.pSpace.x, 2)]
            for j in axes(w, 2), i in axes(w, 1)
                if j > KS.pSpace.ny ÷ 2
                    f[i, j] .= deepcopy(fR)
                end
            end

            #--- interface ---#
            faceff1 = [zeros(axes(fL)) for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
            faceff2 = [zeros(axes(fL)) for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]

            sol = Solution2D1F(w, prim, f)
            flux = Flux2D1F(n1, facew1, facefw1, faceff1, n2, facew2, facefw2, faceff2)

        elseif KS.set.space[3:4] == "2f"

            #--- cell ---#
            hL, bL = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primL, uq, KS.gas.K)
            hR, bR = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primR, uq, KS.gas.K)

            h = [deepcopy(hL) for i in axes(KS.pSpace.x, 1), j in axes(KS.pSpace.x, 2)]
            b = [deepcopy(bL) for i in axes(KS.pSpace.x, 1), j in axes(KS.pSpace.x, 2)]

            for j in axes(w, 2), i in axes(w, 1)
                if j > KS.pSpace.ny ÷ 2
                    w[i, j] .= deepcopy(wR)
                    prim[i, j] .= deepcopy(primR)
                    h[i, j] .= deepcopy(hR)
                    b[i, j] .= deepcopy(bR)
                end
            end

            #--- interface ---#
            facefh1 = [zeros(axes(hL)) for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
            facefb1 = [zeros(axes(bL)) for i = 1:KS.pSpace.nx+1, j = 1:KS.pSpace.ny]
            facefh2 = [zeros(axes(hL)) for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]
            facefb2 = [zeros(axes(bL)) for i = 1:KS.pSpace.nx, j = 1:KS.pSpace.ny+1]

            sol = Solution2D2F(w, prim, h, b)
            flux = Flux2D2F(
                n1,
                facew1,
                facefw1,
                facefh1,
                facefb1,
                n2,
                facew2,
                facefw2,
                facefh2,
                facefb2,
            )

        end

        return sol, flux

    end
end

"""
Initialize finite volume structs

"""
function init_fvm(KS::AbstractSolverSet, uq::AbstractUQ)

    if ndims(KS.ib.wL) == 1

        ctr, face = pure_fvm(KS, KS.pSpace, uq)

    elseif ndims(KS.ib.wL) == 2

        if KS.set.space[3:4] in ("3f", "4f")
            ctr, face = plasma_fvm(KS, KS.pSpace, uq)
        else
            ctr, face = mixture_fvm(KS, KS.pSpace, uq)
        end

    end

    return ctr, face

end

function pure_fvm(KS::AbstractSolverSet, pSpace::PSpace1D, uq::AbstractUQ)

    #--- setup of control volume ---#
    idx0 = (eachindex(pSpace.x)|>collect)[1]
    idx1 = (eachindex(pSpace.x)|>collect)[end]

    if KS.set.space[3:4] == "1f"
        ctr = OffsetArray{ControlVolume1D1F}(undef, idx0:idx1) # with ghost cells
    elseif KS.set.space[3:4] == "2f"
        ctr = OffsetArray{ControlVolume1D2F}(undef, idx0:idx1)
    elseif KS.set.space[3:4] == "4f"
        ctr = OffsetArray{ControlVolume1D4F}(undef, idx0:idx1)
    end

    if uq.method == "galerkin"

        # upstream
        primL = zeros(axes(ks.ib.bc(ks.ps.x0, ks.ib.p), 1), uq.nr + 1)
        primL[:, 1] .= ks.ib.bc(ks.ps.x0, ks.ib.p)
        wL = uq_prim_conserve(primL, KS.gas.γ, uq)
        fL = uq_maxwellian(KS.vSpace.u, primL, uq)

        # downstream
        primR = zeros(axes(ks.ib.bc(ks.ps.x1, ks.ib.p), 1), uq.nr + 1)
        primR[:, 1] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
        wR = uq_prim_conserve(primR, KS.gas.γ, uq)
        fR = uq_maxwellian(KS.vSpace.u, primR, uq)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D1F(KS.pSpace.x[i], KS.pSpace.dx[i], wL, primL, fL)
            else
                ctr[i] = ControlVolume1D1F(KS.pSpace.x[i], KS.pSpace.dx[i], wR, primR, fR)
            end
        end

        # --- setup of cell interface ---#
        face = Array{Interface1D1F}(undef, KS.pSpace.nx + 1)
        for i in eachindex(face)
            face[i] = Interface1D1F(wL, fL)
        end

    elseif uq.method == "collocation"

        # upstream
        primL = zeros(axes(ks.ib.bc(ks.ps.x0, ks.ib.p), 1), uq.op.quad.Nquad)
        for j in axes(primL, 2)
            primL[:, j] .= ks.ib.bc(ks.ps.x0, ks.ib.p)
        end
        wL = uq_prim_conserve(primL, KS.gas.γ, uq)
        fL = uq_maxwellian(KS.vSpace.u, primL, uq)

        # downstream
        primR = zeros(axes(ks.ib.bc(ks.ps.x0, ks.ib.p), 1), uq.op.quad.Nquad)
        for j in axes(primL, 2)
            primR[:, j] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
        end
        wR = uq_prim_conserve(primR, KS.gas.γ, uq)
        fR = uq_maxwellian(KS.vSpace.u, primR, uq)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D1F(KS.pSpace.x[i], KS.pSpace.dx[i], wL, primL, fL)
            else
                ctr[i] = ControlVolume1D1F(KS.pSpace.x[i], KS.pSpace.dx[i], wR, primR, fR)
            end
        end

        # --- setup of cell interface ---#
        face = Array{Interface1D1F}(undef, KS.pSpace.nx + 1)
        for i in eachindex(face)
            face[i] = Interface1D1F(wL, fL)
        end

    end

    return ctr, face

end

function mixture_fvm(KS::AbstractSolverSet, pSpace::PSpace1D, uq::AbstractUQ)

    idx0 = (eachindex(pSpace.x)|>collect)[1]
    idx1 = (eachindex(pSpace.x)|>collect)[end]

    ctr = OffsetArray{ControlVolume1D4F}(undef, idx0:idx1)

    if uq.uqMethod == "galerkin"

        # upstream
        primL = zeros(5, uq.nr + 1, 2)
        primL[:, 1, :] .= ks.ib.bc(ks.ps.x0, ks.ib.p)

        wL = get_conserved(primL, KS.gas.γ, uq)
        h0L, h1L, h2L, h3L = get_maxwell(KS.vSpace.u, primL, uq)

        EL = zeros(3, uq.nr + 1)
        BL = zeros(3, uq.nr + 1)
        BL[:, 1] .= KS.ib.BL

        # downstream
        primR = zeros(5, uq.nr + 1, 2)
        primR[:, 1, :] .= ks.ib.bc(ks.ps.x1, ks.ib.p)

        wR = get_conserved(primR, KS.gas.γ, uq)
        h0R, h1R, h2R, h3R = get_maxwell(KS.vSpace.u, primR, uq)

        ER = zeros(3, uq.nr + 1)
        BR = zeros(3, uq.nr + 1)
        BR[:, 1] .= KS.ib.BR

        lorenz = zeros(3, uq.nr + 1, 2)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wL,
                    primL,
                    h0L,
                    h1L,
                    h2L,
                    h3L,
                    EL,
                    BL,
                    lorenz,
                )
            else
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wR,
                    primR,
                    h0R,
                    h1R,
                    h2R,
                    h3R,
                    ER,
                    BR,
                    lorenz,
                )
            end
        end

        face = Array{Interface1D4F}(undef, KS.pSpace.nx + 1)
        for i = 1:KS.pSpace.nx+1
            face[i] = Interface1D4F(wL, h0L, EL)
        end

    elseif uq.uqMethod == "collocation"

        # upstream
        primL = zeros(5, uq.op.quad.Nquad, 2)
        for j in axes(primL, 2)
            primL[:, j, :] .= ks.ib.bc(ks.ps.x0, ks.ib.p)
        end

        wL = get_conserved(primL, KS.gas.γ)
        h0L, h1L, h2L, h3L = get_maxwell(KS.vSpace.u, primL)

        EL = zeros(3, uq.op.quad.Nquad)
        BL = zeros(3, uq.op.quad.Nquad)
        for j in axes(BL, 2)
            BL[:, j] .= KS.ib.BL
        end

        # downstream
        primR = zeros(5, uq.op.quad.Nquad, 2)
        for j in axes(primL, 2)
            primR[:, j, :] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
        end

        wR = get_conserved(primR, KS.gas.γ)
        h0R, h1R, h2R, h3R = get_maxwell(KS.vSpace.u, primR)

        ER = zeros(3, uq.op.quad.Nquad)
        BR = zeros(3, uq.op.quad.Nquad)
        for j in axes(BR, 2)
            BR[:, j] .= KS.ib.BR
        end

        lorenz = zeros(3, uq.op.quad.Nquad, 2)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wL,
                    primL,
                    h0L,
                    h1L,
                    h2L,
                    h3L,
                    EL,
                    BL,
                    lorenz,
                )
            else
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wR,
                    primR,
                    h0R,
                    h1R,
                    h2R,
                    h3R,
                    ER,
                    BR,
                    lorenz,
                )
            end
        end

        #--- setup of cell interface ---#
        face = Array{Interface1D4F}(undef, KS.pSpace.nx + 1)
        for i = 1:KS.pSpace.nx+1
            face[i] = Interface1D4F(wL, h0L, EL)
        end

    end

    return ctr, face

end

function plasma_fvm(KS::AbstractSolverSet, pSpace::PSpace1D, uq::AbstractUQ)

    idx0 = (eachindex(pSpace.x)|>collect)[1]
    idx1 = (eachindex(pSpace.x)|>collect)[end]

    ctr = OffsetArray{ControlVolume1D4F}(undef, idx0:idx1)

    if uq.method == "galerkin"

        # upstream
        primL = zeros(5, uq.nr + 1, 2)
        primL[:, 1, :] .= ks.ib.bc(ks.ps.x0, ks.ib.p)

        wL = uq_prim_conserve(primL, KS.gas.γ, uq)
        h0L, h1L, h2L, h3L = uq_maxwellian(KS.vSpace.u, primL, uq)

        EL = zeros(3, uq.nr + 1)
        BL = zeros(3, uq.nr + 1)
        BL[:, 1] .= KS.ib.BL

        # downstream
        primR = zeros(5, uq.nr + 1, 2)
        primR[:, 1, :] .= ks.ib.bc(ks.ps.x1, ks.ib.p)

        wR = uq_prim_conserve(primR, KS.gas.γ, uq)
        h0R, h1R, h2R, h3R = uq_maxwellian(KS.vSpace.u, primR, uq)

        ER = zeros(3, uq.nr + 1)
        BR = zeros(3, uq.nr + 1)
        BR[:, 1] .= KS.ib.BR

        lorenz = zeros(3, uq.nr + 1, 2)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wL,
                    primL,
                    h0L,
                    h1L,
                    h2L,
                    h3L,
                    EL,
                    BL,
                    lorenz,
                )
            else
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wR,
                    primR,
                    h0R,
                    h1R,
                    h2R,
                    h3R,
                    ER,
                    BR,
                    lorenz,
                )
            end
        end

        face = Array{Interface1D4F}(undef, KS.pSpace.nx + 1)
        for i = 1:KS.pSpace.nx+1
            face[i] = Interface1D4F(wL, h0L, EL)
        end

    elseif uq.method == "collocation"

        # upstream
        primL = zeros(5, uq.op.quad.Nquad, 2)
        for j in axes(primL, 2)
            primL[:, j, :] .= ks.ib.bc(ks.ps.x0, ks.ib.p)
        end

        wL = uq_prim_conserve(primL, KS.gas.γ, uq)
        h0L, h1L, h2L, h3L = uq_maxwellian(KS.vSpace.u, primL, uq)

        EL = zeros(3, uq.op.quad.Nquad)
        BL = zeros(3, uq.op.quad.Nquad)
        for j in axes(BL, 2)
            BL[:, j] .= KS.ib.BL
        end

        # downstream
        primR = zeros(5, uq.op.quad.Nquad, 2)
        for j in axes(primL, 2)
            primR[:, j, :] .= ks.ib.bc(ks.ps.x1, ks.ib.p)
        end

        wR = uq_prim_conserve(primR, KS.gas.γ, uq)
        h0R, h1R, h2R, h3R = uq_maxwellian(KS.vSpace.u, primR, uq)

        ER = zeros(3, uq.op.quad.Nquad)
        BR = zeros(3, uq.op.quad.Nquad)
        for j in axes(BR, 2)
            BR[:, j] .= KS.ib.BR
        end

        lorenz = zeros(3, uq.op.quad.Nquad, 2)

        for i in eachindex(ctr)
            if i <= KS.pSpace.nx ÷ 2
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wL,
                    primL,
                    h0L,
                    h1L,
                    h2L,
                    h3L,
                    EL,
                    BL,
                    lorenz,
                )
            else
                ctr[i] = ControlVolume1D4F(
                    KS.pSpace.x[i],
                    KS.pSpace.dx[i],
                    wR,
                    primR,
                    h0R,
                    h1R,
                    h2R,
                    h3R,
                    ER,
                    BR,
                    lorenz,
                )
            end
        end

        #--- setup of cell interface ---#
        face = Array{Interface1D4F}(undef, KS.pSpace.nx + 1)
        for i = 1:KS.pSpace.nx+1
            face[i] = Interface1D4F(wL, h0L, EL)
        end

    end

    return ctr, face

end


text:
using Pkg
Pkg.activate(@__DIR__)

using ProgressiveHedging
import JuMP
import Ipopt

function create_model(scenario_id::ScenarioID)

    model = JuMP.Model(()->Ipopt.Optimizer())
    JuMP.set_optimizer_attribute(model, "print_level", 0)
    JuMP.set_optimizer_attribute(model, "tol", 1e-12)

    c = [1.0, 10.0, 0.01]
    d = 7.0
    a = 16.0

    α = 1.0
    β = 1.0
    γ = 1.0
    δ = 1.0
    ϵ = 1.0

    s1 = 8.0
    s2 = 4.0
    s11 = 9.0
    s12 = 16.0
    s21 = 5.0
    s22 = 18.0

    stage1 = JuMP.@variable(model, x[1:3] >= 0.0)
    JuMP.@constraint(model, x[3] <= 1.0)
    obj = zero(JuMP.GenericQuadExpr{Float64,JuMP.VariableRef})
    JuMP.add_to_expression!(obj, sum(c.*x))

    # Second stage
    stage2 = Vector{JuMP.VariableRef}()
    if scenario_id < scid(2)
        vref = JuMP.@variable(model, y >= 0.0)
        JuMP.@constraint(model, α*sum(x) + β*y >= s1)
        JuMP.add_to_expression!(obj, d*y)
    else
        vref = JuMP.@variable(model, y >= 0.0)
        JuMP.@constraint(model, α*sum(x) + β*y >= s2)
        JuMP.add_to_expression!(obj, d*y)
    end
    push!(stage2, vref)

    # Third stage
    stage3 = Vector{JuMP.VariableRef}()
    if scenario_id == scid(0)
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s11)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    elseif scenario_id == scid(1)
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s12)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    elseif scenario_id == scid(2)
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s21)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))

    else
        vref = JuMP.@variable(model, z[1:2])
        JuMP.@constraint(model, ϵ*sum(x) + γ*y + δ*sum(z) == s22)
        JuMP.add_to_expression!(obj, a*sum(z[i]^2 for i in 1:2))
    end
    append!(stage3, vref)

    JuMP.@objective(model, Min, obj)

    vdict = Dict{StageID, Vector{JuMP.VariableRef}}([stid(1) => stage1,
                                                        stid(2) => stage2,
                                                        stid(3) => stage3,
                                                        ])

    return JuMPSubproblem(model, scenario_id, vdict)
end

scen_tree = ScenarioTree()
branch_node_1 = add_node(scen_tree, root(scen_tree))
branch_node_2 = add_node(scen_tree, root(scen_tree))

add_leaf(scen_tree, branch_node_1, 0.375)
add_leaf(scen_tree, branch_node_1, 0.125)

add_leaf(scen_tree, branch_node_2, 0.375)
add_leaf(scen_tree, branch_node_2, 0.125)

(niter, abs_res, rel_res, obj, soln_df, phd) = solve(scen_tree,
                                                     create_model,
                                                     ScalarPenaltyParameter(25.0);
                                                     atol=1e-8, rtol=1e-12, max_iter=500)
@show niter
@show abs_res
@show rel_res
@show obj
@show soln_df


text:
using BenchmarkTools

function reconstruct(a)
    b = similar(a)
    for i in eachindex(a)
        @inbounds b[i] = a[i]
    end
    return b
end

nn = 200000
trn = rand(Int, nn)

@btime deepcopy($trn)
@btime copy($trn)
@btime reconstruct($trn)
@assert copy(trn) == deepcopy(trn) == reconstruct(trn)

nn = 200000
xrand = rand()
trn = collect(ceil(Int, xrand / 2) * nn:ceil(Int, (1 - xrand) / 20) * nn)

@btime deepcopy($trn)
@btime copy($trn)
@btime reconstruct($trn)
@assert copy(trn) == deepcopy(trn) == reconstruct(trn)