text:
export BrickletRGBLEDMatrixSPITFPErrorCount
struct BrickletRGBLEDMatrixSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

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

export BrickletRGBLEDMatrix
"""
RGB LED Matrix with 8x8 pixel
"""
mutable struct BrickletRGBLEDMatrix <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export set_red
"""
    $(SIGNATURES)

Sets the 64 red LED values of the matrix.
"""
function set_red(device::BrickletRGBLEDMatrix, red)
    device.deviceInternal.set_red(red)
end

export get_red
"""
    $(SIGNATURES)

Returns the red LED values as set by :func:`Set Red`.
"""
function get_red(device::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_red()
end

export set_green
"""
    $(SIGNATURES)

Sets the 64 green LED values of the matrix.
"""
function set_green(device::BrickletRGBLEDMatrix, green)
    device.deviceInternal.set_green(green)
end

export get_green
"""
    $(SIGNATURES)

Returns the green LED values as set by :func:`Set Green`.
"""
function get_green(device::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_green()
end

export set_blue
"""
    $(SIGNATURES)

Sets the 64 blue LED values of the matrix.
"""
function set_blue(device::BrickletRGBLEDMatrix, blue)
    device.deviceInternal.set_blue(blue)
end

export get_blue
"""
    $(SIGNATURES)

Returns the blue LED values as set by :func:`Set Blue`.
"""
function get_blue(device::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_blue()
end

export set_frame_duration
"""
    $(SIGNATURES)

Sets the frame duration.

Example: If you want to achieve 20 frames per second, you should
set the frame duration to 50ms (50ms * 20 = 1 second).

Set this value to 0 to turn the automatic frame write mechanism off.

Approach:

* Call :func:`Set Frame Duration` with value > 0.
* Set LED values for first frame with :func:`Set Red`, :func:`Set Green`, :func:`Set Blue`.
* Wait for :cb:`Frame Started` callback.
* Set LED values for second frame with :func:`Set Red`, :func:`Set Green`, :func:`Set Blue`.
* Wait for :cb:`Frame Started` callback.
* and so on.

For frame duration of 0 see :func:`Draw Frame`.
"""
function set_frame_duration(device::BrickletRGBLEDMatrix, frame_duration)
    device.deviceInternal.set_frame_duration(frame_duration)
end

export get_frame_duration
"""
    $(SIGNATURES)

Returns the frame duration as set by :func:`Set Frame Duration`.
"""
function get_frame_duration(device::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_frame_duration()
end

export draw_frame
"""
    $(SIGNATURES)

If you set the frame duration to 0 (see :func:`Set Frame Duration`), you can use this
function to transfer the frame to the matrix.

Approach:

* Call :func:`Set Frame Duration` with 0.
* Set LED values for first frame with :func:`Set Red`, :func:`Set Green`, :func:`Set Blue`.
* Call :func:`Draw Frame`.
* Wait for :cb:`Frame Started` callback.
* Set LED values for second frame with :func:`Set Red`, :func:`Set Green`, :func:`Set Blue`.
* Call :func:`Draw Frame`.
* Wait for :cb:`Frame Started` callback.
* and so on.
"""
function draw_frame(device::BrickletRGBLEDMatrix)
    device.deviceInternal.draw_frame()
end

export get_supply_voltage
"""
    $(SIGNATURES)

Returns the current supply voltage of the Bricklet.
"""
function get_supply_voltage(device::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_supply_voltage()
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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix, 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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix, 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::BrickletRGBLEDMatrix, 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::BrickletRGBLEDMatrix, 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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix, 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::BrickletRGBLEDMatrix)
    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::BrickletRGBLEDMatrix)
    return device.deviceInternal.get_identity()
end

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


text:
module Experiment

using JuMP, COSMO
using Statistics, LinearAlgebra, Random
using Plots
using MATLAB
using FFTW
# Write your package code here.
include("data.jl")
include("position.jl")
include("vector_matrix.jl")
include("solver.jl")
include("basic_function.jl")
include("plot_map.jl")

end

ext: jl


text:
import RDates: WeekendCalendar, @rd_str, setcachedcalendar!, calendar
import Dates: Date

setcachedcalendar!(CALENDARS, "UK/SETTLEMENT", WeekendCalendar() + RuleBasedCalendar([
    # New Year's Day
    PeriodicCalendarRule(rd"1y", rd"1JAN@WEEKEND[NBD]"; cal_mgr=CALENDARS),
    # Good Friday
    PeriodicCalendarRule(rd"1y", rd"0E-1FRI"),
    # Easter Monday
    PeriodicCalendarRule(rd"1y", rd"0E+1MON"),
    # Early May Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+1st MON"),
    # Spring Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+Last MON"),
    # Late Summer Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1AUG+Last MON"),
    # Christmas Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+0b@WEEKEND"; cal_mgr=CALENDARS),
    # Boxing Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+1b@WEEKEND"; cal_mgr=CALENDARS),
    # 50th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(1995,5,8)], exclusions=[Date(1995,5,1)]),
    # Golden Jubilee
    ExplicitDateRule(inclusions=[Date(2002,6,3), Date(2002,6,4)], exclusions=[Date(2002,5,27)]),
    # 2011 Royal Wedding
    ExplicitDateRule(inclusions=[Date(2011,4,29)]),
    # Diamond Jubilee
    ExplicitDateRule(inclusions=[Date(2012,6,4), Date(2012,6,5)], exclusions=[Date(2002,5,28)]),
    # 75th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(2020,5,8)], exclusions=[Date(2020,5,4)]),
    # Platinum Jubilee
    ExplicitDateRule(inclusions=[Date(2022,6,2), Date(2022,6,3)], exclusions=[Date(2022,5,30)])
]))

setcachedcalendar!(CALENDARS, "LONDON", calendar(CALENDARS, "UK/SETTLEMENT"))

setcachedcalendar!(CALENDARS, "UK/STOCK EXCHANGE", WeekendCalendar() + RuleBasedCalendar([
    # New Year's Day
    PeriodicCalendarRule(rd"1y", rd"1JAN@WEEKEND[NBD]"; cal_mgr=CALENDARS),
    # Good Friday
    PeriodicCalendarRule(rd"1y", rd"0E-1FRI"),
    # Easter Monday
    PeriodicCalendarRule(rd"1y", rd"0E+1MON"),
    # Early May Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+1st MON"),
    # Spring Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+Last MON"),
    # Late Summer Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1AUG+Last MON"),
    # Christmas Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+0b@WEEKEND"; cal_mgr=CALENDARS),
    # Boxing Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+1b@WEEKEND"; cal_mgr=CALENDARS),
    # 50th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(1995,5,8)], exclusions=[Date(1995,5,1)]),
    # Golden Jubilee
    ExplicitDateRule(inclusions=[Date(2002,6,3), Date(2002,6,4)], exclusions=[Date(2002,5,27)]),
    # 2011 Royal Wedding
    ExplicitDateRule(inclusions=[Date(2011,4,29)]),
    # Diamond Jubilee
    ExplicitDateRule(inclusions=[Date(2012,6,4), Date(2012,6,5)], exclusions=[Date(2002,5,28)]),
    # 75th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(2020,5,8)], exclusions=[Date(2020,5,4)]),
    # Platinum Jubilee
    ExplicitDateRule(inclusions=[Date(2022,6,2), Date(2022,6,3)], exclusions=[Date(2022,5,30)]),
    # 31Dec1999
    ExplicitDateRule(inclusions=[Date(1999,12,31)])
]))

setcachedcalendar!(CALENDARS, "UK/LSE", calendar(CALENDARS, "UK/STOCK EXCHANGE"))

setcachedcalendar!(CALENDARS, "UK/METAL EXCHANGE", WeekendCalendar() + RuleBasedCalendar([
    # New Year's Day
    PeriodicCalendarRule(rd"1y", rd"1JAN@WEEKEND[NBD]"; cal_mgr=CALENDARS),
    # Good Friday
    PeriodicCalendarRule(rd"1y", rd"0E-1FRI"),
    # Easter Monday
    PeriodicCalendarRule(rd"1y", rd"0E+1MON"),
    # Early May Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+1st MON"),
    # Spring Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1MAY+Last MON"),
    # Late Summer Bank Holiday
    PeriodicCalendarRule(rd"1y", rd"1AUG+Last MON"),
    # Christmas Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+0b@WEEKEND"; cal_mgr=CALENDARS),
    # Boxing Day
    PeriodicCalendarRule(rd"1y", rd"25DEC+1b@WEEKEND"; cal_mgr=CALENDARS),
    # 50th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(1995,5,8)], exclusions=[Date(1995,5,1)]),
    # Golden Jubilee
    ExplicitDateRule(inclusions=[Date(2002,6,3), Date(2002,6,4)], exclusions=[Date(2002,5,27)]),
    # 2011 Royal Wedding
    ExplicitDateRule(inclusions=[Date(2011,4,29)]),
    # Diamond Jubilee
    ExplicitDateRule(inclusions=[Date(2012,6,4), Date(2012,6,5)], exclusions=[Date(2002,5,28)]),
    # 75th VE Day Anniversary
    ExplicitDateRule(inclusions=[Date(2020,5,8)], exclusions=[Date(2020,5,4)]),
    # Platinum Jubilee
    ExplicitDateRule(inclusions=[Date(2022,6,2), Date(2022,6,3)], exclusions=[Date(2022,5,30)]),
    # 31Dec1999
    ExplicitDateRule(inclusions=[Date(1999,12,31)])
]))

setcachedcalendar!(CALENDARS, "UK/LME", calendar(CALENDARS, "UK/METAL EXCHANGE"))


text:
using CuArrays, LazyArrays
CuArrays.allowscalar(false)

include("reduce.jl")

len = 10^7 + 13
input = ones(Int32, len)
output = similar(input, 1024)
cpu_val = reduce(+, input)

gpu_input = CuArray(input)
gpu_output = CuArray(output)

gpu_reduce(+, gpu_input, gpu_output)
gpu_val = Array(gpu_output)[1]
@assert cpu_val == gpu_val

gpu_reduce(+, (@~ gpu_input .+ gpu_input), gpu_output)
gpu_val = Array(gpu_output)[1]
cpu_val = reduce(+, (@~ input .+ input))
@assert cpu_val == gpu_val

a = round.(rand(Float32, (3, 4)) * 100)
b = round.(rand(Float32, (1, 4)) * 100)
d_a = CuArray(a)
d_b = CuArray(b)
d_c = similar(d_a, 1024)

gpu_reduce(+, (@~ d_a .* d_b), d_c)
gpu_val = Array(d_c)[1]
@assert gpu_val == reduce(+, (@~ a .* b))


text:
function histogram2d(dataobject::DataSetType, data1::Array{<:Number,1}, data2::Array{<:Number,1}, xrange::Array{<:Any,1}, yrange::Array{<:Any,1};
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

        return histogram2d(dataobject, [data1, data2];
                            xrange=xrange,
                            yrange=yrange,
                            nbins=nbins,
                            weight=weight,
                            closed=closed,
                            mask=mask,
                            mode=mode,
                            verbose=verbose)
end

function histogram2d(dataobject::DataSetType, data1::Array{<:Number,1}, data2::Array{<:Number,1};
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

        return histogram2d(dataobject, [data1, data2];
                            xrange=xrange,
                            yrange=yrange,
                            nbins=nbins,
                            weight=weight,
                            closed=closed,
                            mask=mask,
                            mode=mode,
                            verbose=verbose)
end

# tested for HydroDataType: todo: test on PartDataType ClumpDataType ...
function histogram2d(dataobject::DataSetType, data::Array{Tuple{Symbol,Symbol,Symbol},1};
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

        if data[1][3] == :linear
            data1 = getvar(dataobject, data[1][1], data[1][2])
        elseif data[1][3] == :log10
            data1 = log10.( getvar(dataobject, data[1][1], data[1][2]) )
        end

        if data[2][3] == :linear
            data2 = getvar(dataobject, data[2][1], data[2][2])
        elseif data[2][3] == :log10
            data2 = log10.( getvar(dataobject, data[2][1], data[2][2]) )
        end

        return histogram2d(dataobject, [data1, data2];
                            xrange=xrange,
                            yrange=yrange,
                            nbins=nbins,
                            weight=weight,
                            closed=closed,
                            mask=mask,
                            mode=mode,
                            verbose=verbose)
end

function histogram2d(dataobject::DataSetType, data::Array{Tuple{Symbol,Symbol},1};
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

        # default given variable tuple: linear scale
        data1 = getvar(dataobject, data[1][1], data[1][2])
        data2 = getvar(dataobject, data[2][1], data[2][2])

        return histogram2d(dataobject, [data1, data2];
                            xrange=xrange,
                            yrange=yrange,
                            nbins=nbins,
                            weight=weight,
                            closed=closed,
                            mask=mask,
                            mode=mode,
                            verbose=verbose)
end

"""
### 2d histogram of 2 given arrays or derived quantities (calculations based on StatsBase.jl)
- pass two 1d data arrays
- or give the quantity with units
- costum/automatic data ranges (left-/right-closed)
- binsizes
- give quantity for weighting with units
- pass a mask to exclude elements (cells/particles/...) from the calculation
- toggle verbose mode

```julia
histogram2d(dataobject::DataSetType, data::Array{Array{Float64,1},1};
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

return Histogram2DMapType
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
- **`data`:** provide two 1d data arrays for axis [x,y]
##### Predefined/Optional Keywords:
- **`xrange`:** the considered interval [xmin, xmax] in units of the given x-axis data
- **`yrange`:** the considered interval [ymin, ymax] in units of the given y-axis data
- **`closed`:** bin intervals :left [a,b)  or :right (a,b] closed
- **`nbins`:** number of bins for x- and y-axis
- **`weight`:** give a tuple with the quantity und unit that is used for weighting (see getvar() for possible derived quantities)
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
***
**`mode`:** see Histogram from StatsBase.jl:
*  `:pdf`: Normalize by sum of weights and bin sizes. Resulting histogram
   has norm 1 and represents a PDF.
* `:density`: Normalize by bin sizes only. Resulting histogram represents
   count density of input and does not have norm 1. Will not modify the
   histogram if it already represents a density.
* `:probability`: Normalize by sum of weights only. Resulting histogram
   represents the fraction of probability mass for each bin and does not have
   norm 1.
*  `:none`: Leaves histogram unchanged. Useful to simplify code that has to
   conditionally apply different modes of normalization.

***

### Defined Methods - function defined for different arguments

- histogram2d(dataobject::DataSetType, data::Array{Array{Float64,1},1}; ...) # enclose two 1d arrays in an array -> [xdata, ydata], ...
- histogram2d(dataobject::DataSetType, data1::Array{<:Number,1}, data2::Array{<:Number,1}; ...) # give two 1d arrays separately without keywords -> xdata, ydata,...
- histogram2d(dataobject::DataSetType, data1::Array{<:Number,1}, data2::Array{<:Number,1}, xrange::Array{<:Any,1}, yrange::Array{<:Any,1}; ...) # give two 1d arrays separately and the intervals without keywords -> xdata, ydata, [xmin,xmax], [ymin,ymax],...
- histogram2d(dataobject::DataSetType, data::Array{Tuple{Symbol,Symbol},1}; ...) # give the quantity names with units for the 2 data components in a tuple (default: assumes linear scale) -> ((:rho, :Msol_pc3), (:cs, :km_s)), ....
- histogram2d(dataobject::DataSetType, data::Array{Tuple{Symbol,Symbol,Symbol},1}; ...) # give the quantity names with units and scaling (:linear or :log10) for the 2 data components in a tuple -> ((:rho, :Msol_pc3, :log10), (:cs, :km_s, :log10)), ....

### Short summary how to give the data and ranges (keep order):
- histogram2d( dataobject, [xdata, ydata], ...)
- histogram2d( dataobject,  xdata, ydata,  ...)
- histogram2d( dataobject,  xdata, ydata, [xmin,xmax], [ymin,ymax], ...)
- histogram2d( dataobject, ((:rho, :nH), (:cs, :km_s)), ....)
- histogram2d( dataobject, ((:rho, :nH, :log10), (:cs, :km_s, :log10)), ....)

"""
function histogram2d(dataobject::DataSetType, data::Array{Array{Float64,1},1};
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        nbins::Array{Int,1}=[100,100],
                        weight::Tuple{Symbol,Symbol}=(:mass,:Msol),
                        closed::Symbol=:left,
                        mask=[false],
                        mode::Symbol=:none,
                        verbose::Bool=verbose_mode)

    printtime("", verbose)

    binning1, binning2, extent, ratio = histogram_binning(data, xrange, yrange, nbins, closed, verbose)

    rows = length(dataobject.data)
    if length(mask) > 1
        if length(mask) !== rows
            error("[Mera] ",now()," : array-mask length: $(length(mask)) does not match with data-table length: $(rows)")

        else

            h_var = fit(Histogram, (data[1] , data[2]),
                        weights( getvar(dataobject, weight[1], weight[2]) .* mask  ),
                        closed=closed,
                        (binning1,binning2) )
        end

    else
        h_var = fit(Histogram, (data[1], data[2]),
                    weights( getvar(dataobject, weight[1], weight[2])  ),
                    closed=closed,
                    (binning1,binning2) )
    end

    h_var = StatsBase.normalize(h_var; mode=mode)

    return Histogram2DMapType(h_var.weights,
                             closed,
                             weight,
                             nbins,
                             binning1,
                             binning2,
                             extent,
                             ratio,
                             dataobject.scale,
                             dataobject.info)
end

function histogram_binning(data::Array{Array{Float64,1},1},
                        xrange::Array{<:Any,1},
                        yrange::Array{<:Any,1},
                        nbins::Array{Int,1},
                        closed::Symbol,
                        verbose::Bool)

    extent = zeros(Float64, 4)
    if xrange[1] === missing
        xmin = minimum(data[1])
    else
        xmin = xrange[1]
    end

    if xrange[2] === missing
        xmax = maximum(data[1])
    else
        xmax = xrange[2]
    end

    if yrange[1] === missing
        ymin = minimum(data[2])
    else
        ymin = yrange[1]
    end

    if yrange[2] === missing
        ymax = maximum(data[2])
    else
        ymax = yrange[2]
    end
    extent[1] = xmin
    extent[2] = xmax
    extent[3] = ymin
    extent[4] = ymax

    ratio = (extent[2]-extent[1]) / (extent[4]-extent[3])

    binning1 = range(xmin, stop=xmax, length=nbins[1])
    binning2 = range(ymin, stop=ymax, length=nbins[2])

    if verbose
        if closed == :left
            println("bin intervals left-closed [a,b)  (default):")
        elseif closed == :right
            println("bin intervals right-closed (a,b]:")
        end
        println("xrange: ", binning1, "\t ->bins=",nbins[1])
        println("yrange: ", binning2, "\t ->bins=",nbins[2])

        println()
    end

    return binning1, binning2, extent, ratio
end


text:
const MIMEVec = AbstractVector{MIME}
const MIMEOrVec = Union{MIME, MIMEVec}

"""
    extramimes(x)

Get a vector of extra (non-default) MIME types to include in a MIMEBundle.
"""
extramimes(x::Any) = MIME[]

"""
Default MIME types to show.

Every MIME type listed here will be included in the default MIMEBundle for an
object if it has a corresponding `show` method. If a nested vector of MIMEs is
specified, the first applicable MIME will be included (and the result will not).

This list is as long as it is for legacy reasons. Types should define an
`extramimes` methods to include extra mimes in their bundle (recommended), or
define a custom `MIMEBundle` method for complete control over exactly what MIME
types are included.
"""
const DEFAULT_MIMES = Vector{MIMEOrVec}([
    MIME("text/plain"),
    MIME("text/latex"),
    MIME("image/svg+xml"),
    [
        MIME("image/png"),
        MIME("image/jpeg"),
    ],
    [
        MIME("text/markdown"),
        MIME("text/html"),
    ],
    [
        MIME("application/vnd.vegalite.v3+json"),
        MIME("application/vnd.vegalite.v2+json"),
        MIME("application/vnd.vega.v5+json"),
        MIME("application/vnd.vega.v4+json"),
        MIME("application/vnd.vega.v3+json"),
    ],
    [
        MIME("application/vnd.dataresource+json"),
        MIME("application/vnd.plotly.v1+json"),
    ],
])


text:
using Documenter, TriMatrices

makedocs(;
	modules=[TriMatrices],
	authors="Jared Lumpe",
	repo="https://github.com/jlumpe/TriMatrices.jl/blob/{commit}{path}#L{line}",
	sitename="TriMatrices.jl",
	pages=[
		"index.md",
		"api.md",
	],
)

deploydocs(
	repo="github.com/jlumpe/TriMatrices.jl.git",
)


text:
export ImmutablePolynomial

"""
    ImmutablePolynomial{T<:Number, N}(coeffs::AbstractVector{T}, var=:x)

Construct an immutable (static) polynomial from its coefficients
`a₀, a₁, …, aₙ`,
lowest order first, optionally in terms of the given variable `x`
where `x` can be a character, symbol, or string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct
this through `ImmutablePolynomial((a_0, a_1, ..., a_n))` (assuming
`a_n ≠ 0`). As well, a vector or number can be used for construction.

The usual arithmetic operators are overloaded to work with polynomials
as well as with combinations of polynomials and scalars. However,
operations involving two non-constant polynomials of different variables causes an
error. Unlike other polynomials, `setindex!` is not defined for `ImmutablePolynomials`.

As the coefficient size is a compile-time constant, several performance
improvements are possible. For example, immutable polynomials can take advantage of
faster polynomial evaluation provided by `evalpoly` from Julia 1.4.

!!! note
    `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
    index always corresponding to the constant term.

# Examples

```jldoctest
julia> using  Polynomials

julia> ImmutablePolynomial((1, 0, 3, 4))
ImmutablePolynomial(1 + 3*x^2 + 4*x^3)

julia> ImmutablePolynomial((1, 2, 3), :s)
ImmutablePolynomial(1 + 2*s + 3*s^2)

julia> one(ImmutablePolynomial)
ImmutablePolynomial(1.0)
```

!!! note
    This was modeled after https://github.com/tkoolen/StaticUnivariatePolynomials.jl by @tkoolen.

"""
struct ImmutablePolynomial{T <: Number,  N} <: StandardBasisPolynomial{T}
    coeffs::NTuple{N, T}
    var::Symbol

    function ImmutablePolynomial{T,N}(coeffs::NTuple{N,T}, var::SymbolLike=:x) where {T <: Number,N}
        N == 0 && return new{T,0}(coeffs, Symbol(var))
        iszero(coeffs[end]) &&  throw(ArgumentError("Leading  term must  be  non-zero"))
        new{T,N}(coeffs, Symbol(var))
    end

end

@register ImmutablePolynomial

## Various interfaces
function ImmutablePolynomial{T,N}(coeffs::Tuple, var::SymbolLike=:x)  where {T,N}
    ImmutablePolynomial{T,N}(NTuple{N,T}(T.(coeffs)), var)
end

function ImmutablePolynomial{T,N}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {T <: Number, N, S}
    if Base.has_offset_axes(coeffs)
      @warn "ignoring the axis offset of the coefficient vector"
    end
    ImmutablePolynomial{T,N}(NTuple{N,T}(tuple(coeffs...)), var)
end

## --
function ImmutablePolynomial{T}(coeffs::NTuple{M,S}, var::SymbolLike=:x) where {T, S<: Number, M}
    N = findlast(!iszero, coeffs)
    if N == nothing
        return zero(ImmutablePolynomial{T}, var)
    else
        cs = NTuple{N,T}(coeffs[i] for  i in  1:N)
        return ImmutablePolynomial{T,N}(cs, var)
    end
end

function ImmutablePolynomial{T}(coeffs::Tuple, var::SymbolLike=:x)  where {T}
    ImmutablePolynomial{T}(T.(coeffs), var)
end

# entry point from abstract.jl; note T <: Number
function ImmutablePolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
    if Base.has_offset_axes(coeffs)
      @warn "ignoring the axes of the coefficient vector and treating it as a list"
    end
    M = length(coeffs)
    ImmutablePolynomial{T}(NTuple{M,T}(tuple(coeffs...)), var)
end

## --

function ImmutablePolynomial(coeffs::Tuple, var::SymbolLike=:x)
    cs = NTuple(promote(coeffs...))
    T = eltype(cs)
    ImmutablePolynomial{T}(cs, var)
end

# Convenience; pass tuple to Polynomial
# Not documented, not sure this is a good idea as P(...)::P is not true...
# Deprecated
function Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T}
    Base.depwarn("Use of `Polynomial(NTuple, var)` is deprecated. Use the `ImmutablePolynomial` constructor",
                 :Polynomial)
    ImmutablePolynomial(coeffs, var)
end
function Polynomial{T}(coeffs::NTuple{N,S}, var::SymbolLike = :x) where{N,T,S}
    Base.depwarn("Use of `Polynomial(NTuple, var)` is deprecated. Use the `ImmutablePolynomial` constructor",
                 :Polynomial)
    ImmutablePolynomial{N,T}(T.(coeffs), var)
end

##
## ----
##
# overrides from common.jl due to  coeffs being non mutable, N in type parameters
Base.collect(p::P) where {P <: ImmutablePolynomial} = [pᵢ for pᵢ ∈ p]

Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p), p.var)

## defining these speeds things up
function Base.zero(P::Type{<:ImmutablePolynomial}, var::SymbolLike=:x)
    R = eltype(P)
    ImmutablePolynomial{R,0}(NTuple{0,R}(),var)
end

function  Base.one(P::Type{<:ImmutablePolynomial}, var::SymbolLike=:x)
    R = eltype(P)
    ImmutablePolynomial{R,1}(NTuple{1,R}(1),var)
end
function variable(P::Type{<:ImmutablePolynomial}, var::SymbolLike=:x)
    R  = eltype(P)
    ImmutablePolynomial{R,2}(NTuple{2,R}((0,1)),var)
end

# degree, isconstant
degree(p::ImmutablePolynomial{T,N}) where {T,N} = N - 1 # no trailing zeros
isconstant(p::ImmutablePolynomial{T,N}) where {T,N}  = N <= 1

function Base.getindex(p::ImmutablePolynomial{T,N}, idx::Int) where {T <: Number,N}
    (idx <  0 || idx > N-1) && return zero(T)
    return p.coeffs[idx + 1]
end

Base.setindex!(p::ImmutablePolynomial, val::Number,  idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable"))

for op in [:isequal, :(==)]
    @eval function Base.$op(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
        check_same_variable(p1,p2) || return false
        p1s, p2s = coeffs(p1), coeffs(p2)
        (N == M  && $op(p1s,p2s)) &&  return  true
        n1 = findlast(!iszero, p1s) # now trim out zeros
        n2 = findlast(!iszero, p2s)
        (n1 == nothing && n2 == nothing) && return true
        (n1 == nothing || n2  == nothing) && return false
        $op(p1s[1:n1],p2s[1:n2]) &&  return true
        false
    end
end

# in common.jl these call chop! and truncate!
function Base.chop(p::ImmutablePolynomial{T,N};
              rtol::Real = Base.rtoldefault(real(T)),
              atol::Real = 0)  where {T,N}
    cs = coeffs(p)
    for i in N:-1:1
        if !isapprox(cs[i], zero(T), rtol=rtol, atol=atol)
            return ImmutablePolynomial{T,i}(cs[1:i], p.var)
        end
    end
    zero(ImmutablePolynomial{T}, p.var)
end

function Base.truncate(p::ImmutablePolynomial{T,N};
                       rtol::Real = Base.rtoldefault(real(T)),
                       atol::Real = 0)  where {T,N}
    q = chop(p, rtol=rtol, atol=atol)
    iszero(q) && return  q
    cs = coeffs(q)
    thresh = maximum(abs,cs) * rtol + atol
    cs′ = map(c->abs(c) <= thresh ? zero(T) : c, cs)
    ImmutablePolynomial{T}(tuple(cs′...), p.var)
end

# no in-place chop! and truncate!
chop!(p::ImmutablePolynomial; kwargs...) =  chop(p; kwargs...)
truncate!(p::ImmutablePolynomial; kwargs...) =  truncate(p; kwargs...)

##
## --------------------
##

(p::ImmutablePolynomial{T,N})(x::S) where {T,N,S} = evalpoly(x, p.coeffs)

function Base.:+(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}

    R = promote_type(S,T)
    iszero(N) && return ImmutablePolynomial{R}(coeffs(p2), p2.var)
    iszero(M) && return ImmutablePolynomial{R}(coeffs(p1), p1.var)

    isconstant(p1) && p1.var != p2.var && return p2 + p1[0]*one(ImmutablePolynomial{R}, p2.var)
    isconstant(p2) && p1.var != p2.var && return p1 + p2[0]*one(ImmutablePolynomial{R}, p1.var)

    p1.var != p2.var && error("Polynomials must have same variable")

    if  N == M
        cs = NTuple{N,R}(p1[i] + p2[i] for i in 0:N-1)
        ImmutablePolynomial{R}(cs, p1.var)
    elseif N < M
        cs = (p2.coeffs) ⊕ (p1.coeffs)
        ImmutablePolynomial{R,M}(cs, p1.var)
    else
        cs = (p1.coeffs) ⊕ (p2.coeffs)
        ImmutablePolynomial{R,N}(cs, p1.var)
    end

end

# not type stable!!!
function Base.:*(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
    isconstant(p1) && return p2 * p1[0]
    isconstant(p2) && return p1 * p2[0]
    p1.var != p2.var && error("Polynomials must have same variable")
    R = promote_type(S,T)
    cs = (p1.coeffs) ⊗ (p2.coeffs)
    if !iszero(cs[end])
        return ImmutablePolynomial{R, N+M-1}(cs, p1.var)
    else
        n = findlast(!iszero, cs)
        return ImmutablePolynomial{R, n}(cs[1:n], p1.var)
    end
end

# Padded vector sum of two tuples assuming N > M
# assume N > M.
# As N ≠ M, we are assured of size of output (max(N,M)), so we generate the function
@generated function ⊕(p1::NTuple{N,T}, p2::NTuple{M,S}) where {T,N,S,M}

    R = promote_type(T,S)

    exprs = Any[nothing for i = 1:N]
    for i in  1:M
        exprs[i] = :(p1[$i] + p2[$i])
    end
    for i in M+1:N
        exprs[i] =:(p1[$i])
    end

    return quote
        Base.@_inline_meta
        tuple($(exprs...))
    end

end

## Static size of product makes generated functions  a good choice
## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl
## convolution of two tuples
@generated function ⊗(p1::NTuple{N,T}, p2::NTuple{M,S}) where {T,N,S,M}
    P = M + N - 1
    R = promote_type(T,S)
    exprs = Any[nothing for i = 1 : P]
    for i in 1 : N
        for j in 1 : M
            k = i + j - 1
            if exprs[k] === nothing
                exprs[k] = :(p1[$i] * p2[$j])
            else
                exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k])))
            end
        end
    end

    return quote
        Base.@_inline_meta
        tuple($(exprs...))
    end

end

# scalar ops
function Base.:+(p::ImmutablePolynomial{T,N}, c::S) where {T, N, S<:Number}
    R = promote_type(T,S)

    iszero(c) && return ImmutablePolynomial{R,N}(p.coeffs, p.var)
    N == 0 && return ImmutablePolynomial{R,1}((c,), p.var)
    N == 1 && return ImmutablePolynomial((p[0]+c,), p.var)

    q = ImmutablePolynomial{R,1}((c,), p.var)
    return p + q

end

function Base.:*(p::ImmutablePolynomial{T,N}, c::S) where {T, N, S <: Number}
    R = promote_type(T,S)
    iszero(c) && return zero(ImmutablePolynomial{R}, p.var)
    ImmutablePolynomial{R,N}(p.coeffs .* c, p.var)
end

function Base.:/(p::ImmutablePolynomial{T,N}, c::S) where {T,N,S <: Number}
    R = eltype(one(T)/one(S))
    isinf(c)  && return zero(ImmutablePolynomial{R}, p.var)
    ImmutablePolynomial{R,N}(p.coeffs ./ c, p.var)
end

Base.:-(p::ImmutablePolynomial{T,N}) where {T,N} = ImmutablePolynomial{T,N}(.-p.coeffs, p.var)

Base.to_power_type(p::ImmutablePolynomial{T,N}) where {T,N} = p

## more performant versions borrowed from StaticArrays
## https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/linalg.jl
LinearAlgebra.norm(q::ImmutablePolynomial{T,0}) where {T} = zero(real(float(T)))
LinearAlgebra.norm(q::ImmutablePolynomial) = _norm(q.coeffs)
LinearAlgebra.norm(q::ImmutablePolynomial, p::Real) = _norm(q.coeffs, p)

@generated function _norm(a::NTuple{N,T}) where {T, N}

    expr = :(abs2(a[1]))
    for j = 2:N
        expr = :($expr + abs2(a[$j]))
    end

    return quote
        $(Expr(:meta, :inline))
        @inbounds return sqrt($expr)
    end

end

_norm_p0(x) = iszero(x) ? zero(x) : one(x)
@generated function _norm(a::NTuple{N,T}, p::Real) where {T, N}
    expr = :(abs(a[1])^p)
    for j = 2:N
        expr = :($expr + abs(a[$j])^p)
    end

    expr_p1 = :(abs(a[1]))
    for j = 2:N
        expr_p1 = :($expr_p1 + abs(a[$j]))
    end

    return quote
        $(Expr(:meta, :inline))
        if p == Inf
            return mapreduce(abs, max, a)
        elseif p == 1
            @inbounds return $expr_p1
        elseif p == 2
            return norm(a)
        elseif p == 0
            return mapreduce(_norm_p0, +, a)
        else
            @inbounds return ($expr)^(inv(p))
        end
    end

end


text:
using JSON
Nproc=4
rmprocs(workers())
addprocs(Nproc)

println("workers: ",length(workers()))

@everywhere using Base.Test
# @everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim
@everywhere gc(true)

# cluster idibaps
@everywhere path="/home/emma/github/fit_behaviour_julia/"
# cluster
# @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
#
@everywhere include(path*"functions/functions_julia.jl")
@everywhere include(path*"functions/parameter_functions.jl")
@everywhere include(path*"functions/utils.jl")

# rat_num = ARGS[1]
rat_num = "1"
path_load=path*"scripts_rats/data/"
filename="processed_data_rat"*rat_num*".json"
filename_1 = "set_trials_10.json"
data = Dict()
open(path_load*filename, "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end
data_1 = Dict()
open(path_load*filename_1, "r") do f
    global data_1
    dicttxt = readstring(f)  # file information to string
    data_1=JSON.parse(dicttxt)  # parse and transform data
end

dataset = data_1[parse(Int64,rat_num)]
n_stim = size(dataset)[1]
stim_10 = data["stim_10"]
stim_20 = data["stim_20"]
choices_10 = data["choices_10"]
choices_20 = data["choices_20"]

n_stim_10 = data["n_stim_10"]
n_stim_20 = data["n_stim_20"]

# if you want to try with less n_stim
# n_stim_10 = 10000
# n_stim_20 = 10000

model_num="10"

args=[("c2",[1]),("c4",0),("st",[1,2]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",[1])]
consts=[("c2",[2]),("c4",[1]),("st",[3,4]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",0)]

z = [0,1,0,0,0,0,0,0]

lb=[-10.,0.0,0.0,0.1]
ub=[20.,20.,20.,20.]

model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>0.0)

aux=2*(rand(4)-0.5)
x_ini=[ aux[1]*10.,abs(aux[2]*10.),abs(aux[3]*10.),rand()*9.9+0.1]

rats_sec = ["2","3","5"]
rats_half = ["1","4"]
tic()
mu_vec = [-1,-0.4,-0.1,0.1,0.4,1]
if n_stim_10 != 0
    history_bias=zeros(n_stim)
    Tframe = 0.05
    stim_10 = permutedims(reshape(stim_10,(10,n_stim_10)))[1:10,:]
    choices_10 = data["choices_10"][1:10,:]
    # mu_10 = data["coherences_10"][1:10,:]
    # sigma_10 = data["sigma_10"][1:10,:]
    sigma_10 = [0,1,2,3,4,5,6,7,8,9]
    fit = ModelFitting(stim_10,sigma_10,choices_10,history_bias, Tframe, args,consts,ub,lb, x_ini,z,model_f)
end
if n_stim_20 != 0
    history_bias=zeros(n_stim)
    Tframe = 0.05
    stim_20 = permutedims(reshape(stim_20,(20,n_stim_20)))[1:10,:]
    choices_20 = data["choices_20"][1:10,:]
    mu_20 = data["coherences_20"][1:10,:]
    sigma_20 = data["sigma_20"][1:10,:]
    fit = ModelFitting_rats(stim_20,sigma_20,mu_20,mu_vec,choices_20,history_bias, Tframe, args,consts,ub,lb, x_ini,z,model_f)
end
# else # la rata 2 té trials tant de 10 com de 20 frames COMPROBAR SI ALGUNA RATA MÉS
#     if ARGS[2] == "10"
#         Nstim_train = trunc(Int,n_stim_10*0.8)
#         stim_train = stim_10[1:Nstim_train,:]
#         stim_test = stim_10[Nstim_train+1:n_stim_10,:]
#         choices_train = choices_10[1:Nstim_train,:]
#         history_bias=zeros(Nstim_train)
#         history_bias_test=zeros(n_stim_10-Nstim_train)
#         Tframe = 10.
#         fit = ModelFitting(stim_train,choices_train,history_bias, Tframe, args,consts,ub,lb, x_ini,z,model_f)
#     elseif ARGS[2] == "20"
#         Nstim_train = trunc(Int,n_stim_10*0.8)
#         stim_train = stim_20[1:Nstim_train,:]
#         stim_test = stim_20[Nstim_train+1:n_stim_20,:]
#         choices_train = choices_20[1:Nstim_train,:]
#         history_bias=zeros(Nstim_train)
#         history_bias_test=zeros(n_stim_20-Nstim_train)
#         Tframe = 20.
#         fit = ModelFitting(stim_train,choices_train,history_bias, Tframe, args,consts,ub,lb, x_ini,z,model_f)
#     end
# end
time_fit = toc()
println("time to fit: "*string(time_fit))
param=make_dict(args,fit["x_bf"],consts,z)

## SAVE THE RESULTS ##
j=1
path_save_res=path*"results/2model10C2St1St2Sig/Rat"*rat_num*"/"
mkpath(path_save_res)
filename_save=path_save_res*"icond"*string(j)*".json"

while isfile(filename_save)
        j=j+1
        filename_save=path_save_res*"icond"*string(j)*".json"
end
tic()
results = Dict("res"=>fit,"icond"=>j,"param"=>param)
time_results = toc()
println("time to dict: "*string(time_results))

tic()
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(results)

# write the file with the stringdata variable information
open(filename_save, "w") do f
        write(f, stringdata)
     end
time_json = toc()
println("time_to_json: "*string(time_json))

## FINISHED ##

println("done")

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
tests = [
  # "util",
  # "variance",
  # "invHessian",
  # "skron",
  # "ustats",
  "reduced"
]

for t in tests
	f = "$t.jl"
	println("* running $f ...")
    t = @elapsed include(f)
    println("done (took $t seconds).")
end


text:
mutable struct PlanktonDiagnostics
    plankton::NamedTuple       # for each species
    tracer::NamedTuple         # for tracers
    time_interval::Int64       # time interval that the diagnostics is time averaged
end

"""
    PlanktonDiagnostics(model; tracer=(:PAR, :NH4, :NO3, :DOC),
                        plankton=(:num, :graz, :mort, :dvid),
                        time_interval = 3600)

Generate a `PlanktonDiagnostics` structure.

Keyword Arguments (Optional)
============================
- `tracer` : a `Tuple` containing the names of nutrient fields to be diagnosed.
- `plankton` : a `Tuple` containing the names of physiological processes of plankton individuals to be diagnosed.
- `time_interval` : The time interval that diagnostics is averaged, an hour (3600 seconds) by default.
"""
function PlanktonDiagnostics(model; tracer=(:PAR,),
                            plankton=(:num, :graz, :mort, :dvid),
                            time_interval::Int64 = 3600)

    @assert isa(tracer, Tuple)
    @assert isa(plankton, Tuple)

    diag_avail(tracer, plankton, model.mode)

    ntr   = length(tracer)
    nproc = length(plankton)
    trs   = []
    procs = []

    total_size = (model.grid.Nx+model.grid.Hx*2, model.grid.Ny+model.grid.Hy*2, model.grid.Nz+model.grid.Hz*2)

    for i in 1:ntr
        tr = zeros(total_size) |> array_type(model.arch)
        push!(trs, tr)
    end
    tr_d = zeros(total_size) |> array_type(model.arch)
    tr_default = (PAR = tr_d,)

    diag_tr = NamedTuple{tracer}(trs)
    diag_tr = merge(diag_tr, tr_default) # add PAR as default diagnostic

    plank_name = keys(model.individuals.phytos)
    Nsp = length(plank_name)

    for j in 1:Nsp
        procs_sp = []
        for k in 1:nproc
            proc = zeros(total_size) |> array_type(model.arch)
            push!(procs_sp, proc)
        end
        diag_proc = NamedTuple{plankton}(procs_sp)

        procs_sp_d = []
        for l in 1:4
            proc = zeros(total_size) |> array_type(model.arch)
            push!(procs_sp_d, proc)
        end
        diag_proc_default = NamedTuple{(:num, :graz, :mort, :dvid)}(procs_sp_d)

        diag_proc = merge(diag_proc, diag_proc_default) # add num, graz, mort, and dvid as default diagnostics

        push!(procs, diag_proc)
    end
    diag_sp = NamedTuple{plank_name}(procs)

    diagnostics = PlanktonDiagnostics(diag_sp, diag_tr, time_interval)

    return diagnostics
end

function show(io::IO, diags::PlanktonDiagnostics)
    print(io, "PlanktonDiagnostics:\n",
              "├── diagnostics of tracers: $(keys(diags.tracer))\n",
              "├── diagnostics of individuals: $(keys(diags.plankton.sp1))\n",
              "└── save averaged diagnostics every $(diags.time_interval) seconds")
end

function diag_avail(tracer, plank, mode::AbstractMode)
    tracer_avail = tracer_avail_diags()
    plank_avail  = plank_avail_diags(mode)
    for i in 1:length(tracer)
        if length(findall(x->x==tracer[i], tracer_avail)) == 0
            throw(ArgumentError("$(tracer[i]) is not one of the diagnostics"))
        end
    end

    for i in 1:length(plank)
        if length(findall(x->x==plank[i], plank_avail)) == 0
            throw(ArgumentError("$(plank[i]) is not one of the diagnostics"))
        end
    end
end

function tracer_avail_diags()
    return (:PAR, :DIC, :DOC, :POC, :NH4, :NO3, :DON, :PON, :PO4, :DOP, :POP)
end

function plank_avail_diags(mode::AbstractMode)
    plank_avail = (:num, :graz, :mort, :dvid, :PS, :resp, :Bm, :Chl, :Th)
    if isa(mode, CarbonMode)
        nothing
    elseif isa(mode, QuotaMode)
        plank_avail = (:num, :graz, :mort, :dvid, :PS, :BS, :VDOC, :VNH4, :VNO3, :VPO4, :resp, :exu, :Bm, :Cq, :Nq, :Pq, :Chl)
    elseif isa(mode, MacroMolecularMode)
        plank_avail = (:num, :graz, :mort, :dvid, :PS, :VDOC, :VNH4, :VNO3, :VPO4, :resp, :ρChl, :S_PRO, :S_DNA, :S_RNA, :exu, :CH, :NST, :PST, :PRO, :DNA, :RNA, :Chl)
    end
    return plank_avail
end


text:
# contact point velocity
function contact_point_velocity(model::Contact, x, q, v, ϕ, c)
    return v + skew(vector_rotate(ϕ, q)) * (c - x)
end

function ∂contact_point_velocity∂x(model::Contact, x, q, v, ϕ, c)
    return -skew(vector_rotate(ϕ, q))
end

function ∂contact_point_velocity∂q(model::Contact, x, q, v, ϕ, c)
    return -skew(c - x) * ∂vector_rotate∂q(ϕ, q)
end

function ∂contact_point_velocity∂v(model::Contact, x, q, v, ϕ, c)
    return 1.0 * I(3)
end

function ∂contact_point_velocity∂ϕ(model::Contact, x, q, v, ϕ, c)
    return -skew(c - x) * rotation_matrix(q)
end

function ∂contact_point_velocity∂c(model::Contact, x, q, v, ϕ, c)
    return skew(vector_rotate(ϕ, q))
end

# relative tangential velocity
function relative_tangential_velocity(model::Contact, xp, qp, vp, ϕp, xc, qc, vc, ϕc)
    # contact points
    cp = contact_point(:parent, model.collision, xp, qp, xc, qc)
    cc = contact_point(:child, model.collision, xp, qp, xc, qc)

    # contact point velocities
    vp = contact_point_velocity(model, xp, qp, vp, ϕp, cp)
    vc = contact_point_velocity(model, xc, qc, vc, ϕc, cc)

    # tangent projection
    return contact_tangent(model.collision, xp, qp, xc, qc) * (vp - vc)
end

function ∂relative_tangential_velocity∂x(jacobian::Symbol, model::Contact, xp, qp, vp, ϕp, xc, qc, vc, ϕc)
    # contact points
    cp = contact_point(:parent, model.collision, xp, qp, xc, qc)
    cc = contact_point(:child, model.collision, xp, qp, xc, qc)

    # contact point velocities
    vp = contact_point_velocity(model, xp, qp, vp, ϕp, cp)
    vc = contact_point_velocity(model, xc, qc, vc, ϕc, cc)

    Δv = vp - vc

    # Jacobians
    if jacobian == :parent
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂x(model, xp, qp, vp, ϕp, cp)
        X += contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xp, qp, vp, ϕp, cp) * ∂contact_point∂x(:parent, jacobian, model.collision, xp, qp, xc, qc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xc, qc, vc, ϕc, cc) * ∂contact_point∂x(:child, jacobian, model.collision, xp, qp, xc, qc)
    elseif jacobian == :child
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xp, qp, vp, ϕp, cp) * ∂contact_point∂x(:parent, jacobian, model.collision, xp, qp, xc, qc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂x(model, xc, qc, vc, ϕc, cc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xc, qc, vc, ϕc, cc) * ∂contact_point∂x(:child, jacobian, model.collision, xp, qp, xc, qc)
    end

    X += [
        Δv' * ∂contact_tangent_one_transpose∂x(jacobian, model.collision, xp, qp, xc, qc);
        Δv' * ∂contact_tangent_two_transpose∂x(jacobian, model.collision, xp, qp, xc, qc);
    ]

    return X
end

function ∂relative_tangential_velocity∂q(jacobian::Symbol, model::Contact, xp, qp, vp, ϕp, xc, qc, vc, ϕc)
    # contact points
    cp = contact_point(:parent, model.collision, xp, qp, xc, qc)
    cc = contact_point(:child, model.collision, xp, qp, xc, qc)

    # contact point velocities
    vp = contact_point_velocity(model, xp, qp, vp, ϕp, cp)
    vc = contact_point_velocity(model, xc, qc, vc, ϕc, cc)

    Δv = vp - vc

    # Jacobians
    if jacobian == :parent
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂q(model, xp, qp, vp, ϕp, cp)
        X += contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xp, qp, vp, ϕp, cp) * ∂contact_point∂q(:parent, jacobian, model.collision, xp, qp, xc, qc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xc, qc, vc, ϕc, cc) * ∂contact_point∂q(:child, jacobian, model.collision, xp, qp, xc, qc)
    elseif jacobian == :child
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xp, qp, vp, ϕp, cp) * ∂contact_point∂q(:parent, jacobian, model.collision, xp, qp, xc, qc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂q(model, xc, qc, vc, ϕc, cc)
        X -= contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂c(model, xc, qc, vc, ϕc, cc) * ∂contact_point∂q(:child, jacobian, model.collision, xp, qp, xc, qc)
    end

    # X += ∂contact_tangent_jvp∂q(jacobian, model.collision, xp, qp, xc, qc, Δv)
    X += [
        Δv' * ∂contact_tangent_one_transpose∂q(jacobian, model.collision, xp, qp, xc, qc);
        Δv' * ∂contact_tangent_two_transpose∂q(jacobian, model.collision, xp, qp, xc, qc);
    ]

    return X
end

function ∂relative_tangential_velocity∂v(jacobian::Symbol, model::Contact, xp, qp, vp, ϕp, xc, qc, vc, ϕc)
    # contact points
    cp = contact_point(:parent, model.collision, xp, qp, xc, qc)
    cc = contact_point(:child, model.collision, xp, qp, xc, qc)

    # contact point velocities
    vp = contact_point_velocity(model, xp, qp, vp, ϕp, cp)
    vc = contact_point_velocity(model, xc, qc, vc, ϕc, cc)

    # Δv = vp - vc

    # Jacobians
    if jacobian == :parent
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂v(model, xp, qp, vp, ϕp, cp)
    elseif jacobian == :child
        X = -1.0 * contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂v(model, xc, qc, vc, ϕc, cc)
    end

    return X
end

function ∂relative_tangential_velocity∂ϕ(jacobian::Symbol, model::Contact, xp, qp, vp, ϕp, xc, qc, vc, ϕc)
    # contact points
    cp = contact_point(:parent, model.collision, xp, qp, xc, qc)
    cc = contact_point(:child, model.collision, xp, qp, xc, qc)

    # contact point velocities
    vp = contact_point_velocity(model, xp, qp, vp, ϕp, cp)
    vc = contact_point_velocity(model, xc, qc, vc, ϕc, cc)

    # Δv = vp - vc

    if jacobian == :parent
        X = contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂ϕ(model, xp, qp, vp, ϕp, cp)
    elseif jacobian == :child
        X = -1.0 * contact_tangent(model.collision, xp, qp, xc, qc) * ∂contact_point_velocity∂ϕ(model, xc, qc, vc, ϕc, cc)
    end

    return X
end


text:
using PyCall
using PyPlot

function plot2Dclassifier(X,y,model;Xtest=[],ytest=[],biasIncluded=false,k=2)

	if biasIncluded
		f1 = 2
		f2 = 3
	else
		f1 = 1
		f2 = 2
	end

	increment = 100

	# Pick some colors for the classes (bor binary it uses red/blue)
	colours = [(0,1,0)
			(1,0,0)
			(0,0,1)
			(1,1,1)
			(1,0,1)
			(0,1,1)
			(1,1,0)
			(.1,.1,.1)
			(1,.5,0)
			(0,.5,0)
			(.5,.5,.5)
			(.5,.25,0)
			(.5,0,.5)
			(0,.5,1)]

	figure()
	if k == 2
		plot(X[y.==1,f1],X[y.==1,f2],"b+")
		plot(X[y.==2,f1],X[y.==2,f2],"ro")
		if !isempty(Xtest)
			plot(Xtest[ytest.==1,f1],Xtest[ytest.==1,f2],"bx");
			plot(Xtest[ytest.==2,f1],Xtest[ytest.==2,f2],"rs");
		end
	else

		for c in 1:k
			plot(X[y.==c,f1],X[y.==c,f2],"+",color=colours[c])
			if !isempty(Xtest)
				plot(Xtest[ytest.==c,f1],Xtest[ytest.==c,f2],"x",color=colours[c])
			end
		end
	end

	(xmin,xmax) = xlim()
	xDomain = linspace(xmin,xmax,increment)
	(ymin,ymax) = ylim()
	yDomain = linspace(ymin,ymax,increment)

	xValues = repmat(xDomain,1,length(xDomain))
	yValues = repmat(yDomain',length(yDomain),1)

	if biasIncluded
		t = length(xValues)
		z = model.predict([ones(t,1) xValues[:] yValues[:]])
	else
		z = model.predict([xValues[:] yValues[:]])
	end
	@assert(length(z) == length(xValues),"Size of model function's output is wrong");

	zValues = reshape(z,size(xValues))

	@pyimport matplotlib.colors as matcolors
	if k == 2
		if all(zValues[:] == 1)
    		cm = [(0,0,.5)];
		elseif all(zValues[:] == 2)
    		cm = [(.5,0,0)];
		else
    		cm = [(0,0,.5);(.5,0,0)];
		end
		cmap = matcolors.ListedColormap(cm,"A")

		contourf(xValues,yValues,zValues,cmap=cmap)
	else
		cm = []
		nColours = 0
		for c in 1:k
			if any(zValues[:] .== c)
				nColours += 1
				colour = (.5colours[c][1],.5colours[c][2],.5colours[c][3])
				zValues[zValues .== c] = nColours
				push!(cm,colour)
			end
		end
		levels = 0:nColours
		cmap = matcolors.ListedColormap(cm,"A")
		contourf(xValues,yValues,zValues,levels=levels,cmap=cmap)
	end
end


text:
push!(LOAD_PATH, string(pwd(),"/src"))
using ModifiedStokesSolver
using PyPlot
using LaTeXStrings

n = 500
zmax = 5
z = linspace(1e-8, zmax, n)
T1S = zeros(n)
T2S = zeros(n)
T3S = zeros(n)
T1L = zeros(n)
T2L = zeros(n)
T3L = zeros(n)

for i=1:n
    T1S[i], T2S[i], T3S[i], T1L[i], T2L[i], T3L[i] = ModifiedStokesSolver.Ti_split(z[i])
end

figure(1)
clf()
plot(z, T1S, label="T1S(z)")
plot(z, T2S, label="T2S(z)")
plot(z, T3S, label="T3S(z)")
plot(z, T1L, label="T1L(z)")
plot(z, T2L, label="T2L(z)")
plot(z, T3L, label="T3L(z)")

xlim(0, zmax)
legend()
xlabel("z")
grid("on")
show()

figure(2)
clf()
f = [1.0, 0.0]
alpha = 1.0
x = linspace(-3.0, 3.0, 10000)
u = map( x -> stokeslet([x, 0.0], f, alpha)[1], x)
plot(x, u, label=L"S(\alpha|x-y|)")
axis([-3,3,0,0.5])
xlabel(L"\alpha |x-y|")
legend()


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_script/script.jl
# A module used to manipulate script function inputs.
#############################################################################

module Script

    using MathOptInterface: AbstractOptimizer, features_available, initialize, NLPBlockData
    using SparseArrays: rowvals, nzrange, SparseMatrixCSC, spzeros
    import Base: afoldl, getindex, iterate
    import JuMP: _NLPData, _NonlinearExprData, NLPEvaluator
    import JuMP._Derivatives: univariate_operator_to_id, operator_to_id,
                              comparison_operator_to_id, CALL, CALLUNIVAR,
                              COMPARISON, VARIABLE, MOIVARIABLE, VALUE,
                              SUBEXPRESSION, PARAMETER, NodeType, NodeData,
                              USER_OPERATOR_ID_START, USER_UNIVAR_OPERATOR_ID_START,
                              UserOperatorRegistry, adjmat
    import Cassette: @context, Cassette.overdub, Cassette.prehook
    import Cassette
    #import CodeTransformation: addmethod! TODO: add this later
    import ForwardDiff: derivative, gradient!
    export dag_flattening!

    include("codetransformation.jl") # TODO: delete this later
    #include("scrubber.jl")
    include("substitute.jl")
    include("patterns.jl")
    include("tracer.jl")
    include("loader.jl")
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
"""
    testing the mean free path calculation
"""

import Particula as P
using Unitful, Test

@testset "λg" begin
    @test P.λg() |> upreferred |> unit == u"m"
    @test isapprox(P.λg() |> ustrip, 65e-9, atol=2e-9)
    @test P.λg(T = 300) |> unit == u"m"
    @test_throws ArgumentError P.λg(T = 300 * u"C")
end


text:
module OpenTelemetryInstrumentationDownloads

using OpenTelemetryAPI

using ArgTools
using Downloads
using Downloads.Curl

using TOML

const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
const INSTRUMENTATION_INFO = InstrumentationInfo(
    name = string(@__MODULE__),
    version = PKG_VERSION,
    schema_url = "https://oolong.dev/OpenTelemetry.jl/dev/OpenTelemetryInstrumentationDownloads/",
)

const download = Downloads.download

const DOWNLOAD_METRICS = Ref{Counter{UInt}}()
const DOWNLOAD_TRACER = Ref{Tracer}()

if VERSION <= v"1.7.2"
    function Downloads.download(
        url::AbstractString,
        output::Union{ArgWrite,Nothing} = nothing;
        method::Union{AbstractString,Nothing} = nothing,
        headers::Union{AbstractVector,AbstractDict} = Pair{String,String}[],
        timeout::Real = Inf,
        progress::Union{Function,Nothing} = nothing,
        verbose::Bool = false,
        downloader::Union{Downloader,Nothing} = nothing,
    )::ArgWrite
        with_span(
            "Downloads $method",
            DOWNLOAD_TRACER[];
            kind = SPAN_KIND_CLIENT,
            attributes = Dict{String,TAttrVal}(
                "http.method" => isnothing(method) ? "" : method,
                "http.url" => url,
            ),
        ) do
            s = current_span()
            arg_write(output) do output
                response = request(
                    url,
                    output = output,
                    method = method,
                    headers = inject!(headers),
                    timeout = timeout,
                    progress = progress,
                    verbose = verbose,
                    downloader = downloader,
                )::Response
                DOWNLOAD_METRICS[](; proto = response.proto, status = response.status)
                s["http.scheme"] = response.proto
                s["http.status_code"] = response.status
                status_ok(response.proto, response.status) && return output
                throw(RequestError(url, Curl.CURLE_OK, "", response))
            end
        end
    end
end

"""
## Metrics

| Meter Name  | Instrument Name | Instrument Type | Unit | Description          |
|:----------- |:--------------- |:--------------- |:---- |:-------------------- |
| `Downloads` | `download`      | `Counter{UInt}` |      | Number of downloads. |

## Spans

We try to follow [semantic conventions for HTTP
spans](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/trace/semantic_conventions/http.md)
here.

The following attributes are added on span creation:

  - `http.method`
  - `http.url`

The following attributes are added when response is received:

  - `http.scheme`
  - `http.status_code`
"""
function init(;
    meter_provider = global_meter_provider(),
    tracer_provider = global_tracer_provider(),
)
    DOWNLOAD_METRICS[] = Counter{UInt}(
        "download",
        Meter(
            "Downloads";
            provider = meter_provider,
            instrumentation_info = INSTRUMENTATION_INFO,
        );
        unit = "",
        description = "Number of downloads.",
    )

    DOWNLOAD_TRACER[] =
        Tracer(; provider = tracer_provider, instrumentation_info = INSTRUMENTATION_INFO)
end

function __init__()
    init()
end

end # module


text:
# This file contains functions for the Hölder DEA model
"""
    AbstractHolderDEAModel
An abstract type representing a Holder DEA model.
"""
abstract type AbstractHolderDEAModel <: AbstractTechnicalDEAModel end

"""
    HolderL1DEAModel
An data structure representing a Höler L1 DEA model.
"""
struct HolderL1DEAModel <: AbstractHolderDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    isweighted::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    effmin::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    HolderL2DEAModel
An data structure representing a Höler L2 DEA model.
"""
struct HolderL2DEAModel <: AbstractHolderDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    isweighted::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    HolderLInfDEAModel
An data structure representing a Höler LInf DEA model.
"""
struct HolderLInfDEAModel <: AbstractHolderDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    isweighted::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaholder(X, Y; l)
Compute the Hölder distance function model using data envelopment analysis for inputs `X` and outputs `Y`,
using Hölder norm `l`.

# Hölder norm `l` specification
- `1`.
- `2`.
- `Inf`.

# Optional Arguments
- `weigt=false`:  set to `true` for weighted (weakly) Hölder distance function.
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

julia> deaholder(X, Y, l = 1, rts = :VRS)
Hölder L1 DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Orientation = Graph; Returns to Scale = VRS
────────────────────────────────────────────────
   efficiency  minimum      slackX1      slackY1
────────────────────────────────────────────────
1         0.0       X1  0.0          0.0
2         0.0       X1  0.0          0.0
3         0.0       X1  0.0          0.0
4         0.0       X1  0.0          0.0
5         3.0       X1  0.0          1.01506e-15
6         2.0       Y1  2.0          0.0
7         0.0       Y1  2.0          0.0
8         6.0       Y1  1.77636e-15  0.0
────────────────────────────────────────────────
```
"""
function deaholder(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    l::Union{Int64,Float64}, weight::Bool = false,
    orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::AbstractHolderDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    if (orient != :Input) && (orient != :Output) && (orient != :Graph)
        throw(ArgumentError("`orient` must be :Input or :Output"));
    end

    # Check l
    if l == 1
        return deaholderl1(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
    elseif l == 2
        return deaholderl2(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
    elseif l == Inf
        return deaholderlinf(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
    else
        throw(ArgumentError("`l` must be 1, 2 or Inf"));
    end

end

function deaholderl1(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    weight::Bool = false,
    orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderL1DEAModel

    # Get parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    # Compute DDF efficiency with 1 for each input and output
    n = nx
    nref = nrefx

    effXY = fill(Inf, (n,m+s))
    lambdaXY = Array{SparseMatrixCSC{Float64, Int64}}(undef, m + s)

    if slack == true
        slackXms = Array{Matrix}(undef, m + s)
        slackYms = Array{Matrix}(undef, m + s)
    end

    if orient == :Graph || orient == :Input
        for i = 1:m
            Gxi = zeros(n,m)
            if weight
                Gxi[:,i] .= X[:,i]
            else
                Gxi[:,i] .= 1
            end

            tempdea = deaddf(X, Y, Gx = Gxi, Gy = :Zeros, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
            effXY[:,i] = efficiency(tempdea)
            lambdaXY[i] = tempdea.lambda

            if slack == true
                slackXms[i] = slacks(tempdea, :X)
                slackYms[i] = slacks(tempdea, :Y)
            end
        end
    end

    if orient == :Graph || orient == :Output
        for i = 1:s
            Gyi = zeros(n,s)
            if weight
                Gyi[:,i] .= Y[:,i]
            else
                Gyi[:,i] .= 1
            end

            tempdea = deaddf(X, Y, Gx = :Zeros, Gy = Gyi, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
            effXY[:,i + m] = efficiency(tempdea)
            lambdaXY[i + m] = tempdea.lambda

            if slack == true
                slackXms[i + m] = slacks(tempdea, :X)
                slackYms[i + m] = slacks(tempdea, :Y)
            end
        end
    end

    # Obtain the minimum efficiency scores
    effi = zeros(n)
    effmin = fill(0, n)
    lambda = spzeros(n, nref)

    if slack == true
        slackX = zeros(n, m)
        slackY = zeros(n, s)
    else
        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    for i = 1:n
        effi[i], effmin[i] = findmin(effXY[i,:])
        lambda[i, :] = lambdaXY[effmin[i]][i,:]

        if slack == true
            slackX[i,:] = slackXms[effmin[i]][i,:]
            slackY[i,:] = slackYms[effmin[i]][i,:]
        end
    end

    # Get X and Y targets
    Xtarget = convert.(Float64, X[:,:])
    Ytarget = convert.(Float64, Y[:,:])
    for i = 1:n
        if effmin[i] <= m
            if weight
                Xtarget[i,effmin[i]] = X[i,effmin[i]] - effi[i] * X[i,effmin[i]]
            else
                Xtarget[i,effmin[i]] = X[i,effmin[i]] - effi[i]
            end
        else
            if weight
                Ytarget[i,effmin[i] - m] = Y[i,effmin[i] - m] + effi[i] * Y[i,effmin[i] - m]
            else
                Ytarget[i,effmin[i] - m] = Y[i,effmin[i] - m] + effi[i]
            end
        end
    end

    if slack == true
        Xtarget = Xtarget - slackX
        Ytarget = Ytarget + slackY
    end

    return HolderL1DEAModel(n, m, s, orient, rts, weight, names,  effi, effmin, slackX, slackY, lambda, Xtarget, Ytarget)

end

function deaholderl2(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    weight::Bool = false,
    orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderL2DEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if rts != :VRS
        throw(ArgumentError("`rts` must be :VRS if l = 2"));
    end

    # Default optimizer
    if optimizer === nothing
        throw(ArgumentError("An optimizer that supports SOS constraints must be specified for Hölder l = 2"));
    end

    noSOS = false

    # Compute efficiency for each DMU
    n = nx
    nref = n

    effi = zeros(n)
    lambdaeff = spzeros(n, nref)
    Xtarget = convert.(Float64, X[:,:])
    Ytarget = convert.(Float64, Y[:,:])
    slackX = zeros(n, m)
    slackY = zeros(n, s)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        if orient == :Input || orient == :Graph
            @variable(deamodel, Xeff[1:m] >= 0)
            set_start_value.(Xeff, x0)
        end
        if orient == :Output || orient == :Graph
            @variable(deamodel, Yeff[1:s] >= 0)
            set_start_value.(Yeff, y0)
        end

        if weight
            if orient == :Graph
                @objective(deamodel, Min, sum(((x0[j] - Xeff[j]) / x0[j])^2 for j in 1:m) + sum(((Yeff[j] - y0[j]) / y0[j])^2 for j in 1:s)  )
            elseif orient == :Input
                @objective(deamodel, Min, sum(((x0[j] - Xeff[j]) / x0[j])^2 for j in 1:m) )
            elseif orient == :Output
                @objective(deamodel, Min, sum(((Yeff[j] - y0[j]) / y0[j])^2 for j in 1:s)  )
            end
        else
            if orient == :Graph
                @objective(deamodel, Min, sum((x0[j] - Xeff[j])^2 for j in 1:m) + sum((Yeff[j] - y0[j])^2 for j in 1:s)  )
            elseif orient == :Input
                @objective(deamodel, Min, sum((x0[j] - Xeff[j])^2 for j in 1:m) )
            elseif orient == :Output
                @objective(deamodel, Min, sum((Yeff[j] - y0[j])^2 for j in 1:s)  )
            end
        end

        @variable(deamodel, gammatau[1:2, 1:nref] >= 0)
        @variable(deamodel, vsX[1:2,1:m] >= 0)
        @variable(deamodel, musY[1:2,1:s] >= 0)
        @variable(deamodel, psi)

        set_start_value.(gammatau[1,i], 1)

        if orient == :Input || orient == :Graph
            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * gammatau[1,t] for t in 1:nref) == Xeff[j] - vsX[2,j])
        elseif orient == :Output
            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * gammatau[1,t] for t in 1:nref) == x0[j] - vsX[2,j])
        end
        if orient == :Output || orient == :Graph
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * gammatau[1,t] for t in 1:nref) == Yeff[j] + musY[2,j])
        elseif orient == :Input
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * gammatau[1,t] for t in 1:nref) == y0[j] + musY[2,j])
        end

        @constraint(deamodel, sum(gammatau[1,:]) == 1)

        if orient == :Graph
            @constraint(deamodel, sum(vsX[1,t] for t in 1:m) + sum(musY[1,t] for t in 1:s) == 1)
        elseif orient == :Input
            @constraint(deamodel, sum(vsX[1,t] for t in 1:m) == 1)
        elseif orient == :Output
            @constraint(deamodel, sum(musY[1,t] for t in 1:s) == 1)
        end

        @constraint(deamodel, [j in 1:nref], sum(vsX[1,t] * Xref[j,t] for t in 1:m) - sum(musY[1,t] * Yref[j,t] for t in 1:s) + psi - gammatau[2,j] == 0 )

        # SOS constraints
        try
            @constraint(deamodel, [j in 1:nref], gammatau[:,j] in SOS1())
        catch
            @constraint(deamodel, [j in 1:nref], gammatau[1,j] * gammatau[2,j] == 0)
            noSOS = true
        end

        try
            @constraint(deamodel, [j in 1:m], vsX[:,j] in SOS1())
        catch
            @constraint(deamodel, [j in 1:m], vsX[1,j] * vsX[2,j] == 0)
            noSOS = true
        end

        try
            @constraint(deamodel, [j in 1:s], musY[:,j] in SOS1())
        catch
            @constraint(deamodel, [j in 1:s], musY[1,j] * musY[2,j] == 0)
            noSOS = true
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i] = JuMP.objective_value(deamodel)
        if effi[i] < 0
            # If objective function is negative change to 0, to avoid error when calculatin the square root.
            @warn ("DMU $i negative value in the objective function: $(effi[i]). Replaced to 0")
            effi[i] = 0
        end
        effi[i] = sqrt(effi[i] )
        lambdaeff[i,:] = JuMP.value.(gammatau[1,:])

        if orient == :Input || orient == :Graph
            Xtarget[i,:]  = JuMP.value.(Xeff)
        end
        if orient == :Output || orient == :Graph
            Ytarget[i,:]  = JuMP.value.(Yeff)
        end

        slackX[i,:] = JuMP.value.(vsX[2,:])
        slackY[i,:] = JuMP.value.(musY[2,:])

        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    # Add slacks to X and Y targets
    Xtarget = Xtarget - slackX
    Ytarget = Ytarget + slackY

    # Display warning if model not solved with SOS constraints
    if noSOS
        @warn ("Model solved with $(optimizer.optimizer) is innacuate. Use a solver that supports SOS1 constraints.")
    end

    return HolderL2DEAModel(n, m, s, orient, rts, weight, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function deaholderlinf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    weight::Bool = true,
    orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderLInfDEAModel

    # Get parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    n = nx

    # Buil direction based on orientation
    if orient == :Input
        if weight
            Gx = X
        else
            Gx = :Ones
        end
        Gy = :Zeros
    elseif orient == :Output
        Gx = :Zeros
        if weight
            Gy = Y
        else
            Gy = :Ones
        end
    elseif orient == :Graph
        if weight
            Gx = X
            Gy = Y
        else
            Gx = :Ones
            Gy = :Ones
        end
    end

    # Get DDF technical efficiency
    tempdea = deaddf(X, Y, Gx = Gx, Gy = Gy, rts = rts, slack = slack, Xref = Xref, Yref = Yref, optimizer = optimizer)
    effi = efficiency(tempdea)
    lambda = peersmatrix(tempdea)

    slackX = slacks(tempdea, :X)
    slackY = slacks(tempdea, :Y)

    Xtarget = targets(tempdea, :X)
    Ytarget = targets(tempdea, :Y)

    return HolderLInfDEAModel(n, m, s, orient, rts, weight, names, effi, slackX, slackY, lambda, Xtarget, Ytarget)

end

function efficiency(model::HolderL1DEAModel, type::Symbol)::Vector

    if (type == :min) return model.effmin end

    throw(ArgumentError("$(typeof(model)) with orienation has no efficiency $(type)"));

end

function Base.show(io::IO, x::HolderL1DEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    isweighted = x.isweighted
    dmunames = names(x)

    eff = efficiency(x)
    effmin = x.effmin
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    effmincol = Array{String}(undef,n)
    for i = 1:n
        if effmin[i] <= m
            effmincol[i] = "X" * string(Int(effmin[i]))
        else
            effmincol[i] = "Y" * string(Int(effmin[i] - m))
        end
    end

    if !compact
        print(io, "Hölder L1 DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if isweighted
            print(io, "Weighted (weakly) Hölder distance function \n")
        end

        if hasslacks == true
            show(io, CoefTable(hcat(eff, effmincol, slackX, slackY), ["efficiency"; "minimum"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff, effmincol), ["efficiency"; "minimum"], dmunames))
        end
    end

end

function Base.show(io::IO, x::HolderL2DEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    isweighted = x.isweighted
    dmunames = names(x)

    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Hölder L2 DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if isweighted
            print(io, "Weighted (weakly) Hölder distance function \n")
        end

        show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ; ["slackY$i" for i in 1:s ]], dmunames))
    end

end

function Base.show(io::IO, x::HolderLInfDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    isweighted = x.isweighted
    dmunames = names(x)

    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Hölder LInf DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if isweighted
            print(io, "Weighted (weakly) Hölder distance function \n")
        end

        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], dmunames))
        end
    end

end


text:
g = Multigraph([0 2 0; 2 0 1; 0 1 0])
collect(edges(g))
ps = [i//4 for i = 1:3]
v_t = [SpiderType.Z, SpiderType.Z, SpiderType.X]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:f}(), zxd)
rewrite!(Rule{:f}(), zxd, matches)
@test sort!(spiders(zxd)) == [1, 3]
@test phase(zxd, 1) == phase(zxd, 3) == 3//4

g = Multigraph(path_graph(5))
add_edge!(g, 1, 2)
ps = [1, 0//1, 0, 0, 1]
v_t = [SpiderType.X, SpiderType.X, SpiderType.Z, SpiderType.Z, SpiderType.Z]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:i1}(), zxd)
rewrite!(Rule{:i1}(), zxd, matches)
@test nv(zxd) == 3 && ne(zxd, count_mul = true) == 3 && ne(zxd) == 2

g = Multigraph([0 2 0; 2 0 1; 0 1 0])
ps = [i//4 for i = 1:3]
v_t = [SpiderType.X, SpiderType.X, SpiderType.Z]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:h}(), zxd)
rewrite!(Rule{:h}(), zxd, matches)
@test nv(zxd) == 8 && ne(zxd) == 8

matches = match(Rule{:i2}(), zxd)
rewrite!(Rule{:i2}(), zxd, matches)
@test nv(zxd) == 4 && ne(zxd, count_mul = true) == 4 && ne(zxd) == 3

g = Multigraph(6)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
add_edge!(g, 3, 4)
add_edge!(g, 3, 5)
add_edge!(g, 3, 6)
ps = [0, 1, 1//2, 0, 0, 0]
v_t = [SpiderType.In, SpiderType.X, SpiderType.Z, SpiderType.Out, SpiderType.Out, SpiderType.Out]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:pi}(), zxd)
rewrite!(Rule{:pi}(), zxd, matches)
@test nv(zxd) == 8 && ne(zxd) == 7
@test zxd.scalar == Scalar(0, 1//2)

g = Multigraph([0 2 0; 2 0 1; 0 1 0])
ps = [1, 1//2, 0]
v_t = [SpiderType.X, SpiderType.Z, SpiderType.In]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:pi}(), zxd)
rewrite!(Rule{:pi}(), zxd, matches)
@test nv(zxd) == 4 && ne(zxd) == 3 && ne(zxd, count_mul = true) == 4
@test zxd.scalar == Scalar(0, 1//2)

g = Multigraph(5)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3, 2)
add_edge!(g, 2, 4)
add_edge!(g, 2, 5)
ps = [0, 1//2, 0, 0, 0]
v_t = [SpiderType.X, SpiderType.Z, SpiderType.Out, SpiderType.Out, SpiderType.Out]
zxd = ZXDiagram(g, v_t, ps)
matches = match(Rule{:c}(), zxd)
rewrite!(Rule{:c}(), zxd, matches)
@test nv(zxd) == 7 && ne(zxd) == 4
@test zxd.scalar == Scalar(-3, 0//1)

g = Multigraph(6)
add_edge!(g, 1, 3)
add_edge!(g, 2, 4)
add_edge!(g, 3, 4)
add_edge!(g, 3, 5)
add_edge!(g, 4, 6)
ps = [0//1 for i = 1:6]
v_t = [SpiderType.In, SpiderType.In, SpiderType.X, SpiderType.Z, SpiderType.Out, SpiderType.Out]
layout = ZXCalculus.ZXLayout(2, Dict(zip(1:6, [1//1, 2, 1, 2, 1, 2])), Dict(zip(1:6, [1//1, 1, 2, 2, 3, 3])))
zxd = ZXDiagram(g, v_t, ps, layout)
matches = match(Rule{:b}(), zxd)
rewrite!(Rule{:b}(), zxd, matches)
@test nv(zxd) == 8 && ne(zxd) == 8
@test zxd.scalar == Scalar(1, 0//1)

g = Multigraph(9)
for e in [[2,6],[3,7],[4,8],[5,9]]
    add_edge!(g, e[1], e[2])
end
ps = [1//2, 0, 1//4, 1//2, 3//4, 0, 0, 0, 0]
st = [SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.In, SpiderType.In, SpiderType.Out, SpiderType.Out]
zxg = ZXGraph(ZXDiagram(g, st, ps))
for e in [[1,2],[1,3],[1,4],[1,5],[2,3]]
    add_edge!(zxg, e[1], e[2])
end
replace!(Rule{:lc}(), zxg)
@test !has_edge(zxg, 2, 3) && ne(zxg) == 9
@test phase(zxg, 2) == 3//2 && phase(zxg, 3) == 7//4 && phase(zxg, 4) == 0//1 && phase(zxg, 5) == 1//4

g = Multigraph(14)
for e in [[3,9],[4,10],[5,11],[6,12],[7,13],[8,14]]
    add_edge!(g, e[1], e[2])
end
ps = [1//1, 0, 1//4, 1//2, 3//4, 1, 5//4, 3//2, 0, 0, 0, 0, 0, 0]
st = [SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.In, SpiderType.Out, SpiderType.In, SpiderType.Out, SpiderType.In, SpiderType.Out]
zxg = ZXGraph(ZXDiagram(g, st, ps))
for e in [[1,2],[1,3],[1,4],[1,5],[1,6],[2,5],[2,6],[2,7],[2,8]]
    add_edge!(zxg, e[1], e[2])
end

replace!(Rule{:p1}(), zxg)
@test !has_edge(zxg, 3, 4) && !has_edge(zxg, 5, 6) && !has_edge(zxg, 7, 8)
@test nv(zxg) == 12 && ne(zxg) == 18
@test phase(zxg, 3) == 1//4 && phase(zxg, 4) == 1//2 && phase(zxg, 5) == 3//4 && phase(zxg, 6) == 1//1 && phase(zxg, 7) == 1//4 && phase(zxg, 8) == 1//2

g = Multigraph(6)
for e in [[2,6]]
    add_edge!(g, e[1], e[2])
end
ps = [1//1, 1//4, 1//2, 3//4, 1, 0]
st = [SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.In]
zxg = ZXGraph(ZXDiagram(g, st, ps))
for e in [[1,2],[2,3],[1,4],[1,5]]
    add_edge!(zxg, e[1], e[2])
end

@test length(match(Rule{:p1}(), zxg)) == 1
replace!(Rule{:pab}(), zxg)
@test nv(zxg) == 7 && ne(zxg) == 6

g = Multigraph(14)
for e in [[3,9],[4,10],[5,11],[6,12],[7,13],[8,14]]
    add_edge!(g, e[1], e[2])
end
ps = [1//1, 1//4, 1//4, 1//2, 3//4, 1, 5//4, 3//2, 0, 0, 0, 0, 0, 0]
st = [SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z,
    SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.In, SpiderType.Out,
    SpiderType.In, SpiderType.Out, SpiderType.In, SpiderType.Out]
zxg = ZXGraph(ZXDiagram(g, st, ps))
for e in [[1,2],[1,3],[1,4],[1,5],[1,6],[2,5],[2,6],[2,7],[2,8]]
    add_edge!(zxg, e[1], e[2])
end
match(Rule{:p2}(), zxg)
replace!(Rule{:p2}(), zxg)
@test zxg.phase_ids[15] == (2, -1)

g = Multigraph(15)
for e in [[3,9],[4,10],[5,11],[6,12],[7,13],[8,14],[2,15]]
    add_edge!(g, e[1], e[2])
end
ps = [1//1, 1//4, 1//2, 1//2, 3//2, 1, 1//2, 3//2, 0, 0, 0, 0, 0, 0, 0, 0]
st = [SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.Z,
    SpiderType.Z, SpiderType.Z, SpiderType.Z, SpiderType.In, SpiderType.Out,
    SpiderType.In, SpiderType.Out, SpiderType.In, SpiderType.Out, SpiderType.Out]
zxg = ZXGraph(ZXDiagram(g, st, ps))
for e in [[1,2],[1,3],[1,4],[1,5],[1,6],[2,5],[2,6],[2,7],[2,8]]
    add_edge!(zxg, e[1], e[2])
end
replace!(Rule{:p3}(), zxg)

@test nv(zxg) == 16 && ne(zxg) == 28
@test ZXCalculus.is_hadamard(zxg, 2, 15) && ZXCalculus.is_hadamard(zxg, 1, 16)


text:
module MoreAnalysis

export bound_variable_connections_graph

import ..Pluto
import ..Pluto: Cell, Notebook, NotebookTopology, ExpressionExplorer

"Find all subexpressions of the form `@bind symbol something`, and extract the `symbol`s."
function find_bound_variables(expr)
	found = Set{Symbol}()
	_find_bound_variables!(found, ExpressionExplorer.maybe_macroexpand(expr; recursive=true, expand_bind=false))
	found
end

function _find_bound_variables!(found::Set{Symbol}, expr::Expr)
	if expr.head === :macrocall && expr.args[1] === Symbol("@bind") && length(expr.args) == 4 && expr.args[3] isa Symbol
		push!(found, expr.args[3])
		_find_bound_variables!(found, expr.args[4])
    elseif expr.args === :quote
        found
    else
		for a in expr.args
			_find_bound_variables!(found, a)
		end
	end
end

function _find_bound_variables!(found::Set{Symbol}, expr::Any) end

"Return the given cells, and all cells that depend on them (recursively)."
function downstream_recursive(notebook::Notebook, topology::NotebookTopology, from::Union{Vector{Cell},Set{Cell}})::Set{Cell}
    found = Set{Cell}(copy(from))
    _downstream_recursive!(found, notebook, topology, from)
    found
end

function _downstream_recursive!(found::Set{Cell}, notebook::Notebook, topology::NotebookTopology, from::Vector{Cell})::Nothing
    for cell in from
        one_down = Pluto.where_referenced(notebook, topology, cell)
        for next in one_down
            if next ∉ found
                push!(found, next)
                _downstream_recursive!(found, notebook, topology, Cell[next])
            end
        end
    end
end

"Return all cells that are depended upon by any of the given cells."
function upstream_recursive(notebook::Notebook, topology::NotebookTopology, from::Union{Vector{Cell},Set{Cell}})::Set{Cell}
    found = Set{Cell}(copy(from))
    _upstream_recursive!(found, notebook, topology, from)
    found
end

function _upstream_recursive!(found::Set{Cell}, notebook::Notebook, topology::NotebookTopology, from::Vector{Cell})::Nothing
    for cell in from
        references = topology.nodes[cell].references
        for upstream in Pluto.where_assigned(notebook, topology, references)
            if upstream ∉ found
                push!(found, upstream)
                _upstream_recursive!(found, notebook, topology, Cell[upstream])
            end
        end
    end
end

"All cells that can affect the outcome of changing the given variable."
function codependents(notebook::Notebook, topology::NotebookTopology, var::Symbol)::Set{Cell}
    assigned_in = filter(notebook.cells) do cell
        var ∈ topology.nodes[cell].definitions
    end

    downstream = collect(downstream_recursive(notebook, topology, assigned_in))

    downupstream = upstream_recursive(notebook, topology, downstream)
end

"Return a `Dict{Symbol,Vector{Symbol}}` where the _keys_ are the bound variables of the notebook.

For each key (a bound symbol), the value is the list of (other) bound variables whose values need to be known to compute the result of setting the bond."
function bound_variable_connections_graph(notebook::Notebook)::Dict{Symbol,Vector{Symbol}}
    topology = notebook.topology
    bound_variables = union(map(notebook.cells) do cell
        find_bound_variables(topology.codes[cell].parsedcode)
    end...)
    Dict{Symbol,Vector{Symbol}}(
        var => let
            cells = codependents(notebook, topology, var)
            defined_there = union!(Set{Symbol}(), (topology.nodes[c].definitions for c in cells)...)
            # Set([var]) ∪
            collect((defined_there ∩ bound_variables))
        end
        for var in bound_variables
    )
end
end


text:
using IterativeSolvers
using FactCheck

facts("Randomized linear algebra") do
    srand(1)
    m = n = 100
    B = randn(m, n)
    nB = norm(B)
    p = 1e-5 #Probability of failure

    for j=1:n
        @fact rnorm(B, 2j, p) --> greater_than_or_equal(nB)
        @fact rnorms(B, j, p) --> greater_than_or_equal(nB)
    end

    A = B*B'
    k = 1
    l, u = reigmin(A, k, p)
    @fact l <= eigmin(A) <= u --> true

    l, u = reigmax(A, k, p)
    @fact l <= eigmax(A) <= u --> true

    l, u = rcond(A, k, p)
    @fact l <= cond(A) <= u --> true
end


text:
struct MetaBall{T}
    i_center::T
    j_center::T
    radius::T
end

function generate_field!(field::AbstractMatrix, metaballs::Array{MetaBall{T}, 1}) where {T}
    height, width = size(field)
    for j in 1:width
        for i in 1:height
            for metaball in metaballs
                field[i, j] += (metaball.radius ^ 2) / ((i - metaball.i_center) ^ 2 + (j - metaball.j_center) ^ 2)
            end
        end
    end

    return nothing
end