text:
export BrickletOutdoorWeatherStationIdentifiersLowLevel
struct BrickletOutdoorWeatherStationIdentifiersLowLevel
    identifiers_length::Integer
    identifiers_chunk_offset::Integer
    identifiers_chunk_data::Vector{Integer}
end

export BrickletOutdoorWeatherSensorIdentifiersLowLevel
struct BrickletOutdoorWeatherSensorIdentifiersLowLevel
    identifiers_length::Integer
    identifiers_chunk_offset::Integer
    identifiers_chunk_data::Vector{Integer}
end

export BrickletOutdoorWeatherStationData
struct BrickletOutdoorWeatherStationData
    temperature::Integer
    humidity::Integer
    wind_speed::Integer
    gust_speed::Integer
    rain::Integer
    wind_direction::Integer
    battery_low::Bool
    last_change::Integer
end

export BrickletOutdoorWeatherSensorData
struct BrickletOutdoorWeatherSensorData
    temperature::Integer
    humidity::Integer
    last_change::Integer
end

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

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

export BrickletOutdoorWeather
"""
433MHz receiver for outdoor weather station
"""
mutable struct BrickletOutdoorWeather <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_station_identifiers_low_level
"""
    $(SIGNATURES)

Returns the identifiers (number between 0 and 255) of all `stations
<https://www.tinkerforge.com/en/shop/accessories/sensors/outdoor-weather-station-ws-6147.html>`__
that have been seen since the startup of the Bricklet.

Each station gives itself a random identifier on first startup.

Since firmware version 2.0.2 a station is removed from the list if no data was received for
12 hours.
"""
function get_station_identifiers_low_level(device::BrickletOutdoorWeather)
    return device.deviceInternal.get_station_identifiers_low_level()
end

export get_sensor_identifiers_low_level
"""
    $(SIGNATURES)

Returns the identifiers (number between 0 and 255) of all `sensors
<https://www.tinkerforge.com/en/shop/accessories/sensors/temperature-humidity-sensor-th-6148.html>`__
that have been seen since the startup of the Bricklet.

Each sensor gives itself a random identifier on first startup.

Since firmware version 2.0.2 a sensor is removed from the list if no data was received for
12 hours.
"""
function get_sensor_identifiers_low_level(device::BrickletOutdoorWeather)
    return device.deviceInternal.get_sensor_identifiers_low_level()
end

export get_station_data
"""
    $(SIGNATURES)

Returns the last received data for a station with the given identifier.
Call :func:`Get Station Identifiers` for a list of all available identifiers.

The return values are:

* Temperature,
* Humidity,
* Wind Speed,
* Gust Speed,
* Rain Fall (accumulated since station power-up),
* Wind Direction,
* Battery Low (true if battery is low) and
* Last Change (seconds since the reception of this data).
"""
function get_station_data(device::BrickletOutdoorWeather, identifier)
    return device.deviceInternal.get_station_data(identifier)
end

export get_sensor_data
"""
    $(SIGNATURES)

Returns the last measured data for a sensor with the given identifier.
Call :func:`Get Sensor Identifiers` for a list of all available identifiers.

The return values are:

* Temperature,
* Humidity and
* Last Change (seconds since the last reception of data).
"""
function get_sensor_data(device::BrickletOutdoorWeather, identifier)
    return device.deviceInternal.get_sensor_data(identifier)
end

export set_station_callback_configuration
"""
    $(SIGNATURES)

Turns callback for station data on or off.
"""
function set_station_callback_configuration(device::BrickletOutdoorWeather, enable_callback)
    device.deviceInternal.set_station_callback_configuration(enable_callback)
end

export get_station_callback_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Station Callback Configuration`.
"""
function get_station_callback_configuration(device::BrickletOutdoorWeather)
    return device.deviceInternal.get_station_callback_configuration()
end

export set_sensor_callback_configuration
"""
    $(SIGNATURES)

Turns callback for sensor data on or off.
"""
function set_sensor_callback_configuration(device::BrickletOutdoorWeather, enable_callback)
    device.deviceInternal.set_sensor_callback_configuration(enable_callback)
end

export get_sensor_callback_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Sensor Callback Configuration`.
"""
function get_sensor_callback_configuration(device::BrickletOutdoorWeather)
    return device.deviceInternal.get_sensor_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

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

The errors are divided into

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

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather, 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::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather, 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::BrickletOutdoorWeather, 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::BrickletOutdoorWeather, 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::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather, 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::BrickletOutdoorWeather)
    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::BrickletOutdoorWeather)
    return device.deviceInternal.get_identity()
end

export get_station_identifiers
"""
Returns the identifiers (number between 0 and 255) of all `stations
<https://www.tinkerforge.com/en/shop/accessories/sensors/outdoor-weather-station-ws-6147.html>`__
that have been seen since the startup of the Bricklet.

Each station gives itself a random identifier on first startup.

Since firmware version 2.0.2 a station is removed from the list if no data was received for
12 hours.
"""
function get_station_identifiers(device::BrickletOutdoorWeather)

    lock(device.stream_lock) do
        ret = get_station_identifiers_low_level(device, )
            identifiers_length = ret.identifiers_length
        identifiers_out_of_sync = ret.identifiers_chunk_offset != 0

        identifiers_data = ret.identifiers_chunk_data

        while !identifiers_out_of_sync && length(identifiers_data) < identifiers_length
            ret = get_station_identifiers_low_level(device, )
                identifiers_length = ret.identifiers_length
            identifiers_out_of_sync = ret.identifiers_chunk_offset != length(identifiers_data)
            identifiers_data += ret.identifiers_chunk_data
        end

        if identifiers_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.identifiers_chunk_offset + 60 < identifiers_length
                ret = get_station_identifiers_low_level(device, )
                    identifiers_length = ret.identifiers_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Identifiers stream is out-of-sync"))
        end
    end

        return identifiers_data[:identifiers_length]
end

export get_sensor_identifiers
"""
Returns the identifiers (number between 0 and 255) of all `sensors
<https://www.tinkerforge.com/en/shop/accessories/sensors/temperature-humidity-sensor-th-6148.html>`__
that have been seen since the startup of the Bricklet.

Each sensor gives itself a random identifier on first startup.

Since firmware version 2.0.2 a sensor is removed from the list if no data was received for
12 hours.
"""
function get_sensor_identifiers(device::BrickletOutdoorWeather)

    lock(device.stream_lock) do
        ret = get_sensor_identifiers_low_level(device, )
            identifiers_length = ret.identifiers_length
        identifiers_out_of_sync = ret.identifiers_chunk_offset != 0

        identifiers_data = ret.identifiers_chunk_data

        while !identifiers_out_of_sync && length(identifiers_data) < identifiers_length
            ret = get_sensor_identifiers_low_level(device, )
                identifiers_length = ret.identifiers_length
            identifiers_out_of_sync = ret.identifiers_chunk_offset != length(identifiers_data)
            identifiers_data += ret.identifiers_chunk_data
        end

        if identifiers_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.identifiers_chunk_offset + 60 < identifiers_length
                ret = get_sensor_identifiers_low_level(device, )
                    identifiers_length = ret.identifiers_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Identifiers stream is out-of-sync"))
        end
    end

        return identifiers_data[:identifiers_length]
end

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


text:
module Remove_response_obspy
export remove_response_obspy, remove_response_obspy!
using PyCall, SeisIO

# Suppressing warning associated with PyCall
using Suppressor

"""
remove_response_obspy(S::SeisDatat)
remove Instrumental response with obspy.
Currently available only one channel.
This function will be replaced by julia original function
"""
function remove_response_obspy!(S::SeisData, xf::String; pre_filt::NTuple{4,Float64}=(0.001, 0.002, 10.0, 20.0), output::String="VEL")

    @suppress begin
    Stream          = pyimport("obspy.core.stream")["Stream"]
    Trace           = pyimport("obspy.core.trace")["Trace"]
    read_inventory  = pyimport("obspy")["read_inventory"]
    UTCDateTime     = pyimport("obspy.core.utcdatetime")["UTCDateTime"]

    trace=Trace()
    trace.data = S.x[1]
    trace.stats.sampling_rate = S.fs
    trace.stats.delta = 1.0./ S.fs
    trace.stats.starttime = UTCDateTime(string(u2d( S.t[1][1,2]*1e-6)))
    str1 = split(S.id[1], ".")
    trace.stats.network = str1[1]
    trace.stats.station=  str1[2]
    trace.stats.location= str1[3]
    trace.stats.channel=  str1[4]

    stream = Stream(traces=trace)
    #read stationXML
    inv = read_inventory(xf)
    #stresm.remove_sensitivity(inventory=inv)
    stream.remove_response(inventory=inv, output=output, pre_filt=pre_filt, plot=false)
    S.x[1] = stream.traces[1].data
    #add note
    SeisIO.note!(S, "remove_response_obspy!, pre_filt=$pre_filt")
    end
    return nothing
end

function remove_response_obspy(S::SeisData; pre_filt::NTuple{4,Float64}=(0.001, 0.002, 10.0, 20.0), output::String="VEL")
  U = deepcopy(S)
  remove_response_obspy!(U, pre_filt=pre_filt, output=output)
  return U
end

end

ext: jl

max_stars_repo_name: kura-okubo/SeisDownloadwJulia

max_issues_repo_name: kura-okubo/SeisDownloadwJulia

max_forks_repo_name: kura-okubo/SeisDownloadwJulia


text:
module PCCA

include("pccap.jl")
include("schur.jl")

export pcca

end


text:
function holtslag_psis(xi)
if xi>0
    y=1+5*xi;
else
    y=(1-16*xi)^(-1/2);
end
return y
end


text:
include("helpers.jl")
include("../src/UnityGymWrapper.jl")
using .UnityGymWrapper
using Plots
using Flux
using Flux: Optimise.update!

function run!(model, opt, env;
        epochs=30, steps=300, maxt=100, batchsize=128,
        trainiters=100, ϵ=0.07, γ=99f-2, τ=1f-2)
    rewards = zeros(Float32, epochs)
    target = deepcopy(model)
    memory = ReplayMemory{eltype(env.observationspace), eltype(env.actionspace), Float32}(
        length(env.observationspace), length(env.actionspace), 5*steps, batchsize
    )
    ps = params(model)

    for i in 1:epochs
        j = 0
        while j < steps
            s1 = reset!(env)
            for _ in 1:(min(maxt, steps-j))
                a = ifelse(rand() <= ϵ, sample(env.actionspace), Flux.onecold(model(s1)))
                s2, r, done, _ = step!(env, a)
                append!(memory, s1, a, r, done)
                s1 = s2
                j += 1
                done && break
            end
        end

        for _ in 1:trainiters
            s1, a, r, s2, d = sample(memory)
            Qtar = r .+ (1 .- d) .* γ .* vec(maximum(target(s2), dims=1))
            gs = gradient(ps) do
                Qs1a = gather(model(s1), a)
                return Flux.mse(Qs1a, Qtar)
            end
            update!(opt, ps, gs)
            softupdate!(target, model, τ)
        end

        ϵ = max(ϵ * 0.99, 0.01)
        rewards[i] = test(s -> Flux.onecold(model(s)), env)
        println("$(i): r = $(rewards[i])")
    end
    return rewards
end

env = Environment("envs/Basic", logFile=pwd()*"/envs/logs/Basic.log")

model = Chain(
    Dense(length(env.observationspace), 64, relu),
    Dense(64, 64, relu),
    Dense(64, env.actionspace.actions[1])
)

r = run!(model, ADAM(), env)
close!(env)

plt = plot(r, labels="Reward")
display(plt)


text:
"""
link(net,rp)
link(net,s,rp)
link_f(net,rp)
link_f(net,s,rp)

Determine link id from relative *rp*  position and node state *s*.

The link function provides a way to access links form the set of nodes in a certain link state,
or from the set of all links. The two-argument version returns the id of the
link at poition *rp* in network *net*. The three-argument version returns the id of the link at
poition *rp* within the set of links that are in state *s*.

All version of this function run in constant time, but fast (_f) verions sacrifice some safty
checks for better performance. See [basic concepts](concepts.md) for details.

See also [adjacent](#Fastnet.adjacent)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[LinkType(1,2),LinkType(1,1),LinkType(2,2)])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,2);

julia> makelink!(net,n1,n2);

julia> makelink!(net,n2,n3);

julia> makelink!(net,n3,n1);

julia> lid=link(net,2,1);

julia> linkdst(net,lid)==n2
true

julia> linksrc(net,lid)==n1
true
```
"""
function link(net::FastNet,cls,rp)
    task="Trying to retrieve id of link number $rp in state $cls"
    checklinkstate(net::FastNet,cls,task)
    if rp<1 || rp>net.kclen[cls]
        throw(ArgumentError("$task, but position $rp does not exist in state $cls"))
    end
    return link_f(net,cls,rp)
end

function link(net::FastNet,rp)
    task="Trying to retrieve id of link number $rp"
    if rp<1 || rp>countlinks_f(net)
        throw(ArgumentError("$task, but position $rp does not exist"))
    end
    return link_f(net,rp)
end

"""
    linkcounts(net)
    linkcounts_f(net)

Return an Array containing the number of link in the vairous link states.

The elements of the array will show the counts in the same order in which the link types
were passed to the FastNet Constructor.

The time required for this function scales only with the number of link states
(it is independent of the number of links).

The alternative (_f) version of this function is identical only provided for convenience.

See also [countlinks](#Fastnet.countlinks)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[LinkType(1,2),LinkType(1,1),LinkType(2,2)])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,2);

julia> makelink!(net,n1,n2);

julia> makelink!(net,n2,n3);

julia> makelink!(net,n3,n1);

julia> linkcounts(net)
3-element Vector{Int64}:
 2
 1
 0
```
"""
function linkcounts(net::FastNet)
    linkcounts_f(net)
end

"""
    linkstate(net,kid)
    linkstate_f(net,kid)

Return the state of the link with id *kid* in network *net*.

Note that the link states are numbered in the order in which they were passed to
the FastNet Constructor.

All version of this function run in constant time, but fast (_f) verion sacrifices some safty
checks for better performance. See [basic concepts](concepts.md) for details.

See also [nodestate!](#Fastnet.nodestate!), [FastNet](#FastNet.FastNet)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[LinkType(1,2),LinkType(1,1),LinkType(2,2)])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> lnk=makelink!(net,n1,n2);

julia> linkstate(net,lnk)
2
```
"""
function linkstate(net::FastNet,kid)
    task="Trying to determine state of link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    linkstate_f(net,kid)
end

"""
    makelink!(net,src,dst)
    makelink_f!(net,src,dst)

Create a new link from node *src* to node*dst* in the network *net* and return it's id.

Worst-case performance of both versions of this function scales only with the number of tracked link states.

The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [destroylink!](#Fastnet.destrolink!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> makelink!(net,n1,n2);

julia> net
Network of 2 nodes and 1 links
```
"""
function makelink!(net::FastNet,src,dst)
    task="Trying to create link from node $src to node $dst"
    checknodeid(net::FastNet,src,task)
    checknodeexists(net::FastNet,src,task)
    checknodeid(net::FastNet,dst,task)
    checknodeexists(net::FastNet,dst,task)
    if net.kclen[net.L]<1
        throw(ArgumentError("Out of links -- Try creating a larger net"))
    end
    makelink_f!(net,src,dst)
end

"""
    destroylink!(net,kid)
    destroylink_f!(net,kid)

Destroy the link with id *kid* in network *net*.

All versions of this function run in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makelink!](#Fastnet.makelink!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,2);

julia> makelink!(net,n1,n2);

julia> net
Network of 2 nodes and 1 links

julia> lnk=randomlink(net);

julia> destroylink!(net,lnk)

julia> net
Network of 2 nodes and 0 links
```
"""
function destroylink!(net::FastNet,kid)
    task="Trying to destroy link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    destroylink_f!(net,kid)
end

"""
    linksrc!(net,kid)
    linksrc_f!(net,kid)

Return the id of the node at the source of link *kid* in *net*.

All versions of this function run in constant time.
The fast (_f) verion sacrifices some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makelink!](#Fastnet.makelink!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,2);

julia> k1=makelink!(net,n1,n2);

julia> linksrc(net,k1)==n1
true
```
"""
function linksrc(net::FastNet,kid)
    task="Trying to retrieve source of link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    linksrc_f(net,kid)
end

"""
    linkdst!(net,kid)
    linkdst_f!(net,kid)

Return the id of the node at the destination of link *kid* in *net*.

All versions of this function run in constant time.
The fast (_f) verion sacrifices some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makelink!](#Fastnet.makelink!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,2);

julia> k1=makelink!(net,n1,n2);

julia> linkdst(net,k1)==n2
true
```
"""
function linkdst(net::FastNet,kid)
    task="Trying to retrieve destination of link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    linkdst_f(net,kid)
end

"""
    nextlinkout!(net,kid)
    nextlinkout_f!(net,kid)

Get the id of the next outgoing link from the node at the source of link *kid* in *net*.

This function can be used to iterate over the outgoing links of a node. If *kid* is the nodes
last link the return value is zero.

All versions of this function run in constant time.
The fast (_f) verion sacrifices some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [firstlinkout](#Fastnet.firstlinkout), [nextlinkin](#Fastnet.nextlinkin)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,1);

julia> n4=makenode!(net,1);

julia> makelink!(net,n1,n2)
1

julia> makelink!(net,n1,n3)
2

julia> makelink!(net,n1,n4)
3

julia> k=firstlinkout(net,n1);

julia> while k!=0
           println(k)
           k=nextlinkout(net,k)
           end
3
2
1
```
"""
function nextlinkout(net::FastNet,kid)
    task="Trying to retrive next outgoing link from the same node as link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    nextlinkout_f(net,kid)
end

"""
    nextlinkin!(net,kid)
    nextlinkin_f!(net,kid)

Get the id of the next incoming link to the node at the destination of link *kid* in *net*.

This function can be used to iterate over the incoming links of a node. If *kid* is the nodes
last link the return value is zero.

All versions of this function run in constant time.
The fast (_f) verion sacrifices some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [firstlinkin](#Fastnet.firstlinkin), [nextlinkout](#Fastnet.nextlinkout)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> n3=makenode!(net,1);

julia> n4=makenode!(net,1);

julia> makelink!(net,n2,n1)
1

julia> makelink!(net,n3,n1)
2

julia> makelink!(net,n4,n1)
3

julia> k=firstlinkin(net,n1);

julia> while k!=0
           println(k)
           k=nextlinkin(net,k)
           end
3
2
1
```
"""
function nextlinkin(net::FastNet,kid)
    task="Trying to retrive next incoming link to the same node as link $kid"
    checklinkid(net::FastNet,kid,task)
    checklinkexists(net::FastNet,kid,task)
    nextlinkin_f(net,kid)
end

"""
    countlinks(net)
    countlinks(net,s)
    countlinks_f(net)
    countlinks_f(net,s)

Count the links in state *s*, or, if no state is provided, in the entire network.

Instead of the state *s* also an Array or Tuple of states can be passed.
In this case the total number of nodes in all of the listed states is returned.

The links in a sincle class or the entire network are counted in constant time.
For the tuples or array arguments the performance scales with the number of elements in the
Tulps/Array.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

If performance is critical use this function rather than linkcounts.

See also [linkcounts](#Fastnet.linkcounts)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[LinkType(1,2,1),LinkType(2,2,2)])
Network of 0 nodes and 0 links

julia> makenodes!(net,20,1)

julia> makenodes!(net,20,2)

julia> for i=1:20
         from=node(net,1,i)
         to=node(net,2,i)
         makelink!(net,from,to)
       end

julia> makelink!(net,node(net,2,1),node(net,2,2));

julia> countlinks(net)
21

julia> countlinks(net,1)
20

julia> countlinks(net,2)
1

julia> countlinks(net,[1,2])
21
```
"""
function countlinks(net::FastNet)
    countlinks_f(net)
end

function countlinks(net::FastNet,cls::Integer)
    task="Trying to count links in state $cls"
    checklinkstate(net::FastNet,cls,task)
    countlinks_f(net,cls)
end

function countlinks(net::FastNet,cls::Union{Vector,Tuple})
    ret=0
    for c in cls
        ret+=countlinks(net,c)
    end
    ret
end

"""
    linkexists(net,kid)
    linkexists_f(net,kid)

Return true node with id *kid* exists in *net*, false otherwise.

This function runs in constant time.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

See also [makelink!](#Fastnet.makelink!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,1);

julia> linkexists(net,7)
false

julia> k=makelink!(net,n1,n2);

julia> linkexists(net,k)
true
```
"""
function linkexists(net::FastNet,kid)
    task="Trying to determine if link $kid is in the network"
    checklinkid(net::FastNet,kid,task)
    linkexists_f(net,kid)
end

"""
    randomlink(net)
    randomlink(net,s)
    randomlink_f(net)
    randomlink_f(net,s)

Return the id of a random link drawn from *net*.

If the second argument *s* is not provided the link will be drawn uniformly from
all links in the network. If *s* is an integer then the link will be drawn uniformly
from the links in state *s*. If *s* is an Array or Tuple of Ints then the link will be
drawn uniformly from the links in the states listed.

This function runs in constant time if *s* is integer or omitted. If *s* is an Array or Tuple the
worst case performance scales only with the number of tracked link states.
The fast (_f) verions sacrifice some safty checks for better performance.
See [basic concepts](concepts.md) for details.

The safe versions of this function will throw an ArgumentError with an informative error message
when trying to pick a link from an empty set. With the fast (_f) version, trying to pick a link from
an empty set will also result in an ArgumentError being thrown, but in this case the message will be
something like "Range must be non-empty".

See also [link](#Fastnet.link)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[LinkType(2,2,2)])
Network of 0 nodes and 0 links

julia> randomgraph!(net);

julia> for i=1:500
         nd=randomnode(net,1)
         nodestate!(net,nd,2)
       end

julia> lnk=randomlink(net,1);

julia> linkstate(net,lnk)
1
```
"""
function randomlink(net::FastNet)
    task="Trying to select random link"
    checknonempty(net,task)
    randomlink_f(net)
end

function randomlink(net::FastNet,cls::Integer)
    task="Trying to select random link from state $cls"
    checklinkstate(net,cls,task)
    checknonempty(net,cls,task)
    randomlink_f(net,cls)
end

function randomlink(net::FastNet,cls::Union{Array,Tuple})
    task="Trying to select a random link from a set of states"
    for c in cls
        checklinkstate(net,c,task)
    end
    tot=countlinks_f(net,cls)
    if tot===0
        throw(ArgumentError("$task, but there isn't any node in any of the states."))
    end
    randomlink_f(net,cls)
end


text:
    using Distributed

    using DifferentialEquations
    using OrdinaryDiffEq
    using DiffEqBase
    using Sundials
    using Plots
    # using DiffEqGPU
    using ODEInterfaceDiffEq
    using DataFrames
    using CSV
    using Statistics
    using LinearAlgebra
    using JLD
    using StatsBase
    using Random
    using SharedArrays
    using BlackBoxOptim
    using BlackBoxOptim: num_func_evals

    include("E:\\UNI\\D_Drive\\PhD\\JULIAlang\\Generalisation_PLacExample\\Examples\\ModelsFunctions\\PLacExample_Model.jl")

    function selectObsSim(simul)

        simObs = zeros(size(simul)[1], length(Obs));
        if typeof(Obs)==Array{Int,1}
            simObs = simul[:,(Obs),1];
        else
            coun = 1;
            for j in 1:length(Obs)
               if !(occursin.("+", Obs[j])) && !(occursin.("-", Obs[j])) &&
               !(occursin.("/", Obs[j])) && !(occursin.("*", Obs[j])) &&
                !(occursin.("^", Obs[j]))
                    tind = findall(x->x==Obs[j], stName)[1];
                    simObs[:,coun] = simul[:,tind,1];
                    coun +=1;
                else #This is the complicated stuff where we need to evaluate the function written as a string
                    tmp1 = Obs[j];
                    for k in 1:length(stName)
                        if occursin(stName[k], Obs[j])
                            tmp1 = replace(tmp1, stName[k]=>string(simul[:,k,1]))
                        end
                    end
                    simObs[:,coun] = eval(Meta.parse(tmp1));
                    coun +=1;
                end
            end
        end
        return(simObs)
    end

    function parOptim(opts)
        ress = convert(SharedArray{Float64,2}, zeros(10,length(opts)));
        for i in 1:length(opts)
            resT = bboptimize(opts[string("Opt_", i)]);
            ress[:,i] = resT.archive_output.best_candidate;
        end;
        return(ress);
    end

    function loglike(dats, mes, errs)

        # Check that the dimensions of all the vectors are the same
        if (size(dats) != size(mes)) || (size(dats) != size(errs))
            println("Please, check the vectors dimensions!")
            return
        end

        # Computation of log likelihood
        nd = length(dats); # Number of time points
        llk_nd = zeros(nd);

        for i in 1:nd
            tmp1 = log(2*pi);
            tmp2 = log(errs[i]^2);
            tmp3 = ((mes[i]-dats[i])^2)/(errs[i]^2);

            llk_nd[i] = (-0.5)*(tmp1+tmp2+tmp3);

            # Check to avoid Nans (should not be required)
            if isnan(llk_nd[i])
                llk_nd[i] = 0;
            end
        end
        llk = sum(llk_nd);

        return(llk);
    end

    function objectiveMLEPLacExample(theta)

        LLK_Exps = zeros(1, Nexp);
        simobs = Array{Any,1}(undef,Nexp);
        for i in 1:Nexp
            simul = PLacExample_SolveAll(tim[i], theta, esp[i], [], ini[i], smp[i], prr[i])
            simobs[i] = selectObsSim(simul);
            llkobs = zeros(1,length(Obs));
            for j in 1:length(Obs)
                llkobs[1,j] = loglike(dataM[i][:,j], simobs[i][:,j], dataE[i][j]);
            end
            LLK_Exps[1,i] = sum(llkobs);
        end
        return(-sum(LLK_Exps))
    end

    function restructInputs(model_def, mle_def, expp)
        if model_def["nInp"] != 0
            r = convert.(Int, 1:model_def["nInp"]:(length(mle_def["uInd"][expp])));
            inputs = zeros(convert.(Int,length(mle_def["uInd"][expp])));
            for j in 1:convert.(Int,length(mle_def["uInd"][expp])/model_def["nInp"])
                for k in 0:(model_def["nInp"]-1)
                    inputs[r[j]+k] = mle_def["uInd"][expp][j,(k+1)];
                end
            end
        else
            inputs = [];
        end
        return(inputs)
    end

    function RunMLEPLacExample(model_def2, mle_def2)

        global Nexp = mle_def2["Nexp"];
        global Obs = mle_def2["Obs"];
        global stName = model_def2["stName"];

        global dataM = mle_def2["DataMean"];
        global dataE = mle_def2["DataError"];

        global ts = Array{Any,1}(undef,mle_def2["Nexp"]);
        global sp = Array{Any,1}(undef,mle_def2["Nexp"]);
        global ivss = Array{Any,1}(undef,mle_def2["Nexp"]);
        global pre = Array{Any,1}(undef,mle_def2["Nexp"]);
        global samps = Array{Any,1}(undef,mle_def2["Nexp"]);
        global inputs = Array{Any,1}(undef,mle_def2["Nexp"]);

        for i in 1:mle_def2["Nexp"]
            ts[i] = collect(0.0:round(mle_def2["finalTime"][i]))';
            sp[i] = [convert(Int,v) for v in (round.(mle_def2["switchT"][i])')];
            ivss[i] = mle_def2["y0"][i];
            if model_def2["Y0eqs"] != [] # ON inputs
                pre[i] = mle_def2["preInd"][i];
            else
                pre[i] = [];
            end
            samps[i] = convert.(Int, round.(mle_def2["tsamps"][i]));
            inputs[i] = restructInputs(model_def2, mle_def2, i);
        end

        global tim = ts;
        global esp = sp;
        global inp = inputs;
        global ini = ivss;
        global prr = pre;
        global smp = samps;

        lb = mle_def2["thetaMIN"];
        up = mle_def2["thetaMAX"];

        rang = [Array{Tuple{Int, Int}}(undef, length(lb))];
        rang = [(lb[i], up[i]) for i in 1:length(lb)];

        opts = Dict()

        convcur = Array{Any,1}(undef, 3);
        for i in 1:3
            convcur[i] = Array{Tuple{Int, Float64},1}();
            callback = oc -> push!(convcur[i], (num_func_evals(oc), best_fitness(oc)))
            opts[string("Opt_", i)] = bbsetup(objectiveMLEPLacExample; SearchRange = rang, NumDimensions = length(lb),
                Method = :adaptive_de_rand_1_bin_radiuslimited, MaxTime = 10, TraceMode = :silent, CallbackFunction = callback, CallbackInterval = 0.0);
        end

        println("----------------------------------------- OPTIMISATION STARTS! -----------------------------------------")

        @time begin
            tet= parOptim(opts);
        end

        println("----------------------------------------- OPTIMISATION ENDED -----------------------------------------")

        names = model_def2["parName"];

        alltog = Array{Dict{String,Any},1}(undef,3)
        for i in 1:(3)
            global tmp = Dict{String, Any}()
            for j in 1:length(names)
                tmp[names[j]] = tet[j,i];
            end
            alltog[i] = tmp
        end

        mle_res = Dict();
        mle_res["Theta"] = convert(Array,tet); # Best thetas results
        mle_res["convCurv"] = convcur; # Convergence curves
        mle_res["StanDict"] = alltog; # Best thetas in the structure Stan needs as initial guess.

        bcfv = zeros(length(convcur))
        for k in 1:length(convcur)
            bcfv[k] = convcur[k][end][2]
        end
        mle_res["BestCFV"] = bcfv; # Best cost function values for each run

        tpin = findfirst(x->x==minimum(bcfv), bcfv);
        mle_res["BestTheta"] = tet[:,tpin]; # Best theta value from all runs

        return(mle_res)

    end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
"""
Module: SymbolicDiff (Symbolic Operation for Arithmetic)
"""

"""
seval(f, dvar, env, cache)
Return the first derivative of expr f with respect to dvar
"""

function seval(f, dvar::Symbol)
    seval(f, dvar, globalenv, SymbolicCache())
end

function seval(f, dvar::Symbol, env::SymbolicEnv)
    seval(f, dvar, env, SymbolicCache())
end

function seval(f, dvar::Symbol, cache::SymbolicCache)
    seval(f, dvar, globalenv, cache)
end

###

function seval(f, dvar::SymbolicVariable{Tv}) where Tv
    seval(f, dvar.var)
end

function seval(f, dvar::SymbolicVariable{Tv}, env::SymbolicEnv) where Tv
    seval(f, dvar.var, env)
end

function seval(f, dvar::SymbolicVariable{Tv}, cache::SymbolicCache) where Tv
    seval(f, dvar.var, cache)
end

function seval(f, dvar::SymbolicVariable{Tv}, env::SymbolicEnv, cache::SymbolicCache) where Tv
    seval(f, dvar.var, env, cache)
end

###

function seval(f::SymbolicValue{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    Tv(0)
end

function seval(f::SymbolicVariable{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    f.var == dvar ? 1 : 0
end

function seval(f::AbstractNumberSymbolic{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    (dvar in f.params) || return 0
    get(cache, (f,dvar)) do
        retval = _eval(Val(f.op), f, dvar, env, cache)
        cache[(f,dvar)] = retval
    end
end

"""
_eval(::Val{xx}, dvar, f, env, cache)

Dispached function to evaluate the first derivative of f
"""

function _eval(::Val{:+}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    +(args...)
end

function _eval(::Val{:-}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    -(args...)
end

function _eval(::Val{:*}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, env, cache) for x = f.args]
    dargs = [seval(x, dvar, env, cache) for x = f.args]
    ret = dargs[1]
    s = args[1]
    for i = 2:length(args)
        ret *= args[i]
        ret += s * dargs[i]
        (i == length(args)) && break
        s *= args[i]
    end
    ret
end

function _eval(::Val{:/}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx,dy = [seval(x, dvar, env, cache) for x = f.args]
    (dx * y - x * dy) / y^2
end

function _eval(::Val{:^}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx,dy = [seval(x, dvar, env, cache) for x = f.args]
    x^(y-1) * (x * log(x) * dy + y * dx)
end

function _eval(::Val{:exp}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx, = [seval(x, dvar, env, cache) for x = f.args]
    exp(x) * dx
end

function _eval(::Val{:log}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx, = [seval(x, dvar, env, cache) for x = f.args]
    dx / x
end

function _eval(::Val{:sqrt}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx, = [seval(x, dvar, env, cache) for x = f.args]
    dx /(2 * sqrt(x))
end

function _eval(::Val{:sum}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    dx, = [seval(x, dvar, env, cache) for x = f.args]
    sum(dx)
end

function _eval(::Val{:dot}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx,dy = [seval(x, dvar, env, cache) for x = f.args]
    dot(x,dy) + dot(dx,y)
end


text:
# This is the unit test suite for the Regions.jl package.

using Regions
using Test

@testset "Regions" begin

    include("test_range.jl")
    include("test_run.jl")
    include("test_region.jl")

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ContainerStateWaiting is a waiting state of a container.

    IoK8sApiCoreV1ContainerStateWaiting(;
        message=nothing,
        reason=nothing,
    )

    - message::String : Message regarding why the container is not yet running.
    - reason::String : (brief) reason the container is not yet running.
"""
mutable struct IoK8sApiCoreV1ContainerStateWaiting <: SwaggerModel
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApiCoreV1ContainerStateWaiting(;message=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerStateWaiting, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1ContainerStateWaiting, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApiCoreV1ContainerStateWaiting

const _property_map_IoK8sApiCoreV1ContainerStateWaiting = Dict{Symbol,Symbol}(Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApiCoreV1ContainerStateWaiting = Dict{Symbol,String}(Symbol("message")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerStateWaiting }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerStateWaiting))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerStateWaiting }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ContainerStateWaiting[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerStateWaiting }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerStateWaiting[property_name]

function check_required(o::IoK8sApiCoreV1ContainerStateWaiting)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerStateWaiting }, name::Symbol, val)
end


text:
function build_integrator(obj::ODEIntegratorBuilder, arg0::Orbit, arg1::OrbitType)
    return jcall(obj, "buildIntegrator", AbstractIntegrator, (Orbit, OrbitType), arg0, arg1)
end


text:
# rewrites the expressions generated by clang

# dictionary to map pointer types to desired type

ptr_dict = Dict{Any, Any} (
:(Ptr{PetscInt}) => :(Union(Ptr{PetscInt}, AbstractArray{PetscInt}, Ptr{Void}))
)

ccall_dict = Dict{Any, Any} (
:Mat => :(Ptr{Void}),
)

function petsc_rewriter(obuf)

  for i=1:length(obuf)
    println("i = ", i)
    ex_i = obuf[i]  # obuf is an array of expressions

    println("ex_i = ", ex_i)
    println("typeof(ex_i) = ", typeof(ex_i))

    if typeof(ex_i) == Expr
      head_i = ex_i.head  # function
      # each ex_i has 2 args, a function signature and a function body
      println("head_i = ", head_i)

      # figure out what kind of expression it is, do the right modification
      if head_i == :function
        rewrite_sig(ex_i.args[1])
        rewrite_body(ex_i.args[2])
      elseif ex_i.head == :const
        make_global_const(ex_i)
      elseif ex_i.head == :typealias
        make_typealias_void(ex_i)
      else
        println("not processing expression", ex_i)
      end

    else
      println("not processing ", typeof(ex_i))
    end  # end if Expr

  end

  return obuf  # return modified obuf
end

##### functions to rewrite function signature #####
function rewrite_sig(ex)  # rewrite the function signature

  @assert ex.head == :call  # verify this is a function signature expression

  # ex.args[1] = function name, as a symbol

  for i=2:length(ex.args)  # process all arguments of the function,
                           # each of which is an expression containing arg name, argtype
    println("typeof(ex.args[$i]) = ", typeof(ex.args[i]))
    @assert typeof(ex.args[i]) == Expr  || typeof(ex.args[i]) == Symbol  # verify these are all expressions
    process_sig_arg(ex.args[i])  # process each expression
  end

end

function process_sig_arg(ex)  # take the expression that is an argument to the function and rewrite it

   @assert ex.head == :(::)

   # modify args here
#   arg_name = ex.args[1]  # symbol
#   arg_type = ex.args[2]  # Expr contianing type tag
   println("ex.args[2] = ", ex.args[2])
   println("type = ", typeof(ex.args[2]))
   ex.args[2] = modify_typetag(ex.args[2])

end

function modify_typetag(ex)

#  @assert ex.head == :curly  # verify this is a typetag
  return get(ptr_dict, ex, ex)  # get new typetag if one was given, either an Expr or a Symbol
end

#####  function to  rewrite the body of the function #####
function rewrite_body(ex)  # rewrite body of a function

  @assert ex.head == :block
  # ex has only one argument, the ccall
  # could insert other statements arond the ccall here?

  process_ccall(ex.args[1])
end

function process_ccall(ex)

  @assert ex.head == :ccall  # verify this is the ccall statement
  # args[1] = Expr, tuple of fname, libname
  # args[2] = return type, symbol
  # args[3] = Expr, tuple of types to ccall
  # args[4] = symbol, first argument name,
  # ...
  modify_types(ex.args[3])

end

function modify_types(ex)

  @assert ex.head == :tuple

  for i=1:length(ex.args)
    ex.args[i] = get(ccall_dict, ex.args[i], ex.args[i])  # get the new argument type symbol from dictionary
  end                                                      # use existing type of no key found

  return nothing
end

##### functions to make consts into global consts #####

function make_global_const(ex)

  @assert ex.head == :const

  val = ex.args[1]  # get the assignment
  ex.args[1] = :(global $val)  # make it global

end

function make_typealias_void(ex)

  @assert ex.head == :typealias

  if typeof(ex.args[2]) == Expr  # this a compound type declaration
    ex2 = ex.args[2]
    @assert ex2.head == :curly
    str = string(ex2.args[2])  # get the pointee type name
    println("str = ", str)
    println("str[1] = ", str[1])
    if str[1] == '_'  # if it begin with an underscore
      println("rewriting typealias")
      ex2.args[2] = :Void  # make it a Ptr{Void}
      println("ex2 = ", ex2)
      println("ex = ", ex)
    end
    # else do nothing
  end
end


text:
"""
for case when we want to subtract two masks
"""
module MaskDiffrence
using ModernGL, ..DisplayWords, ..StructsManag, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW,Dictionaries,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic
export undoDiffrence,displayMaskDiffrence
"""
for case when we want to subtract two masks
"""
function processKeysInfo(maskNumbs::Identity{Tuple{Identity{TextureSpec{T}}, Identity{TextureSpec{G}}}},actor::SyncActor{Any, ActorWithOpenGlObjects},keyInfo::KeyboardStruct ) where {T,G}
    textSpecs =  maskNumbs.value
    maskA = textSpecs[1].value
    maskB = textSpecs[2].value

    if(keyInfo.isCtrlPressed)    # when we want to stop displaying diffrence
        undoDiffrence(actor,maskA,maskB )
        addToforUndoVector(actor, ()->displayMaskDiffrence(maskA,maskB,actor ) )

    elseif(keyInfo.isShiftPressed)  # when we want to display diffrence
        displayMaskDiffrence(maskA,maskB,actor )
        addToforUndoVector(actor, ()-> undoDiffrence(actor,maskA,maskB ))

    end #if

end#processKeysInfo
"""
for case  we want to undo subtracting two masks
"""
function undoDiffrence(actor::SyncActor{Any, ActorWithOpenGlObjects},maskA,maskB )
    @uniforms! begin
    actor.actor.mainForDisplayObjects.mainImageUniforms.isMaskDiffrenceVis:=0
           end
setTextureVisibility(true,maskA.uniforms )
setTextureVisibility(true,maskB.uniforms )
basicRender(actor.actor.mainForDisplayObjects.window)

end#undoDiffrence

"""
SUBTRACTING MASKS
used in order to enable subtracting one mask from the other - hence displaying
pixels where value of mask a is present but mask b not (order is important)
automatically both masks will be set to be invisible and only the diffrence displayed

In order to achieve this  we need to have all of the samplers references stored in a list
1) we need to set both masks to invisible - it will be done from outside the shader
2) we set also from outside uniform marking visibility of diffrence to true
3) also from outside we need to set which texture to subtract from which we will achieve this by setting maskAtoSubtr and maskBtoSubtr int ..Uniforms
    those integers will mark which samplers function will use
4) in shader function will be treated as any other mask and will give contribution to output color multiplied by its visibility(0 or 1)
5) inside the function color will be defined as multiplication of two colors of mask A and mask B - colors will be acessed similarly to samplers
6) color will be returned only if value associated with  maskA is greater than mask B and proportional to this difffrence

In order to provide maximum performance and avoid branching inside shader multiple shader programs will be attached and one choosed  that will use diffrence needed
maskToSubtrastFrom,maskWeAreSubtracting - specifications o textures we are operating on
"""
function displayMaskDiffrence(maskA::TextureSpec, maskB::TextureSpec,actor::SyncActor{Any, ActorWithOpenGlObjects})
    if(!maskA.isMainImage && !maskB.isMainImage)
        #defining variables
        dispObj = actor.actor.mainForDisplayObjects
        vertex_shader = dispObj.vertex_shader
        listOfTextSpecsc=  actor.actor.mainForDisplayObjects.listOfTextSpecifications
        fragment_shade,shader_prog= createAndInitShaderProgram(dispObj.vertex_shader,  listOfTextSpecsc,maskA,maskB,dispObj.gslsStr)
        # saving new variables to the actor
        newForDisp = setproperties(dispObj,(shader_program=shader_prog,fragment_shader=fragment_shade ) )
        actor.actor.mainForDisplayObjects=newForDisp
        #making all ready to display
        reactivateMainObj(shader_prog, newForDisp.vbo,actor.actor.calcDimsStruct  )
        activateTextures(listOfTextSpecsc )
        #making diffrence visible
        @uniforms! begin
        dispObj.mainImageUniforms.isMaskDiffrenceVis:=1
                end
        setTextureVisibility(false,maskA.uniforms )
        setTextureVisibility(false,maskB.uniforms )
        basicRender(actor.actor.mainForDisplayObjects.window)
    end#if
end#displayMaskDiffrence

end#MaskDiffrence


text:
using Test, Random, QPDAS
using LinearAlgebra, SparseArrays

import OSQP

Random.seed!(12345)

model = OSQP.Model()

# Test qp problem
me, mi, n = 20, 20, 1000
# Equality
A = randn(me, n)
b = randn(me)
# Inequality
C = randn(mi, n)
d = randn(mi)
# Linear term
z = randn(n)

M = [A;C]
# Ax=b
# Cx≤d
u = [b;d]
l = [b;fill(-Inf, length(d))]

OSQP.setup!(model; P=SparseMatrixCSC{Float64}(I, n, n), l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

QP = QuadraticProgram(A,b,C,d,z, scaling=false)
@time x2, val2 = solve!(QP)
#x2 = solveQP(A,b,C,d,z)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

# Test again with P=2.0I

OSQP.setup!(model; P=SparseMatrixCSC{Float64}(2.0I, n, n), l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

QP = QuadraticProgram(A,b,C,d,z,2.0I, scaling=false)
@time x2, val2 = solve!(QP)
#x2 = solveQP(A,b,C,d,z)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

# Test with semidefinite=false

OSQP.setup!(model; P=SparseMatrixCSC{Float64}(I, n, n), l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

QP = QuadraticProgram(A,b,C,d,z, semidefinite=false, scaling=false)
@time x2, val2 = solve!(QP)
#x2 = solveQP(A,b,C,d,z)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

# # New test
# b = randn(me)
# d = randn(mi)
z = randn(n)

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

update!(QP, z=z)
@time x2, val2 = solve!(QP)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

# Update multiple
b = randn(me)
d = randn(mi)
z = randn(n)

# Rebuild for OSQP
u = [b;d]
l = [b;fill(-Inf, length(d))]

OSQP.setup!(model; P=SparseMatrixCSC{Float64}(I, n, n), l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

update!(QP, z=z, b=b, d=d)
@time x2, val2 = solve!(QP)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

#
# ## Compare gurobi
#
# env = Gurobi.Env()
# model = gurobi_model(env, f = -z, H = sparse(1.0*I, length(z), length(z)),
#     Aeq = A, beq = b, A = -C, b = -d)
#
# optimize(model)

### TEST P SPARSE

me, mi, n = 20, 20, 1000
# Equality
A = randn(me, n)
b = randn(me)
# Inequality
C = randn(mi, n)
d = randn(mi)
# Linear term
z = randn(n)

M = [A;C]
# Ax=b
# Cx≤d
u = [b;d]
l = [b;fill(-Inf, length(d))]

dig = fill(2.0,n)
dig[1:2:end] .= 4.0
dig2 = 0.3*dig
dig2[1:3:end] .= 1
P = spdiagm(-1 => dig2[1:end-1], 0 => dig, 1 => dig2[1:end-1])

OSQP.setup!(model; P=P, l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

QP = QuadraticProgram(A,b,C,d,z,P, scaling=false)
@time x2, val2 = solve!(QP)
#x2 = solveQP(A,b,C,d,z)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

## TEST UPDATE WITH P
b = randn(me)
d = randn(mi)
z = randn(n)

u = [b;d]
l = [b;fill(-Inf, length(d))]
OSQP.setup!(model; P=P, l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())
OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

update!(QP, z=z, b=b, d=d)
@time x2, val2 = solve!(QP)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

### TEST P DENSE

PA = randn(n,n)
P = PA*PA' + I

# Needs to be sparse for OSQP
OSQP.setup!(model; P=SparseMatrixCSC(P), l=l, A=sparse(M), u=u, verbose=false,
    eps_abs=eps(), eps_rel=eps(),
    eps_prim_inf=eps(), eps_dual_inf=eps())

OSQP.update!(model; q=z)
@time results = OSQP.solve!(model)
x1 = results.x
val1 = results.info.obj_val

QP = QuadraticProgram(A,b,C,d,z,P, scaling=false)
@time x2, val2 = solve!(QP)
#x2 = solveQP(A,b,C,d,z)

@test A*x2 ≈ b atol=1e-12 # works up to 1e-14
@test maximum(C*x2 - d) < 1e-12 # works up to 1e-14

@test norm(x1-z) ≈ norm(x2-z) rtol=1e-11 # works up to 1e-12
@test x1 ≈ x2 rtol=1e-10  # works up to 1e-11

@test val1 ≈ val2 rtol=1e-11

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
using Test, SearchLight, Foos

@testset "Foo unit tests" begin

  ### Your tests here
  @test 1 == 1

end;


text:
"""
    mutable struct MatpowerGenerator <: AbstractNodeLabel

Mutable structure descendant of AbstractNodeLabel

# Fields
- `power_min::Complex128`
- `power_max::Complex128`
- `cost_degree::Int64`
- `cost_coeffs::Array{Float64}`
- `gen_status::Bool`: generator on/off

"""
mutable struct MatpowerGenerator <: AbstractNodeLabel
  power_min::Complex128 #Smin = Pmin + i Qmin
  power_max::Complex128 #Smax = Pmax + i Qmax
  cost_degree::Int64 # degree of the cost polynomial
  cost_coeffs::Array{Float64} # cost polynomial coeffs from higher to lower degree
  gen_status::Bool
end

"""
    create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String) where T <: MatpowerGenerator

Create power variables n for generator `elemid` at `bus` in `bus_vars` if `elem_formulation == :NewVar`\n
Return nothing

"""
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String) where T <: MatpowerGenerator
  if elem_formulation == :NewVar
    bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)
  end
  return
end

"""
[cstrname, polynom, lb, ub] = Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MatpowerGenerator

Return the polynomial contribution in power of generator `elemid` at `bus`(name, value, lower bound, upper bound). Will be used to construct power balance constraints in polynomial problem.\n
If `elem_formulation == :NbMinVar`, return generator bounds `["UNIT", Polynomial() ,Smin, Smax]`\n
If `elem_formulation == :NewVar`, return -Sgen `["UNIT", -Sgen, 0, 0]`\n
Return no contribution `["", Polynomial(), 0, 0]` otherwise
"""
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MatpowerGenerator
  cstrname = "UNIT"
  if elem_formulation == :NbMinVar
    lb = element.power_min
    ub = element.power_max
    if !element.gen_status ## generator is off
      lb = ub = 0
    end
    return [cstrname, Polynomial(), lb, ub]
  elseif elem_formulation == :NewVar
    return [cstrname, Polynomial(-bus_vars[elemid]), 0, 0]
  end
  return ["", Polynomial(), 0, 0] # Raise error here ? Keep doing nothing ? One would expect a generator exists if this is called.
end

"""
    constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: MatpowerGenerator

Return all the constraints defined by generator `elemid` at `bus`. Will be used to construct constraints in polynomial problem.\n
If `elem_formulation == :NewVar`, return generator bounds : "Genbounds" => Smin <= Sgen <= Smax\n
Return empty dictionary otherwise
"""
function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: MatpowerGenerator
  if elem_formulation == :NewVar
    lb = element.power_min
    ub = element.power_max
    if !element.gen_status
      ub = lb = 0
    end
    cstrname = get_GenBounds_cstrname()
    p = Polynomial(-bus_vars[elemid])
    return Dict{String, Constraint}(cstrname => lb << p << ub)
  end
  return Dict{String, Constraint}()
end

"""
    cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_vars::Dict{String, Variable}, Snode::Polynomial, lb, ub) where T <: MatpowerGenerator

Return the polynomial contribution in objective generator `elemid` at `bus`.
"""
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_var::Dict{String, Variable}, Snode::Polynomial, lb, ub) where T <: MatpowerGenerator
  Sgen = Polynomial()
  gencost = Polynomial()
  if elem_formulation == :NbMinVar
    _, S, lb_gen, ub_gen = Snodal(element, bus, elemid, elem_formulation, bus_elems_var)
    abs((lb-lb_gen) - (ub-ub_gen)) ≤ 5e-5 || warn("Power missing at $bus... $(abs((lb-lb_gen) - (ub-ub_gen)))")
    Sgen = add(Snode, -S - (ub-ub_gen))
  elseif elem_formulation == :NewVar
    Sgen = Polynomial(bus_elems_var[elemid])
  else
    warn("No cost contribution for gen $elemid, $bus")
  end

  add!(Sgen, conj(Sgen))
  Sactive = (Sgen)*0.5
  Sexp = 1 + Polynomial()
  d = element.cost_degree
  degrees = element.cost_coeffs
  for i=1:2                                               ## NOTE: NOT normal behaviour, required for linear cost...
    add!(gencost, degrees[d-i+1] * Sexp)
    Sexp = Sexp*Sactive
  end
  if !element.gen_status
    gencost = 0
  end
  return gencost
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
"""
Flac stream decoder object.  Contains an opaque pointer only.
"""
type StreamDecoderPtr  # type not immutable so that finalizer can be applied
    v::Ptr{Void}
end

Base.unsafe_convert(::Type{Ptr{Void}},en::StreamDecoderPtr) = en.v

function StreamDecoderPtr()
    en = StreamDecoderPtr(ccall((:FLAC__stream_decoder_new,libflac),Ptr{Void},()))
    finalizer(en,x->ccall((:FLAC__stream_decoder_delete,libflac),Void,(Ptr{Void},),x.v))
    en
end

for (nm,typ) in (("ogg_serial_number",:Clong),
                 ("md5_checking",:Bool),
                 ("metadata_respond",:MetadataType),
                 ("metadata_ignore",:MetadataType))
    @eval begin
        function $(symbol(string("set_",nm)))(dd::StreamDecoderPtr,val)
            ccall(($(string("FLAC__stream_decoder_set_",nm)),libflac),Bool,(Ptr{Void},$typ),dd,val)
        end
    end
end

set_md5_checking(dd::StreamDecoderPtr) = set_md5_checking(dd,true)

"""
Set the stream decoder to respond to all metadata blocks.

Must be called before the decoder is initialized.
"""
set_metadata_respond_all(dd::StreamDecoderPtr) =
    ccall((:FLAC__stream_decoder_set_metadata_respond_all,libflac),Bool,(Ptr{Void},),dd)

"""
Set the stream decoder to ignore all metadata blocks.

Must be called before the decoder is initialized.
"""
set_metadata_ignore_all(dd::StreamDecoderPtr) =
    ccall((:FLAC__stream_decoder_set_metadata_respond_all,libflac),Bool,(Ptr{Void},),dd)

"""
Returns a character string describing the current state of the decoder
"""
get_state_string(dd::StreamDecoderPtr) =
    bytestring(ccall((:FLAC__stream_decoder_get_resolved_state_string,libflac),Ptr{UInt8},
                     (Ptr{Void},),dd))

@enum(StreamDecoderState,
      DecoderMetaDataSearch,
      DecoderMetaDataRead,
      DecoderFrameSyncSearch,
      DecoderFrameRead,
      DecoderEndOfStream,
      DecoderOggError,
      DecoderSeekError,
      DecoderAborted,
      DecoderMemoryAllocationError,
      DecoderUninitialized)

for (nm,typ) in (("state",:DecoderState),
                 ("md5_checking",:Bool),
                 ("total_samples",:Int64),
                 ("channels",:Cint),
                 ("channel_assignment",:ChannelAssignment),
                 ("bits_per_sample",Cint),
                 ("sample_rate",Cint),
                 ("blocksize",Cint))
    @eval begin
        function $(symbol(string("get_",nm)))(dd::StreamDecoderPtr)
            ccall(($(string("FLAC__stream_decoder_get_",nm)),libflac),$typ,(Ptr{Void},),dd)
        end
    end
end

@enum(DecoderInitStatus,
      DecoderInitOK,
      DecoderInitUnsupportedContainer,
      DecoderInitInvalidCallbacks,
      DecoderInitMemoryAllocationError2,
      DecoderInitErrorOpeningFile,
      DecoderInitAlreadyInitialized)

"debugging metadata callback function"
function debug_mcallback(d::Ptr{Void}, mp::Ptr{Void}, client::Ptr{Void})
    typ = unsafe_load(reinterpret(Ptr{MetaDataType},mp))
    println("Metadata callback on typ = ", typ)
    if typ == Info
        show(unsafe_load(reinterpret(Ptr{InfoMetaData}, mp)))
    elseif typ == Padding
        show(unsafe_load(reinterpret(Ptr{PaddingMetaData}, mp)))
    elseif typ == VorbisComment
        show(unsafe_load(reinterpret(Ptr{VorbisComment}, mp)))
    elseif typ == SeekTable
        show(unsafe_load(reinterpret(Ptr{SeekTable}, mp)))
    end
    nothing
end
"silent metadata callback function"
function silent_mcallback(d::Ptr{Void}, mp::Ptr{Void}, client::Ptr{Void})
    nothing
end

"error callback"
function debug_ecallback(d::Ptr{Void}, status::Int32, client::Ptr{Void})
    error("Got error callback with status = $status")
end

"debugging write callback"
function debug_wcallback(dd::Ptr{Void}, hdr::Ptr{FrameHeader},
                   buffer::Ptr{Ptr{Int32}}, client::Ptr{Void})
    fr = unsafe_load(hdr)
    println("Frame")
    println(" blocksize: ", fr.blocksize)
    println(" samplerate: ", fr.sample_rate)
    println(" channels: ", fr.channels)
    println(" channel assignment: ", fr.channel_assignment)
    println(" bits per sample: ", fr.bits_per_sample)
    println(" number type: ", fr.typ)
    println(" frame or sample number: ", fr.num)
    println(" crc: ", fr.crc)

    data = pointer_to_array(unsafe_load(buffer), fr.blocksize)/(2^(fr.bits_per_sample - 1))
    zero(Int32)
end

function initfile!(dd::StreamDecoderPtr, fnm::ByteString; wcallback=debug_wcallback_c,
                   mcallback=debug_mcallback_c, ecallback=debug_ecallback_c, client_data=nothing)
    status = ccall((:FLAC__stream_decoder_init_file,libflac),DecoderInitStatus,
                   (Ptr{Void}, Ptr{UInt8}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
                   dd, fnm, wcallback, mcallback, ecallback, client_data)
    if status != DecoderInitOK
        error("decoder_init_file returned status $status")
    end
    dd
end

process_single(dd::StreamDecoderPtr) = ccall((:FLAC__stream_decoder_process_single,libflac),Bool,(Ptr{Void},),dd)

process_metadata(dd::StreamDecoderPtr) = ccall((:FLAC__stream_decoder_process_until_end_of_metadata,libflac),Bool,(Ptr{Void},),dd)

process_stream(dd::StreamDecoderPtr) = ccall((:FLAC__stream_decoder_process_until_end_of_stream,libflac),Bool,(Ptr{Void},),dd)

function saving_mcallback(d::Ptr{Void}, mp::Ptr{Void}, client::Ptr{Void})
    typ = unsafe_load(reinterpret(Ptr{MetaDataType},mp))
    if typ == Info
        client_data = unsafe_pointer_to_objref(client)
        client_data["metadata"] = InfoMetaData(unsafe_load(reinterpret(Ptr{InfoMetaData}, mp)))
    end
    nothing
end

function buffering_wcallback(dd::Ptr{Void}, hdr::Ptr{FrameHeader},
                   buffer::Ptr{Ptr{Int32}}, client::Ptr{Void})
    fr = unsafe_load(hdr)
    client_data = unsafe_pointer_to_objref(client)

    max_val = Float32(2^(fr.bits_per_sample - 1))
    for chidx = 1:fr.channels
        start_idx = client_data["channel_idxs"][chidx]
        client_data["channels"][start_idx:start_idx + fr.blocksize - 1,chidx] = pointer_to_array(unsafe_load(buffer,chidx), fr.blocksize)/max_val
        client_data["channel_idxs"][chidx] += fr.blocksize
    end
    zero(Int32)
end

function load(f::File{format"FLAC"})
    decoder = StreamDecoderPtr()

    client_data = Dict()
    initfile!(decoder, f.filename, wcallback=buffering_wcallback_c, mcallback=saving_mcallback_c, client_data=pointer_from_objref(client_data))

    # Peek at just the metadata
    process_metadata(decoder)

    # Create a list of Arrays to hold our sampledata for each channel assigned to us
    client_data["channels"] = Array{Float32,2}(client_data["metadata"].totalsamples, client_data["metadata"].channels)
    client_data["channel_idxs"] = ones(UInt64, client_data["metadata"].channels)

    # Process the rest of the stream
    process_stream(decoder)

    # Return data, fs
    return client_data["channels"], client_data["metadata"].samplerate
end

# Calculate cfunction versions of all our callbacks once, at runtime, as is necessary with cfunction's
debug_mcallback_c = Ptr{Void}(C_NULL)
silent_mcallback_c = Ptr{Void}(C_NULL)
saving_mcallback_c = Ptr{Void}(C_NULL)
debug_wcallback_c = Ptr{Void}(C_NULL)
buffering_wcallback_c = Ptr{Void}(C_NULL)
debug_ecallback_c = Ptr{Void}(C_NULL)
function init_decoder_cfunctions()
    global debug_mcallback_c     = cfunction(debug_mcallback, Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    global silent_mcallback_c    = cfunction(silent_mcallback, Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    global saving_mcallback_c    = cfunction(saving_mcallback, Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    global debug_wcallback_c     = cfunction(debug_wcallback, Int32, (Ptr{Void}, Ptr{FrameHeader}, Ptr{Ptr{Int32}}, Ptr{Void}))
    global buffering_wcallback_c = cfunction(buffering_wcallback, Int32, (Ptr{Void}, Ptr{FrameHeader}, Ptr{Ptr{Int32}}, Ptr{Void}))
    global debug_ecallback_c     = cfunction(debug_ecallback, Void, (Ptr{Void}, Int32, Ptr{Void}))
end


text:
# # In-place transforms

# Complex-to-complex and real-to-real transforms can be performed in-place,
# enabling important memory savings.
# The procedure is very similar to that of out-of-place transforms described in
# [the tutorial](@ref Tutorial).
# The differences are illustrated in the sections below.

# ## Creating a domain partition
#
# We start by partitioning a domain of dimensions ``16×32×64`` along all
# available MPI processes.

using PencilFFTs
using MPI
MPI.Init()

dims_global = (16, 32, 64)  # global dimensions

# Such a partitioning is described by a
# [`Pencil`](https://jipolanco.github.io/PencilArrays.jl/dev/Pencils/) object.
# Here we choose to decompose the domain along the last two dimensions.
# In this case, the actual number of processes along each of these dimensions is
# chosen automatically.

decomp_dims = (2, 3)
comm = MPI.COMM_WORLD
pen = Pencil(dims_global, decomp_dims, comm)

# !!! warning "Allowed decompositions"
#
#     Distributed transforms using PencilFFTs.jl require that the first
#     dimension is *not* decomposed.
#     In other words, if one wants to perform transforms, then `decomp_dims`
#     above must *not* contain `1`.

# ## Creating in-place plans

# To create an in-place plan, pass an in-place transform such as
# [`Transforms.FFT!`](@ref) or [`Transforms.R2R!`](@ref) to
# [`PencilFFTPlan`](@ref).
# For instance:

## Perform a 3D in-place complex-to-complex FFT.
transform = Transforms.FFT!()

## Note that one can also combine different types of in-place transforms.
## For instance:
## transform = (
##     Transforms.R2R!(FFTW.REDFT01),
##     Transforms.FFT!(),
##     Transforms.R2R!(FFTW.DHT),
## )

# We can now create a distributed plan from the previously-created domain
# partition and the chosen transform.

plan = PencilFFTPlan(pen, transform)

# Note that in-place real-to-complex transforms are not currently supported.
# (In other words, the `RFFT!` transform type is not defined.)

# ## Allocating data

# As with out-of-place plans, data should be allocated using
# [`allocate_input`](@ref).
# The difference is that, for in-place plans, this function returns
# a [`ManyPencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.ManyPencilArray) object, which is a container holding multiple
# [`PencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.PencilArray) views sharing the same memory space.

## Allocate data for the plan.
## Since `plan` is in-place, this returns a `ManyPencilArray` container.
A = allocate_input(plan)
summary(A)

# Note that [`allocate_output`](@ref) also works for in-place plans.
# In this case, it returns exactly the same thing as `allocate_input`.

# As shown in the next section, in-place plans must be applied on the returned
# `ManyPencilArray`.
# On the other hand, one usually wants to access and modify data, and for this
# one needs the `PencilArray` views contained in the `ManyPencilArray`.
# The input and output array views can be obtained by calling
# [`first(::ManyPencilArray)`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#Base.first-Tuple{ManyPencilArray}) and [`last(::ManyPencilArray)`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#Base.last-Tuple{ManyPencilArray}).

# For instance, we can initialise the input array with some data before
# transforming:

using Random
u_in = first(A)  # input data view
randn!(u_in)
summary(u_in)

# ## Applying plans

# Like in `FFTW.jl`, one can perform in-place transforms using the `*` and
# `\ ` operators.
# As mentioned above, in-place plans must be applied on the `ManyPencilArray`
# containers returned by `allocate_input`.

plan * A;  # performs in-place forward transform

# After performing an in-place transform, data contained in `u_in` has been
# overwritten and has no "physical" meaning.
# In other words, `u_in` should **not** be used at this point.
# To access the transformed data, one should retrieve the output data view using
# `last(A)`.
#
# For instance, to compute the global sum of the transformed data:

u_out = last(A)  # output data view
sum(u_out)       # sum of transformed data (note that `sum` reduces over all processes)

# Finally, we can perform a backward transform and do stuff with the input view:

plan \ A;  # perform in-place backward transform

# At this point, the data can be once again found in the input view `u_in`,
# while `u_out` should not be accessed.


text:
using ValueGraphFromDataFrame
using SimpleValueGraphs
using LightGraphs
using DataFrames
using Test

vertex_attributes = DataFrame(
    name=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"],
    col1=1:10,
    col2=vcat(trues(5), falses(5))
)

edgelist = DataFrame(
    from=["a", "a", "a", "b", "f", "g", "d", "d", "c"],
    to  =["b", "c", "i", "a", "d", "f", "e", "a", "j"],
    #eval=1:9
)

g = ValDiGraph(edgelist, vertex_attributes)
e = LightGraphs.edges(g) |> first

@testset "ValueGraphFromDataFrame.jl" begin
    get_vertexval(g, 4, :name) == "d"
    get_vertexval(g, 9, :col1) == 9
    get_vertexval(g, 3, :col2) == true
    get_vertexval(g, src(e), :name) == "a"
    get_vertexval(g, dst(e), :col1) == 2
    #get_edgeval(g, 1, 3, :eval) == 3
end


text:
using ChaosTools
using Test, StaticArrays
using DynamicalSystemsBase: CDS, DDS
using DynamicalSystemsBase.Systems: hoop, hoop_jac, hiip, hiip_jac
using DynamicalSystemsBase.Systems: loop, loop_jac, liip, liip_jac
using OrdinaryDiffEq

println("\nTesting lyapunov exponents...")
let
u0 = [0, 10.0, 0]
p = [10, 28, 8/3]
u0h = ones(2)
ph = [1.4, 0.3]

FUNCTIONS = [liip, liip_jac, loop, loop_jac, hiip, hiip_jac, hoop, hoop_jac]
INITCOD = [u0, u0h]
PARAMS = [p, ph]
MLE = [[0.75, 0.95], [0.41, 0.43]]
SLE = [[-0.1, 0.1], [-1.63, -1.61]]

@testset "Lyapunovs" begin
for i in 1:8
    @testset "combination $i" begin
        sysindx = i < 5 ? 1 : 2
        if i < 5
            if isodd(i)
                ds = CDS(FUNCTIONS[i], INITCOD[sysindx], PARAMS[sysindx])
            else
                ds = CDS(FUNCTIONS[i-1], INITCOD[sysindx], PARAMS[sysindx], FUNCTIONS[i])
            end
        else
            if isodd(i)
                ds = DDS(FUNCTIONS[i], INITCOD[sysindx], PARAMS[sysindx])
            else
                ds = DDS(FUNCTIONS[i-1], INITCOD[sysindx], PARAMS[sysindx], FUNCTIONS[i])
            end
        end

        λ = lyapunovs(ds, 2000)
        if i < 5
            λ2 = lyapunovs(ds, 2000; dt = 2.0, Ttr = 10.0,
            diff_eq_kwargs = Dict(:abstol=>1e-9, :solver => Tsit5()))
        else
            λ2 = lyapunovs(ds, 1000; dt = 5, Ttr = 20)
        end

        @test MLE[sysindx][1] < λ[1] < MLE[sysindx][2]
        @test SLE[sysindx][1] < λ[2] < SLE[sysindx][2]
        @test MLE[sysindx][1] < λ2[1] < MLE[sysindx][2]
        @test SLE[sysindx][1] < λ2[2] < SLE[sysindx][2]

        if isodd(i) # Jacobians dont matter for MLE
            T = i < 5 ? 10000 : 1000
            λ = lyapunov(ds, T)
            if i < 5
                λ2 = lyapunov(ds, T; dt = 2.0, Ttr = 10.0,
                diff_eq_kwargs = Dict(:abstol=>1e-9, :solver => Tsit5()))
            else
                λ2 = lyapunovs(ds, T; dt = 5, Ttr = 20)
            end
            @test MLE[sysindx][1] < λ < MLE[sysindx][2]
            @test MLE[sysindx][1] < λ2[1] < MLE[sysindx][2]
        end
    end
end
end

@testset "1D Lyapunovs" begin
    ds = Systems.logistic(;r = 4.0)
    λ = lyapunov(ds, 10000; Ttr = 100)
    @test 0.692 < λ < 0.694
end

@testset "Negative λ, continuous" begin
    f(u, p, t) = -0.9u
    g(du, u, p, t) = (du .= -0.9u)

    ds = ContinuousDynamicalSystem(f, rand(SVector{3}), nothing)
    λ1 = lyapunov(ds, 1000)
    @test λ1 < 0
    ds = ContinuousDynamicalSystem(g, rand(3), nothing)
    λ2 = lyapunov(ds, 1000)
    @test λ2 < 0

    @testset "Lorenz stable" begin
        ds = Systems.lorenz(ρ = 20.0)
        @test lyapunov(ds, 2000, Ttr = 100) < 0
    end
end

@testset "Negative λ, discrete" begin
    f(u, p, t) = 0.9u
    g(du, u, p, t) = (du .= 0.9u)

    ds = DiscreteDynamicalSystem(f, rand(SVector{3}), nothing)
    λ1 = lyapunov(ds, 100000)
    @test isapprox(λ1, log(0.9); rtol = 1e-4)
    ds = DiscreteDynamicalSystem(g, rand(3), nothing)
    λ2 = lyapunov(ds, 100000)
    @test isapprox(λ1, log(0.9); rtol = 1e-4)
end

@testset "Lyapunov convergence" begin

ds = Systems.towel()
tinteg = tangent_integrator(ds)
ls, t = ChaosTools.lyapunovs_convergence(tinteg, 20000, 1, 0)
l1 = [x[1] for x in ls]
@test 0.434 < l1[end] < 0.436
end

end


text:
module Bamberg

using StaticArrays, Plots

export fnp, findminnorm, pointsvecnorm, pointsvecnorm_opt
export generateidealdata, plotdist

include("pointdistance.jl")
include(joinpath("../test/idealdata.jl"))
include("plots.jl")

end

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
include("./aa_properties.jl")
include("./which_natural_aminoacid.jl")
include("./selmacros.jl")
include("./select.jl")

export isacidic,
    isaliphatic,
    isaromatic,
    isbackbone,
    isbasic,
    ischarged,
    ishydrophobic,
    isneutral,
    isnonpolar,
    ispolar,
    isprotein,
    issidechain,
    iswater

export select, selindex

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
# Context management

import Base: unsafe_convert, @deprecate_binding

export
    CuContext, destroy, current_context,
    push, pop, synchronize, device

@enum(CUctx_flags, SCHED_AUTO           = 0x00,
                   SCHED_SPIN           = 0x01,
                   SCHED_YIELD          = 0x02,
                   SCHED_BLOCKING_SYNC  = 0x04,
                   MAP_HOST             = 0x08,
                   LMEM_RESIZE_TO_MAX   = 0x10)
@deprecate_binding BLOCKING_SYNC SCHED_BLOCKING_SYNC

typealias CuContext_t Ptr{Void}

"""
Create a CUDA context for device. A context on the GPU is analogous to a process on the
CPU, with its own distinct address space and allocated resources. When a context is
destroyed, the system cleans up the resources allocated to it.
"""
immutable CuContext
    handle::CuContext_t
end

unsafe_convert(::Type{CuContext_t}, ctx::CuContext) = ctx.handle

function CuContext(dev::CuDevice, flags::CUctx_flags)
    handle_ref = Ref{CuContext_t}()
    @apicall(:cuCtxCreate, (Ptr{CuContext_t}, Cuint, Cint),
                           handle_ref, flags, dev)
    CuContext(handle_ref[])
end

CuContext(dev::CuDevice) = CuContext(dev, SCHED_AUTO)

"Destroy the CUDA context, releasing resources allocated to it."
function destroy(ctx::CuContext)
    @apicall(:cuCtxDestroy, (CuContext_t,), ctx)
end

function current_context()
    handle_ref = Ref{CuContext_t}()
    @apicall(:cuCtxGetCurrent, (Ptr{CuContext_t},), handle_ref)
    CuContext(handle_ref[])
end

function push(ctx::CuContext)
    @apicall(:cuCtxPushCurrent, (CuContext_t,), ctx)
end

function pop()
    handle_ref = Ref{CuContext_t}()
    @apicall(:cuCtxPopCurrent, (Ptr{CuContext_t},), handle_ref)
    return nothing
end

function device(ctx::CuContext)
    if current_context() != ctx
        # TODO: we could push and pop here
        throw(ArgumentError("context should be active"))
    end

    # TODO: cuCtxGetDevice returns the device ordinal, but as a CUDevice*?
    #       This can't be right...
    device_ref = Ref{Cint}()
    @apicall(:cuCtxGetDevice, (Ptr{Cint},), device_ref)
    return CuDevice(device_ref[])
end

synchronize(ctx::CuContext) = @apicall(:cuCtxSynchronize, (CuContext_t,), ctx)


text:
macro ArbFloat(x)
    convert(ArbFloat, string(:($x)))
end
macro ArbFloat(p,x)
    convert(ArbFloat{:($p)}, string(:($x)))
end

# interconvert Arb with Arf

function convert{P}(::Type{ArbFloat{P}}, x::ArfFloat{P})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_arf), Void, (Ptr{ArbFloat{P}}, Ptr{ArfFloat{P}}), &z, &x)
    z
end

function convert{P}(::Type{ArfFloat{P}}, x::ArbFloat{P})
    z = initializer(ArfFloat{P})
    z.mid_exp  = x.mid_exp
    z.mid_size = x.mid_size
    z.mid_d1   = x.mid_d1
    z.mid_d2   = x.mid_d2
    z
end

#interconvert ArbFloat{P} with ArbFloat{Q}

function convert{P,Q}(::Type{ArbFloat{Q}}, a::ArbFloat{P})
    if (Q < P)
        a = round(a, Q, 2)
    end

    z = initializer(ArbFloat{Q})
    z.mid_exp  = a.mid_exp
    z.mid_size = a.mid_size
    z.mid_d1   = a.mid_d1
    z.mid_d2   = a.mid_d2
    z.rad_exp  = a.rad_exp
    z.rad_man  = a.rad_man

    z
end

#

function convert{P}(::Type{ArbFloat{P}}, x::UInt)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_ui), Void, (Ptr{ArbFloat{P}}, UInt), &z, x)
    z
end
if sizeof(Int)==sizeof(Int64)
   convert{P}(::Type{ArbFloat{P}}, x::UInt32) = convert(ArbFloat{P}, convert(UInt64,x))
else
   convert{P}(::Type{ArbFloat{P}}, x::UInt64) = convert(ArbFloat{P}, convert(UInt32,x))
end
convert{P}(::Type{ArbFloat{P}}, x::UInt16) = convert(ArbFloat{P}, convert(UInt,x))

function convert{P}(::Type{ArbFloat{P}}, x::Int)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_si), Void, (Ptr{ArbFloat{P}}, Int), &z, x)
    z
end
if sizeof(Int)==sizeof(Int64)
   convert{P}(::Type{ArbFloat{P}}, x::Int32) = convert(ArbFloat{P}, convert(Int64,x))
else
   convert{P}(::Type{ArbFloat{P}}, x::Int64) = convert(ArbFloat{P}, convert(Int32,x))
end
convert{P}(::Type{ArbFloat{P}}, x::Int16) = convert(ArbFloat{P}, convert(Int,x))

function convert{P}(::Type{ArbFloat{P}}, x::Float64)
    z = initializer(ArbFloat{P})
    fp=copy(x)
    ccall(@libarb(arb_set_d), Void, (Ptr{ArbFloat{P}}, Float64), &z, fp)
    z
end
convert{P}(::Type{ArbFloat{P}}, x::Float32) = convert(ArbFloat{P}, convert(Float64,x))
convert{P}(::Type{ArbFloat{P}}, x::Float16) = convert(ArbFloat{P}, convert(Float64,x))

function convert{P}(::Type{ArbFloat{P}}, x::String)
    s = String(x)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_str), Void, (Ptr{ArbFloat}, Ptr{UInt8}, Int), &z, s, P)
    z
end

convert(::Type{BigInt}, x::String) = parse(BigInt,x)
convert(::Type{BigFloat}, x::String) = parse(BigFloat,x)

function convert{P}(::Type{ArbFloat{P}}, x::BigFloat)
     x = round(x,P,2)
     s = string(x)
     ArbFloat{P}(s)
end

function convert{P}(::Type{BigFloat}, x::ArbFloat{P})
     s = string(midpoint(x))
     round(parse(BigFloat, s),P,2)
end

function convert{I<:Integer,P}(::Type{Rational{I}}, x::ArbFloat{P})
    bf = convert(BigFloat, x)
    convert(Rational{I}, bf)
end

convert{P}(::Type{ArbFloat{P}}, x::BigInt)   = convert(ArbFloat{P}, convert(BigFloat,x))
convert{P}(::Type{ArbFloat{P}}, x::Rational) = convert(ArbFloat{P}, convert(BigFloat,x))
convert{P,S}(::Type{ArbFloat{P}}, x::Irrational{S}) =
    convert(ArbFloat{P}, convert(BigFloat,x))

#= returns 256.0 for convert(big(1.5))
function convert{P}(::Type{ArbFloat{P}}, x::BigFloat)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_round_fmpz), Void, (Ptr{ArbFloat}, Ptr{BigFloat}, Int), &z, &x, P)
    z
end
=#

convert{P}(::Type{ArbFloat{P}}, y::ArbFloat{P}) = y

function convert{P}(::Type{Float64}, x::ArbFloat{P})
    s = smartarbstring(x)
    parse(Float64,s)
end
convert{P}(::Type{Float32}, x::ArbFloat{P}) = convert(Float32,convert(Float64,x))

function convert{I<:Integer,P}(::Type{I}, x::ArbFloat{P})
    s = smartarbstring(x)
    parse(I,split(s,".")[1])
end

for T in (:Int128, :Int64, :Int32, :Int16, :Float64, :Float32, :Float16,
          :(Rational{Int64}), :(Rational{Int32}), :(Rational{Int16}),
          :String)
  @eval convert(::Type{ArbFloat}, x::$T) = convert(ArbFloat{precision(ArbFloat)}, x)
end

# Promotion
for T in (:Int128, :Int64, :Int32, :Int16, :Float64, :Float32, :Float16,
          :(Rational{Int64}), :(Rational{Int32}), :(Rational{Int16}),
          :String)
  @eval promote_rule{P}(::Type{ArbFloat{P}}, ::Type{$T}) = ArbFloat{P}
end

promote_rule{P}(::Type{ArbFloat{P}}, ::Type{BigFloat}) = BigFloat
promote_rule{P}(::Type{ArbFloat{P}}, ::Type{BigInt}) = ArbFloat{P}
promote_rule{P}(::Type{ArbFloat{P}}, ::Type{Rational{BigInt}}) = Rational{BigInt}

promote_rule{P,Q}(::Type{ArbFloat{P}}, ::Type{ArbFloat{Q}}) =
    ifelse(P>Q, ArbFloat{P}, ArbFloat{Q})


text:
push!(LOAD_PATH, "../src/")
using Documenter, Limbdark
include("../src/transit_structure.jl")

makedocs(
    sitename="Limbdark",
    pages = [
        "Home" => "index.md",
        "Installation" => "install.md",
        "Quick start" => "quickstart.md",
        "API" => "api.md"
    ],
    format = Documenter.HTML(prettyurls = false)
)