text:
"""
    geneigsolve((A::AbstractMatrix, B::AbstractMatrix), [howmany = 1, which = :LM,
                                    T = promote_type(eltype(A), eltype(B))]; kwargs...)
    geneigsolve(f, n::Int, [howmany = 1, which = :LM, T = Float64]; kwargs...)
    geneigsolve(f, x₀, [howmany = 1, which = :LM]; kwargs...)
    geneigsolve(f, x₀, howmany, which, algorithm)

Compute at least `howmany` generalized eigenvalues ``λ`` and generalized eigenvectors ``x``
of the form ``(A - λB)x = 0``, where `A` and `B` are either instances of `AbstractMatrix`,
or some function that implements the matrix vector product. In case functions are used, one
could either specify the action of `A` and `B` via a tuple of two functions (or a function
and an `AbstractMatrix`), or one could use a single function that takes a single argument
`x` and returns two results, corresponding to `A*x` and `B*x`. Return the computed
eigenvalues, eigenvectors and a `ConvergenceInfo` structure.

### Arguments:
The first argument is either a tuple of two linear maps, so a function or an
`AbstractMatrix` for either of them, representing the action of `A` and `B`. Alternatively,
a single function can be used that takes a single argument `x` and returns the equivalent of
`(A*x, B*x)` as result. This latter form is compatible with the `do` block syntax of Julia.
If an `AbstractMatrix` is used for either `A` or `B`, a starting vector `x₀` does not need
to be provided, it is then chosen as `rand(T, size(A,1))` if `A` is an `AbstractMatrix` (or
similarly if only `B` is an `AbstractMatrix`). Here
`T = promote_type(eltype(A), eltype(B))` if both `A` and `B` are instances of
`AbstractMatrix`, or just the `eltype` of whichever is an `AbstractMatrix`. If both `A` and
`B` are encoded more generally as a callable function or method, the best approach is to
provide an explicit starting guess `x₀`. Note that `x₀` does not need to be of type
`AbstractVector`, any type that behaves as a vector and supports the required methods (see
KrylovKit docs) is accepted. If instead of `x₀` an integer `n` is specified, it is assumed
that `x₀` is a regular vector and it is initialized to `rand(T,n)`, where the default value
of `T` is `Float64`, unless specified differently.

The next arguments are optional, but should typically be specified. `howmany` specifies how
many eigenvalues should be computed; `which` specifies which eigenvalues should be
targeted. Valid specifications of `which` are given by
*   `:LM`: eigenvalues of largest magnitude
*   `:LR`: eigenvalues with largest (most positive) real part
*   `:SR`: eigenvalues with smallest (most negative) real part
*   `:LI`: eigenvalues with largest (most positive) imaginary part, only if `T <: Complex`
*   `:SI`: eigenvalues with smallest (most negative) imaginary part, only if `T <: Complex`
*   [`EigSorter(f; rev = false)`](@ref): eigenvalues `λ` that appear first (or last if
    `rev == true`) when sorted by `f(λ)`
!!! note "Note about selecting `which` eigenvalues"
    Krylov methods work well for extremal eigenvalues, i.e. eigenvalues on the periphery of
    the spectrum of the linear map. Even with `ClosestTo`, no shift and invert is performed.
    This is useful if, e.g., you know the spectrum to be within the unit circle in the
    complex plane, and want to target the eigenvalues closest to the value `λ = 1`.

The argument `T` acts as a hint in which `Number` type the computation should be performed,
but is not restrictive. If the linear map automatically produces complex values, complex
arithmetic will be used even though `T<:Real` was specified.

### Return values:
The return value is always of the form `vals, vecs, info = geneigsolve(...)` with
*   `vals`: a `Vector` containing the eigenvalues, of length at least `howmany`, but could
    be longer if more eigenvalues were converged at the same cost.
*   `vecs`: a `Vector` of corresponding eigenvectors, of the same length as `vals`.
    Note that eigenvectors are not returned as a matrix, as the linear map could act on any
    custom Julia type with vector like behavior, i.e. the elements of the list `vecs` are
    objects that are typically similar to the starting guess `x₀`, up to a possibly
    different `eltype`. When the linear map is a simple `AbstractMatrix`, `vecs` will be
    `Vector{Vector{<:Number}}`.
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: indicates how many eigenvalues and eigenvectors were actually
        converged to the specified tolerance `tol` (see below under keyword arguments)
    -   `info.residual::Vector`: a list of the same length as `vals` containing the
        residuals `info.residual[i] = f(vecs[i]) - vals[i] * vecs[i]`
    -   `info.normres::Vector{<:Real}`: list of the same length as `vals` containing the
        norm of the residual `info.normres[i] = norm(info.residual[i])`
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
    No warning is printed if not all requested eigenvalues were converged, so always check
    if `info.converged >= howmany`.

### Keyword arguments:
Keyword arguments and their default values are given by:
*   `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
*   `tol::Real`: the requested accuracy (corresponding to the 2-norm of the residual for
    Schur vectors, not the eigenvectors). If you work in e.g. single precision (`Float32`),
    you should definitely change the default value.
*   `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
    constructed. Note that the dimension of the vector space is not known or checked, e.g.
    `x₀` should not necessarily support the `Base.length` function. If you know the actual
    problem dimension is smaller than the default value, it is useful to reduce the value
    of `krylovdim`, though in principle this should be detected.
*   `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below
    for further details on the algorithms.
*   `orth::Orthogonalizer`: the orthogonalization method to be used, see
    [`Orthogonalizer`](@ref)
*   `issymmetric::Bool`: if both linear maps `A` and `B` are symmetric, only meaningful if
    `T<:Real`
*   `ishermitian::Bool`: if both linear maps `A` and `B` are hermitian
*   `isposdef::Bool`: if the linear map `B` is positive definite
The default values are given by `tol = KrylovDefaults.tol`, `krylovdim =
KrylovDefaults.krylovdim`, `maxiter = KrylovDefaults.maxiter`, `orth =
KrylovDefaults.orth`; see [`KrylovDefaults`](@ref) for details.

The default value for the last three parameters depends on the method. If an
`AbstractMatrix` is used, `issymmetric`, `ishermitian` and `isposdef` are checked for that
matrix, otherwise the default values are `issymmetric = false` and `ishermitian = T <:
Real && issymmetric`. When values are provided, no checks will be performed even in the
matrix case.

### Algorithm
The last method, without default values and keyword arguments, is the one that is finally
called, and can also be used directly. Here the algorithm is specified, though currently
only [`GolubYe`](@ref) is available. The Golub-Ye algorithm is an algorithm for solving
hermitian (symmetric) generalized eigenvalue problems `A x = λ B x` with positive definite
`B`, without the need for inverting `B`. It builds a Krylov subspace of size `krylovdim`
starting from an estimate `x` by acting with `(A - ρ(x) B)`, where `ρ(x) = dot(x, A*x)/
dot(x, B*x)`, and employing the Lanczos algorithm. This process is repeated at most
`maxiter` times. In every iteration `k>1`, the subspace will also be expanded to size
`krylovdim+1` by adding ``x_k - x_{k-1}``, which is known as the LOPCG correction and was
suggested by Money and Ye. With `krylovdim = 2`, this algorithm becomes equivalent to
`LOPCG`.

!!! warning "Restriction to symmetric definite generalized eigenvalue problems"
    While the only algorithm so far is restricted to symmetric/hermitian generalized
    eigenvalue problems with positive definite `B`, this is not reflected in the default
    values for the keyword arguments `issymmetric` or `ishermitian` and `isposdef`. Make
    sure to set these to true to understand the implications of using this algorithm.
"""
function geneigsolve end

function geneigsolve(AB::Tuple{AbstractMatrix,AbstractMatrix}, howmany::Int = 1,
        which::Selector = :LM, T = promote_type(eltype.(AB)...); kwargs...)
    if !(size(AB[1],1) == size(AB[1],2) == size(AB[2],1) == size(AB[2],2))
        throw(DimensionMismatch("Matrices `A` and `B` should be square and have matching size"))
    end
    geneigsolve(AB, rand(T, size(AB[1],1)), howmany::Int, which; kwargs...)
end
geneigsolve(AB::Tuple{Any,AbstractMatrix}, howmany::Int = 1, which::Selector = :LM,
            T = eltype(AB[2]); kwargs...) =
    geneigsolve(AB, rand(T, size(AB[2],1)), howmany, which; kwargs...)
geneigsolve(AB::Tuple{AbstractMatrix,Any}, howmany::Int = 1, which::Selector = :LM,
            T = eltype(AB[1]); kwargs...) =
    geneigsolve(AB, rand(T, size(AB[1],1)), howmany, which; kwargs...)

geneigsolve(f, n::Int, howmany::Int = 1, which::Selector = :LM, T = Float64; kwargs...) =
    geneigsolve(f, rand(T, n), howmany, which; kwargs...)

function geneigsolve(f, x₀, howmany::Int = 1, which::Selector = :LM; kwargs...)
    Tx = typeof(x₀)
    Tfx = Core.Compiler.return_type(geneigfun(f), Tuple{Tx}) # should be a tuple type
    Tfx1 = Base.tuple_type_head(Tfx)
    Tfx2 = Base.tuple_type_head(Base.tuple_type_tail(Tfx))
    T1 = Core.Compiler.return_type(dot, Tuple{Tx, Tfx1})
    T2 = Core.Compiler.return_type(dot, Tuple{Tx, Tfx2})
    T = promote_type(T1, T2)
    alg = geneigselector(f, T; kwargs...)
    if alg isa GolubYe && (which == :LI || which == :SI)
        error("Eigenvalue selector which = $which invalid: real eigenvalues expected with Lanczos algorithm")
    end
    geneigsolve(f, x₀, howmany, which, alg)
end

function geneigselector(AB::Tuple{AbstractMatrix,AbstractMatrix}, T::Type;
                        issymmetric = T <: Real && all(LinearAlgebra.issymmetric, AB),
                        ishermitian = issymmetric || all(LinearAlgebra.ishermitian, AB),
                        isposdef = ishermitian && LinearAlgebra.isposdef(AB[2]), kwargs...)
    if (issymmetric || ishermitian) && isposdef
        return GolubYe(; kwargs...)
    else
        throw(ArgumentError("Only symmetric or hermitian generalized eigenvalue problems with positive definite `B` matrix are currently supported."))
    end
end
function geneigselector(f, T::Type; issymmetric = false, ishermitian = issymmetric,
                            isposdef = false, kwargs...)
    if (issymmetric || ishermitian) && isposdef
        return GolubYe(; kwargs...)
    else
        throw(ArgumentError("Only symmetric or hermitian generalized eigenvalue problems with positive definite `B` matrix are currently supported."))
    end
end


text:
""" Conventions for serialization of wiring diagrams.

Defines a consistent set of names for boxes, ports, and wires to be used when
serializing wiring diagrams, as well as conventions for serializing box, port,
and wire attributes.
"""
module WiringDiagramSerialization
export box_id, wire_id, port_name,
  convert_from_graph_data, convert_to_graph_data

using ..DirectedWiringDiagrams

# Identifiers
#############

function box_id(path::Vector{Int})::String
  if isempty(path)
    "root"
  else
    prod(["n$v" for v in path])
  end
end
function box_id(diagram::WiringDiagram, path::Vector{Int})
  v = last(path)
  if v == input_id(diagram) || v == output_id(diagram)
    box_id(path[1:end-1])
  else
    box_id(path)
  end
end

function wire_id(path::Vector{Int}, wire::Int)::String
  if isempty(path)
    string("e", wire)
  else
    string(box_id(path), "e", wire)
  end
end

function port_name(kind::PortKind, port::Int)
  string(kind == InputPort ? "in" : "out", port)
end
function port_name(diagram::WiringDiagram, port::Port)
  kind = if port.box == input_id(diagram)
    InputPort
  elseif port.box == output_id(diagram)
    OutputPort
  else
    port.kind
  end
  port_name(kind, port.port)
end

# Attributes
############

convert_from_graph_data(::Type{Nothing}, data::AbstractDict) = nothing
convert_from_graph_data(::Type{T}, data::AbstractDict) where T <: AbstractDict =
  convert(T, data)

function convert_from_graph_data(Value::Type, data::AbstractDict)
  @assert length(data) == 1
  convert(Value, first(values(data)))
end
function convert_from_graph_data(::Type{Symbol}, data::AbstractDict)
  @assert length(data) == 1
  Symbol(first(values(data)))
end

convert_to_graph_data(::Nothing) = Dict()
convert_to_graph_data(value::AbstractDict) = value
convert_to_graph_data(value) = Dict("value" => value)

end

ext: jl


text:
@testset "Median, Mean, and Variance (on Sphere(2))" begin
    M = Sphere(2)
    x = SnPoint([1., 0, 0.])
    ξ = SnTVector([0., 1., 1.])
    η = SnTVector([0., 1., -1.])
    y = exp(M,x,0.2*ξ)
    z = exp(M,x,0.2*η)
    # CPPA not that exact
    @test distance(M, mean(M, [x,y]; method = :GradientDescent), midPoint(M,x,y) ) ≈ 0
    @test distance(M, mean(M, [x,y]; method = :CyclicProximalPoint), midPoint(M,x,y) ) ≈ 0 atol=8*10^(-5)
    @test distance(M,
        mean(M, [x,y,z]; method = :GradientDescent),
        mean(M, [x,y,z]; method = :CyclicProximalPoint)
    ) ≈ 0 atol=10^(-4)
    # Variance
    m = mean(M, [x,y,z])
    v = 1/(2 * manifoldDimension(M)) * sum( distance.(Ref(M),Ref(m),[x,y,z]).^2 )
    @test variance(M,[x,y,z]) ≈ v atol=10^-8
    # Median
    p = SnPoint([0.960266, 0.197344, 0.197344])
    @test distance(M,median(M,[x,y]),p) ≈ 0
end


text:
# --- All libIIO functions wrapped --- #

module AdalmPluto

using Reexport
using Printf

include("libIIO/libIIO.jl");
@reexport using .libIIO_jl;

# --- Functions taken from libad9361 --- #
include("libad9361.jl");
# --- utils --- #
include("utils.jl");

# --- struct exports --- #
export
    ChannelCfg,
    PlutoSDR
;

# --- functions exports --- #
export
    openPluto,
    close,
    recv,
    recv!,
    send,
    refillJuliaBufferRX,
    updateGain!,
    updateGainMode!,
    updateCarrierFreq!,
    updateSamplingRate!,
    updateBandwidth!,
    getGain,
    getBackends
;

# Preferences importation to change libIIO origin
# We do not export set_provider and get_provider
# To change the lib origin, type
# AdalmPluto.set_provider("local")
set_provider = libIIO_jl.set_provider
get_provider = libIIO_jl.get_provider

# constants
const TX_DEVICE_NAME  = "cf-ad9361-dds-core-lpc";
const RX_DEVICE_NAME  = "cf-ad9361-lpc";
const PHY_DEVICE_NAME = "ad9361-phy";

@enum GainMode begin
    MANUAL
    FAST_ATTACK
    SLOW_ATTACK
    HYBRID
    DEFAULT
end

# ------------------ #
# --- Structures --- #
# ------------------ #

# --- TRx structure --- #
mutable struct ChannelCfg
    rfport::String;
    bandwidth::Int64;
    samplingRate::Int64;
    carrierFreq::Int64;
end

mutable struct IIO_Buffer
    C_ptr::Ptr{iio_buffer};
    C_size::Csize_t;
    C_sample_size::Cssize_t;
    C_first::Ptr{Cuchar};
    C_last::Ptr{Cuchar};
    C_step::Cssize_t;
    i_raw_samples::Array{UInt8};
    q_raw_samples::Array{UInt8};
    samples::Array{ComplexF32};
    nb_samples::Int;
end

# --- Rx structures --- #
struct rxWrapper
    rx::Ptr{iio_device};
    rx0_i::Ptr{iio_channel};
    rx0_q::Ptr{iio_channel};
    chn::Ptr{iio_channel};
    chn_lo::Ptr{iio_channel};
end

mutable struct PlutoRx
    iio::rxWrapper;
    cfg::ChannelCfg;
    buf::IIO_Buffer;
    effectiveSamplingRate::Float64;
    effectiveCarrierFreq::Float64;
    released::Bool;
end

# --- Tx structures --- #
struct txWrapper
    tx::Ptr{iio_device};
    tx0_i::Ptr{iio_channel};
    tx0_q::Ptr{iio_channel};
    chn::Ptr{iio_channel};
    chn_lo::Ptr{iio_channel};
end

mutable struct PlutoTx
    iio::txWrapper;
    cfg::ChannelCfg;
    buf::IIO_Buffer;
    effectiveSamplingRate::Float64;
    effectiveCarrierFreq::Float64;
    released::Bool;
end

# --- final wrap --- #
"""
    PlutoSDR
    +-- ctx::Ptr{iio_context}
    |
    +-- tx::PlutoTx
    |   |
    |   +-- iio::txWrapper
    |   |   |
    |   |   +-- tx::Ptr{iio_device};
    |   |   +-- tx0_i::Ptr{iio_channel};
    |   |   +-- tx0_q::Ptr{iio_channel};
    |   |   +-- chn::Ptr{iio_channel};
    |   |   +-- chn_lo::Ptr{iio_channel};
    |   |
    |   +-- cfg::ChannelCfg
    |   |   |
    |   |   +-- rfport::String;
    |   |   +-- bandwidth::Int64;
    |   |   +-- samplingRate::Int64;
    |   |   +-- carrierFreq::Int64;
    |   |
    |   +-- buf::IIO_Buffer
    |   |   |
    |   |   +-- C_ptr::Ptr{iio_buffer};         |< Pointer to the C buffer
    |   |   +-- C_size::Csize_t;                |< Size of the C buffer in samples
    |   |   +-- C_sample_size::Cssize_t;        |< Size of a sample in the C buffer
    |   |   +-- C_first::Ptr{Cuchar};           |< Pointer to the first sample in the C buffer
    |   |   +-- C_last::Ptr{Cuchar};            |< Pointer to the last sample in the C buffer
    |   |   +-- C_step::Cssize_t;               |< Distance between to sample pointers in the C buffer
    |   |   +-- i_raw_samples::Array{UInt8}     |< Temporary buffer to store I bytes
    |   |   +-- q_raw_samples::Array{UInt8}     |< Temporary buffer to store Q bytes
    |   |   +-- samples::Array{ComplexF32}      |< Samples stored in a Julia array
    |   |   +-- nb_samples::Int                 |< Number of samples in the Julia array
    |   |
    |   +-- effectiveSamplingRate::Float64
    |   +-- effectiveCarrierFreq::Float64
    |   +-- released::Bool
    |
    +-- rx::PlutoRx
    |   |
    |   +-- iio::rxWrapper
    |   |   +-- rx::Ptr{iio_device};
    |   |   +-- rx0_i::Ptr{iio_channel};
    |   |   +-- rx0_q::Ptr{iio_channel};
    |   |   +-- chn::Ptr{iio_channel};
    |   |   +-- chn_lo::Ptr{iio_channel};
    |   |
    |   +-- cfg::ChannelCfg
    |   |   |
    |   |   +-- rfport::String;
    |   |   +-- bandwidth::Int64;
    |   |   +-- samplingRate::Int64;
    |   |   +-- carrierFreq::Int64;
    |   |
    |   +-- buf::IIO_Buffer
    |   |   |
    |   |   +-- C_ptr::Ptr{iio_buffer};         |< Pointer to the C buffer
    |   |   +-- C_size::Csize_t;                |< Size of the C buffer in samples
    |   |   +-- C_sample_size::Cssize_t;        |< Size of a sample in the C buffer
    |   |   +-- C_first::Ptr{Cuchar};           |< Pointer to the first sample in the C buffer
    |   |   +-- C_last::Ptr{Cuchar};            |< Pointer to the last sample in the C buffer
    |   |   +-- C_step::Cssize_t;               |< Distance between to sample pointers in the C buffer
    |   |   +-- i_raw_samples::Array{UInt8}     |< Temporary buffer to store I bytes
    |   |   +-- q_raw_samples::Array{UInt8}     |< Temporary buffer to store Q bytes
    |   |   +-- samples::Array{ComplexF32}      |< Samples decoded and stored in a Julia array
    |   |   +-- nb_samples::Int                 |< Number of samples available in the Julia array
    |   |
    |   +-- effectiveSamplingRate::Float64
    |   +-- effectiveCarrierFreq::Float64
    |   +-- released::Bool
    |
    +-- released::Bool
"""
mutable struct PlutoSDR
    ctx::Ptr{iio_context};
    tx::PlutoTx;
    rx::PlutoRx;
    released::Bool;
end

# ------------------------ #
# --- Helper functions --- #
# ------------------------ #
"""
    scan(backend[, infoIndex, doPrint])

Returns a device URI.

# Arguments
- `backend::String` : the backend to scan (local, xml, ip, usb).
- `infoIndex::Integer=1` : the device index.
- `doPrint::Bool=true` : toggles console printing of the uri.

# Returns
- `uri::String` : the device URI.

[C equivalent](https://analogdevicesinc.github.io/libiio/master/libiio/iio-monitor_8c-example.html#_a15)
"""
function scan(backend::String, deviceIndex=1, doPrint=true)
    # Check if backend is available and create scan context
    C_iio_has_backend(backend) || error("Specified backend $backend is not available");
    scan_context = C_iio_create_scan_context(backend);

    # Getting metadata from scan
    info = Ref{Ptr{Ptr{iio_context_info}}}(0);
    ret = C_iio_scan_context_get_info_list(scan_context, info);

    # Get usb address
    uri = "";
    if ret < 0
        C_iio_context_info_list_free(info[]);
        C_iio_scan_context_destroy(scan_context);
        error("iio_scan_context_get_info_list failed with error $ret :\n", C_iio_strerror(ret));
    elseif ret == 0
        (doPrint) && (@info "No $backend device found");
    else
        loaded_info = unsafe_load(info[], deviceIndex);
        description = C_iio_context_info_get_description(loaded_info);
        uri = C_iio_context_info_get_uri(loaded_info);
        (doPrint) && (@info "Found $ret device(s) with $backend backend.\nSelected $description [$uri]");
    end

    C_iio_context_info_list_free(info[]);
    C_iio_scan_context_destroy(scan_context);

    return uri;
end

"""
    createContext(uri)

Returns the `iio_context` corresponding the the provided uri. Throws an error if no devices are found in the context.

# Arguments
- `uri::String` : the radio URI to get the context from.

# Returns
- `context::Ptr{iio_context}` : the context for the given URI.
"""
function createContext(uri::String)
    context = C_iio_create_context_from_uri(uri);
    if(C_iio_context_get_devices_count(context) == 0) error("No device found in context from uri $uri"); end
    return context;
end

"""
    findTRXDevices(context)

Returns `txDevice::Ptr{iio_device}, rxDevice::Ptr{iio_device}` from the given context.

# Arguments
- `context::Ptr{iio_context}` : the IIO context to get the TX and RX devices from.

# Returns
- `txDevice::Ptr{iio_device}` : a pointer to the TX C iio_device structure.
- `rxDevice::Ptr{iio_device}` : a pointer to the RX C iio_device structure.
"""
function findTRXDevices(context::Ptr{iio_context})
    txDevice = C_iio_context_find_device(context, TX_DEVICE_NAME);
    rxDevice = C_iio_context_find_device(context, RX_DEVICE_NAME);
    return txDevice, rxDevice;
end

"""
    findTRXChannels(context[, txID, rxID])

Returns `txChannel::Ptr{iio_channel}, rxChannel::Ptr{iio_channel}` from the given context.

# Arguments
- `context::Ptr{iio_context}` : the context to get the channels from.
- `txID::Integer=0` : the tx channel number (ex : 0 for tx channel "voltage0").
- `rxID::Integer=0` : the rx channel number (ex : 0 for rx channel "voltage0").

# Returns
- `txChannel::Ptr{iio_channel}` : a pointer to the TX C iio_channel structure.
- `rxChannel::Ptr{iio_channel}` : a pointer to the RX C iio_channel structure.
"""
function findTRXChannels(context::Ptr{iio_context}, txID=0, rxID=0)
    device = C_iio_context_find_device(context, PHY_DEVICE_NAME);
    txChannel = C_iio_device_find_channel(device, "voltage" * string(txID), true);
    rxChannel = C_iio_device_find_channel(device, "voltage" * string(rxID), false);

    return txChannel, rxChannel;
end

"""
    findLoChannels(context[, txLoID, rxLoID])

Returns `txLoChannel::Ptr{iio_channel}, rxLoChannel::Ptr{iio_channel}` from the given context.

# Arguments
- `context::Ptr{iio_context}` : the context to get the channels from.
- `txLoID::Integer=1` : the tx lo channel number (ex : 1 for tx lo channel "altvoltage1").
- `rxLoID::Integer=0` : the rx lo channel number (ex : 0 for rx lo channel "altvoltage0").

# Returns
- `txLoChannel::Ptr{iio_channel}` : a pointer to the TX Lo C iio_channel structure.
- `rxLoChannel::Ptr{iio_channel}` : a pointer to the RX Lo C iio_channel structure.
"""
function findLoChannels(context::Ptr{iio_context}, txLoID=1, rxLoID=0)
    device = C_iio_context_find_device(context, PHY_DEVICE_NAME);
    txLoChannel = C_iio_device_find_channel(device, "altvoltage" * string(txLoID), true);
    rxLoChannel = C_iio_device_find_channel(device, "altvoltage" * string(rxLoID), true);

    return txLoChannel, rxLoChannel;
end

"""
    cfgChannels(context, txCfg, rxCfg)

Configures the RX, RX LO, TX, and TX LO channels from the provided configurations.
Returns the configured channels.

# Arguments
- `context::Ptr{iio_context}` : the context to get the channels from.
- `txCfg::ChannelCfg` : the tx configuration (port / bandwidth / sampling rate / frequency).
- `rxCfg::ChannelCfg` : the rx configuration (port / bandwidth / sampling rate / frequency).

# Returns
- `txChannel::Ptr{iio_channel}` : a pointer to the configured TX C iio_channel structure.
- `rxChannel::Ptr{iio_channel}` : a pointer to the configured RX C iio_channel structure.
- `txLoChannel::Ptr{iio_channel}` : a pointer to the configured TX Lo C iio_channel structure.
- `rxLoChannel::Ptr{iio_channel}` : a pointer to the configured RX Lo C iio_channel structure.
"""
function cfgChannels(context::Ptr{iio_context}, txCfg::ChannelCfg, rxCfg::ChannelCfg)
    txChannel, rxChannel = findTRXChannels(context);
    txLoChannel, rxLoChannel = findLoChannels(context);

    return cfgChannels(txChannel, rxChannel, txLoChannel, rxLoChannel, txCfg, rxCfg);
end

"""
    cfgChannels(txChannel, rxChannel, txLoChannel, rxLoChannel, txCfg, rxCfg)

Configures the RX, RX LO, TX, and TX LO channels from the provided configurations.
Returns the configured channels.

# Arguments
- `txChannel::Ptr{iio_channel}` : the tx channel.
- `rxChannel::Ptr{iio_channel}` : the rx channel.
- `txLoChannel::Ptr{iio_channel}` : the tx lo channel.
- `rxLoChannel::Ptr{iio_channel}` : the rx lo channel.
- `txCfg::ChannelCfg` : the tx configuration (port / bandwidth / sampling rate / frequency).
- `rxCfg::ChannelCfg` : the rx configuration (port / bandwidth / sampling rate / frequency).

# Returns
- `txChannel::Ptr{iio_channel}` : a pointer to the configured TX C iio_channel structure.
- `rxChannel::Ptr{iio_channel}` : a pointer to the configured RX C iio_channel structure.
- `txLoChannel::Ptr{iio_channel}` : a pointer to the configured TX Lo C iio_channel structure.
- `rxLoChannel::Ptr{iio_channel}` : a pointer to the configured RX Lo C iio_channel structure.
"""
function cfgChannels(
    txChannel::Ptr{iio_channel},
    rxChannel::Ptr{iio_channel},
    txLoChannel::Ptr{iio_channel},
    rxLoChannel::Ptr{iio_channel},
    txCfg::ChannelCfg,
    rxCfg::ChannelCfg
)
    C_iio_channel_attr_write(txChannel, "rf_port_select", txCfg.rfport);
    C_iio_channel_attr_write_longlong(txChannel, "rf_bandwidth", txCfg.bandwidth);
    C_iio_channel_attr_write_longlong(txChannel, "sampling_frequency", txCfg.samplingRate);
    C_iio_channel_attr_write_longlong(txLoChannel, "frequency", txCfg.carrierFreq);

    C_iio_channel_attr_write(rxChannel, "rf_port_select", rxCfg.rfport);
    C_iio_channel_attr_write_longlong(rxChannel, "rf_bandwidth", rxCfg.bandwidth);
    C_iio_channel_attr_write_longlong(rxChannel, "sampling_frequency", rxCfg.samplingRate);
    C_iio_channel_attr_write_longlong(rxLoChannel, "frequency", rxCfg.carrierFreq);

    return txChannel, rxChannel, txLoChannel, rxLoChannel;
end

"""
    findIQChannels(device, iID, qID, isOutput)

Returns `IChannel::Ptr{iio_channel}, QChannel::Ptr{iio_channel}` from the given device.

# Arguments
- `device::Ptr{iio_device}` : the device to get the channels from.
- `iID::String` : identification string for the I channel (ex : "voltage0").
- `qID::String` : identification string for the Q channel (ex : "voltage1").
- `isOutput::Bool` : whether the IQ channels are outputs.

# Returns
- `IChannel::Ptr{iio_channel}` : the I channel.
- `QChannel::Ptr{iio_channel}` : the Q channel.
"""
function findIQChannels(device::Ptr{iio_device}, iID::String, qID::String, isOutput::Bool)
    IChannel = C_iio_device_find_channel(device, iID, isOutput);
    QChannel = C_iio_device_find_channel(device, qID, isOutput);
    return IChannel, QChannel;
end

"""
    updateEffectiveCfg!(trx)

Updates the stored values of the effective sampling rate and carrier frequency of either a `PlutoTx` or `PlutoRx`.
If the values are different than the ones in the `ChannelCfg` of the structure, a warning is printed.
Returns the effective values values.

# Arguments
- `trx::Union{PlutoTx, PlutoRx}` : the structure containing the channel to read the current configuration from.

# Keywords
- `doLog::Bool` : toggles the display of the new carrier frequency

# Returns
- `effectiveSamplingRate::Int` : the current sampling rate.
- `effectiveCarrierFreq::Int` : the current carrier frequency.
"""
function updateEffectiveCfg!(trx::Union{PlutoTx, PlutoRx}; doLog=true)
    if typeof(trx) == PlutoRx; global type = :RX; else; global type = :TX; end;

    ret, effectiveSamplingRate = C_iio_channel_attr_read_longlong(trx.iio.chn, "sampling_frequency");
    if ret < 0
        @warnPluto type "Could not get effective sampling rate (Error $ret):\n" * C_iio_strerror(ret);
    end
    ret, effectiveCarrierFreq = C_iio_channel_attr_read_longlong(trx.iio.chn_lo, "frequency");
    if ret < 0
        @warnPluto type "Could not get effective carrier frequency (Error $ret):\n" * C_iio_strerror(ret);
    end

    trx.effectiveSamplingRate = effectiveSamplingRate;
    trx.effectiveCarrierFreq  = effectiveCarrierFreq;

    if effectiveSamplingRate != trx.cfg.samplingRate && doLog
        @warnPluto type "Effective sampling rate ($effectiveSamplingRate) ≠ Requested sampling rate ($(trx.cfg.samplingRate))";
    end
    if effectiveCarrierFreq != trx.cfg.carrierFreq && doLog
        @warnPluto type "Effective carrier frequency ($effectiveCarrierFreq) ≠ Requested carrier frequency ($(trx.cfg.carrierFreq))";
    end

    return effectiveSamplingRate, effectiveCarrierFreq;
end

"""
    createBuffer(device, samplesCount)

Creates a buffer for the provided device. Returns a wrapper around the buffer with basic info
(C pointer, sample size, first sample, last sample, steps between samples, Julia complex samples, number of non-queried Julia complex samples).

# Arguments
- `device::Ptr{iio_device}` : the device for which the buffer is created.
- `samplesCound::UInt` : the size of the buffer in samples.

# Returns
- `buffer::IIO_Buffer` : a buffer for the given channel with space for samplesCount samples.
"""
function createBuffer(device::Ptr{iio_device}, samplesCount::UInt)
    sampleSize = C_iio_device_get_sample_size(device);
    buf        = C_iio_device_create_buffer(device, samplesCount, false);
    first      = C_iio_buffer_start(buf);
    last       = C_iio_buffer_end(buf);
    step       = C_iio_buffer_step(buf);

    return IIO_Buffer(
        buf, samplesCount, sampleSize, first, last, step,   # values regarding the C buffers
        zeros(UInt8, (samplesCount * sampleSize) ÷ 2),      # array to store raw I samples
        zeros(UInt8, (samplesCount * sampleSize) ÷ 2),      # array to store raw Q samples
        zeros(ComplexF32, samplesCount),                    # array to store decoded complex samples
        0);                                                 # number of samples not yet queried
end

# ------------------------ #
# --- Module functions --- #
# ------------------------ #

"""
    openPluto(carrierFreq, samplingRate, bandwidth[, uri, backend])

Creates a PlutoSDR struct and configures the radio to stream the samples.

# Arguments
- `carrierFreq::Int` : the carrier frequency for both tx and rx.
- `samplingRate::Int` : the sampling rate for both tx and rx.
- `gain::Int` : the analog RX gain.

# Keywords
- `addr::String="auto"` : the radio address (ex: "usb:1.3.5"). "auto" takes the first uri found for the given backend.
- `backend::Union{nothing,String}` : the backend to scan for the auto uri. If not specified, all backend are scanned and the first radio found is used.
- `bufferSize::UInt=1024*1024` : the buffer size in samples.
- `bandwidth::Int` : the bandwidth for both tx and rx.

# Returns
- `radio::PlutoSDR` : a fully initialized PlutoSDR structure.
"""
function openPluto(
    carrierFreq::Int, samplingRate::Int, gain::Int, antenna="A;A_BALANCED";
    addr::String="auto", backend=nothing, bufferSize::UInt=UInt64(1024*1024), bandwidth::Int=Int(20e6)
)
    antenna = split(antenna, ";");
    if !(length(antenna) == 2)
        throw(ArgumentError("Couldn't parse antenna ports. Expected syntax: TXport;RXport"));
    end
    radio = openPluto(
        ChannelCfg(antenna[1], bandwidth, samplingRate, carrierFreq),
        ChannelCfg(antenna[2], bandwidth, samplingRate, carrierFreq),
        bufferSize,
        addr,
        backend
    );
    updateGain!(radio, gain);
    return radio;
end

"""
    openPluto(txCfg, rxCfg[, uri, backend])

Creates a PlutoSDR struct and configures the radio to stream the samples.

# Arguments
- `txCfg::ChannelCfg` : the port / bandwidth / sampling rate / carrier frequency for the tx channels.
- `rxCfg::ChannelCfg` : the port / bandwidth / sampling rate / carrier frequency for the rx channels.
- `bufferSize::UInt=1024*1024` : the buffer size in samples.
- `uri::String="auto"` : the radio uri (ex : "usb:1.3.5"). "auto" takes the first uri found for the given backend.
- `backend::Union{nothing,String}` : the backend to scan for the auto uri. If not specified, all backend are scanned and the first radio found is used.
# Returns
- `radio::PlutoSDR` : a fully initialized PlutoSDR structure.
"""
function openPluto(txCfg::ChannelCfg, rxCfg::ChannelCfg, bufferSize::UInt=UInt64(1024*1024), uri="auto", backend=nothing)
    # --- Address is not specified, we try to auto-discover the device
    if uri == "auto"
        # We do not specify the backend, let's try all
        isnothing(backend) ?  allBackend= getBackends() : allBackend = [backend]
        for b ∈ allBackend
            curr_uri = scan(b);
            if !isempty(curr_uri)
                # Found a device, stop
                uri = curr_uri
                break
            end
        end
    end
    if uri == "auto"
        # If uri is init to auto, we scan and we update uri if we found something. If uri is still auto, we do not have found any radio connected, so try with default address
        @warn "Unable to auto-detect the Pluto device, try with default address ip:192.168.2.1. To force a given different address use addr=\"192.168.X.X\" in openPluto"
        uri = "ip:192.168.2.1"
    end
    context = createContext(uri);
    # printing stuff
    description = C_iio_context_get_description(context);
    println("Description :\n$description\n");

    tx, rx = findTRXDevices(context);

    tx0_i, tx0_q = findIQChannels(tx, "voltage0", "voltage1", true);
    rx0_i, rx0_q = findIQChannels(rx, "voltage0", "voltage1", false);

    C_iio_channel_enable(tx0_i);
    C_iio_channel_enable(tx0_q);
    C_iio_channel_enable(rx0_i);
    C_iio_channel_enable(rx0_q);

    txChannel, rxChannel, txLoChannel, rxLoChannel = cfgChannels(context, txCfg, rxCfg);

    iioTx = txWrapper(
        tx,
        tx0_i,
        tx0_q,
        txChannel,
        txLoChannel
    );
    iioRx = rxWrapper(
        rx,
        rx0_i,
        rx0_q,
        rxChannel,
        rxLoChannel
    );

    # 1 MiS buffers
    txBuffer = createBuffer(tx, bufferSize);
    rxBuffer = createBuffer(rx, bufferSize);

    tx = PlutoTx(
        iioTx,
        txCfg,
        txBuffer,
        -1,
        -1,
        false
    );
    rx = PlutoRx(
        iioRx,
        rxCfg,
        rxBuffer,
        -1,
        -1,
        false
    );

    # effective cfg
    updateEffectiveCfg!(tx);
    updateEffectiveCfg!(rx);

    pluto = PlutoSDR(
        context,
        tx,
        rx,
        false
    );

    return pluto;
end

"""
    close(pluto)

Frees the C allocated memory associated to the PlutoSDR structure.
"""
function Base.close(pluto::PlutoSDR)
    if pluto.released
        @warn "Pluto has already been released, abort call";
    else
        C_iio_channel_disable(pluto.rx.iio.rx0_i);
        C_iio_channel_disable(pluto.rx.iio.rx0_q);
        C_iio_channel_disable(pluto.tx.iio.tx0_i);
        C_iio_channel_disable(pluto.tx.iio.tx0_q);
        C_iio_buffer_destroy(pluto.rx.buf.C_ptr);
        C_iio_buffer_destroy(pluto.tx.buf.C_ptr);
        C_iio_context_destroy(pluto.ctx);
        pluto.tx.released = true;
        pluto.rx.released = true;
        pluto.released = true;
    end
end

"""
    updateGainMode!(pluto[, mode])

Modifies the pluto RX channel gain control mode.
Returns an error code < 0 if it doesn't succeed.

# Arguments
- `pluto::PlutoSDR` : the radio to modify.
- `mode::GainMode=DEFAULT` : the new gain mode. DEFAULT ≡ FAST_ATTACK.

# Returns
- `errno::Int` : 0 or a negative error code.
"""
function updateGainMode!(pluto::PlutoSDR, mode::GainMode=DEFAULT)
    control_mode = "";
    if mode == MANUAL
        control_mode = "manual"
    elseif mode == FAST_ATTACK
        control_mode = "fast_attack"
    elseif mode == SLOW_ATTACK
        control_mode = "slow_attack"
    elseif mode == HYBRID
        control_mode = "hybrid"
    else
        control_mode = "fast_attack"
    end

    return C_iio_channel_attr_write(pluto.rx.iio.chn, "gain_control_mode", control_mode);
end

"""
    updateGain!(pluto, value)

Changes the gain control mode to manual et sets the given value.
Prints a warning and returns the error code if it doesn't succeed.

# Arguments
- `pluto::PlutoSDR` : the radio to modify.
- `value::Int64` : the manual gain value.

# Returns
- `errno::Int` : 0 or a negative error code.
"""
function updateGain!(pluto::PlutoSDR, value::Int64)
    ret = updateGainMode!(pluto, MANUAL);
    if ret < 0
        @warnPluto :RX "Could not set gain_control_mode to manual (Error $ret):\n" * C_iio_strerror(ret);
        return ret;
    end
    ret = C_iio_channel_attr_write_longlong(pluto.rx.iio.chn, "hardwaregain", value);
    if ret < 0
        @warnPluto :RX "Could not set hardwaregain to $value (Error $ret):\n" * C_iio_strerror(Int64(ret));
    end
    return ret;
end

"""
    getGain(pluto)

Returns the gain value
# Arguments
- `pluto::PlutoSDR` : the radio to monitor

# Returns
- `gain` :Current radio gain
"""
function getGain(pluto::PlutoSDR)
    value = Ref{Csize_t}(0)
    ret,gain   = C_iio_channel_attr_read_longlong(pluto.rx.iio.chn, "hardwaregain")
    if ret < 0
        @warnPluto :RX "Could not read hardwaregain"
    end
    return gain;
end

"""
    updateCarrierFreq!(pluto, value; doLog)

Changes the carrier frequency. Prints the new effective frequency.

# Arguments
- `pluto::PlutoSDR` : the radio to modify.
- `value::Int64` : the new carrier frequency.

# Keywords
- `doLog::Bool` : toggles the display of the new carrier frequency

# Returns
- `errno::Int` : 0 or a negative error code.
"""
function updateCarrierFreq!(pluto::PlutoSDR, value::Int64; doLog=true)
    # set the carrier frequencie for RX
    errno = C_iio_channel_attr_write_longlong(pluto.rx.iio.chn_lo, "frequency", value);
    if (errno < 0); return errno; end;

    # store the requested configuration and effective configuration for RX
    pluto.rx.cfg.carrierFreq = value;
    effectiveCarrierFreq = updateEffectiveCfg!(pluto.rx; doLog=doLog)[2];
    if doLog; @infoPluto :RX "New RX carrier frequency : $effectiveCarrierFreq ($value)"; end;

    # set the carrier frequencie for TX
    errno = C_iio_channel_attr_write_longlong(pluto.tx.iio.chn_lo, "frequency", value);
    if (errno < 0); return errno; end;

    # store the requested configuration and effective configuration for TX
    pluto.tx.cfg.carrierFreq = value;
    effectiveCarrierFreq = updateEffectiveCfg!(pluto.tx; doLog=doLog)[2];
    if doLog; @infoPluto :TX "New TX carrier frequency : $effectiveCarrierFreq ($value)"; end;

    # not actually errno because it's equal to the number of bytes written
    return 0;
end

"""
    updateSamplingRate!(pluto, value; doLog)

Changes the sampling rate. Prints the new effective sampling rate.

# Arguments
- `pluto::PlutoSDR` : the radio to modify.
- `value::Int64` : the new sampling rate.

# Keywords
- `doLog::Bool` : toggles the display of the new carrier frequency

# Returns
- `errno::Int` : 0 or a negative error code.
"""
function updateSamplingRate!(pluto::PlutoSDR, value::Int64; doLog=true)
    # set the sampling rate for RX
    errno = C_iio_channel_attr_write_longlong(pluto.rx.iio.chn, "sampling_frequency", value);
    if (errno < 0); return errno; end;

    # store the requested configuration and effective configuration for RX
    pluto.rx.cfg.samplingRate = value;
    effectiveSamplingRate = updateEffectiveCfg!(pluto.rx; doLog=doLog)[1];
    if doLog; @infoPluto :RX "New RX sampling rate : $effectiveSamplingRate ($value)"; end;

    # set the sampling rate for TX
    errno = C_iio_channel_attr_write_longlong(pluto.tx.iio.chn, "sampling_frequency", value);
    if (errno < 0); return errno; end;

    # store the requested configuration and effective configuration for TX
    pluto.tx.cfg.samplingRate = value;
    effectiveSamplingRate = updateEffectiveCfg!(pluto.tx; doLog=doLog)[1];
    if doLog; @infoPluto :TX "New TX sampling rate : $effectiveSamplingRate ($value)"; end;

    # not actually errno because it's equal to the number of bytes written
    return 0;
end

"""
    updateBandwidth!(pluto, value)

Changes the bandwidth. Prints the new value.

# Arguments
- `pluto::PlutoSDR` : the radio to modify.
- `value::Int64` : the new sampling rate.

# Keywords
- `doLog::Bool` : toggles the display of the new carrier frequency

# Returns
- `errno::Int` : 0 or a negative error code.
"""
function updateBandwidth!(pluto, value::Int64; doLog=true)
    # set the bandwidth for RX
    errno = C_iio_channel_attr_write_longlong(pluto.rx.iio.chn, "rf_bandwidth", value);
    if (errno < 0); return errno; end;

    # store the new configuration for RX
    if doLog; @infoPluto :RX "New Rx bandwidth : $value"; end;
    pluto.rx.cfg.bandwidth = value;

    # set the bandwidth for TX
    errno = C_iio_channel_attr_write_longlong(pluto.tx.iio.chn, "rf_bandwidth", value);
    if (errno < 0); return errno; end;

    # store the new configuration for TX
    if doLog; @infoPluto :TX "New Tx bandwidth : $value"; end;
    pluto.tx.cfg.bandwidth = value;

    return 0;
end

"""
    recv(pluto, nbSamples)

Reads nbSamples from the Julia buffer. If there are less than nbSamples samples in the Julia buffer,
the remaining samples are read, the buffer is refilled, and a total of the nbSamples is read.
Returns a newly allocated `Array{ComplexF32}`.

# Arguments
- `pluto::PlutoSDR` : the radio to get receive the samples from.
- `nbSamples::Int` : the number of samples to receive.

# Returns
- `buffer::Array{ComplexF32}` : an array with nbSamples complex values.
"""
function recv(pluto::PlutoSDR, nbSamples::Int)
    buffer = zeros(ComplexF32, nbSamples);
    recv!(buffer, pluto);
    return buffer;
end

"""
    recv!(sig, pluto)

Fills the `sig` input buffer with samples from the Julia buffer. If there are not enough samples in the Julia buffer,
it is refilled until `sig` is full.
Returns the number of samples filled or a negative error number.

# Arguments
- `sig::Array{ComplexF32}` : the buffer to be filled.
- `pluto::PlutoSDR` : the radio to read the samples from.
"""
function recv!(sig::Array{ComplexF32}, pluto::PlutoSDR)
    # FIXME: make this function somewhat error resilient
    samplesNeeded = length(sig);

    while (samplesNeeded > 0)
        # if needed refill the Julia buffer
        if (pluto.rx.buf.nb_samples == 0)
            nbNewSamples = refillJuliaBufferRX(pluto);
        end
        samplesQueried = min(samplesNeeded, pluto.rx.buf.nb_samples);
        # assuming the Julia buffer is filled until the end (which it should)
        sig[(end - samplesNeeded + 1):(end - samplesNeeded + samplesQueried)] =
            pluto.rx.buf.samples[(end - pluto.rx.buf.nb_samples + 1):(end - pluto.rx.buf.nb_samples + samplesQueried)];
        pluto.rx.buf.nb_samples -= samplesQueried;
        samplesNeeded -= samplesQueried;
    end

    return length(sig);
end

"""
    refillJuliaBufferRX(pluto)

Refills the radio buffer, decode the samples into `ComplexF32` values, and store those values in the `pluto` structure.
To access those samples : `pluto.rx.buf.samples`.

# Arguments
- `pluto::PlutoSDR` : the radio to receive the samples from, and the structure storing those samples.

# Returns
- `nbSamples::Int` : the number of samples stored into the Julia buffer.
"""
function refillJuliaBufferRX(pluto::PlutoSDR)
    nbytes = C_iio_buffer_refill(pluto.rx.buf.C_ptr);
    if (nbytes < 0)
        return nbytes;
    end
    # intermediary buffers from which complex values are computed
    nbytes_i = C_iio_channel_read!(pluto.rx.iio.rx0_i, pluto.rx.buf.C_ptr, pluto.rx.buf.i_raw_samples);
    nbytes_q = C_iio_channel_read!(pluto.rx.iio.rx0_q, pluto.rx.buf.C_ptr, pluto.rx.buf.q_raw_samples);

    # find the appropriate types and their max values for normalization
    bytes_per_value = pluto.rx.buf.C_step ÷ 2;
    utype, type, utype_norm, type_norm = nbytesToType(bytes_per_value);

    pluto.rx.buf.samples = reinterpret(Complex{type}, [pluto.rx.buf.i_raw_samples pluto.rx.buf.q_raw_samples]'[:]) / Float32(type_norm);
    pluto.rx.buf.nb_samples = (nbytes_i + nbytes_q) ÷ pluto.rx.buf.C_step;

    return pluto.rx.buf.nb_samples;
end

"""
   print(pluto)
   Print the current radio configuration
# Arguments
- `pluto::PlutoSDR` : the radio to receive the samples from, and the structure storing those samples.

# Returns
- A string with the different configuration aspects
"""
function Base.print(pluto::PlutoSDR)
	# Print message
	strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n",pluto.rx.effectiveCarrierFreq/1e6,pluto.rx.effectiveSamplingRate/1e6);
    println("Current Pluto Configuration in Rx mode\n$strF");
	strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n",pluto.tx.effectiveCarrierFreq/1e6,pluto.tx.effectiveSamplingRate/1e6);
    println("Current Pluto Configuration in Tx mode\n$strF");

end

"""
 send(pluto,buffer;repeat=false)
Send the buffer `buffer` with the Pluto device. If flag is set to true, repeat the transmission indefinively. If use_internal_buffer is set to true, the buffer used is the one already stored in the Pluto device (it saves buffer re-alloc)
"""
function send(pluto::PlutoSDR,buffer,flag=false;use_internal_buffer=false)
    # We need to populate the internal buffer of the Pluto
    # As AD9361 takes 14 bits, we should somehow convert the input stream to Int, then <<2 to have MSB aligned, and write to the Pluto Tx buffer
    # From external buffer, populate internal Pluto buffer
    if use_internal_buffer == false
        # --- Transfer buffer in iternal IIO_Buffer
        N = populateBuffer!(pluto,buffer)
    else
        # We will use previously transfered buffer, but need the size
        N = pluto.tx.buf.nb_samples
    end
    cnt = 0
    while true
        # C_iio_buffer_push(pluto.tx.buf.C_ptr)
        tmp = C_iio_buffer_push_partial(pluto.tx.buf.C_ptr,N)
        # tmp = C_iio_buffer_push(pluto.tx.buf.C_ptr)
        cnt += tmp
        (flag == false) && break
    end
    return cnt
end

"""
Populate Pluto internal IIO buffer with the desired samples from buffer.
populateBuffer(pluto,buffer)
"""
function populateBuffer!(pluto,buffer::Vector)
    N = length(buffer)
    pluto.tx.buf.nb_samples = N
    @assert N < 1024 * 1024 "Transmitted buffer is too large (should be < 1024 * 1024)"
    # Recast pointer...
    cPos = pluto.tx.buf.C_first
    for n ∈ 1 : N
        ptr = Ptr{Int16}(cPos)
        # Input is ]-1;1[ ==> Convert to Int
        # << 2 to be compliant with Q14 of AD
        tmpI = Int16( round(real(buffer[n]) * 1 << 15)) << 2
        tmpQ = Int16( round(imag(buffer[n]) * 1 << 15)) << 2
        # We need to populate Pluto buffer
        # @Tx side we can directly play with pluto.tx.buff.C_first. Position is stored in ptr, and we navigate with the step from IIO_Buffer
        unsafe_store!(ptr,tmpI,1)
        unsafe_store!(ptr,tmpQ,2)
        cPos += pluto.tx.buf.C_step # As we write I and Q but C_step contains 2 words
    end
    return UInt(N)
end

# This function is significantly slower than the one above despite using less arrays and doing stuff manually.
# It is kept here to remind me why the weird interleaving of arrays is still my best solution.
##
# function refillJuliaBufferRX_ButItsSlower(pluto::PlutoSDR)
#     nbytes_C = C_iio_buffer_refill(pluto.rx.buf.C_ptr);
#     #  nbytes_Julia = C_iio_channel_read!(pluto.rx.iio.chn, pluto.rx.buf.C_ptr, pluto.rx.buf.raw_samples); # doesn't work ?
#
#     pluto.rx.buf.C_first = C_iio_buffer_first(pluto.rx.buf.C_ptr, pluto.rx.iio.rx0_i);
#     pluto.rx.buf.C_last  = C_iio_buffer_end(pluto.rx.buf.C_ptr);
#     pluto.rx.buf.C_step  = C_iio_buffer_step(pluto.rx.buf.C_ptr);
#
#     format = unsafe_load(ccall(
#         (:iio_channel_get_data_format, libIIO_jl.libIIO),
#         Ptr{iio_data_format}, (Ptr{iio_channel},),
#         pluto.rx.iio.rx0_q
#     ));
#
#     src_ptr = C_iio_buffer_first(pluto.rx.buf.C_ptr, pluto.rx.iio.rx0_i);
#     dst_ptr = pointer(pluto.rx.buf.raw_samples, 1);
#     len = Int64(format.length / 8 * format.repeat);
#
#     foreach(src_ptr:pluto.rx.buf.C_step:pluto.rx.buf.C_last-1) do src_ptr
#         ccall(
#             (:iio_channel_convert, libIIO_jl.libIIO),
#             Cvoid, (Ptr{iio_channel}, Ptr{Cvoid}, Ptr{Cvoid}),
#             pluto.rx.iio.chn, dst_ptr, src_ptr
#         );
#         dst_ptr += len;
#     end
#
#     bytes_per_value = pluto.rx.buf.C_step ÷ 2;
#     utype, type, utype_norm, type_norm = nbytesToType(bytes_per_value);
#
#     pluto.rx.buf.samples = reinterpret(Complex{type}, pluto.rx.buf.raw_samples) / Float32(type_norm);
#     pluto.rx.buf.nb_samples = pluto.rx.buf.C_size;
#
#     return pluto.rx.buf.nb_samples;
# end

end


text:
"""
Bequest function after DeNardi (2004)
"""
function bequest(b::Float64, p::Param)
	@assert b >= 0.0  # no negative bequests
	p.ν * u(b + p.bbar,p)
end

function bequest(b::Vector{Float64}, p::Param)
	n = length(b)
	y = similar(b)
	for i in 1:n
		y[i] = bequest(b[i],p)
	end
	y
end

#utility functions

# utility without discrete choice
function u(x::Float64,p::Param)
	if p.gamma == 1.0
		log(x)
	else
		p.oneover_oneminusgamma * (x^p.oneminusgamma)
	end
end
function u(x::Array{T}, p::Param) where T
	n = length(x)
	y = zeros(T,n)
	for i in 1:n
		y[i] = u(x[i],p)
	end
	y
end

# utility with discrete choice
function u(x::Float64,working::Bool,p::Param)
	if p.gamma == 1.0
		log(x) - p.alpha*working
	else
		p.oneover_oneminusgamma * (x^p.oneminusgamma - 1.0) - p.alpha*working
	end
end
function u(x::Array{T}, working::Bool, p::Param) where T
	n = length(x)
	y = zeros(T,n)
	for i in 1:n
 		y[i] = u(x[i],working,p)
	end
	y
end

"bequest with discrete choice"
function bequest(b::Float64, working::Bool, p::Param)
	bequest(b,p) - p.alpha * working
end

function bequest(b::Vector{Float64}, working::Bool, p::Param)
	n = length(b)
	y = similar(b)
	for i in 1:n
		y[i] = bequest(b[i],working,p)
	end
	y
end

# partial derivative of utility wrt c
function up(c::Float64,p::Param)
	if p.gamma == 1.0
		1.0 / c
	else
		c ^ (p.neg_gamma)
	end
end
function up(c::Array{Float64,2},p::Param)
	n = length(c)
	x = zeros(size(c))
	for i in 1:n
		x[i] = up(c[i],p)
	end
	x
end
function up!(c::Array{Float64},p::Param)
	for i in eachindex(c)
		c[i] = up(c[i],p)
	end
end
function up(c::Array{Float64},p::Param)
	n = length(c)
	x = zeros(n)
	for i in 1:n
		x[i] = up(c[i],p)
	end
	x
end

# inverse of partial derivative
function iup(u::Float64,p::Param)
	if p.gamma == 1.0
		1.0 / u
	else
		u ^ p.neg_oneover_gamma
	end
end
function iup(u::Vector{Float64},p::Param)
	n = length(u)
	x = zeros(n)
	for i in 1:n
		x[i] = iup(u[i],p)
	end
	return x
end

NBLt(ylow,t::Int,p::Param) = sum(-p.R^(-(j-1) ) * income(j,p,ylow) for j in t:p.nT-1)

nbl(ylow,it,p::Param) = sum(-p.R^(j) * income(it+j,p,ylow) for j in 1:((p.nT)-1)-it )
nblv(ylow,it,p::Param) = [-p.R^(j) * income(it+j,p,ylow) for j in 1:((p.nT)-1)-it ]

# NBLt(ylow,t::Int,p::Param) = sum(-p.R^(-j ) * income(j,p,ylow) for j in (p.nT-1)-t:-1:0)
NBLs(ylow,p::Param) = [-p.R^(-j) * income(j,p,ylow) for j in (p.nT-1)-1:-1:0]

# nextbound(ylow,it,p::Param) = sum(-income(j,p,ylow)*p.R^(-(j-it)) for j in (it+1):(p.nT-1))
nextbound(ylow,it,p::Param) = -income(it+1,p,ylow)*p.R^(-1)
pnextbound(ylow,it,p::Param) = ["income($j,p,ylow)*p.R^($(-(j-it)))" for j in (it+1):(p.nT-1)]

"""
	abounds(ylow,p::Param)

returns a vector with the lowest possible asset level in each period.
This assumes the worst case of income draw in each period and zero consumption, while ``a_T=0``.
"""
function abounds(ylow,p::Param)
	x = [nextbound(ylow,it,p) for it in 1:p.nT-2]
	prepend!(x, -ylow + x[1]) # adds period 1 bound
	x
	# n = NBLs(ylow,p)  # vector of borrowing limit in period t
	# reverse(cumsum(reverse(n)))
end

allincomes(m::Model,p::Param) = [DCEGM.income(j,p,y) for y in m.yvec, j in 1:p.nT]

"expected borrowing limit"
function EBL(m::Model,p::Param)
	yy = [DCEGM.income(j,p,y) for y in m.yvec, j in 1:p.nT]  # income at each age and each state
	r = similar(yy)
	r[:,p.nT] = yy[:,p.nT]

	for it in p.nT-1:-1:1
		r[:,it] = m.ywgt * yy[:,it+1]
	end
	return r
end

"""
rouwenhorst AR1 approximation

This is taken from [http://karenkopecky.net/RouwenhorstPaperFinal.pdf](Karen Kopecky's paper)

"""
function rouwenhorst(rho::Float64,mu_eps,sigma_eps,n)
	q = (rho+1)/2
	nu = ((n-1)/(1-rho^2))^(1/2) * sigma_eps
	P = [q 1-q ; 1-q q]

	for i=2:n-1

		# P = q * vcat(hcat(P , zeros(i,1)),zeros(1,i+1)) .+ (1-q).* vcat( hcat(zeros(i,1),P), zeros(1,i+1)) .+
		# (1-q) .* vcat(zeros(1,i+1),hcat(P,zeros(i,1))) .+ q .*vcat(zeros(1,i+1),hcat(zeros(i,1),P))
		# P[2:i,:] = P[2:i,:] ./ 2

		P = q .* [P zeros(i,1); zeros(1,i+1)] .+ (1-q) .* [zeros(i,1) P; zeros(1,i+1)] .+
		(1-q) .* [zeros(1,i+1); P zeros(i,1)] .+   q   .* [zeros(1,i+1); zeros(i,1) P]

		P[2:i,:] = P[2:i,:] ./ 2
	end

	z = range(mu_eps/(1-rho)-nu,stop = mu_eps/(1-rho)+nu,length = n);
	return (collect(z),P)
end

# asset grid scaling
function scaleGrid(lb::Float64,ub::Float64,n::Int;logorder::Int=1,offset=0.05)
	if lb >= ub
		@warn("lower bound $lb is greater than upper bound $ub: forcing lb < ub.")
		lb = ub - abs(ub)*offset
	end
	out = zeros(n)
	if logorder==0
		out    = collect(range(lb,stop = ub,length = n))
	elseif logorder==1
		off = 1
		if lb<0
			off = 1 - lb #  adjust in case of neg bound
		end
		out[1] = log(lb + off)
		out[n] = log(ub + off)
		out    = collect(range(out[1],stop = out[n],length = n))
		out    = exp.(out) .- off
	elseif logorder==2
		off = 1
		if lb<0
			off = 1 - lb #  adjust in case of neg bound
		end
		out[1] = log( log(lb + off) + off )
		out[n] = log( log(ub + off) + off )
		out    = collect(range(out[1],stop = out[n],length = n))
		out    = exp.( exp.(out) .- off ) .- off
	elseif logorder==3
		off = 1
		if lb<0
			off = 1 - lb #  adjust in case of neg bound
		end
		out[1] = log( log( log(lb + off) + off ) + off )
		out[n] = log( log( log(ub + off) + off ) + off )
		out    = collect(range(out[1],stop = out[n],length = n))
		out    = exp.( exp.( exp.(out) .- off ) .- off ) .- off
	else
		error("supports only up to tripple log grid")
	end
end

"""
lifecycle profile in income
"""
function income(it::Int,p::Param,shock::Float64)
	age = it + 19
	if it > p.retage
		p.pension
	else
		exp( p.inc0 + p.inc1*age - p.inc2*(age^2) + shock)
	end
end
function income(it::Int,p::Param,shock::Array{Float64})
	x = similar(shock)
	for i=1:length(x)
		x[i] = income(it,p,shock[i])
	end
	return x
end

function quadpoints(n,lbnd,ubnd)

   x2  = ubnd
   x1  = lbnd
   x   = zeros(n)
   w   = zeros(n)
   EPS = 3.e-14
   m   = floor((n+1)/2)
   xm  = (x2+x1)/2
   xl  = (x2-x1)/2
   i = 1
   z1 = 1.e99
   pp = p1 = p2 = p3 = z = 0.0
   while (i <= m)
	   z  = cos(pi*(i-0.25)/(n+0.5))
	   while (abs(z-z1)>EPS)
	       p1 = 1
	       p2 = 0
	       j=1
	       while (j <= n)
				 p3 = copy(p2)
				 p2 = copy(p1)
				 p1 = ((2*j-1)*z*p2-(j-1)*p3)/j
		         j=j+1
	       end
	       pp = n*(z*p1-p2)/(z*z-1)
	       z1 = z
	       z  = z1 - p1/pp
    	end
     x[i]     = xm - xl*z
     x[n+1-i] = xm + xl*z
     w[i]     = 2*xl/((1-z*z)*pp*pp)
     w[n+1-i] =w[i]
     i = i+1
   end
   return (x,w)
end

function linearapprox(x::Vector{Float64},y::Vector{Float64},xi::Float64,lo::Int,hi::Int)
	r = 0.0
	n = length(x)
	@assert n==length(y)

	# determining bounds
	if xi == x[1]
		r = y[1]
		return r
	elseif xi < x[1]
		# get linear approx below
		@inbounds r = y[1] + (y[2] - y[1]) * (xi - x[1])  / (x[2] - x[1])
		return r
	end
	if xi == x[n]
		r = y[n]
		return (r,n)
	elseif xi > x[n]
		# get linear approx above
		@inbounds r = y[n] + (y[n] - y[n-1]) * (xi - x[n])  / (x[n] - x[n-1])
		return r
	end

	# if have to find interval
	if hi - lo > 1
		jinf = searchsortedlast(x,xi,lo,hi,Base.Forward)	# get rid
	# if not, lo is jinf
	else
		jinf = lo
	end
	@inbounds r = (y[jinf] * (x[jinf+1] - xi) + y[jinf+1] * (xi - x[jinf]) ) / (x[jinf+1] - x[jinf])
	return r
end
linearapprox(x::Vector{Float64},y::Vector{Float64},xi::Float64) = linearapprox(x,y,xi,1,length(x))

function linearapprox(x::Vector{Float64},y::Vector{Float64},xi::Vector{Float64})
	n = length(xi)
	z = zeros(n)
	for i in 1:n
		z[i] = linearapprox(x,y,xi[i])
	end
	return z
end


text:
# This file was generated, do not modify it. # hide
cap_gr = groupby(capacity, [:country, :primary_fuel]);


text:
x1, x2, y1, y2 = [parse(Int, m.match) for m in eachmatch(r"-?\d+", readline())]
count = 0
function simulate(yspeed, xspeed)
    x, y, maxy = 0, 0, 0
    while x <= x2 && y >= y1
        y, x = y + (yspeed -= 1), x + (xspeed -= xspeed > 0)
        maxy = max(y, maxy)
        x ∈ x1:x2 && y ∈ y1:y2 && (global count += 1) > 0 && return maxy
    end
    return 0
end
maximum(simulate(ys, xs) for ys=y1:200, xs=1:x2+1) |> m -> println("$m\n$count")

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
# Exercise 11-1
function histogram(s)
    d = Dict()
    for c in s
        d[c] = get(d, c, 0) + 1
    end
    d
end

# println(histogram("brontosaurus"))

# Exercise 11-2
function readerset()
    fin = open("words.txt")
    words = Set{String}()
    for line in eachline("words.txt")
        push!(words, line)
    end
    return words
end

# Exercise 11-3
function invertdict(d)
    inverse = Dict{Int64, Array}()
    for key in keys(d)
        val = d[key]
        push!(get!(inverse, val, []), key)
    end
    return inverse
end

# println(invertdict(histogram("parrot")))

# Exercise 11-4
function ack(m, n)
    if m == 0
        return n + 1
    elseif m > 0 && n == 0
        return ack(m - 1, 1)
    else
        return ack(m - 1, ack(m, n - 1))
    end
end

table = Dict{Pair, Int64}()
function ackmemo(m, n)
    if m == 0
        return n + 1
    elseif m > 0 && n == 0
        return ackmemo(m - 1, 1)
    elseif Pair(m, n) in keys(table)
        return table[Pair(m, n)]
    else
        table[Pair(m, n)] = ackmemo(m - 1, ackmemo(m, n - 1))
        return table[Pair(m, n)]
    end
end

# println(ack(3, 4) == ackmemo(3, 4))

function ex11_4(m, n)
    @time begin
        ackmemo(m, n)
    end
    @time begin
        ack(m, n)
    end
end

# ex11_4(3, 4)

# Exercise 11-5
# Same as 10-7, unique uses Set, same idea.
function hasduplicates(t)
    return unique(t) != t
end

# Exercise 11-6
function rotatechar(c, i)
    if i < 0
        i = mod(i, 26)
    end
    if islowercase(c)
        return Char((c - 'a' + i) % 26 + 'a')
    else
        return Char((c - 'A' + i) % 26 + 'A')
    end
end

function rotateword(str, i)
    result = ""
    for char in str
        result *= rotatechar(char, i)
    end
    return result
end

function rotatepairs(word, set)
    pairs = Array{Pair, 1}();
    for i = 1:25
        rot = rotateword(word, i)
        if rot in set
            push!(pairs, Pair(word, rot))
        end
    end
    return pairs
end

function ex11_6()
    set = readerset()
    for word in set
        pairs = rotatepairs(word, set)
        if !isempty(pairs)
            println(pairs)
        end
    end
end

# ex11_6()


text:
#this file is to randomly selected some individuals to have no omics data.
using Distributions, DelimitedFiles,Random
Random.seed!(123)

cd("/Users/tianjing/Box/omics/omics_missing_pattern_legarra")

pct_all=[0.1,0.3,0.5,0.7,0.8,0.9,0.95,0.99]
repAll=20  #rep1,...,rep20; in each rep, randomly selected some individuals to have no omics data.

for i in 1:length(pct_all)
    pct=pct_all[i]
    nMissingInd=Int(floor(955*pct))
    println(pct*100," % training individuals (i.e., $nMissingInd training individuals) has no missing data.")

    for rep in 1:repAll
        missingIndIndex=sort(sample(1:955,nMissingInd,replace=false))
        println(length(missingIndIndex), " individuals are sampled in rep$rep.")
        writedlm("missing_index/missingIndex.pct$pct.rep$rep.txt",missingIndIndex)
    end
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
""" Define a schema and ACSet which can be used as the base for most individual-based epidemiological models.
    Alone, this module allows for simulation of Markov models.
"""
module SchemaBase

export TheoryIBM, AbstractIBM, IBM,
    npeople, nstate, statelabel, get_index_state,
    queue_state_update, apply_state_updates,
    render_states,
    initialize_states, reset_states,
    simulation_loop

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Present
using Catlab.Theories

""" ACSet definition for a basic individual-based model
    See Catlab.jl documentation for description of the `@present` syntax.
"""
@present TheoryIBM(FreeSchema) begin
    Person::Ob
    State::Ob

    state::Hom(Person, State)
    state_update::Hom(Person, State)

    StateLabel::AttrType
    statelabel::Attr(State, StateLabel)
end

""" An abstract ACSet for a basic Markov individual-based model.
"""
@abstract_acset_type AbstractIBM

"""
A concrete ACSet for a basic Markov individual-based model inheriting from `AbstractIBM`.
"""
@acset_type IBM(TheoryIBM, index = [:state, :state_update], unique_index = [:statelabel]) <: AbstractIBM

"""
    npeople(model::AbstractIBM, states)

Return the number of people in some set of `states` (an element of the State Ob).
If called without the argument `states`, simply return the total population size.
"""
npeople(model::AbstractIBM) = nparts(model, :Person)
npeople(model::AbstractIBM, states) = length(incident(model, states, [:state, :statelabel]))

"""
    nstate(model::AbstractIBM)

Return the size of the finite state space.
"""
nstate(model::AbstractIBM) = nparts(model, :State)

"""
    statelabel(model::AbstractIBM)

Return the labels (names) of the states in the finite state space.
"""
statelabel(model::AbstractIBM) = subpart(model, :statelabel)

"""
    get_index_state(model::AbstractIBM, states)

Return an integer vector giving the persons who are in the states specified in `states`.
If called without the argument `states`, simply return everyone's index.
"""
get_index_state(model::AbstractIBM, states) = incident(model, states, [:state, :statelabel])
get_index_state(model::AbstractIBM) = parts(model, :Person)

"""
    queue_state_update(model::AbstractIBM, persons, state)

For persons specified in `persons`, queue a state update to `state`, which will be applied at the
end of the time step.
"""
function queue_state_update(model::AbstractIBM, persons, state)
    if length(persons) > 0
        state_index = incident(model, state, :statelabel)
        state_index > 0 || throw(ArgumentError("state $(state) is not is the set of state labels"))
        set_subpart!(model, persons, :state_update, state_index)
    end
end

"""
    apply_state_updates(model::AbstractIBM)

Apply all queued state updates.
"""
function apply_state_updates(model::AbstractIBM)
    for state = parts(model, :State)
        people_to_update = incident(model, state, :state_update)
        if length(people_to_update) > 0
            set_subpart!(model, people_to_update, :state, state)
        end
    end
    set_subpart!(model, :state_update, 0)
end

"""
    render_states(model::AbstractIBM, steps::Integer)

Return a tuple whose first element is a matrix containing counts of
states (columns) by time step (rows), and whose second element is a _process_
function which can be used in the simulation loop.
"""
function render_states(model::AbstractIBM, steps::Integer)
    out = Array{Int64}(undef, steps, nstate(model))

    output_states(t::Int) = begin
        out[t, :] = [length(incident(model, i, :state)) for i = parts(model, :State)]
    end

    return (out, output_states)
end

"""
    initialize_states(model::AbstractIBM, initial_states, state_labels::Vector{String})

Initialize the categorical states of a model. The argument `initial_states` can either
be provided as a vector of integers, corresponding to the internal storage of the ACSet,
or as a vector of strings. It should be equal in length to the population which is to be
simulated.
"""
function initialize_states(model::AbstractIBM, initial_states::Vector{T}, state_labels::Vector{String}) where {T <: Integer}
    length(unique(initial_states)) <= length(state_labels) || throw(ArgumentError("'initial_states' has more unique values than 'state_labels', please fix"))
    if nparts(model, :State) > 0
        reset_states(model, initial_states)
    else
        add_parts!(model, :State, length(state_labels), statelabel = state_labels)
        people = add_parts!(model, :Person, length(initial_states))
        set_subpart!(model, people, :state, initial_states);
    end
end

function initialize_states(model::AbstractIBM, initial_states::Vector{String}, state_labels::Vector{String})
    length(unique(initial_states)) <= length(state_labels) || throw(ArgumentError("'initial_states' has more unique values than 'state_labels', please fix"))
    if nparts(model, :State) > 0
        reset_states(model, initial_states)
    else
        add_parts!(model, :State, length(state_labels), statelabel = state_labels)
        people = add_parts!(model, :Person, length(initial_states))
        set_subpart!(model, people, :state, indexin(initial_states, state_labels));
    end
end

"""
    initialize_states(model::AbstractIBM, initial_states)

Reset a model's categorical states.
"""
function reset_states(model::AbstractIBM, initial_states::Vector{T}) where {T <: Integer}
    nparts(model, :Person) == length(initial_states) || throw(ArgumentError("'initial_states' must be equal to the number of persons in the model"))
    set_subpart!(model, :state_update, 0)
    set_subpart!(model, :state, initial_states);
end

function reset_states(model::AbstractIBM, initial_states::Vector{String})
    nparts(model, :Person) == length(initial_states) || throw(ArgumentError("'initial_states' must be equal to the number of persons in the model"))
    set_subpart!(model, :state_update, 0)
    set_subpart!(model, :state, indexin(initial_states, subpart(model, :statelabel)));
end

"""
    simulation_loop(model::AbstractIBM, processes::Union{Function, AbstractVector{Function}}, steps::Integer)

A simple predefined simulation loop for basic (no events) individual based models. Processes are called first,
followed by state updates.
"""
function simulation_loop(model::AbstractIBM, processes::Union{Function, AbstractVector{Function}}, steps::Integer)
    if processes isa Function
        processes = [processes]
    end
    for t = 1:steps
        for p = processes
            p(t)
        end
        apply_state_updates(model)
    end
end

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    month_cnv(number[, shor=true, up=true, low=true]) -> month_name
    month_cnv(name) -> number

### Purpose ###

Convert between a month English name and  the equivalent number.

### Explanation ###

For example, converts from "January" to 1  or vice-versa.

### Arguments ###

The functions has two methods, one with numeric input (and three possible
boolean keywords) and the other one with string input.

Numeric input arguments:

* `number`: the number of the month to be converted to month name.
* `short` (optional boolean keyword): if true, the abbreviated (3-character)
  name of the month will be returned, e.g. "Apr" or "Oct".  Default is false.
* `up` (optional boolean keyword): if true, the name of the month will be all in
  upper case, e.g. "APRIL" or "OCTOBER".  Default is false.
* `low` (optional boolean keyword): if true, the name of the month will be all
  in lower case, e.g. "april" or "october".  Default is false.

String input argument:

* `name`: month name to be converted to month number.

All mandatory arguments can be provided either as a single element or as an
array.

### Output ###

The month name or month number, depending on the input.  For numeric input, the
format of the month name is influenced by the optional keywords.

### Example ###

``` julia
month_cnv(["janua", "SEP", "aUgUsT"])
# => 3-element Array{Integer,1}:
#     1
#     9
#     8
month_cnv([2, 12, 6], short=true, low=true)
# => 3-element Array{UTF8String,1}:
#     "feb"
#     "dec"
#     "jun"
```

"""
##### Numeric input
function month_cnv(number::Integer; short::Bool=false,
                   up::Bool=false, low::Bool=false)
    if short
        name = Dates.VALUETOMONTHABBR["english"][number]
    else
        name = Dates.VALUETOMONTH["english"][number]
    end
    if up
        name = uppercase(name)
    elseif low
        name = lowercase(name)
    end
    return name
end

function month_cnv{N<:Integer}(numbers::AbstractArray{N}; short::Bool=false,
                               up::Bool=false, low::Bool=false)
    # Use same type as "Dates.VALUETOMONTH".
    months = similar(numbers, String)
    for i in eachindex(numbers)
        months[i] = month_cnv(numbers[i], short=short, up=up, low=low)
    end
    return months
end

##### String input
function month_cnv(name::AbstractString)
    name = strip(name)
    if length(name) >= 3
        output::Integer = eval(parse("Dates.$(ucfirst(lowercase(name[1:3])))"))
        return output
    else
        # Do the same as original AstroLib: return -1 for input string shorter
        # than 3.
        return -1
    end
end

@vectorize_1arg AbstractString month_cnv


text:
struct BatesEngine{C <: ComplexLogFormula} <: AbstractHestonEngine{HestonGaussLaguerre}
  model::BatesModel
  evaluations::Int
  cpxLog::C
  integration::HestonGaussLaguerre
end

function BatesEngine(batesModel::BatesModel)
  evals = 0
  cpxLog = Gatheral()
  integration = HestonGaussLaguerre(144)

  return BatesEngine(batesModel, evals, cpxLog, integration)
end

function add_on_term(pe::BatesEngine, phi::Float64, t::Float64, j::Int)
  batesModel = pe.model
  nu = get_nu(batesModel)
  delta2 = 0.5 * get_delta(batesModel) * get_delta(batesModel)
  lambda = get_lambda(batesModel)
  i = j == 1 ? 1.0 : 0.0
  g = complex(i, phi)

  return t * lambda * (expm1(nu * g + delta2 * g * g) - g * expm1(nu * delta2))
end


text:
module CryptoGroups

using Random: AbstractRNG

abstract type AbstractGroup end
abstract type CyclicGroup <: AbstractGroup end

import Base.^
^(G::AbstractGroup,n) = typeof(G)(G.G^n,G.q,G.t)

"""
The value of the group element.
"""
value(G::AbstractGroup) = error("Must be implemented by Group.")

"""
The number of bits one should use for the secret a and b in Diffie-Hellman key exchange.
"""
security(G::AbstractGroup) = error("Must be implemented by Group.")
order(G::AbstractGroup) = error("Must be implemented by Group.")
binary(G::AbstractGroup) = value(G) ### Default

include("utils.jl")
include("primegroup.jl")
include("rfc.jl")
include("ellipticgroup.jl")

export PrimeGroup, EllipticGroup, AbstractGroup, CyclicGroup, value, security, order, binary

end # module


text:
using Test
using DiffinDiffsBase

using DataAPI: refarray, refvalue, refpool, invrefpool
using DataFrames
using Dates: Date, Year
using DiffinDiffsBase: @fieldequal, unpack, @unpack, checktable, hastreat, parse_treat,
    isintercept, isomitsintercept, parse_intercept!,
    ncol, nrow, _mult!,
    checkdata!, grouptreatintterms, groupxterms, groupcontrasts,
    checkvars!, groupsample, makeweights,
    _totermset!, parse_didargs!, _treatnames, _parse_bycells!, _parse_subset, _nselected,
    treatindex, checktreatindex, DefaultExportFormat
using LinearAlgebra: Diagonal
using Missings: allowmissing, disallowmissing
using PooledArrays: PooledArray
using StatsBase: Weights, UnitWeights
using StatsModels: termvars
using StatsProcedures: _f, _byid, groupargs, copyargs, pool
using Tables: table

import Base: ==, show
import DiffinDiffsBase: valid_didargs
import StatsProcedures: required, result

include("testutils.jl")

const tests = [
    "tables",
    "utils",
    "time",
    "ScaledArrays",
    "terms",
    "treatments",
    "parallels",
    "operations",
    "procedures",
    "did"
]

printstyled("Running tests:\n", color=:blue, bold=true)

@time for test in tests
    include("$test.jl")
    println("\033[1m\033[32mPASSED\033[0m: $(test)")
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################
"""
Example considering a multi-state Markov model considering a subsystem of three
identical and independent channels taken from:

> Reliability of Safety-Critical System - Theory and Applications by Marvin
  Rausand (2014)

Correspondig slides can be found at:
https://www.ntnu.edu/documents/624876/1277046207/SIS+book+-+chapter+05+-+Introduction+to+Markov+methods/d98c074f-a48e-45b1-a016-668350efbbe7

States:
- 1 : Three channels are functioning
- 2 : Two channels are functioning and one is failed
- 3 : One channel is functioning and two are failed
- 4 : Three channels are failed
"""

# load pkgs
using Unitful
using MultiStateSystems

# initialize the state-transition diagram
std²³ = STD()

# add the states to the std
add_states!(std²³, init = [1.0, 0.0, 0.0, 0.0])

# add the transitions to the std
add_transition!(std²³,  states = (1,2),
                        rate   = 3λ)
add_transitions!(std²³, states = [(2,3), (2,1)],
                        rate   = [2λ, μ])
add_transitions!(std²³, states = [(3,4), (3,1)],
                        rate   = [λ, μ])
add_transition!(std²³,  states = (4,1),
                        rate   = μ)

# return the std
return std²³


text:
Base.@deprecate min(V1::ElementarySpace, V2::ElementarySpace) infimum(V1, V2)
Base.@deprecate max(V1::ElementarySpace, V2::ElementarySpace) supremum(V1, V2)
Base.@deprecate(infinum(V1, V2...), infimum(V1, V2...))

Base.@deprecate(fusiontreetype(::Type{I}, ::StaticLength{N}) where {I<:Sector, N},
    fusiontreetype(I, N))

abstract type RepresentationSpace{I<:Sector} end
Base.@deprecate(RepresentationSpace(args...), GradedSpace(args...))
Base.@deprecate(RepresentationSpace{I}(args...) where {I}, Vect[I](args...))

Base.@deprecate(×(a::Sector, b::Sector), ⊠(a,b))

Base.@deprecate(
    permuteind(t::TensorMap, p1::IndexTuple, p2::IndexTuple=(); copy::Bool = false),
    permute(t, p1, p2; copy = copy))

@noinline function Base.getindex(::ComplexNumbers, I::Type{<:Group})
    S = Rep[I]
    if repr(S) == "GradedSpace[Irrep[$I]]"
        Base.depwarn("`ℂ[$I]` is deprecated because $I represents a group rather than its representations, use `ℂ[Irrep[$I]]`, `Rep[$I]` or `GradedSpace[Irrep[$I]]` instead.",
        ((Base.Core).Typeof(Base.getindex)).name.mt.name)
    else
        Base.depwarn("`ℂ[$I]` is deprecated because $I represents a group rather than its representations, use `ℂ[Irrep[$I]]`, `Rep[$I]`, `$S` or `GradedSpace[Irrep[$I]]` instead.",
        ((Base.Core).Typeof(Base.getindex)).name.mt.name)
    end
    return Rep[I]
end

@noinline function ℤ₂(args...)
    Base.depwarn("`ℤ₂(args...)` is deprecated, use `Z2Irrep(args...)` or ``Irrep[ℤ₂](args...)` instead.", ((Base.Core).Typeof(ℤ₂)).name.mt.name)
    Irrep[ℤ₂](args...)
end
@noinline function ℤ₃(args...)
    Base.depwarn("`ℤ₃(args...)` is deprecated, use `Z3Irrep(args...)` or ``Irrep[ℤ₃](args...)` instead.", ((Base.Core).Typeof(ℤ₃)).name.mt.name)
    Irrep[ℤ₃](args...)
end
@noinline function ℤ₄(args...)
    Base.depwarn("`ℤ₄(args...)` is deprecated, use `Z4Irrep(args...)` or ``Irrep[ℤ₄](args...)` instead.", ((Base.Core).Typeof(ℤ₄)).name.mt.name)
    Irrep[ℤ₄](args...)
end
@noinline function U₁(args...)
    Base.depwarn("`U₁(args...)` is deprecated, use `U1Irrep(args...)` or ``Irrep[U₁](args...)` instead.", ((Base.Core).Typeof(U₁)).name.mt.name)
    Irrep[U₁](args...)
end
@noinline function CU₁(args...)
    Base.depwarn("`CU₁(args...)` is deprecated, use `CU1Irrep(args...)` or ``Irrep[CU₁](args...)` instead.", ((Base.Core).Typeof(CU₁)).name.mt.name)
    Irrep[CU₁](args...)
end
@noinline function SU₂(args...)
    Base.depwarn("`SU₂(args...)` is deprecated, use `SU2Irrep(args...)` or ``Irrep[SU₂](args...)` instead.", ((Base.Core).Typeof(SU₂)).name.mt.name)
    Irrep[SU₂](args...)
end

Base.@deprecate(permuteind!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, p1, p2),
                permute!(tdst, tsrc, p1, p2))

function Base.getindex(::Type{GradedSpace}, ::Type{I}) where {I<:Sector}
    @warn "`getindex(::Type{GradedSpace}, I::Type{<:Sector})` is deprecated, use `ℂ[I]`, `Vect[I]`, or, if `I == Irrep[G]`, `Rep[G]` instead." maxlog = 1
    return Vect[I]
end

function Base.getindex(::ComplexNumbers, d1::Pair{I, Int}, dims::Pair{I, Int}...) where {I<:Sector}
    @warn "`ℂ[s1=>n1, s2=>n2, ...]` is deprecated, use `ℂ[I](s1=>n1, s2=>n2, ...)`, `Vect[I](s1=>n1, s2=>n2, ...)` instead with `I = typeof(s1)`." maxlog = 1
    return Vect[I](d1, dims...)
end

function Base.getindex(::RealNumbers, d::Int)
    @warn "`ℝ[d]` is deprecated, use `ℝ^d` or `CartesianSpace(d)`." maxlog = 1
    return ℝ^d
end

function Base.getindex(::ComplexNumbers, d::Int)
    @warn "`ℂ[d]` is deprecated, use `ℂ^d` or `ComplexSpace(d)`." maxlog = 1
    return ℂ^d
end


text:
# extra Monoid instances

# Note, we do not implement neutral for functions like + or *, as this is already handled by Base.reduce_empty


text:
"""
A log conconcave function is majorized with a piecewise envelop, which on the original scale is piecewise
exponential. As the resulting extremely precise envelop adapts, the rejection rate dramatically decreases.
"""
module AdaptiveRejectionSampling
# ------------------------------
using ForwardDiff # For automatic differentiation, no user nor approximate derivatives
using StatsBase # To include the basic sample from array function
# ------------------------------
export Line, Objective, Envelop, RejectionSampler # Structures/classes
export run_sampler!, sample_envelop, eval_envelop, add_segment! # Methods
# ------------------------------

"""
    Line(slope::Float64, intercept::Float64)
Basic ensamble-unit for an envelop.
"""
mutable struct Line
    slope::Float64
    intercept::Float64
end

"""
    intersection(l1::Line, l2::Line)
Finds the horizontal coordinate of the intersection between lines
"""
function intersection(l1::Line, l2::Line)
    @assert l1.slope != l2.slope "slopes should be different"
    - (l2.intercept - l1.intercept) / (l2.slope - l1.slope)
end

"""
    exp_integral(l::Line, x1::Float64, x2::Float64)
Computes the integral
    ``LaTeX \\int_{x_1} ^ {x_2} \\exp\\{ax + b\\} dx. ``
The resulting value is the weight assigned to the segment [x1, x2] in the envelop
"""
function exp_integral(l::Line, x1::Float64, x2::Float64)
    a, b = l.slope, l.intercept
    exp(b) * (exp(a * x2) - exp(a * x1)) / a
end

"""
    Envelop(lines::Vector{Line}, support::Tuple{Float64, Float64})
A piecewise linear function with k segments defined by the lines `L_1, ..., L_k` and cutpoints
`c_1, ..., c_k+1` with `c1 = support[1]` and `c2 = support[2]`. A line L_k is active in the segment
[c_k, c_k+1], and it's assigned a weight w_k based on [exp_integral](@exp_integral). The weighted integral
over c_1 to c_k+1 is one, so that the envelop is interpreted as a density.
"""
mutable struct Envelop
    lines::Vector{Line}
    cutpoints::AbstractVector{Float64}
    weights::AbstractVector{Float64}
    size::Int

    Envelop(lines::Vector{Line}, support::Tuple{Float64, Float64}) = begin
        @assert issorted([l.slope for l in lines], rev = true) "line slopes must be decreasing"
        intersections = [intersection(lines[i], lines[i + 1]) for i in 1:(length(lines) - 1)]
        cutpoints = [support[1]; intersections; support[2]]
        @assert issorted(cutpoints) "cutpoints must be ordered"
        @assert length(unique(cutpoints)) == length(cutpoints) "cutpoints can't have duplicates"
        weights = [exp_integral(l, cutpoints[i], cutpoints[i + 1]) for (i, l) in enumerate(lines)]
        @assert Inf ∉ weights "Overflow in assigning weights"
        new(lines, cutpoints, weights, length(lines))
    end
end

"""
    add_segment!(e::Envelop, l::Line)
Adds a new line segment to an envelop based on the value of its slope (slopes must be decreasing
always in the envelop). The cutpoints are automatically determined by intersecting the line with
the adjacent lines.
"""
function add_segment!(e::Envelop, l::Line)
    # Find the position in sorted array with binary search
    pos = searchsortedfirst([-line.slope for line in e.lines], -l.slope)
    # Find the new cutpoints
    if pos == 1
        new_cut = intersection(l, e.lines[pos])
        # Insert in second position, first one is the support bound
        insert!(e.cutpoints, pos + 1, new_cut)
    elseif pos == e.size + 1
        new_cut = intersection(l, e.lines[pos - 1])
        insert!(e.cutpoints, pos, new_cut)
    else
        new_cut1 = intersection(l, e.lines[pos - 1])
        new_cut2 = intersection(l, e.lines[pos])
        splice!(e.cutpoints, pos, [new_cut1, new_cut2])
        @assert issorted(e.cutpoints)  "incompatible line: resulting intersection points aren't sorted"
    end
    # Insert the new line
    insert!(e.lines, pos, l)
    e.size += 1
    # Recompute weights (this could be done more efficiently in the future by updating the neccesary ones only)
    e.weights = [exp_integral(line, e.cutpoints[i], e.cutpoints[i + 1]) for (i, line) in enumerate(e.lines)]
end

"""
    sample_envelop(p::Envelop)
Samples an element from the density defined by the envelop `e` with it's exponential weights.
See [`Envelop`](@Envelop) for details.
"""
function sample_envelop(e::Envelop)
    # Randomly select lines based on envelop weights
    i = sample(1:e.size, weights(e.weights))
    a, b = e.lines[i].slope, e.lines[i].intercept
    # Use the inverse CDF method for sampling
    log(exp(-b) * rand() * e.weights[i] * a + exp(a * e.cutpoints[i])) / a
end

"""
    eval_envelop(e::Envelop, x::Float64)
Eval point a point `x` in the piecewise linear function defined by `e`. Necessary for evaluating
the density assigned to the point `x`.
"""
function eval_envelop(e::Envelop, x::Float64)
    # searchsortedfirst is the proper method for and ordered list
    pos = searchsortedfirst(e.cutpoints, x)
    if pos == 1 || pos == length(e.cutpoints) + 1
        return 0.0
    else
        a, b = e.lines[pos - 1].slope, e.lines[pos - 1].intercept
        return exp(a * x + b)
    end
end

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

"""
    Objective(logf::Function, support:)
    Objective(logf::Function, grad::Function)
Convenient structure to store the objective function to be sampled. It must receive the logarithm of
f and not f directly. It uses automatic differentiation by default, but the user can provide the derivative optionally.
"""
struct Objective
    logf::Function
    grad::Function
    Objective(logf::Function) = begin
        # Automatic differentiation
        grad(x) = ForwardDiff.derivative(logf, x)
        new(logf, grad)
    end
    Objective(logf::Function, grad::Function) = new(logf, grad)
end

"""
    RejectionSampler(f::Function, support::Tuple{Float64, Float64}[ ,δ::Float64])
    RejectionSampler(f::Function, support::Tuple{Float64, Float64}, init::Tuple{Float64, Float64})
An adaptive rejection sampler to obtain iid samples from a logconcave function `f`, supported in the
domain `support` = (support[1], support[2]). To create the object, two initial points `init = init[1], init[2]`
such that `loff'(init[1]) > 0` and `logf'(init[2]) < 0` are necessary. If they are not provided, the constructor
will perform a greedy search based on `δ`.

The argument `support` must be of the form `(-Inf, Inf), (-Inf, a), (b, Inf), (a,b)`, and it represent the
interval in which f has positive value, and zero elsewhere.

## Keyword arguments
- `max_segments::Int = 10` : max size of envelop, the rejection-rate is usually slow with a small number of segments
- `max_failed_factor::Float64 = 0.001`: level at which throw an error if one single sample has a rejection rate
    exceeding this value
"""
struct RejectionSampler
    objective::Objective
    envelop::Envelop
    max_segments::Int
    max_failed_rate::Float64
    # Constructor when initial points are provided
    RejectionSampler(
            f::Function,
            support::Tuple{Float64, Float64},
            init::Tuple{Float64, Float64};
            max_segments::Int = 25,
            max_failed_rate::Float64 = 0.001
    ) = begin
        @assert support[1] < support[2] "invalid support, not an interval"
        logf(x) = log(f(x))
        objective = Objective(logf)
        x1, x2 = init
        @assert x1 < x2 "cutpoints must be ordered"
        a1, a2 = objective.grad(x1), objective.grad(x2)
        @assert a1 >= 0 "logf must have positive slope at initial cutpoint 1"
        @assert a2 <= 0 "logf must have negative slope at initial cutpoint 2"
        b1, b2 = objective.logf(x1) - a1 * x1, objective.logf(x2) - a2 * x2
        line1, line2 = Line(a1, b1), Line(a2, b2)
        envelop = Envelop([line1, line2], support)
        new(objective, envelop, max_segments, max_failed_rate)
    end

    # Constructor for greedy search of starting points
    RejectionSampler(
            f::Function,
            support::Tuple{Float64, Float64},
            δ::Float64 = 0.5;
            search_range::Tuple{Float64, Float64} = (-10.0,10.0),
            kwargs...
    ) = begin
        logf(x) = log(f(x))
        grad(x) = ForwardDiff.derivative(logf, x)
        grid_lims = max(search_range[1], support[1]), min(search_range[2], support[2])
        grid = grid_lims[1]:δ:grid_lims[2]
        i1, i2 = findfirst(grad.(grid) .> 0.), findfirst(grad.(grid) .< 0.)
        @assert (i1 != nothing) &&  (i2 != nothing) "couldn't find initial points, please provide them or change `search_range`"
        x1, x2 = grid[i1], grid[i2]
        RejectionSampler(f, support, (x1, x2); kwargs...)
    end
end

"""
    run_sampler!(sampler::RejectionSampler, n::Int)
It draws `n` iid samples of the objective function of `sampler`, and at each iteration it adapts the envelop
of `sampler` by adding new segments to its envelop.
"""
function run_sampler!(s::RejectionSampler, n::Int)
    i = 0
    failed, max_failed = 0, trunc(Int, n / s.max_failed_rate)
    out = zeros(n)
    while i < n
        candidate = sample_envelop(s.envelop)
        acceptance_ratio = exp(s.objective.logf(candidate)) / eval_envelop(s.envelop, candidate)
        if rand() < acceptance_ratio
            i += 1
            out[i] = candidate
        else
            if length(s.envelop.lines) <= s.max_segments
                a = s.objective.grad(candidate)
                b = s.objective.logf(candidate) - a * candidate
                add_segment!(s.envelop, Line(a, b))
            end
            failed += 1
            @assert failed < max_failed "max_failed_factor reached"
        end
    end
    out
end
end #


text:
# generates the Lagrange basis for an array of roots at point x
function Lagrange(roots,x)
    nr    = length(roots)
    li    = Array(typeof(x),nr)

    for j = 1:nr
        li[j] = one(x)
        for i = 1:nr
            if i == j
                continue
            else
                li[j] *= (x-roots[i])/(roots[j]-roots[i])
            end
        end
    end

    return li
end

LegendreP(n::Int,x) = LegendrePgen(n::Int,x)[1]

function LegendrePgen{T}(n::Int,x::T)
    if n == 0
        return one(T), zero(T)
    else
        Pn1, Pn2 = LegendrePgen(n-1,x)
        return x*Pn1*(2n-1)/n-Pn2*(n-1)/n, Pn1
    end
end

# Generates hermite interpolation polynomials for the interval s∈[0,1]
# and valid for the first N derivatives (i.e. N=0 generates
# ps[1,1]=s-1 and ps[1,1]=s) etc.
function hermite01(n::Int)

    T = Rational{Int}
    e = eye(T,2n)
    M = zeros(T,2n,2n)

    for j = 1:2n
        p = Poly(e[:,j])
        for i1 = 0:n-1
            M[i1+1,  j] = polyval(polyder(p,i1),zero(T))
            M[i1+n+1,j] = polyval(polyder(p,i1), one(T))
        end
    end

    M⁻¹ = inv(M)                # here M⁻¹ is exact,
                                # typeof(M⁻¹)=Array{Rational{Int},2}

    ps = [ Poly(M⁻¹*e[:,i]) for i = 1:2n]
    return reshape(ps,n,2)'

end

# generates Gauss points and weights on the interval [0,1], this
# function comes from the RosettaCode page
function gauss(n)
    λ, Q = eig(SymTridiagonal(zeros(n), [ i / sqrt(4i^2 - 1) for i = 1:n-1 ]))
    x = (λ + 1) * 1/2
    w = [ 2*Q[1,i]^2 for i = 1:n ]/2
    return x, w
end

# generates Lobatto points and weights on the interval s∈[0,1]
function lobatto(n)

    if n < 2
        error("Number of points should be at least 2.")
    end

    #1 Compute the exact Legendre polynomial with rational coefficients
    p = LegendreP(n-1, Poly(Rational{Int}[0,1]))
    dp = polyder(p,1)

    #2 get the collocation points by solving Pₙ₋₁'(xᵢ)=0.  The roots
    #function returns Float64.
    x = [-1, sort(roots(dp)), 1]

    #3 compute weights as 2/(n(n-1)(Pₙ₋₁(xᵢ))^2)
    px = polyval(p,x)
    w = 2./px.^2/(n*(n-1))

    # we need to rescale the results from the interval [-1,1] to [0,1]
    return (x.+1)/2, w/2

end


text:
"""
    RMS(x[n], y[n-1], n)

`y[n] = sqrt((x[n]^2 +(n-1)*(y[n-1]^2))/n)`

Returns RMS value of `x[n]` when added to `y[n-1]` which will then be of length `n`.

### Examples

```jldoctest
julia> RMS(32, 19.2267175, 4)
23.0922064638398

julia> RMS(100, 10, 1_000)
10.483320084782301
```
"""
function RMS(xn, yn_, n)
    return sqrt((xn^2+(n-1)*(yn_^2))/n)
end

"""
    Decibels(a...)

`return = 20*log10(a...)`

Takes one to three values.

Returns (one to three values of) decibel level(s) as a function of 20 times the log10 of value `a` (, `b` ...).

### Examples

```jldoctest
julia> Decibels(2)
6.020599913279624

julia> Decibels(10, 100, 1_000)
(20.0, 40.0, 60.0)
```
"""
Decibels(x) = 20*log10(x)
Decibels(a, b) = (20*log10(a), 20*log10(b))
Decibels(a, b, c) = (20*log10(a), 20*log10(b), 20*log10(c))


text:
using Base, Logging

using Test, Distributed

if nworkers() > 1
    rmprocs(workers())
end
processors = 4 # Sys.CPU_CORES

if Base.JLOptions().code_coverage == 1
    addprocs(processors, exeflags = ["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
else
    addprocs(processors, exeflags = "--check-bounds=yes")
end

println("Workers:", nworkers())

using LinearAlgebra
using Statistics
using Random

using JuMP

using Ipopt
using Cbc
using GLPK
# using PowerModels

using MathOptInterface

const MOI = MathOptInterface
const MOIU = MOI.Utilities

using Juniper
include("util.jl")

opt_rtol = 1e-6
opt_atol = 1e-6

sol_rtol = 1e-3
sol_atol = 1e-3

function solve(m::Model)
    JuMP.optimize!(m)
    bm = JuMP.backend(m)
    return MOI.get(bm, MOI.TerminationStatus())
end

function getsolvetime(m::Model)
    bm = JuMP.backend(m)
    return MOI.get(bm, MOI.SolveTime())
end

function internalmodel(m::Model)
    bm = JuMP.backend(m)
    return bm.optimizer.model.inner
end

function getobjgap(m::Model)
    bm = JuMP.backend(m)
    return MOI.get(bm, MOI.RelativeGap())
end

juniper_strong_restart_2 = DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            strong_branching_perc = 25,
            strong_branching_nsteps = 2,
            strong_restart = true
            )

juniper_reliable_restart = DefaultTestSolver(
            branch_strategy=:Reliability,
            reliability_branching_perc = 25,
            reliability_branching_threshold = 2,
            strong_restart = true
            )

juniper_strong_restart = DefaultTestSolver(
                branch_strategy=:StrongPseudoCost,
                strong_branching_perc = 25,
                strong_restart = true
            )
juniper_strong_no_restart = DefaultTestSolver(
                branch_strategy=:StrongPseudoCost,
                strong_branching_perc = 25,
                strong_restart = false
            )

juniper_mosti = DefaultTestSolver(
                branch_strategy=:MostInfeasible,
            )

juniper_pseudo = DefaultTestSolver(
                branch_strategy=:PseudoCost,
            )

start = time()

@testset "Juniper" begin
    include("debug.jl")
    include("functions.jl")
    include("basic.jl")
    include("user_limits.jl")
    include("parallel.jl")
    include("fpump.jl")
    include("pod.jl")
    include("MOI_wrapper.jl")
    # include("power_models_acp.jl")
    # include("power_models_socwr.jl")
end
println("Time for all tests: ", time()-start)


text:
function foo()
    lines = readlines(open("data/7.txt", "r"))
    #= lines = readlines(open("h7", "r")) =#
    deps = Dict{Char, Array{Char, 1}}()
    for line in lines
        a = line[6]
        b = line[37]
        if !(a in keys(deps))
            deps[a] = []
        end
        push!(deps[a], b)
    end

    order = Array{Char, 1}()
    total = Set(vcat(values(deps)...))
    while !isempty(deps)

        s = Set(keys(deps))
        j = Set(vcat(values(deps)...))
        x = sort(collect(setdiff(s, j)))

        push!(order, x[1])
        delete!(deps, x[1])
        delete!(total, x[1])
    end
    x = vcat(order, sort(collect(total)))
    return join(x)
end

println(foo())


text:
using Plots

theme(:mydark)   # or any other dark theme

f(x) = sin(x) * cos(x/10) * exp(-x/10)
f1(x) = 0.7*sin(1.3x) + 0.4*sin(1.4x+5) * exp(x/20) - 0.27
f2(x) = 1.27*sinc(1.3x-6.5) - 0.42

let p = plot(f, 0, 10; size=(1300,500), color=RGBA(1,0,0,0.5), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f(x), 0, 10; color=RGBA(1,0,0,0.03), legend=false)
  end
  plot!(p, f2, 0, 10; color=RGBA(0,1,0,0.25), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f2(x), 0, 10; color=RGBA(0,1,0,0.01), legend=false)
  end
  plot!(p, f1, 0, 10; color=RGBA(0,0,1,0.5), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f1(x), 0, 10; color=RGBA(0,0,1,0.03), legend=false)
  end
  p
end

max_stars_repo_name: org-arl/unet-blog

max_issues_repo_name: org-arl/unet-blog

max_forks_repo_name: org-arl/unet-blog


text:
using Dolo

module Temp
    using Dolo

    struct DisG{ID, n_m, n_s, n_x, G_x, G_e}

        model::AbstractModel
        grid_endo
        grid_exo
        dprocess
        s0::Dolo.ListOfPoints{n_s}

        function DisG(model;
            discretization=Dict())

            grid, dprocess = Dolo.discretize(model; discretization...)
            grid_endo = grid.endo
            grid_exo = grid.exo

            n_m = length(model.symbols[:exogenous])
            n_s = length(model.symbols[:states])
            n_x = length(model.symbols[:controls])

            s0 = Dolo.nodes(grid_endo)

            Gx = typeof(grid_exo)
            Ge = typeof(grid_endo)

            ID = id(model)

            new{ID, n_m, n_s, n_x, Gx, Ge}(model, grid_endo, grid_exo, dprocess, s0)

        end

    end

    function size_mu(disg::DisG)

        return length(Dolo.nodes(disg.grid_endo)) * length(Dolo.nodes(disg.grid_exo))

    end

    function ergodic_distribution(sol, disg::DisG)

        dr = sol.dr

        return Dolo.ergodic_distribution(disg.model, dr, disg.grid_exo, disg.grid_endo, disg.dprocess)[2]

    end

    function (G::DisG)(μ0,x0; exo = nothing)

        ∂G_∂μ = Dolo.new_transition(G.model, G.dprocess, x0, G.grid_exo, G.grid_endo; exo=exo)

        μ1 = Dolo.new_distribution(∂G_∂μ, μ0)

        ∂G_∂x = (x) -> μ0'* (ForwardDiff.gradient(new_transition, x))[3]

        return μ1, ∂G_∂μ, ∂G_∂x
    end

end

filename = "./examples/models/rbc_mc.yaml"
model = Model(filename)
sol = time_iteration(model)
dr = sol.dr
grid_endo, grid_exo = dr.grid_endo, dr.grid_exo
length(Dolo.nodes(grid_endo))
length(Dolo.nodes(grid_exo))
F = Dolo.Euler(model)
x0 = F.x0
μ0 = zeros(40)

G = Main.Temp.DisG(model)

Main.Temp.ergodic_distribution(sol, G)

Main.Temp.size_mu(G)

μ1, ∂G_∂μ, ∂G_∂x = G(μ0,x0)
∂G_∂μ
∂G_∂x


text:
module LogNumbers
using Base

export LogNumber, to_number

struct LogNumber{T}
    sign::Int8
    exponent::T
end
LogNumber(s, e) = LogNumber(Int8(s), e)

LogNumber(x) = LogNumber(sign(x), log(abs(x)))
to_number(x::LogNumber) = x.sign*exp(x.exponent)

Base.sign(x::LogNumber) = x.sign

function mul_div(pm::Union{typeof(+), typeof(-)}, x::LogNumber{T1}, y::LogNumber{T2}) where {T1, T2}
    if x.sign == 0 || y.sign == 0
        return LogNumber(0, oneunit(promote_type(T1, T2)))
    end
    LogNumber(x.sign*y.sign, pm(x.exponent, y.exponent))
end

Base.:(*)(x::LogNumber, y::LogNumber) = mul_div(+, x, y)
Base.:(*)(x::LogNumber, y) = x*LogNumber(y)
Base.:(*)(x, y::LogNumber) = LogNumber(x)*y

Base.:(/)(x::LogNumber, y::LogNumber) = mul_div(-, x, y)
Base.:(/)(x::LogNumber, y) = x/LogNumber(y)
Base.:(/)(x, y::LogNumber) = LogNumber(x)/y

function ordered_plus(x::LogNumber, y::LogNumber)
    (x.sign*exp(x.exponent-y.exponent)+y.sign) * abs(y)
end
Base.:(+)(x::LogNumber, y::LogNumber) = x.exponent > y.exponent ? ordered_plus(y, x) : ordered_plus(x, y)
Base.:(+)(x, y::LogNumber) = LogNumber(x) + y
Base.:(+)(x::LogNumber, y) = x + LogNumber(y)

Base.:(-)(x::LogNumber) = LogNumber(-x.sign, x.exponent)
Base.:(-)(x::LogNumber, y::LogNumber) = x + (-y)
Base.:(-)(x, y::LogNumber) = x + (-y)
Base.:(-)(x::LogNumber, y) = x + (-y)

Base.abs(x::LogNumber) = LogNumber(1, x.exponent)

Base.sqrt(x::LogNumber) = LogNumber(sqrt(x.sign), x.exponent/2)

Base.:(^)(x::LogNumber, n) = LogNumber(x.sign^n, x.exponent*n)

end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/11/2014

module CASCoordination

export
    AbstractCASCoord,

    initialize,
    update,

    setRecord,
    getRecord,
    getAll,

    GenericCoord

using AbstractCASCoordImpl

using GenericCoordImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES