text:
POMDPs.solve(sol::DESPOTSolver, p::POMDP) = DESPOTPlanner(sol, p)

function POMDPModelTools.action_info(p::DESPOTPlanner, b)
    info = Dict{Symbol, Any}()
    try
        Random.seed!(p.rs, rand(p.rng, UInt32))

        D = build_despot(p, b)

        if p.sol.tree_in_info
            info[:tree] = D
        end

        check_consistency(p.rs)

        if isempty(D.children[1]) && D.U[1] <= D.l_0[1]
            throw(NoGap(D.l_0[1]))
        end

        best_l = -Inf
        best_as = actiontype(p.pomdp)[]
        for ba in D.children[1]
            l = ba_l(D, ba)
            if l > best_l
                best_l = l
                best_as = [D.ba_action[ba]]
            elseif l == best_l
                push!(best_as, D.ba_action[ba])
            end
        end

        return rand(p.rng, best_as)::actiontype(p.pomdp), info # best_as will usually only have one entry, but we want to break the tie randomly
    catch ex
        return default_action(p.sol.default_action, p.pomdp, b, ex)::actiontype(p.pomdp), info
    end
end

POMDPs.action(p::DESPOTPlanner, b) = first(action_info(p, b))::actiontype(p.pomdp)

ba_l(D::DESPOT, ba::Int) = D.ba_rho[ba] + sum(D.l[bnode] for bnode in D.ba_children[ba])

POMDPs.updater(p::DESPOTPlanner) = SIRParticleFilter(p.pomdp, p.sol.K, rng=p.rng)

function Random.seed!(p::DESPOTPlanner, seed)
    Random.seed!(p.rng, seed)
    return p
end


text:
"""
    compute(eph,jd0,time,target,center)

Compute position and velocity of target with respect to center at epoch
jd0+time. This method does not support the NAIF numbering scheme.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

This method does not support the NAIF body identification scheme.

Output units are:
* AU and AU/day for position and velocity
* rad and rad/day for librations
* second and unitless for time ephemeris and time ephemeris rate

# Arguments
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required.
- `center::Integer`: The origin of the coordinate system.

The possible values for target and center are :
* 1 : Mercury Barycenter
* 2 : Venus Barycenter
* 3 : Earth
* 4 : Mars Barycenter
* 5 : Jupiter Barycenter
* 6 : Saturn Barycenter
* 7 : Uranus Barycenter
* 8 : Neptune Barycenter
* 9 : Pluto Barycenter
* 10    : Moon
* 11    : Sun
* 12    : Solar Sytem barycenter
* 13    : Earth-moon barycenter
* 14    : Nutation angles
* 15    : Librations
* 16    : TT-TDB
* 17    : TCG-TCB
* asteroid number + 2000000    : asteroid

"""
function compute(eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    result = Array{Float64,1}(undef,6)
    stat = unsafe_compute!(result,eph,jd0,time,target,center)
    @_checkStatus stat "Unable to compute ephemeris"
    return result
end

"""
    unsafe_compute!(eph,jd0,time,target,center)

In place version of the compute function. Does not perform any checks!
Compute position and velocity of target with respect to center at epoch
jd0+time. This method does not support the NAIF numbering scheme.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

This method does not support the NAIF body identification scheme.

Output units are:
* AU and AU/day for position and velocity
* rad and rad/day for librations
* second and unitless for time ephemeris and time ephemeris rate

# Arguments
- `result`: container for result. It is not checked if it is sufficiently large enough!
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required.
- `center::Integer`: The origin of the coordinate system.

# Return:
- status integer from CALCEPH: 0 if an error occured

The possible values for target and center are :
* 1 : Mercury Barycenter
* 2 : Venus Barycenter
* 3 : Earth
* 4 : Mars Barycenter
* 5 : Jupiter Barycenter
* 6 : Saturn Barycenter
* 7 : Uranus Barycenter
* 8 : Neptune Barycenter
* 9 : Pluto Barycenter
* 10    : Moon
* 11    : Sun
* 12    : Solar Sytem barycenter
* 13    : Earth-moon barycenter
* 14    : Nutation angles
* 15    : Librations
* 16    : TT-TDB
* 17    : TCG-TCB
* asteroid number + 2000000    : asteroid

"""
function unsafe_compute!(result,eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer)
    stat = ccall((:calceph_compute, libcalceph), Cint,
    (Ptr{Cvoid},Cdouble,Cdouble,Cint,Cint,Ref{Cdouble}),
    eph.data,jd0,time,target,center,result)
    return stat
end

"""
    compute(eph,jd0,time,target,center,unit)

Compute position and velocity of target with respect to center
at epoch jd0+time.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

# Arguments
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required. The numbering system depends on the parameter unit.
- `center::Integer`: The origin of the coordinate system. The numbering system depends on the parameter unit.
- `unit::Integer` : The units of the result. This integer is a sum of some unit constants (unit*) and/or the constant useNaifId. If the unit contains useNaifId, the NAIF identification numbering system is used for the target and the center. If the unit does not contain useNaifId, the old number system is used for the target and the center.

"""
function compute(eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer,unit::Integer)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    result = Array{Float64,1}(undef,6)
    stat = unsafe_compute!(result,eph,jd0,time,target,center,unit)
    @_checkStatus stat "Unable to compute ephemeris"
    return result
end

"""
    unsafe_compute!(result,eph,jd0,time,target,center,unit)

In place version of the compute function. Does not perform any checks!
Compute position and velocity of target with respect to center
at epoch jd0+time.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

# Arguments
- `result`: container for result. It is not checked if it is sufficiently large enough!
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required. The numbering system depends on the parameter unit.
- `center::Integer`: The origin of the coordinate system. The numbering system depends on the parameter unit.
- `unit::Integer` : The units of the result. This integer is a sum of some unit constants (unit*) and/or the constant useNaifId. If the unit contains useNaifId, the NAIF identification numbering system is used for the target and the center. If the unit does not contain useNaifId, the old number system is used for the target and the center.

# Return:
- status integer from CALCEPH: 0 if an error occured

"""
function unsafe_compute!(result,eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer,unit::Integer)
    stat = ccall((:calceph_compute_unit, libcalceph), Cint,
    (Ptr{Cvoid},Cdouble,Cdouble,Cint,Cint,Cint,Ref{Cdouble}),
    eph.data,jd0,time,target,center,unit,result)
    return stat
end

"""
    compute(eph,jd0,time,target,center,unit,order)

Compute position and derivatives up to order of target with respect to center
at epoch jd0+time.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

# Arguments
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required. The numbering system depends on the parameter unit.
- `center::Integer`: The origin of the coordinate system. The numbering system depends on the parameter unit.
- `unit::Integer` : The units of the result. This integer is a sum of some unit constants (unit*) and/or the constant useNaifId. If the unit contains useNaifId, the NAIF identification numbering system is used for the target and the center. If the unit does not contain useNaifId, the old number system is used for the target and the center.
- `order::Integer` : The order of derivatives
    * 0: only the position is computed.
    * 1: only the position and velocity are computed.
    * 2: only the position, velocity and acceleration are computed.
    * 3: the position, velocity and acceleration and jerk are computed.

"""
function compute(eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer,unit::Integer,order::Integer)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    @_checkOrder order
    result = Array{Float64,1}(undef,3+3order)
    stat = unsafe_compute!(result,eph,jd0,time,target,center,unit,order)
    @_checkStatus stat "Unable to compute ephemeris"
    return result
end

"""
    unsafe_compute!(result,eph,jd0,time,target,center,unit,order)

In place version of the compute function. Does not perform any checks!
Compute position and derivatives up to order of target with respect to center
at epoch jd0+time.
To get the best precision for the interpolation, the time is splitted in two
floating-point numbers. The argument jd0 should be an integer and time should
be a fraction of the day. But you may call this function with time=0 and jd0,
the desired time, if you don't care about precision.

# Arguments
- `result`: container for result. It is not checked if it is sufficiently large enough!
- `eph`: ephemeris
- `jd0::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `time::Float64`: jd0+time must be equal to the Julian Day for the time coordinate corresponding to the ephemeris (usually TDB or TCB)
- `target::Integer`: The body or reference point whose coordinates are required. The numbering system depends on the parameter unit.
- `center::Integer`: The origin of the coordinate system. The numbering system depends on the parameter unit.
- `unit::Integer` : The units of the result. This integer is a sum of some unit constants (unit*) and/or the constant useNaifId. If the unit contains useNaifId, the NAIF identification numbering system is used for the target and the center. If the unit does not contain useNaifId, the old number system is used for the target and the center.
- `order::Integer` : The order of derivatives
    * 0: only the position is computed.
    * 1: only the position and velocity are computed.
    * 2: only the position, velocity and acceleration are computed.
    * 3: the position, velocity and acceleration and jerk are computed.

# Return:
- status integer from CALCEPH: 0 if an error occured

"""
function unsafe_compute!(result,eph::Ephem,jd0::Float64,time::Float64,
   target::Integer,center::Integer,unit::Integer,order::Integer)
    stat = ccall((:calceph_compute_order, libcalceph), Cint,
    (Ptr{Cvoid},Cdouble,Cdouble,Cint,Cint,Cint,Cint,Ref{Cdouble}),
    eph.data,jd0,time,target,center,unit,order,result)
    return stat
end

ext: jl


text:
struct PolicyRewardTest <: TrainingTest
    idx::Integer
    state::AbstractState
    image::AbstractImage
    actions::Vector{<:Integer}
    action::Integer
    reward::Float32
end

mutable struct PolicyRewardTests <: TrainingTests
    tests::Vector{PolicyRewardTest}
    images::AbstractArray{Float32, 4}
    actual_reward::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function PolicyRewardTests(tests; period=1, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            Flux.unsqueeze(get_reward.(tests), 1),
            period,
            warmup,
        )
    end
end

function gettestfn(tests::PolicyRewardTests)
    function fn(n::Network)
        initial_output = initialinferencefn(tests)(n)
        acc, fails = getpolicyacc(tests, initial_output)
        recurrent_output = recurrentinferencefn(tests, initial_output)(n)
        acc, fails, getrewardmse(tests, recurrent_output)
    end
end

struct PolicyRewardTestsSVG
    tests::PolicyRewardTests
end
svg(x::PolicyRewardTests) = PolicyRewardTestsSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::PolicyRewardTestsSVG)
    num_tests = length(x.tests)
    range = 1:(min(10, num_tests))
    subset = 1:num_tests
    subset_indices = subset[range]
    svgbatch = [gridsvg(x.tests[i]) for i in subset_indices]
    Base.show(io, m, GridSVG(svgbatch, main_title="PolicyReward Tests"))
end


text:
@inline UNITLESS_ABS2(x::Number) = abs2(x)
@inline UNITLESS_ABS2(x::AbstractArray) = sum(UNITLESS_ABS2, x)
@inline UNITLESS_ABS2(x::RecursiveArrayTools.ArrayPartition) = sum(UNITLESS_ABS2, x.x)

@inline recursive_length(u::AbstractArray{<:Number}) = length(u)
@inline recursive_length(u::Number) = length(u)
@inline recursive_length(u::AbstractArray{<:AbstractArray}) = sum(recursive_length, u)
@inline recursive_length(u::RecursiveArrayTools.ArrayPartition) = sum(recursive_length, u.x)
@inline recursive_length(u::RecursiveArrayTools.VectorOfArray) = sum(recursive_length, u.u)

@inline ODE_DEFAULT_NORM(u::Union{AbstractFloat,Complex},t) = @fastmath abs(u)
@inline ODE_DEFAULT_NORM(u::Array{T},t) where T<:Union{AbstractFloat,Complex} =
                                         sqrt(real(sum(abs2,u)) / length(u))
@inline ODE_DEFAULT_NORM(u::StaticArrays.StaticArray{T},t) where T<:Union{AbstractFloat,Complex} =
                                            sqrt(real(sum(abs2,u)) / length(u))

@inline ODE_DEFAULT_NORM(u::AbstractArray,t) = sqrt(UNITLESS_ABS2(u)/recursive_length(u))
@inline ODE_DEFAULT_NORM(u,t) = norm(u)

@inline ODE_DEFAULT_ISOUTOFDOMAIN(u,p,t) = false
@inline ODE_DEFAULT_PROG_MESSAGE(dt,u,p,t) =
           "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(maximum(abs.(u)))

@inline NAN_CHECK(x::Number) = isnan(x)
@inline NAN_CHECK(x::Float64) = isnan(x) || (x>1e50)
@inline NAN_CHECK(x::Enum) = false
@inline NAN_CHECK(x::AbstractArray) = any(NAN_CHECK, x)
@inline NAN_CHECK(x::RecursiveArrayTools.ArrayPartition) = any(NAN_CHECK, x.x)

@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u,p,t) = false
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::Union{Number,AbstractArray},p,t) = NAN_CHECK(u)


text:
# Autogenerated wrapper script for eudev_jll for i686-linux-gnu
export libudev, udevadm, udevd

using gperf_jll
JLLWrappers.@generate_wrapper_header("eudev")
JLLWrappers.@declare_library_product(libudev, "libudev.so.1")
JLLWrappers.@declare_executable_product(udevadm)
JLLWrappers.@declare_executable_product(udevd)
function __init__()
    JLLWrappers.@generate_init_header(gperf_jll)
    JLLWrappers.@init_library_product(
        libudev,
        "lib/libudev.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        udevadm,
        "bin/udevadm",
    )

    JLLWrappers.@init_executable_product(
        udevd,
        "sbin/udevd",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JetPack, Jets, Test

@testset "JopPermutedims, 2D" begin
    A = JopPermutedims(JetSpace(Float64,3,2),[2;1])
    m = rand(domain(A))
    d = A*m
    @test d ≈ transpose(m)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopPermutedims, 3D" begin
    A = JopPermutedims(JetSpace(Float64,4,3,2),[3;1;2])
    m = rand(domain(A))
    d = A*m
    @test d ≈ permutedims(m,[3;1;2])
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end


text:
export BrickHATSleepMode
struct BrickHATSleepMode
    power_off_delay::Integer
    power_off_duration::Integer
    raspberry_pi_off::Bool
    bricklets_off::Bool
    enable_sleep_indicator::Bool
end

export BrickHATVoltages
struct BrickHATVoltages
    voltage_usb::Integer
    voltage_dc::Integer
end

export BrickHATVoltagesCallbackConfiguration
struct BrickHATVoltagesCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

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

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

export BrickHAT
"""
HAT for Raspberry Pi with 8 Bricklets ports and real-time clock
"""
mutable struct BrickHAT <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export set_sleep_mode
"""
    $(SIGNATURES)

Sets the sleep mode.

.. note::
 Calling this function will cut the Raspberry Pi's power after Power Off Delay seconds.
 You have to shut down the operating system yourself, e.g. by calling 'sudo shutdown -h now'.

Parameters:

* Power Off Delay: Time before the RPi/Bricklets are powered off.
* Power Off Duration: Duration that the RPi/Bricklets stay powered off.
* Raspberry Pi Off: RPi is powered off if set to true.
* Bricklets Off: Bricklets are powered off if set to true.
* Enable Sleep Indicator: If set to true, the status LED will blink in a 1s interval
  during the whole power off duration. This will draw additional 0.3mA.

Example: To turn RPi and Bricklets off in 5 seconds for 10 minutes with sleep
indicator enabled, call (5, 60*10, true, true, true).

This function can also be used to implement a watchdog. To do this you can
write a program that calls this function once per second in a loop with
(10, 2, true, false, false). If the RPi crashes or gets stuck
the HAT will reset the RPi after 10 seconds.
"""
function set_sleep_mode(device::BrickHAT, power_off_delay, power_off_duration, raspberry_pi_off, bricklets_off, enable_sleep_indicator)
    device.deviceInternal.set_sleep_mode(power_off_delay, power_off_duration, raspberry_pi_off, bricklets_off, enable_sleep_indicator)
end

export get_sleep_mode
"""
    $(SIGNATURES)

Returns the sleep mode settings as set by :func:`Set Sleep Mode`.
"""
function get_sleep_mode(device::BrickHAT)
    return device.deviceInternal.get_sleep_mode()
end

export set_bricklet_power
"""
    $(SIGNATURES)

Set to true/false to turn the power supply of the connected Bricklets on/off.
"""
function set_bricklet_power(device::BrickHAT, bricklet_power)
    device.deviceInternal.set_bricklet_power(bricklet_power)
end

export get_bricklet_power
"""
    $(SIGNATURES)

Returns the power status of the connected Bricklets as set by :func:`Set Bricklet Power`.
"""
function get_bricklet_power(device::BrickHAT)
    return device.deviceInternal.get_bricklet_power()
end

export get_voltages
"""
    $(SIGNATURES)

Returns the USB supply voltage and the DC input supply voltage.

There are three possible combinations:

* Only USB connected: The USB supply voltage will be fed back to the
  DC input connector. You will read the USB voltage and a slightly lower
  voltage on the DC input.
* Only DC input connected: The DC voltage will not be fed back to the
  USB connector. You will read the DC input voltage and the USB voltage
  will be 0.
* USB and DC input connected: You will read both voltages. In this case
  the USB supply will be without load, but it will work as backup if you
  disconnect the DC input (or if the DC input voltage falls below the
  USB voltage).
"""
function get_voltages(device::BrickHAT)
    return device.deviceInternal.get_voltages()
end

export set_voltages_callback_configuration
"""
    $(SIGNATURES)

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

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

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

.. versionadded:: 2.0.1\$nbsp;(Firmware)
"""
function set_voltages_callback_configuration(device::BrickHAT, period, value_has_to_change)
    device.deviceInternal.set_voltages_callback_configuration(period, value_has_to_change)
end

export get_voltages_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set Voltages Callback Configuration`.

.. versionadded:: 2.0.1\$nbsp;(Firmware)
"""
function get_voltages_callback_configuration(device::BrickHAT)
    return device.deviceInternal.get_voltages_callback_configuration()
end

export set_rtc_driver
"""
    $(SIGNATURES)

Configures the RTC driver that is given to the Raspberry Pi to be used.
Currently there are two different RTCs used:

* Hardware version <= 1.5: PCF8523
* Hardware version 1.6: DS1338

The correct driver will be set during factory flashing by Tinkerforge.

.. versionadded:: 2.0.3\$nbsp;(Firmware)
"""
function set_rtc_driver(device::BrickHAT, rtc_driver)
    device.deviceInternal.set_rtc_driver(rtc_driver)
end

export get_rtc_driver
"""
    $(SIGNATURES)

Returns the RTC driver as set by :func:`Set RTC Driver`.

.. versionadded:: 2.0.3\$nbsp;(Firmware)
"""
function get_rtc_driver(device::BrickHAT)
    return device.deviceInternal.get_rtc_driver()
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::BrickHAT)
    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::BrickHAT, 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::BrickHAT)
    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::BrickHAT, 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::BrickHAT, 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::BrickHAT, 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::BrickHAT)
    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::BrickHAT)
    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::BrickHAT)
    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::BrickHAT, 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::BrickHAT)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the HAT is connected to
(typically '0' as the HAT is the root device in the topology),
the position, the hardware and firmware version as well as the
device identifier.

The HAT (Zero) Brick is always at position 'i'.

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

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


text:
path = ["../examples/opf/14bus/case.raw"; "../examples/opf/14bus/case.rop"]

M = create_PowersenseData(path)

run_opf!(M, solver = Ipopt.Optimizer, obj_type = "linear", formulation = PNPAPVmodel);

@test isapprox(M.cost, 1.467e+02, rtol=1e-1)

run_opf!(M, solver = Ipopt.Optimizer, obj_type = "linear", box_constraints = true, formulation = PBRARVmodel);

@test isapprox(M.cost, 1.467e+02, rtol=1e-1)

run_opf!(M, solver = Ipopt.Optimizer, obj_type = "linear", box_constraints = false, formulation = CBRARVmodel);

@test isapprox(M.cost, 1.467e+02, rtol=1e-1)

run_opf!(M, solver = Ipopt.Optimizer, obj_type = "linear", box_constraints = true, formulation = CBRARVmodel);

@test isapprox(M.cost, 1.467e+02, rtol=1e-1)

path2 = "../examples/opf/case300.m"

M2 = create_PowersenseData(path2)

run_opf!(M2, solver = Ipopt.Optimizer, obj_type = "quadratic", formulation = PNPAPVmodel);

@test isapprox(M.cost, 7.1973e+05, rtol=1e+3)


text:
module GXBeam

using ArnoldiMethod
using ForwardDiff
using LinearAlgebra
using LinearMaps
using NLsolve
using SparseArrays
using StaticArrays
using WriteVTK
using DiffEqBase

export curve_length, discretize_beam
export System, Assembly, PrescribedConditions, DistributedLoads
export combine_loads
export system_state, reset_state!, set_state!
export set_element_deflections!, set_element_rotations!
export set_element_forces!, set_element_moments!
export set_element_linear_momenta!, set_element_angular_momenta!
export set_point_deflections!, set_point_rotations!
export set_point_forces!, set_point_moments!
export static_analysis, static_analysis!
export steady_state_analysis, steady_state_analysis!
export eigenvalue_analysis, eigenvalue_analysis!
export initial_condition_analysis, initial_condition_analysis!
export time_domain_analysis, time_domain_analysis!
export AssemblyState
export PointState, extract_point_state, extract_point_states, extract_point_states!
export ElementState, extract_element_state, extract_element_states, extract_element_states!
export wiener_milenkovic
export rotate, rotate!
export translate, translate!
export cross_section_velocities, cross_section_velocities!
export deform_cross_section, deform_cross_section!

export write_vtk
export left_eigenvectors, correlate_eigenmodes

const GAUSS_NODES = SVector(-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526)
const GAUSS_WEIGHTS = SVector(0.34785484513745385, 0.6521451548625462, 0.6521451548625462, 0.34785484513745385)

include("math.jl")
include("element.jl")
include("point.jl")
include("loads.jl")
include("assembly.jl")
include("system.jl")
include("analyses.jl")
include("postprocess.jl")

include("interfaces/differentialequations.jl")
include("interfaces/forwarddiff.jl")

end


text:
"""
Subtypes of `AbstractRule` are types which represent the primitive derivative
propagation "rules" that can be composed to implement forward- and reverse-mode
automatic differentiation.

More specifically, a `rule::AbstractRule` is a callable Julia object generally
obtained via calling [`frule`](@ref) or [`rrule`](@ref). Such rules accept
differential values as input, evaluate the chain rule using internally stored/
computed partial derivatives to produce a single differential value, then
return that calculated differential value.

For example:

```julia-repl
julia> using ChainRules: frule, rrule, AbstractRule

julia> x, y = rand(2);

julia> h, dh = frule(hypot, x, y);

julia> h == hypot(x, y)
true

julia> isa(dh, AbstractRule)
true

julia> Δx, Δy = rand(2);

julia> dh(Δx, Δy) == ((y / h) * Δx + (x / h) * Δy)
true

julia> h, (dx, dy) = rrule(hypot, x, y);

julia> h == hypot(x, y)
true

julia> isa(dx, AbstractRule) && isa(dy, AbstractRule)
true

julia> Δh = rand();

julia> dx(Δh) == (y / h) * Δh
true

julia> dy(Δh) == (x / h) * Δh
true
```

See also: [`frule`](@ref), [`rrule`](@ref), [`Rule`](@ref), [`DNERule`](@ref), [`WirtingerRule`](@ref)
"""
abstract type AbstractRule end

# this ensures that consumers don't have to special-case rule destructuring
Base.iterate(rule::AbstractRule) = (rule, nothing)
Base.iterate(::AbstractRule, ::Any) = nothing

"""
    accumulate(Δ, rule::AbstractRule, args...)

Return `Δ + rule(args...)` evaluated in a manner that supports ChainRules'
various `AbstractDifferential` types.

This method intended to be customizable for specific rules/input types. For
example, here is pseudocode to overload `accumulate` w.r.t. a specific forward
differentiation rule for a given function `f`:

```
df(x) = # forward differentiation primitive implementation

frule(::typeof(f), x) = (f(x), Rule(df))

accumulate(Δ, rule::Rule{typeof(df)}, x) = # customized `accumulate` implementation
```

See also: [`accumulate!`](@ref), [`store!`](@ref), [`AbstractRule`](@ref)
"""
accumulate(Δ, rule::AbstractRule, args...) = add(Δ, rule(args...))

"""
    accumulate!(Δ, rule::AbstractRule, args...)

Similar to [`accumulate`](@ref), but compute `Δ + rule(args...)` in-place,
storing the result in `Δ`.

Note that this function internally calls `Base.Broadcast.materialize!(Δ, ...)`.

See also: [`accumulate`](@ref), [`store!`](@ref), [`AbstractRule`](@ref)
"""
function accumulate!(Δ, rule::AbstractRule, args...)
    return materialize!(Δ, broadcastable(add(cast(Δ), rule(args...))))
end

accumulate!(Δ::Number, rule::AbstractRule, args...) = accumulate(Δ, rule, args...)

"""
    store!(Δ, rule::AbstractRule, args...)

Compute `rule(args...)` and store the result in `Δ`, potentially avoiding
intermediate temporary allocations that might be necessary for alternative
approaches (e.g. `copyto!(Δ, extern(rule(args...)))`)

Note that this function internally calls `Base.Broadcast.materialize!(Δ, ...)`.

Like [`accumulate`](@ref) and [`accumulate!`](@ref), this function is intended
to be customizable for specific rules/input types.

See also: [`accumulate`](@ref), [`accumulate!`](@ref), [`AbstractRule`](@ref)
"""
store!(Δ, rule::AbstractRule, args...) = materialize!(Δ, broadcastable(rule(args...)))

# Special purpose updating for operations which can be done in-place. This function is
# just internal and free-form; it is not a method of `accumulate!` directly as it does
# not adhere to the expected method signature form, i.e. `accumulate!(value, rule, args)`.
# Instead it's `_update!(old, new, extrastuff...)` and is not specific to any particular
# rule.

_update!(x, y) = x + y
_update!(x::Array{T,N}, y::AbstractArray{T,N}) where {T,N} = x .+= y

_update!(x, ::Zero) = x
_update!(::Zero, y) = y
_update!(::Zero, ::Zero) = Zero()

function _update!(x::NamedTuple{Ns}, y::NamedTuple{Ns}) where Ns
    return NamedTuple{Ns}(map(p->_update!(getproperty(x, p), getproperty(y, p)), Ns))
end

function _update!(x::NamedTuple, y, p::Symbol)
    new = NamedTuple{(p,)}((_update!(getproperty(x, p), y),))
    return merge(x, new)
end

function _update!(x::NamedTuple{Ns}, y::NamedTuple{Ns}, p::Symbol) where Ns
    return _update!(x, getproperty(y, p), p)
end

#####
##### `Rule`
#####

Cassette.@context RuleContext

const RULE_CONTEXT = Cassette.disablehooks(RuleContext())

Cassette.overdub(::RuleContext, ::typeof(+), a, b) = add(a, b)
Cassette.overdub(::RuleContext, ::typeof(*), a, b) = mul(a, b)

Cassette.overdub(::RuleContext, ::typeof(add), a, b) = add(a, b)
Cassette.overdub(::RuleContext, ::typeof(mul), a, b) = mul(a, b)

"""
    Rule(propation_function[, updating_function])

Return a `Rule` that wraps the given `propation_function`. It is assumed that
`propation_function` is a callable object whose arguments are differential
values, and whose output is a single differential value calculated by applying
internally stored/computed partial derivatives to the input differential
values.

If an updating function is provided, it is assumed to have the signature `u(Δ, xs...)`
and to store the result of the propagation function applied to the arguments `xs` into
`Δ` in-place, returning `Δ`.

For example:

```
frule(::typeof(*), x, y) = x * y, Rule((Δx, Δy) -> Δx * y + x * Δy)

rrule(::typeof(*), x, y) = x * y, (Rule(ΔΩ -> ΔΩ * y'), Rule(ΔΩ -> x' * ΔΩ))
```

See also: [`frule`](@ref), [`rrule`](@ref), [`accumulate`](@ref), [`accumulate!`](@ref), [`store!`](@ref)
"""
struct Rule{F,U<:Union{Function,Nothing}} <: AbstractRule
    f::F
    u::U
end

# NOTE: Using `Core.Typeof` instead of `typeof` here so that if we define a rule for some
# constructor based on a `UnionAll`, we get `Rule{Type{Thing}}` instead of `Rule{UnionAll}`
Rule(f) = Rule{Core.Typeof(f),Nothing}(f, nothing)

(rule::Rule{F})(args...) where {F} = Cassette.overdub(RULE_CONTEXT, rule.f, args...)

# Specialized accumulation
# TODO: Does this need to be overdubbed in the rule context?
accumulate!(Δ, rule::Rule{F,U}, args...) where {F,U<:Function} = rule.u(Δ, args...)

#####
##### `DNERule`
#####

"""
    DNERule(args...)

Construct a `DNERule` object, which is an `AbstractRule` that signifies that the
current function is not differentiable with respect to a particular parameter.
**DNE** is an abbreviation for Does Not Exist.
"""
struct DNERule <: AbstractRule end

DNERule(args...) = DNE()

#####
##### `WirtingerRule`
#####

"""
TODO
"""
struct WirtingerRule{P<:AbstractRule,C<:AbstractRule} <: AbstractRule
    primal::P
    conjugate::C
end

function (rule::WirtingerRule)(args...)
    return Wirtinger(rule.primal(args...), rule.conjugate(args...))
end

#####
##### `frule`/`rrule`
#####

#=
In some weird ideal sense, the fallback for e.g. `frule` should actually be "get
the derivative via forward-mode AD". This is necessary to enable mixed-mode
rules, where e.g. `frule` is used within a `rrule` definition. For example,
broadcasted functions may not themselves be forward-mode *primitives*, but are
often forward-mode *differentiable*.

ChainRules, by design, is decoupled from any specific AD implementation. How,
then, do we know which AD to fall back to when there isn't a primitive defined?

Well, if you're a greedy AD implementation, you can just overload `frule` and/or
`rrule` to use your AD directly. However, this won't play nice with other AD
packages doing the same thing, and thus could cause load-order-dependent
problems for downstream users.

It turns out, Cassette solves this problem nicely by allowing AD authors to
overload the fallbacks w.r.t. their own context. Example using ForwardDiff:

```
using ChainRules, ForwardDiff, Cassette

Cassette.@context MyChainRuleCtx

# ForwardDiff, itself, can call `my_frule` instead of
# `frule` to utilize the ForwardDiff-injected ChainRules
# infrastructure
my_frule(args...) = Cassette.overdub(MyChainRuleCtx(), frule, args...)

function Cassette.execute(::MyChainRuleCtx, ::typeof(frule), f, x::Number)
    r = frule(f, x)
    if isa(r, Nothing)
        fx, df = (f(x), Rule(Δx -> ForwardDiff.derivative(f, x) * Δx))
    else
        fx, df = r
    end
    return fx, df
end
```
=#

"""
    frule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_ΔΩ₁::AbstractRule, rule_for_ΔΩ₂::AbstractRule, ...))

where each returned propagation rule `rule_for_ΔΩᵢ` can be invoked as

    rule_for_ΔΩᵢ(Δx₁, Δx₂, ...)

to yield `Ωᵢ`'s corresponding differential `ΔΩᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    ΔΩᵢ = ∂Ωᵢ_∂x₁ * Δx₁ + ∂Ωᵢ_∂x₂ * Δx₂ + ...

If no method matching `frule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```julia-repl
julia> x = rand();

julia> sinx, dsin = frule(sin, x);

julia> sinx == sin(x)
true

julia> dsin(1) == cos(x)
true
```

unary input, binary output scalar function:

```julia-repl
julia> x = rand();

julia> sincosx, (dsin, dcos) = frule(sincos, x);

julia> sincosx == sincos(x)
true

julia> dsin(1) == cos(x)
true

julia> dcos(1) == -sin(x)
true
```

See also: [`rrule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
frule(::Any, ::Vararg{Any}) = nothing

"""
    rrule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_Δx₁::AbstractRule, rule_for_Δx₂::AbstractRule, ...))

where each returned propagation rule `rule_for_Δxᵢ` can be invoked as

    rule_for_Δxᵢ(ΔΩ₁, ΔΩ₂, ...)

to yield `xᵢ`'s corresponding differential `Δxᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    Δxᵢ = ∂Ω₁_∂xᵢ * ΔΩ₁ + ∂Ω₂_∂xᵢ * ΔΩ₂ + ...

If no method matching `rrule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```julia-repl
julia> x = rand();

julia> sinx, dx = rrule(sin, x);

julia> sinx == sin(x)
true

julia> dx(1) == cos(x)
true
```

binary input, unary output scalar function:

```julia-repl
julia> x, y = rand(2);

julia> hypotxy, (dx, dy) = rrule(hypot, x, y);

julia> hypotxy == hypot(x, y)
true

julia> dx(1) == (y / hypot(x, y))
true

julia> dy(1) == (x / hypot(x, y))
true
```

See also: [`frule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
rrule(::Any, ::Vararg{Any}) = nothing

#####
##### macros
#####

"""
    @scalar_rule(f(x₁, x₂, ...),
                 @setup(statement₁, statement₂, ...),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

A convenience macro that generates simple scalar forward or reverse rules using
the provided partial derivatives. Specifically, generates the corresponding
methods for `frule` and `rrule`:

    function ChainRules.frule(::typeof(f), x₁, x₂, ...)
        Ω = f(x₁, x₂, ...)
        \$(statement₁, statement₂, ...)
        return Ω, (Rule((Δx₁, Δx₂, ...) -> ∂f₁_∂x₁ * Δx₁ + ∂f₁_∂x₂ * Δx₂ + ...),
                   Rule((Δx₁, Δx₂, ...) -> ∂f₂_∂x₁ * Δx₁ + ∂f₂_∂x₂ * Δx₂ + ...),
                   ...)
    end

    function ChainRules.rrule(::typeof(f), x₁, x₂, ...)
        Ω = f(x₁, x₂, ...)
        \$(statement₁, statement₂, ...)
        return Ω, (Rule((ΔΩ₁, ΔΩ₂, ...) -> ∂f₁_∂x₁ * ΔΩ₁ + ∂f₂_∂x₁ * ΔΩ₂ + ...),
                   Rule((ΔΩ₁, ΔΩ₂, ...) -> ∂f₁_∂x₂ * ΔΩ₁ + ∂f₂_∂x₂ * ΔΩ₂ + ...),
                   ...)
    end

Note that the result of `f(x₁, x₂, ...)` is automatically bound to `Ω`. This
allows the primal result to be conveniently referenced (as `Ω`) within the
derivative/setup expressions.

Note that the `@setup` argument can be elided if no setup code is need. In other
words:

    @scalar_rule(f(x₁, x₂, ...),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

is equivalent to:

    @scalar_rule(f(x₁, x₂, ...),
                 @setup(nothing),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

For examples, see ChainRules' `rules` directory.

See also: [`frule`](@ref), [`rrule`](@ref), [`AbstractRule`](@ref)
"""
macro scalar_rule(call, maybe_setup, partials...)
    if Meta.isexpr(maybe_setup, :macrocall) && maybe_setup.args[1] == Symbol("@setup")
        setup_stmts = map(esc, maybe_setup.args[3:end])
    else
        setup_stmts = (nothing,)
        partials = (maybe_setup, partials...)
    end
    @assert Meta.isexpr(call, :call)
    f, inputs = esc(call.args[1]), esc.(call.args[2:end])
    if all(Meta.isexpr(partial, :tuple) for partial in partials)
        forward_rules = Any[rule_from_partials(partial.args...) for partial in partials]
        reverse_rules = Any[]
        for i in 1:length(inputs)
            reverse_partials = [partial.args[i] for partial in partials]
            push!(reverse_rules, rule_from_partials(reverse_partials...))
        end
    else
        @assert length(inputs) == 1 && all(!Meta.isexpr(partial, :tuple) for partial in partials)
        forward_rules = Any[rule_from_partials(partial) for partial in partials]
        reverse_rules = Any[rule_from_partials(partials...)]
    end
    forward_rules = length(forward_rules) == 1 ? forward_rules[1] : Expr(:tuple, forward_rules...)
    reverse_rules = length(reverse_rules) == 1 ? reverse_rules[1] : Expr(:tuple, reverse_rules...)
    return quote
        function ChainRules.frule(::typeof($f), $(inputs...))
            $(esc(:Ω)) = $call
            $(setup_stmts...)
            return $(esc(:Ω)), $forward_rules
        end
        function ChainRules.rrule(::typeof($f), $(inputs...))
            $(esc(:Ω)) = $call
            $(setup_stmts...)
            return $(esc(:Ω)), $reverse_rules
        end
    end
end

function rule_from_partials(∂s...)
    wirtinger_indices = findall(x -> Meta.isexpr(x, :call) && x.args[1] === :Wirtinger,  ∂s)
    ∂s = map(esc, ∂s)
    Δs = [Symbol(string(:Δ, i)) for i in 1:length(∂s)]
    Δs_tuple = Expr(:tuple, Δs...)
    if isempty(wirtinger_indices)
        ∂_mul_Δs = [:(mul(@thunk($(∂s[i])), $(Δs[i]))) for i in 1:length(∂s)]
        return :(Rule($Δs_tuple -> add($(∂_mul_Δs...))))
    else
        ∂_mul_Δs_primal = Any[]
        ∂_mul_Δs_conjugate = Any[]
        ∂_wirtinger_defs = Any[]
        for i in 1:length(∂s)
            if i in wirtinger_indices
                Δi = Δs[i]
                ∂i = Symbol(string(:∂, i))
                push!(∂_wirtinger_defs, :($∂i = $(∂s[i])))
                ∂f∂i_mul_Δ = :(mul(wirtinger_primal($∂i), wirtinger_primal($Δi)))
                ∂f∂ī_mul_Δ̄ = :(mul(conj(wirtinger_conjugate($∂i)), wirtinger_conjugate($Δi)))
                ∂f̄∂i_mul_Δ = :(mul(wirtinger_conjugate($∂i), wirtinger_primal($Δi)))
                ∂f̄∂ī_mul_Δ̄ = :(mul(conj(wirtinger_primal($∂i)), wirtinger_conjugate($Δi)))
                push!(∂_mul_Δs_primal, :(add($∂f∂i_mul_Δ, $∂f∂ī_mul_Δ̄)))
                push!(∂_mul_Δs_conjugate, :(add($∂f̄∂i_mul_Δ, $∂f̄∂ī_mul_Δ̄)))
            else
                ∂_mul_Δ = :(mul(@thunk($(∂s[i])), $(Δs[i])))
                push!(∂_mul_Δs_primal, ∂_mul_Δ)
                push!(∂_mul_Δs_conjugate, ∂_mul_Δ)
            end
        end
        primal_rule = :(Rule($Δs_tuple -> add($(∂_mul_Δs_primal...))))
        conjugate_rule = :(Rule($Δs_tuple -> add($(∂_mul_Δs_conjugate...))))
        return quote
            $(∂_wirtinger_defs...)
            WirtingerRule($primal_rule, $conjugate_rule)
        end
    end
end


text:
using TensND, LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations
sympy.init_printing(use_unicode=true)

Polar = coorsys_polar() ; r, θ = getcoords(Polar) ; 𝐞ʳ, 𝐞ᶿ = unitvec(Polar)
@set_coorsys Polar
ℬˢ = get_normalized_basis(Polar)
𝐱 = getOM(Polar)
Cartesian = coorsys_cartesian(symbols("x y", real = true))
𝐞₁, 𝐞₂ = unitvec(Cartesian)
x₁, x₂ = getcoords(Cartesian)
𝕀, 𝕁, 𝕂 = ISO(Val(2),Val(Sym))
𝟏 = tensId2(Val(2),Val(Sym))

E, k, μ = symbols("E k μ", positive = true)
ν, κ = symbols("ν κ", real = true)
k = E / (3(1-2ν)) ; μ = E / (2(1+ν))
λ = k -2μ/3

𝐆 =simplify(1/(8 * PI * μ * (1-ν)) * (𝐞ʳ ⊗ 𝐞ʳ -(3-4ν) * ln(r) * 𝟏))
HG = -simplify(HESS(𝐆))
aHG = getarray(HG)
𝕄 = SymmetricTensor{4,2}((i,j,k,l)->(aHG[i,k,j,l]+aHG[j,k,i,l]+aHG[i,l,j,k]+aHG[j,l,i,k])/4)
ℾ = simplify(Tens(𝕄,ℬˢ))
ℾ₂ = simplify(1/(8PI * μ * (1-ν) * r^2) * (-2𝕁 +2(1-2ν)*𝕀 + 2(𝟏⊗𝐞ʳ⊗𝐞ʳ + 𝐞ʳ⊗𝐞ʳ⊗𝟏) + 8ν*𝐞ʳ⊗ˢ𝟏⊗ˢ𝐞ʳ -8𝐞ʳ⊗𝐞ʳ⊗𝐞ʳ⊗𝐞ʳ))
simplify(ℾ-ℾ₂)

ℂ = 2λ * 𝕁 + 2μ * 𝕀

𝕜 = simplify(ℾ ⊡ ℂ)
d = Dict(r => sqrt(x₁^2+x₂^2), sin(θ) => x₂/sqrt(x₁^2+x₂^2), cos(θ) => x₁/sqrt(x₁^2+x₂^2), ν => (3-κ)/4)

Spherical = coorsys_spherical() ; θ, ϕ, r = getcoords(Spherical) ; 𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical) ;
ℬˢ = get_normalized_basis(Spherical)
@set_coorsys Spherical
𝕀, 𝕁, 𝕂 = ISO(Val(3),Val(Sym))
𝟏 = tensId2(Val(3),Val(Sym))

E, k, μ = symbols("E k μ", positive = true)
ν = symbols("ν", real = true)
k = E / (3(1-2ν)) ; μ = E / (2(1+ν))
λ = k -2μ/3

𝐆 = 1/ (8PI * μ * (3k+4μ) * r) * ( (3k+7μ) * 𝟏 + (3k+μ) * 𝐞ʳ⊗𝐞ʳ)
𝐆₂ = 1 / (16PI * μ * (1-ν) * r) * ( (3-4ν) * 𝟏 + 𝐞ʳ⊗𝐞ʳ)
simplify(𝐆-𝐆₂)

HG = -simplify(HESS(𝐆))
aHG = getarray(HG)
𝕄 = SymmetricTensor{4,3}((i,j,k,l)->(aHG[i,k,j,l]+aHG[j,k,i,l]+aHG[i,l,j,k]+aHG[j,l,i,k])/4)
ℾ = simplify(Tens(𝕄,ℬˢ))
ℾ₂ = simplify(1/(16PI * μ * (1-ν) * r^3) * (-3𝕁 +2(1-2ν)*𝕀 + 3(𝟏⊗𝐞ʳ⊗𝐞ʳ + 𝐞ʳ⊗𝐞ʳ⊗𝟏) + 12ν*𝐞ʳ⊗ˢ𝟏⊗ˢ𝐞ʳ -15𝐞ʳ⊗𝐞ʳ⊗𝐞ʳ⊗𝐞ʳ))
simplify(ℾ-ℾ₂)

Cartesian = coorsys_cartesian(symbols("x y z", real = true))
𝐞₁, 𝐞₂, 𝐞3 = unitvec(Cartesian)
F = symbols("F", real = true)
J = simplify(det(𝟏 + F*GRAD(𝐆⋅𝐞₁)))
factor(simplify(subs(J, θ => PI/2, ϕ => 0)))


text:
function [J,T1,T2]=JacobiSphere2TT(ksi1,ksi2,F,Grid)
Rad=Grid.Rad;

X1=0.25*(F.P(1,1).*(1-ksi1)*(1-ksi2)+...
  F.P(1,2).*(1+ksi1)*(1-ksi2)+...
  F.P(1,3).*(1+ksi1)*(1+ksi2)+...
  F.P(1,4).*(1-ksi1)*(1+ksi2));
X2=0.25*(F.P(2,1).*(1-ksi1)*(1-ksi2)+...
  F.P(2,2).*(1+ksi1)*(1-ksi2)+...
  F.P(2,3).*(1+ksi1)*(1+ksi2)+...
  F.P(2,4).*(1-ksi1)*(1+ksi2));
X3=0.25*(F.P(3,1).*(1-ksi1)*(1-ksi2)+...
  F.P(3,2).*(1+ksi1)*(1-ksi2)+...
  F.P(3,3).*(1+ksi1)*(1+ksi2)+...
  F.P(3,4).*(1-ksi1)*(1+ksi2));

r=sqrt(X1^2+X2^2+X3^2);
f=Rad/r;
X1=X1/r;
X2=X2/r;
X3=X3/r;
[lam,theta]=cart2sphere(X1,X2,X3);

DD=[-sin(lam) cos(lam) 0
        0        0     1];
sinlam=sin(lam);
coslam=cos(lam);
sinth=sin(theta);
costh=cos(theta);
a11=sinlam*sinlam*costh*costh+sinth*sinth;
a12=-sinlam*coslam*costh*costh;
a13=-coslam*sinth*costh;
a21=a12;
a22=coslam*coslam*costh*costh+sinth*sinth;
a23=-sinlam*sinth*costh;
a31=-coslam*sinth;
a32=-sinlam*sinth;
a33=costh;
A=[a11 a12 a13
   a21 a22 a23
   a31 a32 a33];

B=[Grid.Nodes(F.N(1)).P(1) Grid.Nodes(F.N(2)).P(1) Grid.Nodes(F.N(3)).P(1) Grid.Nodes(F.N(4)).P(1)
   Grid.Nodes(F.N(1)).P(2) Grid.Nodes(F.N(2)).P(2) Grid.Nodes(F.N(3)).P(2) Grid.Nodes(F.N(4)).P(2)
   Grid.Nodes(F.N(1)).P(3) Grid.Nodes(F.N(2)).P(3) Grid.Nodes(F.N(3)).P(3) Grid.Nodes(F.N(4)).P(3)];
C=0.25*[-1+ksi2  -1+ksi1
         1-ksi2  -1-ksi1
         1+ksi2   1+ksi1
        -1-ksi2   1-ksi1];
J=f*DD*A*B*C;
D=abs(det(J));
X=[X1 X2 X3]*Rad;
T1=DD*A;
T2=B*C;
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license
# David Gleich modified it to support parsing Array{UInt8} as well

# Except now I'm not really using this code. See the myparse function
# below for a simple solution.
# Will delete soon!

import Base.Checked: add_with_overflow, mul_with_overflow
import Base.next

## string to integer functions ##

function parse(::Type{T}, c::UInt8, base::Integer=36) where T<:Integer
    a::Int = (base <= 36 ? 10 : 36)
    2 <= base <= 62 || throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
    d = UInt8('0') <= c <= UInt8('9') ? c-UInt8('0')    :
        UInt8('A') <= c <= UInt8('Z') ? c-UInt8('A')+10 :
        UInt8('a') <= c <= UInt8('z') ? c-UInt8('a')+a  : throw(ArgumentError("invalid digit: $(repr(c))"))
    d < base || throw(ArgumentError("invalid base $base digit $(repr(c))"))
    convert(T, d)
end

@inline isspacecode(c::UInt8) = c == UInt8(' ') || UInt8('\t') <= c <= UInt8('\r') || c == UInt8('\u85')

function parseint_next(s::Array{UInt8}, startpos::Int, endpos::Int)
    (0 < startpos <= endpos) || (return UInt8(0), 0, 0)
    j = startpos
    c, startpos = next(s,startpos)
    c, startpos, j
end

function parseint_preamble(signed::Bool, base::Int, s::Array{UInt8}, startpos::Int, endpos::Int)
    c, i, j = parseint_next(s, startpos, endpos)

    while isspacecode(c)
        c, i, j = parseint_next(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    sgn = 1
    if signed
        if c == UInt8('-') || c == UInt8('+')
            (c == UInt8('-')) && (sgn = -1)
            c, i, j = parseint_next(s,i,endpos)
        end
    end

    while isspacecode(c)
        c, i, j = parseint_next(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    if base == 0
        if c == UInt8('0') && !done(s,i)
            c, i = next(s,i)
            base = c==UInt8('b') ? 2 : c==UInt8('o') ? 8 : c==UInt8('x') ? 16 : 10
            if base != 10
                c, i, j = parseint_next(s,i,endpos)
            end
        else
            base = 10
        end
    end
    return sgn, base, j
end

function tryparse_internal(::Type{T}, s::Array{UInt8}, startpos::Int, endpos::Int, base_::Integer, raise::Bool) where T<:Integer
    _n = Nullable{T}()
    sgn, base, i = parseint_preamble(T<:Signed, Int(base_), s, startpos, endpos)
    if sgn == 0 && base == 0 && i == 0
        raise && throw(ArgumentError("input string is empty or only contains whitespace"))
        return _n
    end
    if !(2 <= base <= 62)
        raise && throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
        return _n
    end
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(SubString(s,startpos,endpos)))"))
        return _n
    end
    c, i = parseint_next(s,i,endpos)
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(SubString(s,startpos,endpos)))"))
        return _n
    end

    base = convert(T,base)
    m::T = div(typemax(T)-base+1,base)
    n::T = 0
    a::Int = base <= 36 ? 10 : 36
    while n <= m
        d::T = UInt8('0') <= c <= UInt8('9') ? c-UInt8('0')    :
               UInt8('A') <= c <= UInt8('Z') ? c-UInt8('A')+10 :
               UInt8('a') <= c <= UInt8('z') ? c-UInt8('a')+a  : base
        if d >= base
            raise && throw(ArgumentError("invalid base $base digit $(repr(c)) in $(repr(String(s[startpos:endpos])))"))
            return _n
        end
        n *= base
        n += d
        if i > endpos
            n *= sgn
            return Nullable{T}(n)
        end
        c, i = next(s,i)
        isspacecode(c) && break
    end
    (T <: Signed) && (n *= sgn)
    while !isspacecode(c)
        d::T = UInt8('0') <= c <= UInt8('9') ? c-UInt8('0')    :
               UInt8('A') <= c <= UInt8('Z') ? c-UInt8('A')+10 :
               UInt8('a') <= c <= UInt8('z') ? c-UInt8('a')+a  : base

        if d >= base
            raise && throw(ArgumentError("invalid base $base digit $(repr(c)) in $(repr(SubString(s,startpos,endpos)))"))
            return _n
        end
        (T <: Signed) && (d *= sgn)

        n, ov_mul = mul_with_overflow(n, base)
        n, ov_add = add_with_overflow(n, d)
        if ov_mul | ov_add
            raise && throw(OverflowError())
            return _n
        end
        (i > endpos) && return Nullable{T}(n)
        c, i = next(s,i)
    end
    while i <= endpos
        c, i = next(s,i)
        if !isspacecode(c)
            raise && throw(ArgumentError("extra characters after whitespace in $(repr(SubString(s,startpos,endpos)))"))
            return _n
        end
    end
    return Nullable{T}(n)
end

@inline function check_valid_base(base)
    if 2 <= base <= 62
        return base
    end
    throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
end

"""
    tryparse(type, str, [base])

Like [`parse`](@ref), but returns a [`Nullable`](@ref) of the requested type. The result
will be null if the string does not contain a valid number.
"""
tryparse(::Type{T}, s::Array{UInt8}, base::Integer) where {T<:Integer} =
    tryparse_internal(T, s, start(s), endof(s), check_valid_base(base), false)
tryparse(::Type{T}, s::Array{UInt8}) where {T<:Integer} =
    tryparse_internal(T, s, start(s), endof(s), 0, false)

#=
function myparse(::Type{T}, s::Array{UInt8}, base::Integer) where T<:Integer
    get(tryparse_internal(T, s, start(s), endof(s), check_valid_base(base), true))
end

function myparse(::Type{T}, s::Array{UInt8}) where T<:Integer
    get(tryparse_internal(T, s, start(s), endof(s), 0, true)) # Zero means, "figure it out"
end

function myparse(::Type{T}, s::Array{UInt8}, pos::Int64, len::Int64) where T<:Integer
    get(tryparse_internal(T, s, pos, len, 0, true)) # Zero means, "figure it out"
end
=#

function myparse(::Type{T}, a::Vector{UInt8}, start::Integer, len::Integer) where {T <: Integer}
  i = start
  n::T = zero(T)
  s::T = one(T)
  if a[1] == UInt8('-')
    s = -s
    i += 1
  elseif a[1] == UInt8('+')
    i += 1
  end
  max_without_overflow = div(typemax(T)-9,10) # the larg
  while i <= len && n <= max_without_overflow
    if UInt8('0') <= a[i] <= UInt8('9')
      n *= T(10)
      n += T(a[i]-UInt8('0'))
    else
      throw(ArgumentError("$(String(a[start:len])) is not a valid base 10 integer"))
    end
    i += 1
  end
  (T <: Signed) && (n *= s)
  f::Bool = false
  while i <= len
    if UInt8('0') <= a[i] <= UInt8('9')
      d = T(a[i]-UInt8('0'))
      (T <: Signed) && (d *= s)
      n10,f1 = mul_with_overflow(T(10), n)
      f = f|f1
      n,f1 = add_with_overflow(n10,d)
      f = f|f1
    else
      throw(ArgumentError("$(String(a[start:len])) is not a valid base 10 integer"))
    end
    i += 1
  end
  if f
    throw(ArgumentError("$(String(a[start:len])) caused overflow for type $(T)"))
  end
  return n
end
myparse(::Type{T}, a::Vector{UInt8}) where {T <: Integer}  = parse_int(T, a, 1, length(a))

@inline tryparse(::Type{Float64}, s::Vector{UInt8}, pos::Int64, len::Int64) =
            ccall(:jl_try_substrtod, Nullable{Float64}, (Ptr{UInt8},Csize_t,Csize_t), s, pos, len)

if VERSION < v"0.6"
    @inline function myparse(::Type{Float64}, s::Vector{UInt8}, pos::Int64, last::Int64)
        result = tryparse(Float64, s, pos-1, last-pos+1)
        if isnull(result)
            throw(ArgumentError("cannot parse $(repr(s)) as $Float64"))
        end
        return result.value
    end
else
    @inline function myparse(::Type{Float64}, s::Vector{UInt8}, pos::Int64, last::Int64)
        result = tryparse(Float64, s, pos-1, last-pos+1)
        if isnull(result)
            throw(ArgumentError("cannot parse $(repr(s)) as $Float64"))
        end
        return unsafe_get(result)
    end
end

ParsableNumbers = Union{Float16,Float32,Float64,Integer}


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using DataStructures

using PyCall
p = PyCall.pyimport("pyplasm")

using Plasm

include("./views.jl")
include("./graphic_text.jl")


text:
using Gnuplot
x = -2π:0.05:2π
cbwidth = 0.02
rightmargin = 0.875
@gp x sin.(x) "with linespoints ls 1 t 'sin' "
@gp :- "set style line 1 lc rgb 'black' lt 1 lw 1 pt 7 pi -4 ps 1.0"
@gp :- "set grid" "set auto fix"
@gp :- "set offsets graph .05, graph .05, graph .05, graph .05"
@gp :- "set key bottom left font ',10' box opaque"
# end of general propierties for first line and some general styling
@gp :- x[1:5:end] 0.3cos.(x)[1:5:end] "w p t '0.3cos' pt 6 ps 1 lc rgb 'black'"
@gp :- x -0.5cos.(x) "w l t '-0.5cos' lw 1.5 dashtype 1 lc rgb 'red' "
@gp :- x -0.65sin.(x) x  "w l t '-0.65sin' lw 5 lc palette " palette(:ice)
#colorbar style
@gp :- "set palette rgb '7,5,15' " "set cbtics out nomirror"
@gp :- "set colorbox user origin graph 1.01, graph 0 size $cbwidth, graph 1"
@gp :- "set rmargin at screen $rightmargin"
@gp :- "set colorbox noborder"
save(term="pngcairo size 700,400", output="mixlines.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#!/usr/bin/env julia

# Julia script to run SNaQ! on a test data set with 3 input parameters:
# - a given number of hybridizations
# - a given number of independent runs
# - a integer that will be mapped to a combination of other parameters.

# This is to be used with slurm, which defines a shell variable SLURM_ARRAY_TASK_ID.
# This variable will serve as the 3rd integer argument to the script,
# to determine the combination of parameters that each slurm task will run.

#---------------------------------------#
#         read script arguments         #
#---------------------------------------#

if length(ARGS)<3
  error("need 3 parameters: h (# hybrids), r (# runs) and combination ID")
end
h = parse(Int, ARGS[1]);
runs = parse(Int, ARGS[2]);
id = parse(Int, ARGS[3]);
println("h: ",h,", number of runs: ",runs,", combination id: ",id)

using PhyloNetworks # package that has the function snaq!

#---------------------------------------#
#         set function parameters       #
#         from last integer argument    #
#---------------------------------------#

# lists of parameter values for all parameters to be varied:
lFTA = [0.000001, 0.00001, 0.0001, 0.001, 0.01]
lNF    = [100, 75, 50, 25]
lRatio = [1, 100, 10000]   # controld LTA: Ratio=LTA/FTA
lXTR   = [0.001,    0.01]
lXTA   = [0.000001, 0.001]
# 5 parameters, with a total of 5*4*3*2*2 = 240 combinations
nparams = 5
nlevels = [length(lFTA),length(lNF),length(lRatio),length(lXTR),length(lXTA)]

# parameter fixed to a single value
FTR = 0.00001

"""
comb(index of parameter combination)

Take an integer as input, return a tuple of parameter values.
External objects are used: nparams, nlevels, and lFTA etc.
The integer input should be between 0 and 239, or
between 0 and the total # combinations -1 in general.
index 0 -> first values of all parameters
index 239 -> last values of all parameters
"""
function comb(combID)
  paramID = Vector{Int}(nparams)
  d = combID
  for par in 1:nparams
    d,r = fldmod(d, nlevels[par]) # combid = d * nlevels + r
    paramID[par] = r+1 # indexing in parameter list starts at 1, not 0
  end
  println("parameter levels: ",paramID)
  return lFTA[paramID[1]], lNF[paramID[2]], lRatio[paramID[3]], lXTR[paramID[4]], lXTA[paramID[5]]
end

FTA, NF, Ratio, XTR, XTA = comb(id)
LTA = FTA*Ratio;

#---------------------------------------#
#         read input, set seed          #
#---------------------------------------#

## modify the 2 lines below for the real run, which requires input files:
## replace dummy strings by what was commented out
tableCF = "data"  # readTableCF("tableCF.txt")
startingNet = "net" # readTopology(string("h",h,"BestStartingTree.out"))

# next: control the seed for the generation of random numbers, to make
#       the results reproducible (or to restart the analysis if some fail)
srand(h+1);
for i in 0:id
  global s = round(Int,floor(rand()*100000))
  # seed used below, each combination of parameter will have its own
end

#---------------------------------------#
#    set output file name and           #
#    run main command (or pretend)      #
#---------------------------------------#

rootname =  string("snaq/h",h,"nf",NF,"xta",XTA,"xtr",XTR,"fta",FTA,"ftr",FTR,"lta",LTA);
snaq_cmd = parse( # remove this line and triple quotes to run the real stuff
"""
snaq!(startingNet, tableCF, hmax=h, Nfail=NF, ftolAbs=FTA, ftolRel=FTR, xtolRel=XTR,xtolAbs=XTA, liktolAbs=LTA, runs=runs, seed=s, filename=rootname)
""")
println("real stuff: unquote and run this SNaQ command:")
println("       ", snaq_cmd)

println("trial: arguments to SNaQ command written to output file:")
println("       ", rootname)
open(rootname,"w") do fout # fout for File OUTput
  for arg in snaq_cmd.args[4:end]
    write(fout, arg.args[1],"=", string(eval(arg.args[2])),"\n")
  end
end # closes fout

info("all done, end of STDOUT")

max_stars_repo_name: morrowcj/computingtools

max_issues_repo_name: morrowcj/computingtools

max_forks_repo_name: morrowcj/computingtools


text:
# This script is modified from:
#   `https://github.com/evetion/GeoArrays.jl/blob/master/src/io.jl`
# Copyright (c) 2018 Maarten Pronk, MIT license

"""
    st_read(fn::AbstractString, bands = nothing)

# Return
A AbstractGeoArray object
"""
function st_read(fn::AbstractString, bands = nothing)
    isfile(fn) || error("File not found.")
    dataset = ArchGDAL.unsafe_read(fn)

    if bands === nothing
        nbands = ArchGDAL.nraster(dataset)
        bands = 1:nbands
    end
    A = ArchGDAL.read(dataset, bands)
    am = get_affine_map(dataset)

    # nodata masking
    # A = Array{Union{Missing, eltype(A)}}(A)
    mask = falses(size(A))

    for i = bands
        band = ArchGDAL.getband(dataset, i)
        maskflags = mask_flags(band)

        # All values are valid, skip masking
        if :GMF_ALL_VALID in maskflags
            @debug "No masking"
            continue
        # Mask is valid for all bands
        elseif :GMF_PER_DATASET in maskflags
            @debug "Mask for each band"
            maskband = ArchGDAL.getmaskband(band)
            m = ArchGDAL.read(maskband) .== 0
            mask[:,:,i] = m
        # Alpha layer
        elseif :GMF_ALPHA in maskflags
            @warn "Dataset has band $i with an Alpha band, which is unsupported for now."
            continue
        # Nodata values
        elseif :GMF_NODATA in maskflags
            @debug "Flag NODATA"
            nodata = get_nodata(band)
            mask[:, :, i] = A[:,:,i] .== nodata
        else
            @warn "Unknown/unsupported mask."
        end
    end

    if any(mask)
        A = Array{Union{Missing,eltype(A)}}(A)
        A[mask] .= missing
    end

    wkt = ArchGDAL.getproj(dataset) |> crs2wkt
    ArchGDAL.destroy(dataset)
    names = bandnames(fn)
    # @show names
    GeoArray(A, am, wkt, names)
end

"""
    readGDAL(file::String, options...)
    readGDAL(files::Array{String,1}, options)

# Arguments:
- `options`: other parameters to `ArchGDAL.read(dataset, options...)`.

# Return
"""
function readGDAL(file::AbstractString, options...)
    ArchGDAL.read(file) do dataset
        ArchGDAL.read(dataset, options...)
    end
end

# read multiple tiff files and cbind
function readGDAL(files::Vector{<:AbstractString}, options...)
    # bands = collect(bands)
    # bands = collect(Int32, bands)
    res = map(file -> readGDAL(file, options...), files)
    res
    # vcat(res...)
end

const read = st_read
export st_read, readGDAL


text:
"""
Returns a handle (>=0) for the device at the I2C bus address.

i2c_bus:= >=0.
i2c_address:= 0-0x7F.
i2c_flags:= 0, no flags are currently defined.

Normally you would only use the [*i2c_**] functions if
you are or will be connecting to the Pi over a network.  If
you will always run on the local Pi use the standard SMBus
module instead.

Physically buses 0 and 1 are available on the Pi.  Higher
numbered buses will be available if a kernel supported bus
multiplexor is being used.

For the SMBus commands the low level transactions are shown
at the end of the function description.  The following
abbreviations are used.

. .
S     (1 bit) : Start bit
P     (1 bit) : Stop bit
Rd/Wr (1 bit) : Read/Write bit. Rd equals 1, Wr equals 0.
A, NA (1 bit) : Accept and not accept bit.
Addr  (7 bits): I2C 7 bit address.
reg   (8 bits): Command byte, which often selects a register.
Data  (8 bits): A data byte.
Count (8 bits): A byte defining the length of a block operation.

[..]: Data sent by the device.
. .

...
h = i2c_open(pi, 1, 0x53) # open device at address 0x53 on bus 1
...
"""
function i2c_open(self::Pi, i2c_bus, i2c_address, i2c_flags=0)
    # I p1 i2c_bus
    # I p2 i2c_addr
    # I p3 4
    ## extension ##
    # I i2c_flags
    extents = [pack("I", i2c_flags)]
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_I2CO, i2c_bus, i2c_address, 4, extents))
end

"""
Closes the I2C device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).

...
i2c_close(pi, h)
...
"""
function i2c_close(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_I2CC, handle, 0))
end

"""
Sends a single bit to the device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
bit:= 0 or 1, the value to write.

SMBus 2.0 5.5.1 - Quick command.
. .
S Addr bit [A] P
. .

...
i2c_write_quick(pi, 0, 1) # send 1 to device 0
i2c_write_quick(pi, 3, 0) # send 0 to device 3
...
"""
function i2c_write_quick(self::Pi, handle, bit)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_I2CWQ, handle, bit))
end

"""
Sends a single byte to the device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
byte_val:= 0-255, the value to write.

SMBus 2.0 5.5.2 - Send byte.
. .
S Addr Wr [A] byte_val [A] P
. .

...
i2c_write_byte(pi, 1, 17)   # send byte   17 to device 1
i2c_write_byte(pi, 2, 0x23) # send byte 0x23 to device 2
...
"""
function i2c_write_byte(self::Pi, handle, byte_val)
    return _u2i(
        _pigpio_command(self.sl, _PI_CMD_I2CWS, handle, byte_val))
end

"""
Reads a single byte from the device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).

SMBus 2.0 5.5.3 - Receive byte.
. .
S Addr Rd [A] [Data] NA P
. .

...
b = i2c_read_byte(pi, 2) # read a byte from device 2
...
"""
function i2c_read_byte(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_I2CRS, handle, 0))
end

"""
Writes a single byte to the specified register of the device
associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
byte_val:= 0-255, the value to write.

SMBus 2.0 5.5.4 - Write byte.
. .
S Addr Wr [A] reg [A] byte_val [A] P
. .

...
# send byte 0xC5 to reg 2 of device 1
i2c_write_byte_data(pi, 1, 2, 0xC5)

# send byte 9 to reg 4 of device 2
i2c_write_byte_data(pi, 2, 4, 9)
...
"""
function i2c_write_byte_data(self::Pi, handle, reg, byte_val)
    # I p1 handle
    # I p2 reg
    # I p3 4
    ## extension ##
    # I byte_val
    extents = IOBuffer
    Base.write(extents, byte_val)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_I2CWB, handle, reg, 4, extents))
end

"""
Writes a single 16 bit word to the specified register of the
device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
word_val:= 0-65535, the value to write.

SMBus 2.0 5.5.4 - Write word.
. .
S Addr Wr [A] reg [A] word_val_Low [A] word_val_High [A] P
. .

...
# send word 0xA0C5 to reg 5 of device 4
i2c_write_word_data(pi, 4, 5, 0xA0C5)

# send word 2 to reg 2 of device 5
i2c_write_word_data(pi, 5, 2, 23)
...
"""
function i2c_write_word_data(self::Pi, handle, reg, word_val)
    # I p1 handle
    # I p2 reg
    # I p3 4
    ## extension ##
    # I word_val
    extents = IOBuffer
    write(extents, word_val)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_I2CWW, handle, reg, 4, extents))
end

"""
Reads a single byte from the specified register of the device
associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.

SMBus 2.0 5.5.5 - Read byte.
. .
S Addr Wr [A] reg [A] S Addr Rd [A] [Data] NA P
. .

...
# read byte from reg 17 of device 2
b = i2c_read_byte_data(pi, 2, 17)

# read byte from reg  1 of device 0
b = i2c_read_byte_data(pi, 0, 1)
...
"""
function i2c_read_byte_data(self, handle, reg)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_I2CRB, handle, reg))
end

"""
Reads a single 16 bit word from the specified register of the
device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.

SMBus 2.0 5.5.5 - Read word.
. .
S Addr Wr [A] reg [A] S Addr Rd [A] [DataLow] A [DataHigh] NA P
. .

...
# read word from reg 2 of device 3
w = i2c_read_word_data(pi, 3, 2)

# read word from reg 7 of device 2
w = i2c_read_word_data(pi, 2, 7)
...
"""
function i2c_read_word_data(self, handle, reg)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_I2CRW, handle, reg))
end

"""
Writes 16 bits of data to the specified register of the device
associated with handle and reads 16 bits of data in return.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
word_val:= 0-65535, the value to write.

SMBus 2.0 5.5.6 - Process call.
. .
S Addr Wr [A] reg [A] word_val_Low [A] word_val_High [A]
S Addr Rd [A] [DataLow] A [DataHigh] NA P
. .

...
r = i2c_process_call(pi, h, 4, 0x1231)
r = i2c_process_call(pi, h, 6, 0)
...
"""
function i2c_process_call(self, handle, reg, word_val)
    # I p1 handle
    # I p2 reg
    # I p3 4
    ## extension ##
    # I word_val
    extents = IOBuffer
    write(extents, word_val)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_I2CPC, handle, reg, 4, extents))
end

"""
Writes up to 32 bytes to the specified register of the device
associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
data:= the bytes to write.

SMBus 2.0 5.5.7 - Block write.
. .
S Addr Wr [A] reg [A] length(data) [A] data0 [A] data1 [A] ... [A]
datan [A] P
. .

...
i2c_write_block_data(pi, 4, 5, b'hello')

i2c_write_block_data(pi, 4, 5, "data bytes")

i2c_write_block_data(pi, 5, 0, b'\\x00\\x01\\x22')

i2c_write_block_data(pi, 6, 2, [0, 1, 0x22])
...
"""
function i2c_write_block_data(self, handle, reg, data)
# I p1 handle
# I p2 reg
# I p3 len
## extension ##
# s len data bytes
    if length(data)
        return _u2i(_pigpio_command_ext(
            self.sl, _PI_CMD_I2CWK, handle, reg, length(data), data))
    else
        return 0
    end
end

"""
Reads a block of up to 32 bytes from the specified register of
the device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.

SMBus 2.0 5.5.7 - Block read.
. .
S Addr Wr [A] reg [A]
S Addr Rd [A] [Count] A [Data] A [Data] A ... A [Data] NA P
. .

The amount of returned data is set by the device.

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(b, d) = i2c_read_block_data(pi, h, 10)
if b >= 0
# process data
else
# process read failure
...
"""
function i2c_read_block_data(self::Pi, handle, reg)
    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command(self.sl, _PI_CMD_I2CRK, handle, reg, false))
    if bytes > 0
        data = rxbuf(bytes)
    else
        data = ""
    end

    unlock(self.sl.l)
    return bytes, data
end

"""
Writes data bytes to the specified register of the device
associated with handle and reads a device specified number
of bytes of data in return.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
data:= the bytes to write.

The SMBus 2.0 documentation states that a minimum of 1 byte may
be sent and a minimum of 1 byte may be received.  The total
number of bytes sent/received must be 32 or less.

SMBus 2.0 5.5.8 - Block write-block read.
. .
S Addr Wr [A] reg [A] length(data) [A] data0 [A] ... datan [A]
S Addr Rd [A] [Count] A [Data] ... A P
. .

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(b, d) = i2c_block_process_call(pi, h, 10, b'\\x02\\x05\\x00')

(b, d) = i2c_block_process_call(pi, h, 10, b'abcdr')

(b, d) = i2c_block_process_call(pi, h, 10, "abracad")

(b, d) = i2c_block_process_call(pi, h, 10, [2, 5, 16])
...
"""
function i2c_block_process_call(self::Pi, handle, reg, data)
    # I p1 handle
    # I p2 reg
    # I p3 len
    ## extension ##
    # s len data bytes

    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_I2CPK, handle, reg, length(data), data, false))
    if bytes > 0
        data = rxbuf(self, bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Writes data bytes to the specified register of the device
associated with handle .  1-32 bytes may be written.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
data:= the bytes to write.

. .
S Addr Wr [A] reg [A] data0 [A] data1 [A] ... [A] datan [NA] P
. .

...
i2c_write_i2c_block_data(pi, 4, 5, 'hello')

i2c_write_i2c_block_data(pi, 4, 5, b'hello')

i2c_write_i2c_block_data(pi, 5, 0, b'\\x00\\x01\\x22')

i2c_write_i2c_block_data(pi, 6, 2, [0, 1, 0x22])
...
"""
function i2c_write_i2c_block_data(self::Pi, handle, reg, data)
    # I p1 handle
    # I p2 reg
    # I p3 len
    ## extension ##
    # s len data bytes
    if length(data) > 0
        return _u2i(_pigpio_command_ext(
            self.sl, _PI_CMD_I2CWI, handle, reg, length(data), [data]))
    else
        return 0
    end
end

"""
Reads count bytes from the specified register of the device
associated with handle .  The count may be 1-32.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
reg:= >=0, the device register.
count:= >0, the number of bytes to read.

. .
S Addr Wr [A] reg [A]
S Addr Rd [A] [Data] A [Data] A ... A [Data] NA P
. .

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(b, d) = i2c_read_i2c_block_data(pi, h, 4, 32)
if b >= 0
# process data
else
# process read failure
...
"""
function i2c_read_i2c_block_data(self::Pi, handle, reg, count)
    # I p1 handle
    # I p2 reg
    # I p3 4
    ## extension ##
    # I count
    extents = IOBuffer()
    write(extents, count)
    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_I2CRI, handle, reg, 4, extents, false))
    if bytes > 0
        data = rxbuf(self, bytes)
    else
        data = ""
    end

    unlock(self.sl.l)
    return bytes, data
end

"""
Returns count bytes read from the raw device associated
with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
count:= >0, the number of bytes to read.

. .
S Addr Rd [A] [Data] A [Data] A ... A [Data] NA P
. .

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(count, data) = i2c_read_device(pi, h, 12)
...
"""
function i2c_read_device(self::Pi, handle, count)
    # Don't raise exception.  Must release lock.
    bytes = u2i(
    _pigpio_command(self.sl, _PI_CMD_I2CRD, handle, count, false))
    if bytes > 0
        data = rxbuf(self, bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Writes the data bytes to the raw device associated with handle.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
data:= the bytes to write.

. .
S Addr Wr [A] data0 [A] data1 [A] ... [A] datan [A] P
. .

...
i2c_write_device(pi, h, b"\\x12\\x34\\xA8")

i2c_write_device(pi, h, b"help")

i2c_write_device(pi, h, 'help')

i2c_write_device(pi, h, [23, 56, 231])
...
"""
function i2c_write_device(self::Pi, handle, data)
    # I p1 handle
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes
    if length(data)
        return _u2i(_pigpio_command_ext(
            self.sl, _PI_CMD_I2CWD, handle, 0, length(data), data))
    else
        return 0
    end
end

"""
This function executes a sequence of I2C operations.  The
operations to be performed are specified by the contents of data
which contains the concatenated command codes and associated data.

handle:= >=0 (as returned by a prior call to [*i2c_open*]).
data:= the concatenated I2C commands, see below

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(count, data) = i2c_zip(pi, h, [4, 0x53, 7, 1, 0x32, 6, 6, 0])
...

The following command codes are supported

Name    @ Cmd & Data @ Meaning
End     @ 0          @ No more commands
Escape  @ 1          @ Next P is two bytes
On      @ 2          @ Switch combined flag on
Off     @ 3          @ Switch combined flag off
Address @ 4 P        @ Set I2C address to P
Flags   @ 5 lsb msb  @ Set I2C flags to lsb + (msb << 8)
Read    @ 6 P        @ Read P bytes of data
Write   @ 7 P ...    @ Write P bytes of data

The address, read, and write commands take a parameter P.
Normally P is one byte (0-255).  If the command is preceded by
the Escape command then P is two bytes (0-65535, least significant
byte first).

The address defaults to that associated with the handle.
The flags default to 0.  The address and flags maintain their
previous value until updated.

Any read I2C data is concatenated in the returned bytearray.

...
Set address 0x53, write 0x32, read 6 bytes
Set address 0x1E, write 0x03, read 6 bytes
Set address 0x68, write 0x1B, read 8 bytes
End

0x04 0x53   0x07 0x01 0x32   0x06 0x06
0x04 0x1E   0x07 0x01 0x03   0x06 0x06
0x04 0x68   0x07 0x01 0x1B   0x06 0x08
0x00
...
"""
function i2c_zip(self::Pi, handle, data)
    # I p1 handle
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes

    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_I2CZ, handle, 0, length(data), data, false))
    if bytes > 0
        data = self._rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
This function selects a pair of GPIO for bit banging I2C at a
specified baud rate.

Bit banging I2C allows for certain operations which are not possible
with the standard I2C driver.

o baud rates as low as 50
o repeated starts
o clock stretching
o I2C on any pair of spare GPIO

SDA:= 0-31
SCL:= 0-31
baud:= 50-500000

Returns 0 if OK, otherwise PI_BAD_USER_GPIO, PI_BAD_I2C_BAUD, or
PI_GPIO_IN_USE.

NOTE

The GPIO used for SDA and SCL must have pull-ups to 3V3 connected.
As a guide the hardware pull-ups on pins 3 and 5 are 1k8 in value.

...
h = bb_i2c_open(pi, 4, 5, 50000) # bit bang on GPIO 4/5 at 50kbps
...
"""
function bb_i2c_open(self::Pi, SDA, SCL, baud=100000)
    # I p1 SDA
    # I p2 SCL
    # I p3 4
    ## extension ##
    # I baud
    extents = IOBuffer()
    write(extents, baud)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_BI2CO, SDA, SCL, 4, extents))
end

"""
This function stops bit banging I2C on a pair of GPIO
previously opened with [*bb_i2c_open*].

SDA:= 0-31, the SDA GPIO used in a prior call to [*bb_i2c_open*]

Returns 0 if OK, otherwise PI_BAD_USER_GPIO, or PI_NOT_I2C_GPIO.

...
bb_i2c_close(pi, SDA)
...
"""
function bb_i2c_close(self::Pi, SDA)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_BI2CC, SDA, 0))
end

"""
This function executes a sequence of bit banged I2C operations.
The operations to be performed are specified by the contents
of data which contains the concatenated command codes and
associated data.

SDA:= 0-31 (as used in a prior call to [*bb_i2c_open*])
data:= the concatenated I2C commands, see below

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(count, data) = pi.bb_i2c_zip(
             h, [4, 0x53, 2, 7, 1, 0x32, 2, 6, 6, 3, 0])
...

The following command codes are supported

Name    @ Cmd & Data   @ Meaning
End     @ 0            @ No more commands
Escape  @ 1            @ Next P is two bytes
Start   @ 2            @ Start condition
Stop    @ 3            @ Stop condition
Address @ 4 P          @ Set I2C address to P
Flags   @ 5 lsb msb    @ Set I2C flags to lsb + (msb << 8)
Read    @ 6 P          @ Read P bytes of data
Write   @ 7 P ...      @ Write P bytes of data

The address, read, and write commands take a parameter P.
Normally P is one byte (0-255).  If the command is preceded by
the Escape command then P is two bytes (0-65535, least significant
byte first).

The address and flags default to 0.  The address and flags maintain
their previous value until updated.

No flags are currently defined.

Any read I2C data is concatenated in the returned bytearray.

...
Set address 0x53
start, write 0x32, (re)start, read 6 bytes, stop
Set address 0x1E
start, write 0x03, (re)start, read 6 bytes, stop
Set address 0x68
start, write 0x1B, (re)start, read 8 bytes, stop
End

0x04 0x53
0x02 0x07 0x01 0x32   0x02 0x06 0x06 0x03

0x04 0x1E
0x02 0x07 0x01 0x03   0x02 0x06 0x06 0x03

0x04 0x68
0x02 0x07 0x01 0x1B   0x02 0x06 0x08 0x03

0x00
...
"""
function bb_i2c_zip(self::Pi, SDA, data)
    # I p1 SDA
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes

    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_BI2CZ, SDA, 0, length(data), [data], false))
    if bytes > 0
        data = self._rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end


text:
@testset "l2inner" begin
    ndim = 10
    a = rand(Float,ndim)
    b = rand(Float,ndim)
    ca = rand(Complex{Float},ndim)
    cb = rand(Complex{Float},ndim)
    @test a⋅b ≈ VectorDataUtils.l2inner(a,b)
    @test a⋅cb ≈ VectorDataUtils.l2inner(a,cb)
    @test ca⋅b ≈ VectorDataUtils.l2inner(ca,b)
    @test ca⋅cb ≈ VectorDataUtils.l2inner(ca,cb)
    @test a⋅a ≈ VectorDataUtils.l2inner(a)
    @test a⋅a ≈ VectorDataUtils.sql2norm(a)
    @test sqrt(a⋅a) ≈ VectorDataUtils.l2norm(a)
end

@testset "l2dist" begin
    ndim = 10
    a = rand(Float,ndim)
    b = rand(Float,ndim)
    @test (a-b)⋅(a-b) ≈ VectorDataUtils.sql2dist(a,b)
    @test sqrt(sum(abs2,a-b)) ≈ VectorDataUtils.l2dist(a,b)
end

@testset "lpnorms" begin
    ndim = 10
    a = rand(Float,ndim)
    @test sum(abs.(a)) ≈ VectorDataUtils.l1norm(a)
end

@testset "lpdist" begin
    ndim = 10
    a = rand(Float,ndim)
    b = rand(Float,ndim)
    @test sum(abs.(a-b)) ≈ VectorDataUtils.l1dist(a,b)
end


text:
"""
    Atlas <: Environment

    humanoid robot designed by Boston Dynamics (v5)
"""
struct Atlas end

function atlas(;
    representation=:minimal,
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    damper=0.0,
    spring=0.0,
    parse_damper=true,
    info=nothing,
    model_type=:simple,
    seed=1,
    contact_feet=true,
    contact_body=false,
    vis=Visualizer(),
    name=:robot,
    infeasible_control=false,
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_mechanism(:atlas;
        timestep,
        gravity,
        friction_coefficient,
        damper,
        spring,
        parse_damper,
        contact_feet,
        contact_body,
        model_type)

    initialize!(mechanism, :atlas)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu_inf = input_dimension(mechanism)
    nu = infeasible_control ? nu_inf : nu_inf - 6 # remove first 6 controls
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)
    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z
    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = infeasible_control ? I(nu) : [zeros(nu, 6) I(nu)]
    control_scaling = Diagonal(ones(nu))

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Atlas, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)
end


text:
using OrdinaryDiffEq, DifferentialEquations

tspan = (0.0, 200.0)

"""
Rossler system with chaotic parameter values
"""
function rossler_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    a, b, c = p
    x, y, z = u

    Ft = U_func_time(t)
    Fs = U_func_space(u)

    du[1] = -y - z
    du[2] = x + a*y
    du[3] = b + x*z - c*z
    du .+= Ft .+ Fs
end

# Generate data
p = [0.1, 0.1, 14]
u0 = [1.0, 1.0, 0.0]
ts = range(tspan[1], tspan[2], length=5000)

function solve_rossler_system(;U_func_time=U_func_time_trivial,
                              U_func_space=U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->rossler_system(du, u, p, t,
                                U_func_time=U_func_time,
                                U_func_space=U_func_space),
                    u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=ts);
    return Array(sol)
end

# True model in SINDy syntax
#     x y z c xx xy xz yy yz zz
 A = [[0 -1 -1 0   0 0 0 0 0 0];
      [1 0.1 0 0   0 0 0 0 0 0];
      [0 0 -14 0.1 0 0 1 0 0 0]]
n = size(A, 1)
sindy_library = Dict(
    "cross_terms"=>2,
    "constant"=>nothing
);
core_dyn_true = sindyModel(ts, A,
                            convert_string2function(sindy_library),
                            ["x", "y", "z"])

# Actually export
export ts, solve_rossler_system, rossler_system, core_dyn_true


text:
using BinDeps

@BinDeps.setup

deps = [ arduino = library_dependency("Arduino", aliases = ["libarduino"]) ]

@linux_only begin
    prefix = joinpath(BinDeps.depsdir(arduino), "usr")
    ardsrcdir = joinpath(BinDeps.depsdir(arduino), "src", "arduino-serial-lib")
    ardbuilddir = joinpath(BinDeps.depsdir(arduino), "builds", "arduino-serial-lib")
    provides(BuildProcess,
            (@build_steps begin
                CreateDirectory(ardbuilddir)
                @build_steps begin
                    ChangeDirectory(ardbuilddir)
                    FileRule(joinpath(prefix, "lib", "libarduino.so"), @build_steps begin
                        `g++ -c -O -Wall -fpic -std=gnu99 $ardsrcdir/arduino-serial-lib.c`
                        `g++ -shared -o libarduino.so arduino-serial-lib.o`
                        `cp libarduino.so $prefix/lib`
                        `cp arduino-serial-lib.h $prefix/include`
                    end)
                end
             end), arduino, os = :Linux)
    run()
end

@osx_only begin
    prefix = joinpath(BinDeps.depsdir(arduino), "usr")
    ardsrcdir = joinpath(BinDeps.depsdir(arduino), "src", "arduino-serial-lib")
    ardbuilddir = joinpath(BinDeps.depsdir(arduino), "builds", "arduino-serial-lib")
    provides(BuildProcess,
            (@build_steps begin
                CreateDirectory(ardbuilddir)
                @build_steps begin
                    ChangeDirectory(ardbuilddir)
                    FileRule(joinpath(prefix, "lib", "libarduino.dylib"), @build_steps begin
                        `g++ -c -arch i386 -arch x86_64 -mmacosx-version-min=10.6 -O -Wall -fpic -std=gnu99 $ardsrcdir/arduino-serial-lib.c`
                        `g++ -dynamiclib -undefined suppress -flat_namespace arduino-serial-lib.o -o libarduino.dylib`
                        `cp libarduino.dylib $prefix/lib`
                        `cp arduino-serial-lib.h $prefix/include`
                    end)
                end
             end), arduino, os = :Darwin)
    run()
end

@windows_only begin
    prefix = joinpath(BinDeps.depsdir(arduino), "usr")
    ardsrcdir = joinpath(BinDeps.depsdir(arduino), "src", "arduino-serial-lib")
    ardbuilddir = joinpath(BinDeps.depsdir(arduino), "builds", "arduino-serial-lib")
    provides(BuildProcess,
            (@build_steps begin
                CreateDirectory(ardbuilddir)
                @build_steps begin
                    ChangeDirectory(ardbuilddir)
                    FileRule(joinpath(prefix, "lib", "libarduino.dll"), @build_steps begin
                        `g++ -c -O -Wall -std=gnu99 $ardsrcdir/arduino-serial-lib.c`
                        `g++ -shared -o libarduino.dll arduino-serial-lib.o -Wl,--out-implib,libarduino.a
                        `cp libarduino.dll $prefix/lib`
                        `cp libarduino.a $prefix/lib`
                        `cp arduino-serial-lib.h $prefix/include`
                    end)
                end
             end), arduino, os = :Windows)
    run()
end

@BinDeps.install


text:
module PackageCompiler

using Pkg, Serialization, Libdl, UUIDs
using Pkg: TOML,  Operations, Types

include("compiler_flags.jl")
include("static_julia.jl")
include("api.jl")
include("snooping.jl")
include("system_image.jl")
include("pkg.jl")
include("incremental.jl")

const sysimage_binaries = ("sys.$(Libdl.dlext)",)

function copy_system_image(src, dest, ignore_missing = false)
    for file in sysimage_binaries
        # backup
        srcfile = joinpath(src, file)
        destfile = joinpath(dest, file)
        if !isfile(srcfile)
            ignore_missing && continue
            error("No file: $srcfile")
        end
        if isfile(destfile)
            if isfile(destfile * ".backup")
                rm(destfile * ".backup", force = true)
            end
            mv(destfile, destfile * ".backup", force = true)
        end
        @info "Copying system image: $srcfile to $destfile"
        cp(srcfile, destfile, force = true)
    end
end

julia_cpu_target(x) = error("CPU target needs to be a string or `nothing`")
julia_cpu_target(x::String) = x # TODO: match against available targets
function julia_cpu_target(::Nothing)
    replace(Base.julia_cmd().exec[2], "-C" => "")
end

"""
Reverts a forced compilation of the system image.
This will restore any previously backed up system image files, or
build a new, clean system image.
"""
function revert(debug = false)
    syspath = default_sysimg_path(debug)
    sysimg_backup = dirname(get_backup!(debug))
    copy_system_image(sysimg_backup, syspath)
end

function get_root_dir(path)
    path, name = splitdir(path)
    if isempty(name)
        return splitdir(path)[2]
    else
        name
    end
end

function sysimg_folder(files...)
    base_path = normpath(abspath(joinpath(@__DIR__, "..", "sysimg")))
    isdir(base_path) || mkpath(base_path)
    normpath(abspath(joinpath(base_path, files...)))
end

function sysimgbackup_folder(files...)
    backup = sysimg_folder("backup")
    isdir(backup) || mkpath(backup)
    sysimg_folder("backup", files...)
end

function package_folder(package...)
    packages = normpath(abspath(joinpath(@__DIR__, "..", "packages")))
    isdir(packages) || mkpath(packages)
    normpath(abspath(joinpath(packages, package...)))
end

"""
    compile_package(packages...; kw_args...)

with packages being either a string naming a package, or a tuple `(package_name, precompile_file)`.
If no precompile file is given, it will use the packages `runtests.jl`, which is a good canditate
for figuring out what functions to compile!
"""
function compile_package(packages...; kw_args...)
    args = map(packages) do package
        # If no explicit path to a seperate precompile file, use runtests
        isa(package, String) && return (package, "test/runtests.jl")
        isa(package, Tuple{String, String}) && return package
        error("Unrecognized package. Use `packagename::String`, or `(packagename::String, rel_path_to_testfile::String)`. Found: `$package`")
    end
    compile_package(args...; kw_args...)
end

"""
    compile_package(
        packages::Tuple{String, String}...;
        force = false, reuse = false, debug = false, cpu_target = nothing,
        additional_packages = Symbol[]
    )

Compile a list of packages. Each package comes as a tuple of `(package_name, precompile_file)`
where the precompile file should contain all function calls, that should get compiled into the system image.
Usually the `runtests.jl` file is a good candidate, since it should run all important functions of a package.
You can pass `additional_packages` a vector of symbols with package names, to help AOT compiling
uninstalled, recursive dependencies of `packages`. Look at `compile_incremental` to
use a toml instead.
"""
function compile_package(
        packages::Tuple{String, String}...;
        force = false, reuse = false, debug = false,
        cpu_target = nothing, verbose = false
    )
    userimg = sysimg_folder("precompile.jl")
    if !reuse
        # TODO that's a pretty weak way to check that it's not a path...
        ispackage = all(x-> !occursin(Base.Filesystem.path_separator, x), first.(packages))
        isruntests = all(x-> x == "test/runtests.jl", last.(packages))
        if ispackage && isruntests
            snoop_packages(Symbol.(first.(packages))...; file = userimg)
        else
            ispackage || @warn "Giving path to package deprecated. Use Package name!"
            isruntests || @warn "Giving a snoopfile is deprecated. Use runtests from package!"
        end
    end
    !isfile(userimg) && reuse && error("Nothing to reuse. Please run `compile_package(reuse = true)`")
    image_path = sysimg_folder()
    build_sysimg(image_path, userimg, cpu_target=cpu_target, verbose = verbose)
    imgfile = joinpath(image_path, "sys.$(Libdl.dlext)")
    syspath = joinpath(default_sysimg_path(debug), "sys.$(Libdl.dlext)")
    if force
        try
            backup = syspath * ".packagecompiler_backup"
            isfile(backup) || mv(syspath, backup)
            cp(imgfile, syspath)
            @info """
            Replaced system image successfully. Next start of julia will load the newly compiled system image.
            If you encounter any errors with the new julia image, try `PackageCompiler.revert([debug = false])`.
            """
        catch e
            @warn "An error occured while replacing sysimg files:" error = e
            @info "Recovering old system image from backup"
            # if any file is missing in default system image, revert!
            if !isfile(syspath)
                @info "$syspath missing. Reverting!"
                revert(debug)
            end
        end
    else
        @info """
        Not replacing system image.
        You can start julia with $(`julia -J $imgfile`) at a posix shell to load the compiled files.
        """
    end
    imgfile
end

export compile_package, revert, force_native_image!, executable_ext, build_executable, build_shared_lib, static_julia, compile_incremental

end # module


text:
#!/usr/bin/env julia
using CSV, DataFrames, PyCall, PyPlot
PyDict(matplotlib."rcParams")["axes.formatter.use_locale"] = true
infb_sv_1019 = CSV.File("data/infb_sv_1019.csv") |> DataFrame
infb_sv_1019[:b_prop] = infb_sv_1019[:b]./infb_sv_1019[:samples]
scatter(infb_sv_1019[:b_prop_vic],
    infb_sv_1019[:median_age_b], s = infb_sv_1019[:b_prop]*5000,
    color="lightblue")
grid(1)
for (i, season) in enumerate(infb_sv_1019[:season])
    annotate(season, (infb_sv_1019[:b_prop_vic][i],
        infb_sv_1019[:median_age_b][i]))
end
title("Laboratorierapporterad influensa B Sverige 2010–19")
xlabel("Andel B Victoralinjen")
ylabel("Medianålder fall")
savefig("../../images/infb_sv_1019.svg")

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust