text:
mutable struct AdvectionSpeedUpMatrix

    op       :: MatrixOperators

    T_DIVx_U    :: AbstractArray{Float64, 2}
    T_DIVy_V    :: AbstractArray{Float64, 2}
    T_DIVz_W    :: AbstractArray{Float64, 2}

    U_∂x_T      :: AbstractArray{Float64, 2}
    V_∂y_T      :: AbstractArray{Float64, 2}
    W_∂z_T      :: AbstractArray{Float64, 2}

    T_∂x_U      :: AbstractArray{Float64, 2}
    T_∂y_V      :: AbstractArray{Float64, 2}
    T_∂z_W      :: AbstractArray{Float64, 2}

    U_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto U grid
    V_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto V grid
    W_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto W grid

    U_fluxmask_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto U grid with fluxmask applied
    V_fluxmask_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto V grid with fluxmask applied
    W_fluxmask_interp_T :: AbstractArray{Float64, 2}  # interpolation of T grid onto W grid with fluxmask applied

    T_mask_T       :: AbstractArray{Float64, 2}
    U_mask_U       :: AbstractArray{Float64, 2}
    V_mask_V       :: AbstractArray{Float64, 2}
    W_mask_W       :: AbstractArray{Float64, 2}
    T_bordermask_T :: AbstractArray{Float64, 2}

    U_fluxmask_U   :: AbstractArray{Float64, 2}
    V_fluxmask_V   :: AbstractArray{Float64, 2}
    W_fluxmask_W   :: AbstractArray{Float64, 2}

    T_Δvol_T :: AbstractArray{Float64, 2}

    Δx_U :: AbstractArray{Float64, 3}
    Δy_V :: AbstractArray{Float64, 3}
    Δz_W :: AbstractArray{Float64, 3}

    T_invΔz_T :: AbstractArray{Float64}
    sumΔvol_T :: AbstractArray{Float64}

    function AdvectionSpeedUpMatrix(;
        gi             :: DisplacedPoleCoordinate.GridInfo,
        Nz             :: Int64,
        Nz_av          :: AbstractArray{Int64, 2},
        mask3          :: AbstractArray{Float64, 3},
        Δz_T           :: AbstractArray{Float64, 3},      # thickness of T grid
        Δz_W           :: AbstractArray{Float64, 3},      # thickness of W grid
        nomotionmask3  :: Union{Nothing, AbstractArray{Float64, 3}} = nothing,
    )

        if size(Δz_T) != (Nz, gi.Nx, gi.Ny)
            throw(ErrorException("Δz_T is wrong size"))
        end

        if size(Δz_W) != (Nz+1, gi.Nx, gi.Ny)
            throw(ErrorException("Δz_W is wrong size"))
        end

        println("TODO: there should be no horizontal flux in the last z grid (bottom)")
        println("so as to avoid instability caused by diffusion in a thin artial gridbox.")
        println("Possible solution is to do vertical diffusion implicitly.")

        # define a converter to make 2D variable repeat in z direction for Nz times
        cvt23 = (x,) -> repeat(reshape(x, 1, size(x)...), outer=(Nz, 1, 1))
        cvt23_diagm = (x,) -> spdiagm( 0 => view(cvt23(x), :) )
        cvt3_diagm = (x,) -> spdiagm( 0 => view(x, :) )

        Nx = gi.Nx
        Ny = gi.Ny

        @time op = MatrixOperators(Nx=Nx, Ny=Ny, Nz=Nz)

        mask3_flat = view(mask3,  :)

        onV_if_unblocked_north_onT = op.V_S_T  * mask3_flat
        onV_if_unblocked_south_onT = op.V_N_T  * mask3_flat
        onU_if_unblocked_east_onT  = op.U_W_T  * mask3_flat
        onU_if_unblocked_west_onT  = op.U_E_T  * mask3_flat
        onW_if_unblocked_up_onT    = op.W_DN_T * mask3_flat
        onW_if_unblocked_dn_onT    = op.W_UP_T * mask3_flat

        V_mask = onV_if_unblocked_north_onT .* onV_if_unblocked_south_onT
        U_mask = onU_if_unblocked_east_onT  .* onU_if_unblocked_west_onT
        W_mask = onW_if_unblocked_up_onT    .* onW_if_unblocked_dn_onT

        T_mask_T = spdiagm(0 => mask3_flat)
        V_mask_V = spdiagm(0 => V_mask)
        U_mask_U = spdiagm(0 => U_mask)
        W_mask_W = spdiagm(0 => W_mask)

        T_bordermask_T = spdiagm( 0 => (
               (op.T_N_T  * mask3_flat)
            .* (op.T_S_T  * mask3_flat)
            .* (op.T_E_T  * mask3_flat)
            .* (op.T_W_T  * mask3_flat)
            .* (op.T_UP_T * mask3_flat)
            .* (op.T_DN_T * mask3_flat)
        ))

        nomotionmask3_flat = view(nomotionmask3,  :)

        onV_if_unblocked_north_onT = op.V_S_T  * nomotionmask3_flat
        onV_if_unblocked_south_onT = op.V_N_T  * nomotionmask3_flat
        onU_if_unblocked_east_onT  = op.U_W_T  * nomotionmask3_flat
        onU_if_unblocked_west_onT  = op.U_E_T  * nomotionmask3_flat
        onW_if_unblocked_up_onT    = op.W_DN_T * nomotionmask3_flat
        onW_if_unblocked_dn_onT    = op.W_UP_T * nomotionmask3_flat

        V_fluxmask = onV_if_unblocked_north_onT .* onV_if_unblocked_south_onT
        U_fluxmask = onU_if_unblocked_east_onT  .* onU_if_unblocked_west_onT
        W_fluxmask = onW_if_unblocked_up_onT    .* onW_if_unblocked_dn_onT

        V_fluxmask_V = V_mask_V * spdiagm(0 => V_fluxmask)
        U_fluxmask_U = U_mask_U * spdiagm(0 => U_fluxmask)
        W_fluxmask_W = W_mask_W * spdiagm(0 => W_fluxmask)

        # ===== [ BEGIN face area and lengths on U V ] =====

        Δx_U = gi.dx_w |> cvt23
        Δy_U = gi.DY   |> cvt23

        Δx_V = gi.DX   |> cvt23
        Δy_V = (hcat(gi.dy_s, gi.dy_n[:, end])) |> cvt23

        V_Δx_V    = (Δx_V         |>  cvt3_diagm)
        U_Δy_U    = (Δy_U         |>  cvt3_diagm)
        W_Δz_W    = (Δz_W         |>  cvt3_diagm)

        U_invΔx_U = (Δx_U.^(-1)   |>  cvt3_diagm)
        U_invΔy_U = (Δy_U.^(-1)   |>  cvt3_diagm)

        V_invΔx_V = (Δx_V.^(-1)   |>  cvt3_diagm)
        V_invΔy_V = (Δy_V.^(-1)   |>  cvt3_diagm)

        W_invΔz_W = (Δz_W.^(-1)   |>  cvt3_diagm)

        Δσ_U = Δx_U .* Δy_U
        Δσ_V = Δx_V .* Δy_V

        U_invΔσ_U = Δσ_U.^(-1) |> cvt3_diagm
        V_invΔσ_V = Δσ_V.^(-1) |> cvt3_diagm

        # ===== [ END face area and lengths on U V ] =====

        # ===== [ BEGIN face area and lengths on T ] =====

        Δx_T = gi.dx_c |> cvt23
        Δy_T = gi.dy_c |> cvt23

        Δσ_T = Δx_T .* Δy_T

        T_Δx_T = (Δx_T |>  cvt3_diagm)
        T_Δy_T = (Δy_T |>  cvt3_diagm)
        T_Δz_T = (Δz_T |>  cvt3_diagm)

        T_invΔx_T = (Δx_T.^(-1) |> cvt3_diagm)
        T_invΔy_T = (Δy_T.^(-1) |> cvt3_diagm)
        T_invΔσ_T = (Δσ_T.^(-1) |> cvt3_diagm)

        T_invΔz_T = (Δz_T.^(-1) |>  cvt3_diagm)

        T_Δvol_T = T_Δx_T * T_Δy_T * T_Δz_T

        # Δz is special. Need to clear NaN
        function clearNaN!(m)
            for i = 1:length(m.nzval)
                if isnan(m.nzval[i])
                    m.nzval[i] = 0
                end
            end
            dropzeros!(m)
        end

        clearNaN!(T_Δz_T)
        clearNaN!(W_Δz_W)
        clearNaN!(T_invΔz_T)
        clearNaN!(W_invΔz_W)

        clearNaN!(T_Δvol_T)

#        any(Δz_T .== 0) && throw(ErrorException("0!"))
#        any(Δz_W .== 0) && throw(ErrorException("00!"))

#        isnan(sum(T_invΔz_T)) && throw(ErrorException("000!"))
#        isnan(sum(W_invΔz_W)) && throw(ErrorException("0000!"))

        # ===== [ END face area and lengths on T ] =====

        # ===== [ BEG making matrix ] =====
        # MAGIC!!

        T_DIVx_U = T_mask_T * T_invΔσ_T * ( op.T_W_U - op.T_E_U    ) * U_Δy_U  ; dropzeros!(T_DIVx_U);
        T_DIVy_V = T_mask_T * T_invΔσ_T * ( op.T_S_V - op.T_N_V    ) * V_Δx_V  ; dropzeros!(T_DIVy_V);
        T_DIVz_W = T_mask_T * T_invΔz_T * ( op.T_DN_W - op.T_UP_W  )           ; dropzeros!(T_DIVz_W);

        U_∂x_T = U_mask_U * U_invΔx_U * (op.U_W_T  - op.U_E_T)                 ; dropzeros!(U_∂x_T);
        V_∂y_T = V_mask_V * V_invΔy_V * (op.V_S_T  - op.V_N_T)                 ; dropzeros!(V_∂y_T);
        W_∂z_T = W_mask_W * W_invΔz_W * (op.W_DN_T - op.W_UP_T)                ; dropzeros!(W_∂z_T);

        T_∂x_U  = T_mask_T * T_invΔx_T * ( op.T_W_U - op.T_E_U )               ; dropzeros!(T_∂x_U);
        T_∂y_V  = T_mask_T * T_invΔy_T * ( op.T_S_V - op.T_N_V )               ; dropzeros!(T_∂y_V);
        T_∂z_W  = T_mask_T * T_invΔz_T * ( op.T_DN_W - op.T_UP_W )             ; dropzeros!(T_∂z_W);

        function selfDivision(m, ones_vec)
            local wgts = m * ones_vec
            m_t = transpose(m) |> sparse
            for (i, wgt) in enumerate(wgts)
                if wgt != 0
                    _beg = m_t.colptr[i]
                    _end = m_t.colptr[i+1]-1
                    m_t.nzval[_beg:_end] ./= wgt
                end
            end

            return transpose(m_t) |> sparse
        end

        ones_T = ones(Float64, op.T_pts)

        U_interp_T = (op.U_W_T + op.U_E_T) * T_mask_T
        U_interp_T = selfDivision(U_interp_T, ones_T)

        V_interp_T = (op.V_S_T + op.V_N_T) * T_mask_T
        V_interp_T = selfDivision(V_interp_T, ones_T)

        W_interp_T = (op.W_DN_T + op.W_UP_T) * T_mask_T
        W_interp_T = selfDivision(W_interp_T, ones_T)

        U_fluxmask_interp_T = U_fluxmask_U * U_interp_T
        V_fluxmask_interp_T = V_fluxmask_V * V_interp_T
        W_fluxmask_interp_T = W_fluxmask_W * W_interp_T

        sumΔvol_T = reshape(ones(Float64, op.T_pts), 1, :) * T_mask_T * T_Δvol_T

        return new(
            op,

            T_DIVx_U,
            T_DIVy_V,
            T_DIVz_W,

            U_∂x_T,
            V_∂y_T,
            W_∂z_T,

            T_∂x_U,
            T_∂y_V,
            T_∂z_W,

            U_interp_T,
            V_interp_T,
            W_interp_T,

            U_fluxmask_interp_T,
            V_fluxmask_interp_T,
            W_fluxmask_interp_T,

            T_mask_T,
            U_mask_U,
            V_mask_V,
            W_mask_W,
            T_bordermask_T,

            U_fluxmask_U,
            V_fluxmask_V,
            W_fluxmask_W,

            T_Δvol_T,

            Δx_U,
            Δy_V,
            Δz_W,

            T_invΔz_T,
            sumΔvol_T,
        )
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
export BrickletBarometerAirPressureCallbackThreshold
struct BrickletBarometerAirPressureCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletBarometerAltitudeCallbackThreshold
struct BrickletBarometerAltitudeCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletBarometerAveraging
struct BrickletBarometerAveraging
    moving_average_pressure::Integer
    average_pressure::Integer
    average_temperature::Integer
end

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

export BrickletBarometer
"""
Measures air pressure and altitude changes
"""
mutable struct BrickletBarometer <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_air_pressure
"""
    $(SIGNATURES)

Returns the air pressure of the air pressure sensor.

If you want to get the air pressure periodically, it is recommended to use the
:cb:`Air Pressure` callback and set the period with
:func:`Set Air Pressure Callback Period`.
"""
function get_air_pressure(device::BrickletBarometer)
    return device.deviceInternal.get_air_pressure()
end

export get_altitude
"""
    $(SIGNATURES)

Returns the relative altitude of the air pressure sensor. The value is
calculated based on the difference between the current air pressure
and the reference air pressure that can be set with :func:`Set Reference Air Pressure`.

If you want to get the altitude periodically, it is recommended to use the
:cb:`Altitude` callback and set the period with
:func:`Set Altitude Callback Period`.
"""
function get_altitude(device::BrickletBarometer)
    return device.deviceInternal.get_altitude()
end

export set_air_pressure_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Air Pressure` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Air Pressure` callback is only triggered if the air pressure has
changed since the last triggering.
"""
function set_air_pressure_callback_period(device::BrickletBarometer, period)
    device.deviceInternal.set_air_pressure_callback_period(period)
end

export get_air_pressure_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Air Pressure Callback Period`.
"""
function get_air_pressure_callback_period(device::BrickletBarometer)
    return device.deviceInternal.get_air_pressure_callback_period()
end

export set_altitude_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Altitude` callback is only triggered if the altitude has changed since
the last triggering.
"""
function set_altitude_callback_period(device::BrickletBarometer, period)
    device.deviceInternal.set_altitude_callback_period(period)
end

export get_altitude_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Altitude Callback Period`.
"""
function get_altitude_callback_period(device::BrickletBarometer)
    return device.deviceInternal.get_altitude_callback_period()
end

export set_air_pressure_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Air Pressure Reached` callback.

The following options are possible:

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

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the air pressure is *outside* the min and max values"
 "'i'",    "Callback is triggered when the air pressure is *inside* the min and max values"
 "'<'",    "Callback is triggered when the air pressure is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the air pressure is greater than the min value (max is ignored)"
"""
function set_air_pressure_callback_threshold(device::BrickletBarometer, option, min, max)
    device.deviceInternal.set_air_pressure_callback_threshold(option, min, max)
end

export get_air_pressure_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Air Pressure Callback Threshold`.
"""
function get_air_pressure_callback_threshold(device::BrickletBarometer)
    return device.deviceInternal.get_air_pressure_callback_threshold()
end

export set_altitude_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Altitude Reached` callback.

The following options are possible:

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

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the altitude is *outside* the min and max values"
 "'i'",    "Callback is triggered when the altitude is *inside* the min and max values"
 "'<'",    "Callback is triggered when the altitude is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the altitude is greater than the min value (max is ignored)"
"""
function set_altitude_callback_threshold(device::BrickletBarometer, option, min, max)
    device.deviceInternal.set_altitude_callback_threshold(option, min, max)
end

export get_altitude_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Altitude Callback Threshold`.
"""
function get_altitude_callback_threshold(device::BrickletBarometer)
    return device.deviceInternal.get_altitude_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Air Pressure Reached`,
* :cb:`Altitude Reached`

are triggered, if the thresholds

* :func:`Set Air Pressure Callback Threshold`,
* :func:`Set Altitude Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletBarometer, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletBarometer)
    return device.deviceInternal.get_debounce_period()
end

export set_reference_air_pressure
"""
    $(SIGNATURES)

Sets the reference air pressure for the altitude calculation.
Setting the reference to the current air pressure results in a calculated
altitude of 0cm. Passing 0 is a shortcut for passing the current air pressure as
reference.

Well known reference values are the Q codes
`QNH <https://en.wikipedia.org/wiki/QNH>`__ and
`QFE <https://en.wikipedia.org/wiki/Mean_sea_level_pressure#Mean_sea_level_pressure>`__
used in aviation.
"""
function set_reference_air_pressure(device::BrickletBarometer, air_pressure)
    device.deviceInternal.set_reference_air_pressure(air_pressure)
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature of the air pressure sensor.

This temperature is used internally for temperature compensation of the air
pressure measurement. It is not as accurate as the temperature measured by the
:ref:`temperature_bricklet` or the :ref:`temperature_ir_bricklet`.
"""
function get_chip_temperature(device::BrickletBarometer)
    return device.deviceInternal.get_chip_temperature()
end

export get_reference_air_pressure
"""
    $(SIGNATURES)

Returns the reference air pressure as set by :func:`Set Reference Air Pressure`.
"""
function get_reference_air_pressure(device::BrickletBarometer)
    return device.deviceInternal.get_reference_air_pressure()
end

export set_averaging
"""
    $(SIGNATURES)

Sets the different averaging parameters. It is possible to set
the length of a normal averaging for the temperature and pressure,
as well as an additional length of a
`moving average <https://en.wikipedia.org/wiki/Moving_average>`__
for the pressure. The moving average is calculated from the normal
averages.  There is no moving average for the temperature.

Setting the all three parameters to 0 will turn the averaging
completely off. If the averaging is off, there is lots of noise
on the data, but the data is without delay. Thus we recommend
to turn the averaging off if the Barometer Bricklet data is
to be used for sensor fusion with other sensors.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_averaging(device::BrickletBarometer, moving_average_pressure, average_pressure, average_temperature)
    device.deviceInternal.set_averaging(moving_average_pressure, average_pressure, average_temperature)
end

export get_averaging
"""
    $(SIGNATURES)

Returns the averaging configuration as set by :func:`Set Averaging`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_averaging(device::BrickletBarometer)
    return device.deviceInternal.get_averaging()
end

export set_i2c_mode
"""
    $(SIGNATURES)

Sets the I2C mode. Possible modes are:

* 0: Fast (400kHz)
* 1: Slow (100kHz)

If you have problems with obvious outliers in the
Barometer Bricklet measurements, they may be caused by EMI issues.
In this case it may be helpful to lower the I2C speed.

It is however not recommended to lower the I2C speed in applications where
a high throughput needs to be achieved.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function set_i2c_mode(device::BrickletBarometer, mode)
    device.deviceInternal.set_i2c_mode(mode)
end

export get_i2c_mode
"""
    $(SIGNATURES)

Returns the I2C mode as set by :func:`Set I2C Mode`.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function get_i2c_mode(device::BrickletBarometer)
    return device.deviceInternal.get_i2c_mode()
end

export get_identity
"""
    $(SIGNATURES)

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

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

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

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

ext: jl


text:
struct LogitCrossEntropy end

@generated function loss(::LogitCrossEntropy, x::AbstractStrideMatrix{N,M,T}, y) where {M,N,T}
    quote
        l = zero($T)
        @avx for n ∈ axes(x,1)
            xₙ_1 = x[n,1]
            Base.Cartesian.@nexprs $(M-1) m -> xₙ_{m+1} = max(x[n,m+1], xₙ_m)
            s_1 = exp(x[n,1] - $(Symbol(:xₙ_,M)))
            Base.Cartesian.@nexprs $(M-1) m -> s_{m+1} = s_m + exp(x[n,m+1] - $(Symbol(:xₙ_,M)))
            logspxmax = log($(Symbol(:s_,M))) + $(Symbol(:xₙ_,M))
            Base.Cartesian.@nexprs $M m -> l -= (y[n] == (m-1)) * (x[n,m] - logspxmax)
        end
        l
    end
end
function logit_cross_entropy_grad_expr(loss, M, ::Type{T}) where {T}
    loopbody = quote
        xₙ_1 = x[n,1]
        Base.Cartesian.@nexprs $(M-1) m -> xₙ_{m+1} = max(x[n,m+1], xₙ_m)
        s_1 = exp(x[n,1] - $(Symbol(:xₙ_,M)))
        Base.Cartesian.@nexprs $(M-1) m -> s_{m+1} = s_m + exp(x[n,m+1] - $(Symbol(:xₙ_,M)))
        logspxmax = log($(Symbol(:s_,M))) + $(Symbol(:xₙ_,M))
        Base.Cartesian.@nexprs $M m -> δx_m = (x[n,m] - logspxmax)
        Base.Cartesian.@nexprs $M m -> y_m = y[n] == (m-1)
    end
    loss && push!(loopbody.args, :(Base.Cartesian.@nexprs $M m -> l -= y_m ? δx_m : zero($T)))
    push!(loopbody.args, :(Base.Cartesian.@nexprs $M m -> ∂x[n,m] = exp(δx_m) - y_m))
    quote
        l = zero($T)
        @avx for n ∈ axes(∂x,1)
            $loopbody
        end
        l
    end
end
@generated function gradient!(∂x, ::LogitCrossEntropy, x::AbstractStrideMatrix{N,M,T}, y) where {M,N,T}
    logit_cross_entropy_grad_expr(false, M, T)
end
@generated function loss_and_gradient!(∂x, ::LogitCrossEntropy, x::AbstractStrideMatrix{N,M,T}, y) where {M,N,T}
    logit_cross_entropy_grad_expr(true, M, T)
end


text:
int_rules_1_2_2_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.2.1*(a+b*x^2+c*x^4)^p =#
    #= Antiderivative(((~a)+(~b')*(~x)^2+(~c')*(~x)^4)^(~p),~x) :=  1/c^p*Antiderivative((b/2+c*x^2)^(2*p),x) <-- FreeQ([a,b,c,p],x) &&  EqQ(b^2-4*a*c,0) && IntegerQ(p) =#
    #= Antiderivative(1/((~a)+(~b')*(~x)^2+(~c')*(~x)^4)^(5/4),~x) :=  2*x/(3*a*(a+b*x^2+c*x^4)^(1/4)) +  x*(2*a+b*x^2)/(6*a*(a+b*x^2+c*x^4)^(5/4)) <-- FreeQ([a,b,c],x) &&  EqQ(b^2-4*a*c,0) =#
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p / (~b + 2 * ~c * (~x) ^ 2) ^ (2 * ~p)) * Antiderivative((~b + 2 * ~c * (~x) ^ 2) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p - 1 / 2))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p)) / (1 + (2 * ~c * (~x) ^ 2) / ~b) ^ (2 * FracPart(~p))) * Antiderivative((1 + (2 * ~c * (~x) ^ 2) / ~b) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p) / (4 * ~p + 1) + ((2 * ~p) / (4 * ~p + 1)) * Antiderivative((2 * ~a + ~b * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~p, 0) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (-(~x) * (((~b) ^ 2 - 2 * ~a * ~c) + ~b * ~c * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((((~b) ^ 2 - 2 * ~a * ~c) + 2 * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) + ~b * ~c * (4 * ~p + 7) * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~c / q) * Antiderivative(1 / ((~b / 2 - q / 2) + ~c * (~x) ^ 2), ~x) - (~c / q) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~a / ~c, 2)], With([r = Rt(2q - ~b / ~c, 2)], (1 / (2 * ~c * q * r)) * Antiderivative((r - ~x) / ((q - r * ~x) + (~x) ^ 2), ~x) + (1 / (2 * ~c * q * r)) * Antiderivative((r + ~x) / (q + r * ~x + (~x) ^ 2), ~x))) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], 2 * sqrt(-(~c)) * Antiderivative(1 / (sqrt(~b + q + 2 * ~c * (~x) ^ 2) * sqrt((-(~b) + q) - 2 * ~c * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~c, 0))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 4)], (((1 + q ^ 2 * (~x) ^ 2) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / (~a * (1 + q ^ 2 * (~x) ^ 2) ^ 2))) / (2 * q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(2 * atan(q * ~x), 1 / 2 - (~b * q ^ 2) / (4 * ~c))) <-- FreeQ([~a, ~b, ~c], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~c / ~a, 0) && LtQ(~b / ~a, 0)))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(-2 * ~a - (~b - q) * (~x) ^ 2) * sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / q)) / (2 * sqrt(-(~a)) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(asin(~x / sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / (2q))), (~b + q) / (2q)) <-- IntegerQ(q)) <-- FreeQ([~a, ~b, ~c], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~a, 0) && GtQ(~c, 0)))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt((2 * ~a + (~b - q) * (~x) ^ 2) / (2 * ~a + (~b + q) * (~x) ^ 2)) * sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / q)) / (2 * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) * sqrt(~a / (2 * ~a + (~b + q) * (~x) ^ 2)))) * Elliptic.F(asin(~x / sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / (2q))), (~b + q) / (2q))) <-- FreeQ([~a, ~b, ~c], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~a, 0) && GtQ(~c, 0)))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (((2 * ~a + (~b + q) * (~x) ^ 2) * sqrt((2 * ~a + (~b - q) * (~x) ^ 2) / (2 * ~a + (~b + q) * (~x) ^ 2))) / (2 * ~a * Rt((~b + q) / (2 * ~a), 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(atan(Rt((~b + q) / (2 * ~a), 2) * ~x), (2q) / (~b + q)) <-- PosQ((~b + q) / ~a) && Not(PosQ((~b - q) / ~a) && SimplerSqrtQ((~b - q) / (2 * ~a), (~b + q) / (2 * ~a)))) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (((2 * ~a + (~b - q) * (~x) ^ 2) * sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / (2 * ~a + (~b - q) * (~x) ^ 2))) / (2 * ~a * Rt((~b - q) / (2 * ~a), 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(atan(Rt((~b - q) / (2 * ~a), 2) * ~x), (-2q) / (~b - q)) <-- PosQ((~b - q) / ~a)) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(1 + ((~b + q) * (~x) ^ 2) / (2 * ~a)) * sqrt(1 + ((~b - q) * (~x) ^ 2) / (2 * ~a))) / (Rt(-((~b + q)) / (2 * ~a), 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(asin(Rt(-((~b + q)) / (2 * ~a), 2) * ~x), (~b - q) / (~b + q)) <-- NegQ((~b + q) / ~a) && Not(NegQ((~b - q) / ~a) && SimplerSqrtQ(-((~b - q)) / (2 * ~a), -((~b + q)) / (2 * ~a)))) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(1 + ((~b - q) * (~x) ^ 2) / (2 * ~a)) * sqrt(1 + ((~b + q) * (~x) ^ 2) / (2 * ~a))) / (Rt(-((~b - q)) / (2 * ~a), 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(asin(Rt(-((~b - q)) / (2 * ~a), 2) * ~x), (~b + q) / (~b - q)) <-- NegQ((~b - q) / ~a)) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 4)], (((1 + q ^ 2 * (~x) ^ 2) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / (~a * (1 + q ^ 2 * (~x) ^ 2) ^ 2))) / (2 * q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.F(2 * atan(q * ~x), 1 / 2 - (~b * q ^ 2) / (4 * ~c))) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ(~c / ~a))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Antiderivative(1 / (sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ(~c / ~a))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (((~a) ^ IntPart(~p) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p)) / ((1 + (2 * ~c * (~x) ^ 2) / (~b + q)) ^ FracPart(~p) * (1 + (2 * ~c * (~x) ^ 2) / (~b - q)) ^ FracPart(~p))) * Antiderivative((1 + (2 * ~c * (~x) ^ 2) / (~b + q)) ^ ~p * (1 + (2 * ~c * (~x) ^ 2) / (~b - q)) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~P4) ^ ~p, ~x) => With([a = Coeff(~P4, ~x, 0), b = Coeff(~P4, ~x, 1), c = Coeff(~P4, ~x, 2), d = Coeff(~P4, ~x, 3), e = Coeff(~P4, ~x, 4)], Subst(Antiderivative(SimplifyIntegrand((((a + d ^ 4 / (256 * e ^ 3)) - (b * d) / (8 * e)) + (c - (3 * d ^ 2) / (8 * e)) * (~x) ^ 2 + e * (~x) ^ 4) ^ ~p, ~x), ~x), ~x, d / (4 * e) + ~x) <-- EqQ((d ^ 3 - 4 * c * d * e) + 8 * b * e ^ 2, 0) && NeQ(d, 0)) <-- FreeQ(~p, ~x) && (PolyQ(~P4, ~x, 4) && (NeQ(~p, 2) && NeQ(~p, 3)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.1 (a+b x^2+c x^4)^p.jl


text:
using Gym
using IntervalSets
using Base.Test

function sim(env, nsteps=100, rng=Base.Random.GLOBAL_RNG)
    println("### Playing ", spec(env))
    step = 1
    done = false
    r_tot = 0.0
    obs = reset!(env)
    srand(env, 17)
    while !done && step <= nsteps
        A = action_space(env)
        action = rand(A)
        obs, rew, done, info = step!(env, action)
        r_tot += rew
        step += 1
        # render(env)
        # sleep(0.02)
    end
    # render(env, close=true)

    return r_tot
end

envs = [GymEnv("Breakout-v0"), GymEnv("CartPole-v0"), GymEnv("Pong-v0")]

for env in envs
    r = sim(env)
end

env = GymEnv("CartPole-v0")
@test action_space(env) == 0:1
O = observation_space(env)
@test O.left ≈ [-4.8, -3.40282e38, -0.418879, -3.40282e38] rtol=1e-4
@test O.right ≈ [4.8, 3.40282e38, 0.418879, 3.40282e38]  rtol=1e-4

env = GymEnv("DuplicatedInput-v0")
@test action_space(env) == (0:1,0:1,0:4)
@test observation_space(env) == 0:5


text:
function targetPos(target::AbstractTarget, tank::Tank)
rdoelprimprim = [target.ρ; 0; 0]
Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
RLos = angle_to_dcm(0, -deg2rad(tank.sight.θ), deg2rad(tank.sight.ξ), :XZY)
rdoelprim = Rhull*RLos*rdoelprimprim
rdoel = rdoelprim + [0; tank.rWY; 0]
return rdoel
end

function muzzlePos(tank::Tank)
rprimprim = [tank.canon.lw; 0; 0]
Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
Rloop = angle_to_dcm(0, -deg2rad(tank.canon.θ), 0, :XZY)
Rtoren = angle_to_dcm(0, 0, deg2rad(tank.turret.ξ), :XZY)
rprim = Rhull*Rtoren*Rloop*rprimprim
r = rprim +[0; tank.rWY; 0]
return r
end

function muzzleVel(tank::Tank)
uprimprim = [tank.canon.u₀; 0; 0]
Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
Rloop = angle_to_dcm(0, -deg2rad(tank.canon.θ), 0, :XZY)
Rtoren = angle_to_dcm(0, 0, deg2rad(tank.turret.ξ), :XZY)
uprim = Rhull*Rtoren*Rloop*uprimprim
u =uprim
return u
end

function wind(tank::Tank, w::Wind)
    w = [-w.cross*sind(tank.turret.ξ+tank.hull.χ)-w.range*cosd(tank.turret.ξ+tank.hull.χ);0.0;w.cross*cosd(tank.turret.ξ+tank.hull.χ)-w.range*sind(tank.turret.ξ+tank.hull.χ)]
end


text:
function triangulate(method::DirectLinearTransform, 𝐅::AbstractArray, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ =  𝒟
    𝐏₁, 𝐏₂ = construct(ProjectionMatrix(),𝐅)
    N = length(ℳ)
    𝒴 = Array{Point3D}(undef,N)
    for n = 1:N
        𝐦 = ℳ[n]
        𝐦ʹ = ℳʹ[n]
        eq1 = 𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:]
        eq2 = 𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:]
        eq3 = 𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:]
        eq4 = 𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:]
        𝐀 = SMatrix{4,4}(transpose(hcat(eq1,eq2,eq3,eq4)))
        # 𝐀 = [ (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ;
        #       (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ;
        #       (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ;
        #       (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ]
        # 𝐀₁ = vec2antisym(𝐦)*𝐏₁
        # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂
        # @show typeof(𝐀₁)
        # @show size(𝐀₁)
        # @show size(𝐀₂)
        # 𝐀 = vcat(𝐀₁,𝐀₂)
        # @show size(𝐀)
        # λ, f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀)))
        # @show λ
        # 𝒴[n] = Point3DH(𝑛(f))
        U,S,V = svd(𝐀)
        𝒴[n] = hom⁻¹(V[:,4])
    end
    𝒴
end

function triangulate(method::DirectLinearTransform, 𝐏₁::AbstractArray, 𝐏₂::AbstractArray, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ =  𝒟
    N = length(ℳ)
    𝒴 = Array{Point3D}(undef,N)
    for n = 1:N
        𝐦 = ℳ[n]
        𝐦ʹ = ℳʹ[n]
        eq1 = 𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:]
        eq2 = 𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:]
        eq3 = 𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:]
        eq4 = 𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:]
        𝐀 = SMatrix{4,4}(transpose(hcat(eq1,eq2,eq3,eq4)))
        # 𝐀 = [  (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ;
        #         (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ;
        #         (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ;
        #         (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ]

        # 𝐀₁ = vec2antisym(𝐦)*𝐏₁
        # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂
        # @show typeof(𝐀₁)
        # @show size(𝐀₁)
        # @show size(𝐀₂)
        # 𝐀 = vcat(𝐀₁,𝐀₂)
        # @show size(𝐀)
        # λ, f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀)))
        # @show λ
        # 𝒴[n] = Point3DH(𝑛(f))
        U,S,V = svd(𝐀)
        𝒴[n] = hom⁻¹(V[:,4])
    end
    𝒴
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(centroidal_quadruped)
height = 0.5
q1[3] += height
q1[6 + 3] += height
q1[9 + 3] += height
q1[12 + 3] += height
q1[15 + 3] += height

v1 = zeros(centroidal_quadruped.nq)
# v1= randn(centroidal_quadruped.nq)
# ## Time
h = 0.01
T = 100

# ## Simulator
s = Simulator(centroidal_quadruped, T,
    h=h)

# ## Simulate
simulate!(s, q1, v1)

using BenchmarkTools
@benchmark simulate!(s, q1, v1) setup=(s=s, q1=q1, v1=v1)

# ## Visualizer
vis = Visualizer()
open(vis)

# ## Visualize
visualize!(vis, s,
    fixed_camera=false)


text:
# TABLE LOOKUPS
const Float8toFloat32 = Float32[0.0, 0.015625, 0.03125, 0.046875, 0.0625, 0.078125, 0.09375,
    0.109375, 0.125, 0.140625, 0.15625, 0.171875, 0.1875, 0.203125, 0.21875, 0.234375, 0.25,
    0.265625, 0.28125, 0.296875, 0.3125, 0.328125, 0.34375, 0.359375, 0.375, 0.390625, 0.40625,
    0.421875, 0.4375, 0.453125, 0.46875, 0.484375, 0.5, 0.53125, 0.5625, 0.59375, 0.625, 0.65625,
    0.6875, 0.71875, 0.75, 0.78125, 0.8125, 0.84375, 0.875, 0.90625, 0.9375, 0.96875, 1.0, 1.0625,
    1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375,
    2.0, 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3.0, 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875,
    4.0, 4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25, 6.5, 6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.5,
    9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, Inf, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, -0.0, -0.015625, -0.03125, -0.046875,
    -0.0625, -0.078125, -0.09375, -0.109375, -0.125, -0.140625, -0.15625, -0.171875, -0.1875, -0.203125,
    -0.21875, -0.234375, -0.25, -0.265625, -0.28125, -0.296875, -0.3125, -0.328125, -0.34375, -0.359375,
    -0.375, -0.390625, -0.40625, -0.421875, -0.4375, -0.453125, -0.46875, -0.484375, -0.5, -0.53125,
    -0.5625, -0.59375, -0.625, -0.65625, -0.6875, -0.71875, -0.75, -0.78125, -0.8125, -0.84375, -0.875,
    -0.90625, -0.9375, -0.96875, -1.0, -1.0625, -1.125, -1.1875, -1.25, -1.3125, -1.375, -1.4375, -1.5,
    -1.5625, -1.625, -1.6875, -1.75, -1.8125, -1.875, -1.9375, -2.0, -2.125, -2.25, -2.375, -2.5, -2.625,
    -2.75, -2.875, -3.0, -3.125, -3.25, -3.375, -3.5, -3.625, -3.75, -3.875, -4.0, -4.25, -4.5, -4.75, -5.0,
    -5.25, -5.5, -5.75, -6.0, -6.25, -6.5, -6.75, -7.0, -7.25, -7.5, -7.75, -8.0, -8.5, -9.0, -9.5, -10.0,
    -10.5, -11.0, -11.5, -12.0, -12.5, -13.0, -13.5, -14.0, -14.5, -15.0, -15.5, -Inf, -NaN, -NaN, -NaN,
    -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN, -NaN]

const Float8_4toFloat32 = Float32[0.0, 0.001953125, 0.00390625, 0.005859375, 0.0078125, 0.009765625,
    0.01171875, 0.013671875, 0.015625, 0.017578125, 0.01953125, 0.021484375, 0.0234375, 0.025390625,
    0.02734375, 0.029296875, 0.03125, 0.03515625, 0.0390625, 0.04296875, 0.046875, 0.05078125, 0.0546875,
    0.05859375, 0.0625, 0.0703125, 0.078125, 0.0859375, 0.09375, 0.1015625, 0.109375, 0.1171875, 0.125,
    0.140625, 0.15625, 0.171875, 0.1875, 0.203125, 0.21875, 0.234375, 0.25, 0.28125, 0.3125, 0.34375,
    0.375, 0.40625, 0.4375, 0.46875, 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1.0,
    1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0,
    4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0,
    22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 36.0, 40.0, 44.0, 48.0, 52.0, 56.0, 60.0, 64.0, 72.0, 80.0, 88.0,
    96.0, 104.0, 112.0, 120.0, 128.0, 144.0, 160.0, 176.0, 192.0, 208.0, 224.0, 240.0, Inf, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, -0.0, -0.001953125, -0.00390625, -0.005859375, -0.0078125, -0.009765625, -0.01171875,
    -0.013671875, -0.015625, -0.017578125, -0.01953125, -0.021484375, -0.0234375, -0.025390625, -0.02734375,
    -0.029296875, -0.03125, -0.03515625, -0.0390625, -0.04296875, -0.046875, -0.05078125, -0.0546875, -0.05859375,
    -0.0625, -0.0703125, -0.078125, -0.0859375, -0.09375, -0.1015625, -0.109375, -0.1171875, -0.125, -0.140625,
    -0.15625, -0.171875, -0.1875, -0.203125, -0.21875, -0.234375, -0.25, -0.28125, -0.3125, -0.34375, -0.375,
    -0.40625, -0.4375, -0.46875, -0.5, -0.5625, -0.625, -0.6875, -0.75, -0.8125, -0.875, -0.9375, -1.0, -1.125,
    -1.25, -1.375, -1.5, -1.625, -1.75, -1.875, -2.0, -2.25, -2.5, -2.75, -3.0, -3.25, -3.5, -3.75, -4.0, -4.5,
    -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0, -18.0, -20.0,
    -22.0, -24.0, -26.0, -28.0, -30.0, -32.0, -36.0, -40.0, -44.0, -48.0, -52.0, -56.0, -60.0, -64.0, -72.0, -80.0,
    -88.0, -96.0, -104.0, -112.0, -120.0, -128.0, -144.0, -160.0, -176.0, -192.0, -208.0, -224.0, -240.0, -Inf, -NaN,
    -NaN, -NaN, -NaN, -NaN, -NaN, -NaN]

Base.Float32(x::Float8) = @inbounds Float8toFloat32[reinterpret(UInt8,x) + 0x01]
Base.Float32(x::Float8_4) = @inbounds Float8_4toFloat32[reinterpret(UInt8,x) + 0x01]

# """Old, function-based version."""
# function Float32(val::T) where {T<:AbstractFloat8}
#
#     ival = reinterpret(UInt8, val)
#
#     # seperate into sign, exponent, significand
#     sign = UInt32((ival & sign_mask(T)) >> 7)
#     exp  = UInt32((ival & exponent_mask(T)) >> n_significant_bits(T))
#     sig = UInt32(ival & significand_mask(T))
#
#     if exp == zero(UInt32)
#         if sig == zero(UInt32)          # +-0 case
#             return reinterpret(Float32,sign << 31)
#         else                            # subnormals
#             n_bit = 1
#             # first significand bit set to one, else zero, to check for size of subnormal
#             bit = first_sig_bit_mask(T)
#             while (bit & sig) == 0
#                 n_bit = n_bit + 1
#                 bit = bit >> 1
#             end
#             sign = sign << 31
#
#             # bias = 2^(n_exp-1) - 1, i.e. 127 for Float32, 15 for Float16, 3 for Float8, 7 for Float8_4
#             # difference in bias + 1 has to be added, e.g. 127-3 = 124 = 0x0000007c
#
#             exp = ((bias_difference(T)+1 - n_bit) << 23) % UInt32
#             sig = ((sig & (~bit)) << n_bit) << sig_bit_shift(T)
#             ret = sign | exp | sig
#             reinterpret(Float32,ret)
#         end
#     elseif exp == exp_bits_all_one(T)       # Inf/NaN case
#         if sig == zero(UInt32)              # Infinity
#             if sign == zero(UInt32)
#                 return Inf32
#             else
#                 return -Inf32
#             end
#         else                # NaN, preserve sign and significand (first sig bit always 1)
#             # NaN32 == reinterpret(Flaot32,0x7fc00000)
#             ret = 0x7fc00000 | (sign<<31) | (sig<<sig_bit_shift(T))
#             reinterpret(Float32,ret)
#         end
#     else
#         sign = sign << 31
#
#         # bias = 2^(n_exp-1) - 1, i.e. 127 for Float32, 15 for Float16, 3 for Float8, 7 for Float8_4
#         # difference in bias has to be added, e.g. 127-3 = 124 = 0x0000007c
#
#         exp  = (exp + bias_difference(T)) << 23
#         sig  = sig << sig_bit_shift(T)
#         ret = sign | exp | sig
#         return reinterpret(Float32, ret)
#     end
# end


text:
module SymRCM

# (C) 2020, Petr Krysl

export symrcm

using SparseArrays

"""
    adjgraph(A)

Compute the adjacency graph from a sparse matrix.

The sparse matrix `A` is assumed to be symmetric.
The results will be wrong if it isn't.

- `sortbydeg`: Should the neighbor lists be sorted by column degree? The default is
  `true`, but often results of very similar quality are obtained when this is
  set to `false` and the lists are not sorted. The second option is much
  faster, as the sorting is expensive.
"""
function adjgraph(A::SparseMatrixCSC; sortbydeg = true)
    colptr = A.colptr
    rowval = A.rowval
    ncols = length(colptr)-1
    neighbors = Vector{Vector{eltype(colptr)}}(undef, ncols)
    cdeg = diff(colptr) # the degree is colptr[j+1]-colptr[j]
    for j in 1:ncols
        cstart = colptr[j]
        jdeg = cdeg[j]
        neighbors[j] = [rowval[cstart+m-1] for m in 1:jdeg]
    end
    # All of these sorts can be done in parallel,  they are totally independent.
    # The question is when to switch over to parallel execution so as to
    # amortize the cost of starting up threads.
    if sortbydeg
        for j in 1:ncols
            sort!(neighbors[j], by = j -> cdeg[j])
        end
    end
    return neighbors
end

"""
    nodedegrees(adjgr::Vector{Vector{Int}})

Compute the degrees of the nodes in the adjacency graph.

conn = [9 1 8 4;
       1 3 2 8;
       8 2 7 5;
       2 6 7 7];
nfens = 9;
adjgr = adjgraph(conn, nfens)
nodedegrees(adjgr)

julia> degrees = node_degrees(adjgr)
9-element Array{Int64,1}:
 5
 6
 3
 3
 3
 2
 4
 7
 3
"""
function nodedegrees(adjgr::Vector{Vector{T}}) where {T}
    degrees = fill(0, length(adjgr))
    for k = 1:length(adjgr)
        degrees[k] = length(adjgr[k])
    end
    return degrees
end

"""
    symrcm(adjgr::Vector{Vector{T}}, degrees::Vector{T}) where {T}

Reverse Cuthill-McKee node-renumbering algorithm.
"""
function symrcm(adjgr::Vector{Vector{T}}, degrees::Vector{T}) where {T}
    @assert length(adjgr) == length(degrees)
    # Initialization
    n = length(adjgr)
    ndegperm = sortperm(degrees) # sorted nodal degrees
    inR = fill(false, n) # Is a node in the result list?
    inQ = fill(false, n) # Is a node in the queue?
    R = T[]
    sizehint!(R, n)
    Q = T[] # Node queue
    sizehint!(Q, n)
    while true
        P = zero(T) # Find the next node to start from
        while !isempty(ndegperm)
            i = popfirst!(ndegperm)
            if !inR[i]
                P = i
                break
            end
        end
        if P == zero(T)
            break # That was the last node
        end
        # Now we have a node to start from: put it into the result list
        push!(R, P); inR[P] = true
        empty!(Q) # empty the queue
        append!(Q, adjgr[P]); inQ[adjgr[P]] .= true # put adjacent nodes in queue
        while length(Q) >= 1
            C = popfirst!(Q) # child to put into the result list
            inQ[C] = false # make note: it is not in the queue anymore
            if !inR[C]
                push!(R, C); inR[C] = true
            end
            for i in adjgr[C] # add all adjacent nodes into the queue
                if (!inR[i]) && (!inQ[i]) # contingent on not being in result/queue
                    push!(Q, i); inQ[i] = true
                end
            end
        end
    end
    return reverse(R) # reverse the result list
end

"""
    symrcm(A::SparseMatrixCSC; sortbydeg = true)

Reverse Cuthill-McKee node-renumbering algorithm.

Compute the adjacency graph from a sparse matrix. The sparse matrix `A` is
assumed to be symmetric. The results will be wrong if it isn't.

- `sortbydeg`: Should the neighbor lists be sorted by column degree? The default is
  `true`, but often results of very similar quality are obtained when this is
  set to `false` and the lists are not sorted. The second option can be much
  faster, as the sorting is expensive when the neighbor lists are long.
"""
function symrcm(A::SparseMatrixCSC; sortbydeg = true)
    ag = adjgraph(A; sortbydeg = sortbydeg)
    nd = nodedegrees(ag)
    return symrcm(ag, nd)
end

end # module


text:
struct Categorical <: Distribution{Int} end

"""
    categorical(probs::AbstractArray{U, 1}) where {U <: Real}

Given a vector of probabilities `probs` where `sum(probs) = 1`, sample an `Int` `i` from the set {1, 2, .., `length(probs)`} with probability `probs[i]`.
"""
const categorical = Categorical()

function logpdf(::Categorical, x::Int, probs::AbstractArray{U,1}) where {U <: Real}
    (x > 0 && x <= length(probs)) ? log(probs[x]) : -Inf
end

function logpdf_grad(::Categorical, x::Int, probs::AbstractArray{U,1})  where {U <: Real}
    grad = zeros(length(probs))
    grad[x] = 1.0 / probs[x]
    (nothing, grad)
end

function random(::Categorical, probs::AbstractArray{U,1}) where {U <: Real}
    rand(Distributions.Categorical(probs))
end
is_discrete(::Categorical) = true

(::Categorical)(probs) = random(Categorical(), probs)

has_output_grad(::Categorical) = false
has_argument_grads(::Categorical) = (true,)

export categorical


text:
using Distributions, DataFrames

const θ_dist = Uniform(-1, 1)
const ϕ_dist = Uniform(0, π)

struct ChemicalShift
    σᵢₛₒ::Array{Distribution}
    Δσ::Array{Distribution}
    ησ::Array{Distribution}
    weights::Array{Float64}
end

function ChemicalShift(p::Array{Float64})
    sites = length(p) ÷ 7
    σᵢₛₒ = Array{Distribution}(undef, sites)
    Δσ = Array{Distribution}(undef, sites)
    ησ = Array{Distribution}(undef, sites)
    weights = zeros(sites)

    for i in 1:length(p)
        i % 7 in [0, 2, 4, 5, 6] && p[i] < 0 && return missing
        (isnan(p[i]) || isinf(p[i])) && return missing
        (isnothing(p[i]) || ismissing(p[i])) && return missing
    end

    for i = 1:sites
        σᵢₛₒ[i] = Normal(p[7*(i-1)+1], p[7*(i-1)+2])
        Δσ[i] = Normal(p[7*(i-1)+3], p[7*(i-1)+4])
        ησ[i] = truncated(Normal(p[7*(i-1)+5], p[7*(i-1)+6]), 0.0, 1.0)
        weights[i] = p[7*i]
    end
    weights /= sum(weights)

    return ChemicalShift(σᵢₛₒ, Δσ, ησ, weights)
end

function get_ν(
    θ::Float64,
    ϕ::Float64,
    σᵢₛₒ::Distribution,
    Δσ::Distribution,
    ησ::Distribution,
)
    rand(σᵢₛₒ) + (rand(Δσ) / 2) * (3 * θ^2 - 1 + rand(ησ) * (1-θ^2) * cos(2 * ϕ));
end

function estimate_powder_pattern(p::ChemicalShift, N::Int64)
    powder_pattern = zeros(N)
    i = 1
    for j = 1:(length(p.weights) - 1)
        to_add = floor(Int, p.weights[j] * N)
        θ = rand(θ_dist, to_add)
        ϕ = rand(ϕ_dist, to_add)
        powder_pattern[i:(i + to_add - 1)] = get_ν.(θ, ϕ, p.σᵢₛₒ[j], p.Δσ[j], p.ησ[j])
        i += to_add
    end
    θ = rand(θ_dist, N - i + 1)
    ϕ = rand(ϕ_dist, N - i + 1)
    powder_pattern[i:end] = get_ν.(θ, ϕ, p.σᵢₛₒ[end], p.Δσ[end], p.ησ[end])
    return powder_pattern
end

function get_chemical_shift_starting_values(sites::Int64)
    starting_values = zeros(7 * sites)
    starting_values[1:7:end] = rand(Uniform(-4000, 4000), sites)
    starting_values[2:7:end] = rand(Uniform(0.000001, sqrt(800)), sites)
    starting_values[3:7:end] = rand(Uniform(-4000, 4000), sites)
    starting_values[4:7:end] = rand(Uniform(0.000001, sqrt(400)), sites)
    starting_values[5:7:end] = rand(Uniform(0, 1), sites)
    starting_values[6:7:end] = rand(Uniform(0.000001, 1), sites)
    starting_values[7:7:end] = map(x -> 1 / sites, 1:sites)
    return starting_values
end

function get_output_table(minimizer::Array{Float64})
    minimizer = transform_params(push!(minimizer,1.0), ChemicalShift)
    df = DataFrame(["σᵢₛₒ" minimizer[1] minimizer[2]
                    "Δσ" minimizer[3] minimizer[4]
                    "ησ"  minimizer[5]  minimizer[6]])
    rename!(df, Symbol.(["Parameter","Mean","St. Dev."]))
    return df
end


text:
include("ellpla.jl")
# -------------------------------------------
# Contact between an ellipsoid and a triangle
# ===========================================
# References
# -----------
# [1] http://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates
# -----------------------------------------------------------------------------------------------------------------
# ellipsoid is defined by centre x, radii a, b, c, and orientation R;
# triangle is identified by points p, q, r;
# ----------------------------------------------
# returned (conpnt, depth, A, B, normal), where:
# conpnt - contact point (or [0.0, 0.0, 0.0] if no contact)
# normal - contact normal (or [0.0, 0.0, 0.0] if no contact)
# depth - positive penetration depth (or -1.0 if no contact)
# A, B - contact ellipse radii (or 0.0, 0.0 if no contact)
# ____________________________________________________________
function ellipsoid_triangle_contact(x, a, b, c, R, p, q, r)
  e0 = q-p
  e1 = r-p
  normal = cross(-e1, e0)
  normal = normal/norm(normal)
  (conpnt, depth, A, B) = ellipsoid_plane_contact(x, a, b, c, R, p, normal);
  if depth == -1.0
    normal = [0.0 0.0 0.0];
  else
    vv = conpnt-p;
    aa = dot(e0,e0)
    bb = dot(e0,e1)
    cc = dot(e1,e1)
    dd = dot(e0,vv)
    ee = dot(e1,vv)
    mm = aa*cc-bb*bb
    ss = (cc*dd-bb*ee)/mm
    tt = (aa*ee-bb*dd)/mm
    if !(ss >= 0.0 && tt >= 0.0 && ss + tt <= 1.0)
      conpnt = [0.0; 0.0; 0.0];
      normal = [0.0; 0.0; 0.0];
      depth = -1.0;
      A = 0.0;
      B = 0.0;
    end
  end
  return (conpnt, normal, depth, A, B)
end

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
# Structured matrices
using LinearAlgebra: AbstractTriangular

# Matrix wrapper types that we know are square and are thus potentially invertible. For
# these we can use simpler definitions for `/` and `\`.
const SquareMatrix{T} = Union{Diagonal{T}, AbstractTriangular{T}}

function rrule(::typeof(/), A::AbstractMatrix{<:Real}, B::T) where T<:SquareMatrix{<:Real}
    Y = A / B
    function slash_pullback(Ȳ)
        S = T.name.wrapper
        ∂A = @thunk Ȳ / B'
        ∂B = @thunk S(-Y' * (Ȳ / B'))
        return (NO_FIELDS, ∂A, ∂B)
    end
    return Y, slash_pullback
end

function rrule(::typeof(\), A::T, B::AbstractVecOrMat{<:Real}) where T<:SquareMatrix{<:Real}
    Y = A \ B
    function backslash_pullback(Ȳ)
        S = T.name.wrapper
        ∂A = @thunk S(-(A' \ Ȳ) * Y')
        ∂B = @thunk A' \ Ȳ
        return NO_FIELDS, ∂A, ∂B
    end
    return Y, backslash_pullback
end

#####
##### `Diagonal`
#####

function rrule(::Type{<:Diagonal}, d::AbstractVector)
    function Diagonal_pullback(ȳ::AbstractMatrix)
        return (NO_FIELDS, diag(ȳ))
    end
    function Diagonal_pullback(ȳ::Composite)
        # TODO: Assert about the primal type in the Composite, It should be Diagonal
        # infact it should be exactly the type of `Diagonal(d)`
        # but right now Zygote loses primal type information so we can't use it.
        # See https://github.com/FluxML/Zygote.jl/issues/603
        return (NO_FIELDS, ȳ.diag)
    end
    return Diagonal(d), Diagonal_pullback
end

function rrule(::typeof(diag), A::AbstractMatrix)
    function diag_pullback(ȳ)
        return (NO_FIELDS, Diagonal(ȳ))
    end
    return diag(A), diag_pullback
end
if VERSION ≥ v"1.3"
    function rrule(::typeof(diag), A::AbstractMatrix, k::Integer)
        function diag_pullback(ȳ)
            return (NO_FIELDS, diagm(size(A)..., k => ȳ), DoesNotExist())
        end
        return diag(A, k), diag_pullback
    end

    function rrule(::typeof(diagm), m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)
        function diagm_pullback(ȳ)
            return (NO_FIELDS, DoesNotExist(), DoesNotExist(), _diagm_back.(kv, Ref(ȳ))...)
        end
        return diagm(m, n, kv...), diagm_pullback
    end
end
function rrule(::typeof(diagm), kv::Pair{<:Integer,<:AbstractVector}...)
    function diagm_pullback(ȳ)
        return (NO_FIELDS, _diagm_back.(kv, Ref(ȳ))...)
    end
    return diagm(kv...), diagm_pullback
end

function _diagm_back(p, ȳ)
    k, v = p
    d = diag(ȳ, k)[1:length(v)] # handle if diagonal was smaller than matrix
    return Composite{typeof(p)}(second = d)
end

function rrule(::typeof(*), D::Diagonal{<:Real}, V::AbstractVector{<:Real})
    function times_pullback(Ȳ)
        return (NO_FIELDS, @thunk(Diagonal(Ȳ .* V)), @thunk(D * Ȳ))
    end
    return D * V, times_pullback
end

#####
##### `Symmetric`/`Hermitian`
#####

function frule((_, ΔA, _), T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    return T(A, uplo), T(ΔA, uplo)
end

function rrule(T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    Ω = T(A, uplo)
    function HermOrSym_pullback(ΔΩ)
        return (NO_FIELDS, _symherm_back(T, ΔΩ, Ω.uplo), DoesNotExist())
    end
    return Ω, HermOrSym_pullback
end

function frule((_, ΔA), TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    return TM(A), TM(_symherm_forward(A, ΔA))
end
function frule((_, ΔA), ::Type{Array}, A::LinearAlgebra.HermOrSym)
    return Array(A), Array(_symherm_forward(A, ΔA))
end

function rrule(TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    function Matrix_pullback(ΔΩ)
        TA = _symhermtype(A)
        T∂A = TA{eltype(ΔΩ),typeof(ΔΩ)}
        uplo = A.uplo
        ∂A = T∂A(_symherm_back(A, ΔΩ, uplo), uplo)
        return NO_FIELDS, ∂A
    end
    return TM(A), Matrix_pullback
end
rrule(::Type{Array}, A::LinearAlgebra.HermOrSym) = rrule(Matrix, A)

# Get type (Symmetric or Hermitian) from type or matrix
_symhermtype(::Type{<:Symmetric}) = Symmetric
_symhermtype(::Type{<:Hermitian}) = Hermitian
_symhermtype(A) = _symhermtype(typeof(A))

# for Ω = Matrix(A::HermOrSym), push forward ΔA to get ∂Ω
function _symherm_forward(A, ΔA)
    TA = _symhermtype(A)
    return if ΔA isa TA
        ΔA
    else
        TA{eltype(ΔA),typeof(ΔA)}(ΔA, A.uplo)
    end
end

# for Ω = HermOrSym(A, uplo), pull back ΔΩ to get ∂A
_symherm_back(::Type{<:Symmetric}, ΔΩ, uplo) = _symmetric_back(ΔΩ, uplo)
function _symherm_back(::Type{<:Hermitian}, ΔΩ::AbstractMatrix{<:Real}, uplo)
    return _symmetric_back(ΔΩ, uplo)
end
_symherm_back(::Type{<:Hermitian}, ΔΩ, uplo) = _hermitian_back(ΔΩ, uplo)
_symherm_back(Ω, ΔΩ, uplo) = _symherm_back(typeof(Ω), ΔΩ, uplo)

function _symmetric_back(ΔΩ, uplo)
    L, U, D = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), Diagonal(ΔΩ)
    return uplo == 'U' ? U .+ transpose(L) - D : L .+ transpose(U) - D
end
_symmetric_back(ΔΩ::Diagonal, uplo) = ΔΩ
_symmetric_back(ΔΩ::UpperTriangular, uplo) = Matrix(uplo == 'U' ? ΔΩ : transpose(ΔΩ))
_symmetric_back(ΔΩ::LowerTriangular, uplo) = Matrix(uplo == 'U' ? transpose(ΔΩ) : ΔΩ)

function _hermitian_back(ΔΩ, uplo)
    L, U, rD = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), real.(Diagonal(ΔΩ))
    return uplo == 'U' ? U .+ L' - rD : L .+ U' - rD
end
_hermitian_back(ΔΩ::Diagonal, uplo) = real.(ΔΩ)
function _hermitian_back(ΔΩ::LinearAlgebra.AbstractTriangular, uplo)
    ∂UL = ΔΩ .- Diagonal(_extract_imag.(diag(ΔΩ)))
    return if istriu(ΔΩ)
        return Matrix(uplo == 'U' ? ∂UL : ∂UL')
    else
        return Matrix(uplo == 'U' ? ∂UL' : ∂UL)
    end
end

#####
##### `Adjoint`
#####

# ✖✖✖TODO: Deal with complex-valued arrays as well
function rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, adjoint(ȳ))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::Type{<:Adjoint}, A::AbstractVector{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, vec(adjoint(ȳ)))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractMatrix{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, adjoint(ȳ))
    end
    return adjoint(A), adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractVector{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, vec(adjoint(ȳ)))
    end
    return adjoint(A), adjoint_pullback
end

#####
##### `Transpose`
#####

function rrule(::Type{<:Transpose}, A::AbstractMatrix)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, transpose(ȳ))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::Type{<:Transpose}, A::AbstractVector)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, vec(transpose(ȳ)))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractMatrix)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, transpose(ȳ))
    end
    return transpose(A), transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractVector)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, vec(transpose(ȳ)))
    end
    return transpose(A), transpose_pullback
end

#####
##### Triangular matrices
#####

function rrule(::Type{<:UpperTriangular}, A::AbstractMatrix)
    function UpperTriangular_pullback(ȳ)
        return (NO_FIELDS, Matrix(ȳ))
    end
    return UpperTriangular(A), UpperTriangular_pullback
end

function rrule(::Type{<:LowerTriangular}, A::AbstractMatrix)
    function LowerTriangular_pullback(ȳ)
        return (NO_FIELDS, Matrix(ȳ))
    end
    return LowerTriangular(A), LowerTriangular_pullback
end

function rrule(::typeof(triu), A::AbstractMatrix, k::Integer)
    function triu_pullback(ȳ)
        return (NO_FIELDS, triu(ȳ, k), DoesNotExist())
    end
    return triu(A, k), triu_pullback
end
function rrule(::typeof(triu), A::AbstractMatrix)
    function triu_pullback(ȳ)
        return (NO_FIELDS, triu(ȳ))
    end
    return triu(A), triu_pullback
end

function rrule(::typeof(tril), A::AbstractMatrix, k::Integer)
    function tril_pullback(ȳ)
        return (NO_FIELDS, tril(ȳ, k), DoesNotExist())
    end
    return tril(A, k), tril_pullback
end
function rrule(::typeof(tril), A::AbstractMatrix)
    function tril_pullback(ȳ)
        return (NO_FIELDS, tril(ȳ))
    end
    return tril(A), tril_pullback
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

using IntervalArithmetic
using Test

setprecision(Interval, Float64)

@testset "setdiff" begin
    x = 2..4
    y = 3..5

    d = setdiff(x, y)

    @test typeof(d) ==Vector{Interval{Float64}}
    @test length(d) ==1
    @test d ==[2..3]
    @test setdiff(y, x) ==[4..5]

    x = 2..4
    y = 2..5

    @test typeof(d) ==Vector{Interval{Float64}}
    @test length(setdiff(x, y)) ==0
    @test setdiff(y, x) ==[4..5]

    x = 2..5
    y = 3..4
    @test setdiff(x, y) ==[2..3, 4..5]

#    X = IntervalBox(2..4, 3..5)
#    Y = IntervalBox(3..5, 4..6)

    #@test setdiff(X, Y) ==IntervalBox(2..3, 3..4)
end


text:
"""
ARE fileformat

Newmaxo ASCII Region file

Example:
14 2799 925 0 0 660 0 0 0 0 0 0 0 0 0 LJ
2.799.925
2792 932
2797 943
2785.91 954.62
2780 942
2774 947
2781 957
2784 964
2784 975
2784 977
2790 974
2794 982
2790 992
2799 925

Returns dictionary with airspace description
Key     Name of Volume

Airspace struct
nb_point    Contains the number of vertices of the volume
bottom_fl   Low level of the volume in flight level (FL). Can be 0.
top_fl      High level of the volume in flight level (FL). Can be 0.
surface     Value 2 (can be negative). Can be 0.
sector_num  Value 3 (can be negative). Can be 0.
points      Matrix with all lat/lon [deg] vertices of the volume (first/last the same)

See EUROCONTROL NEST Manual sectino 9.7.6 for Are fileformat description
"""
module Are

export read

# include("utility.jl")
using ..util
using Format
using CSV
using Dates
using Navigation

struct Airspace
    nb_point::Int64
    bottom_fl::Float64
    top_fl::Float64
    surface::Float64
    sector_num::Float64
    points::Matrix{Float64}
    box::Array{Float64, 2}
end

#TODO Remove open. Just use eachline
function read(filename)
    dict = Dict{String, Airspace}()
    open(filename) do file
        j = 1
        nb_point = 3
        name = ""
        points = Matrix{Float64}(undef, 0, 2)
        bottom_fl = 0.0
        top_fl = 0.0
        surface = 0.0
        sector_num = 0.0
        for (i, ln) in enumerate(eachline(file))
            splitline = split(ln, " ")
            if length(splitline) == 2
                lat = parse(Float64, splitline[1]) / 60.0
                lon = parse(Float64, splitline[2]) / 60.0
                points = vcat(points, [lat lon])
                if j == nb_point
                    box = box_spherical_polygon(points)
                    dict[name] = Airspace(nb_point, bottom_fl, top_fl, surface,
                    sector_num, points, box)
                end
                j += 1
            elseif length(splitline) == 15
                nb_point = parse(Int64, splitline[1])
                bottom_fl = parse(Float64, splitline[5])
                top_fl = parse(Float64, splitline[6])
                surface = parse(Float64, splitline[7])
                sector_num = parse(Float64, splitline[8])
                name = splitline[15]
                points = Matrix{Float64}(undef, 0, 2)
                j = 1
            else
            end

        end
    end
    return dict
end

# function isinside(are::Airspace, point::Point)
#     # test if inside box
#     if isinsidebox(are, point)
#         #test if inside polygon
#         return NaN
#     else
#         return false
#     end
# end

function box_spherical_polygon(points::Array{Float64,2})
    box = Array{Float64, 2}(undef, 2, 2)
    lats_max = Array{Float64,1}()
    lats_min = Array{Float64,1}()
    previous_point = Array{Float64,1}()
    for point in eachrow(points)
        if isempty(previous_point)
            previous_point = point
        else
            # pos₁ = Point(previous_point[1], previous_point[2])
            # pos₂ = Point(point[1], point[2])
            pos₁ = Point_deg(previous_point[1], previous_point[2])
            pos₂ = Point_deg(point[1], point[2])
            append!(lats_max, maximum(latitude_options(pos₁, pos₂)))
            append!(lats_min, minimum(latitude_options(pos₁, pos₂)))
            previous_point = point
        end
    end
    box[1,1] = maximum(lats_max)
    box[2,1] = minimum(lats_min)

    # Limiting longitudes
    box[1,2] = maximum(points[:,2])
    box[2,2] = minimum(points[:,2])
    return box
end

function latitude_options(pos₁::Point_deg, pos₂::Point_deg)
    bearing = Navigation.bearing(pos₁, pos₂)
    closest_point_deg = Navigation.closest_point_to_pole(pos₁, bearing)
    if Navigation.distance(pos₁, pos₂) > Navigation.distance(pos₁,
        closest_point_deg)
        return pos₁.ϕ, pos₂.ϕ, closest_point_deg.ϕ
    else
        return pos₁.ϕ, pos₂.ϕ
    end
end

# function isinsidebox(are::Airspace, point::Point)
#     return are.box[2,1] < point.ϕ < are.box[1,1] && are.box[2,2] < point.λ < are.box[1,2]
# end

end # module


text:
type GroupedDataFrameCode
  parent::AbstractDataFrame
  cols::Vector           # columns used for sorting
  idx::Vector{Int}       # indexing vector when sorted by the given columns
  starts::Vector{Int}    # starts of groups
  ends::Vector{Int}      # ends of groups
  levelidx::Vector{Int}  # integer coded groups
end

function _markuniqueofsorted(x)
  n = length(x)
  res = ones(Bool, n)

  xprev = x[1]
  for i=1:(n-1)
    res[i] = x[i+1] != x[i]
  end
  return res
end

pda(da::DataArray, order) = PooledDataArray(da.data, order, da.na, Int)
pda(pda::PooledDataArray, order) = pda

@doc """
taken from the DataFrames package and modified to ensure that groups are encoded in an
order given by the `colorders` argument.
""" ->
function groupby{T}(df::AbstractDataFrame, cols::Vector{T}, colorders::Dict{T,Any})
  ## a subset of Wes McKinney's algorithm here:
  ##     http://wesmckinney.com/blog/?p=489

  ncols = length(cols)
  # use the pool trick to get a set of integer references for each unique item
  da = df[cols[ncols]]

  if haskey(colorders, cols[ncols])
    order = colorders[cols[ncols]]
    dv = pda(da, order)
  else
    dv = PooledDataArray(da)
  end
  # if there are NAs, add 1 to the refs to avoid underflows in x later
  dv_has_nas = (findfirst(dv.refs, 0) > 0 ? 1 : 0)
  x = copy(dv.refs) .+ dv_has_nas
  # also compute the number of groups, which is the product of the set lengths
  ngroups = length(dv.pool) + dv_has_nas
  # if there's more than 1 column, do roughly the same thing repeatedly
  for j = (ncols - 1):-1:1
    da = df[cols[j]]
    if haskey(colorders, cols[j])
      order = colorders[cols[j]]
      dv = pda(da, order)
    else
      dv = PooledDataArray(da)
    end
    dv_has_nas = (findfirst(dv.refs, 0) > 0 ? 1 : 0)
    for i = 1:nrow(df)
      x[i] += (dv.refs[i] + dv_has_nas- 1) * ngroups
    end
    ngroups = ngroups * (length(dv.pool) + dv_has_nas)
    # TODO if ngroups is really big, shrink it
  end
  (idx, starts) = DataArrays.groupsort_indexer(x, ngroups)
  # Remove zero-length groupings
  nonzerogroupings = _markuniqueofsorted(starts)
  starts = starts[nonzerogroupings]

  ends = [starts[2:end] .- 1] # change
  levelidx = Int[x[idx[endi]] for endi in ends]
#   println("idx:\n$idx\nends:\n$ends\nlevelidx:\n$levelidx")
  GroupedDataFrameCode(df, cols, idx, starts[1:end-1], ends, levelidx)
end

Base.start(gd::GroupedDataFrameCode) = 1
Base.next(gd::GroupedDataFrameCode, state::Int) =
    ((gd.levelidx[state], sub(gd.parent, gd.idx[gd.starts[state]:gd.ends[state]])),
     state + 1)
Base.done(gd::GroupedDataFrameCode, state::Int) = state > length(gd.starts)
Base.length(gd::GroupedDataFrameCode) = length(gd.starts)
Base.endof(gd::GroupedDataFrameCode) = length(gd.starts)
Base.first(gd::GroupedDataFrameCode) = gd[1]
Base.last(gd::GroupedDataFrameCode) = gd[end]

max_stars_repo_name: cpritcha/Pivots

max_issues_repo_name: cpritcha/Pivots

max_forks_repo_name: cpritcha/Pivots


text:
module ProjectEulerDownload

using JSON
using Cascadia
using Gumbo
using HTTP

const PED = ProjectEulerDownload

function configure(euler_dir::String)
    dir = pkgdir(PED)
    d = Dict("root" => euler_dir)
    write(joinpath(dir, "config.json"), JSON.json(d))
end

function download(problem_id::Int)
    pdir = pkgdir(PED)
    config_path = joinpath(pdir, "config.json")
    if !isfile(config_path)
        error("Please run PED.configure(\"Path to root directory\") first")
        return
    end
    root_dir = JSON.parsefile(config_path)["root"]
    cd(root_dir)

    url = "https://projecteuler.net/problem=$problem_id";
    res = HTTP.get(url);
    body = String(res.body);
    html = parsehtml(body);

    title = text(eachmatch(sel"h2", html.root)[1][1]);
    folder_title = replace(title, " "=> "_")
    lid = lpad(problem_id, 3, '0')
    folder = "$(lid)_$folder_title"
    # create the new folder and change to there
    mkdir(folder)
    cd(folder)

    # CREATE SEVERAL FILES
    # README.md
    problem_content = eachmatch(sel".problem_content", html.root)[1]
    str_content = string(problem_content)
    str_content = replace(str_content, "</p>" =>"\n")
    str_content = replace(str_content, r"<dfn(>|.*?>)(?<text>.*?)</dfn>" => s"**\g<text>**")
    str_content = replace(str_content, r"<(?<tag>[/]?(sup|sub))>" => s"<_\g<tag>>")
    str_content = replace(str_content, r"<[^_].*?>" => "")
    str_content = replace(str_content, r"<_(?<tag>.*?)>" => s"<\g<tag>>")

    readme = "# $title\n$str_content\n\nSource: [Project Euler #$problem_id]($url)"
    write("README.md", readme)
    println("Checkout the readme!")
    println(joinpath(pwd(), "README.md"))

    # solution.jl
    write("solution.jl", "")

    # runtests.jl
    test_file = "using Test\n\ninclude(\"solution.jl\")\n\n"
    test_file *= "@testset \"example\" begin\n"
    test_file *= "\t@test\n"
    test_file *= "end"
    write("runtests.jl", test_file);

    println("Happy coding!")
end

end


text:
# colors chosen by according to https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106
# as proposed by @tpoisot in https://github.com/JuliaPlots/Plots.jl/issues/1144
const wong_palette = [
    RGB(([230, 159,   0] / 255)...), # orange
    RGB(([ 86, 180, 233] / 255)...), # sky blue
    RGB(([  0, 158, 115] / 255)...), # blueish green
    RGB(([240, 228,  66] / 255)...), # yellow
    RGB(([  0, 114, 178] / 255)...), # blue
    RGB(([213,  94,   0] / 255)...), # vermillion
    RGB(([204, 121, 167] / 255)...), # reddish purple
    ]

_themes[:wong] = PlotTheme(
    palette = expand_palette(colorant"white", wong_palette; lchoices = [57], cchoices = [100]),
    gradient = cgrad(:viridis).colors,
)

_themes[:wong2] = PlotTheme(
    palette = expand_palette(colorant"white", [RGB(0,0,0); wong_palette]; lchoices = [57], cchoices = [100]),
    gradient = cgrad(:viridis).colors,
)


text:
using Base.Threads: lock, unlock, SpinLock
const _lock = SpinLock()

const CacheComplexityThres = 32
const CacheComplexityDepthFactor = 2
const RevShiftBits = 2
const ShiftBits = sizeof(UInt) * 8 - RevShiftBits

const HAS_FLAGS = UInt(0b11) << ShiftBits
struct BitsValue{L}
    value::L
end

const FLAGS = (Any = UInt(0), Expr = UInt(1), Call = UInt(2))

const CompactExpr{Args} = Tuple{Symbol,Args}

const RefValPool = Any[]
const RefValIndex = Base.IdDict{Any,UInt}()

const ExprPool = CompactExpr[]
const ExprIndex = Dict{CompactExpr,UInt}()

struct Call{F,Args<:Tuple}
    f::F
    args::Args
end

const CallPool = Call[]
const CallIndex = Dict{Call,UInt}()

struct ExprMeta
    complexity::Int
    inner_lns::Vector{Tuple{Int,ExprMeta}}
    lns::Vector{Tuple{Int,LineNumberNode}}
end

struct SimpleMeta
    complexity::Int
end

SimpleMeta(meta::ExprMeta) = SimpleMeta(meta.complexity)
SimpleMeta(meta::SimpleMeta) = meta

const leaf_meta = SimpleMeta(0)

@nospecialize

function _assign_id(pool::AbstractVector, flag::UInt)
    i = length(pool)
    if iszero(i & HAS_FLAGS)
        return UInt(i) | (flag << ShiftBits)
    else
        throw(OverflowError("pool too long"))
    end
end

@static if VERSION < v"1.1"
    function _get!(default::Function, d::Base.IdDict{K,V}, @nospecialize(key)) where {K,V}
        val = get(d, key, Base.secret_table_token)
        if val === Base.secret_table_token
            val = default()
            if !isa(val, V)
                val = convert(V, val)::V
            end
            setindex!(d, val, key)
            return val
        else
            return val::V
        end
    end
end

_get!(f, d, key) = get!(f, d, key)

function _get_id!(pool, index, val, flag)
    lock(_lock)
    try
        _get!(index, val) do
            id = _assign_id(pool, flag)
            push!(pool, val)
            index[val] = id
            id
        end
    finally
        unlock(_lock)
    end
end

function _lookup!(pool, idx::Integer)
    lock(_lock)
    try
        pool[idx]
    finally
        unlock(_lock)
    end
end

function compress_impl!(val::Ptr{T}) where {T}
    leaf_meta, Call(Constructor{Ptr{T}}(), tuple(BitsValue(UInt(val))))
end

function compress_impl!(val)
    if isbits(val)
        leaf_meta, BitsValue(val)
    else
        leaf_meta, _get_id!(RefValPool, RefValIndex, val, FLAGS.Any)
    end
end

compress_impl!(val::Symbol) = (leaf_meta, val)
function compress_impl!(ex::Expr)
    args = Any[]
    inner_lns = Tuple{Int,ExprMeta}[]
    lns = Tuple{Int,LineNumberNode}[]
    i = 0
    maxcomplexity = 1
    for each in ex.args
        if each isa LineNumberNode
            push!(lns, (i, each))
        else
            i += 1
            meta, a = compress_impl!(each)
            push!(args, a)
            if meta !== leaf_meta # leaf
                if meta isa SimpleMeta || isempty(meta.inner_lns) && isempty(meta.lns)
                else
                    push!(inner_lns, (i, meta))
                end
                maxcomplexity =
                    max(meta.complexity + CacheComplexityDepthFactor, maxcomplexity)
            end
        end
    end

    base = (ex.head, Tuple(args))
    if maxcomplexity < CacheComplexityThres
        base = _get_id!(ExprPool, ExprIndex, base, FLAGS.Expr)
    end
    ExprMeta(maxcomplexity, inner_lns, lns), base
end

function ln_to_tuple(ln::LineNumberNode)
    (ln.line, Symbol(ln.file))
end

function meta_to_tuple(meta::ExprMeta)
    Tuple((i, meta_to_tuple(m)) for (i, m) in meta.inner_lns),
    Tuple((i, ln_to_tuple(l)) for (i, l) in meta.lns)
end

function compress(val)
    meta, encoded = compress_impl!(val)
    meta isa SimpleMeta && return encoded
    (meta_to_tuple(meta), encoded)
end

decompress_impl(encoded::Symbol, ::Tuple) = encoded
decompress_impl(encoded::BitsValue, ::Tuple) = encoded.value

function decompress_impl(encoded::UInt, meta::Tuple)
    flag = (encoded >> ShiftBits)
    if flag == FLAGS.Any
        return _lookup!(RefValPool, (encoded << RevShiftBits >> RevShiftBits) + 1)
    elseif flag == FLAGS.Expr
        encoded = _lookup!(ExprPool, (encoded << RevShiftBits >> RevShiftBits) + 1)
        return decompress_impl(encoded, meta)
    elseif flag == FLAGS.Call
        encoded = _lookup!(CallPool, (encoded << RevShiftBits >> RevShiftBits) + 1)
        return decompress_impl(encoded, meta)
    else
        error("invalid flag $flag")
    end
end

function get_from_tuple(xs::Tuple, key, default)
    for (k, v) in xs
        k == key && return v
    end
    return default
end

function decompress_impl(encoded::Call, ::Tuple)
    f = encoded.f
    args = decompress.(encoded.args)
    f(args...)
end

function decompress_impl(encoded::Tuple, meta::Tuple)
    inner_lns, lns = meta
    head, args = encoded
    ex = Expr(head)
    j = 1
    for i in eachindex(args)
        while j <= length(lns) && lns[j][1] < i
            line, file = lns[j][2]
            push!(ex.args, LineNumberNode(line, file))
            j += 1
        end
        default_meta = ((), ())
        push!(ex.args, decompress_impl(args[i], get_from_tuple(inner_lns, i, default_meta)))
    end
    ex
end

function decompress(encoded::Tuple)
    meta, encoded = encoded
    decompress_impl(encoded, meta)
end

function decompress(encoded)
    default_meta = ((), ())
    decompress_impl(encoded, default_meta)
end

struct Constructor{A} end
function (::Constructor{A})(args...) where {A}
    A(args...)
end

function compress_impl!(a::QuoteNode)
    meta, sub = compress_impl!(a.value)
    meta = SimpleMeta(meta.complexity)
    encoded = Call(Constructor{QuoteNode}(), tuple(sub))
    if meta.complexity < CacheComplexityThres
        encoded = _get_id!(CallPool, CallIndex, encoded, FLAGS.Call)
    end
    return SimpleMeta(meta), encoded
end

function compress_impl!(a::GlobalRef)
    _, mod_encoded = compress_impl!(a.mod)
    _, name_encoded = compress_impl!(a.mod)
    meta = SimpleMeta(1)
    encoded = Call(Constructor{GlobalRef}(), tuple(mod_encoded, name_encoded))
    return meta, encoded
end

function may_cache_call(meta::Union{ExprMeta,SimpleMeta}, call::Call)
    if meta.complexity < CacheComplexityThres
        SimpleMeta(meta), _get_id!(CallPool, CallIndex, call, FLAGS.Call)
    else
        SimpleMeta(meta), call
    end
end


text:
abstract Projection

type AffineTransformation
    t :: Array{Float64,1}
    m :: Array{Float64,2}
end

function AffineTransformation(x0, x1, y0, y1, dest::BoundingBox)
    sx = width(dest) / (x1 - x0)
    sy = height(dest) / (y1 - y0)
    p = lowerleft(dest)
    tx = p.x - sx * x0
    ty = p.y - sy * y0
    t = [tx, ty]
    m = diagm([sx, sy])
    AffineTransformation(t, m)
end

function project(self::AffineTransformation, x::Real, y::Real)
    #self.m*[x,y] + self.t
    u = self.t[1] + self.m[1,1] * x #+ self.m[1,2] * y
    v = self.t[2] + self.m[2,2] * y #+ self.m[2,1] * x
    u, v
end

project(proj::Projection, p::Point) = Point(project(proj, p.x, p.y)...)

function project(self::AffineTransformation, x::Vector, y::Vector)
    p = self.t[1] + self.m[1,1] * x #+ self.m[1,2] * y
    q = self.t[2] + self.m[2,2] * y #+ self.m[2,1] * x
    return p, q
end

project(self::AffineTransformation, x::AbstractArray, y::AbstractArray) =
    project(self, collect(x), collect(y))

#function compose(self::AffineTransformation, other::AffineTransformation)
#    self.t = call(other.t[1], other.t[2])
#    self.m = self.m * other.m
#end

type PlotGeometry <: Projection
    dest_bbox::BoundingBox
    xlog::Bool
    ylog::Bool
    aff::AffineTransformation
    xflipped::Bool
    yflipped::Bool

    function PlotGeometry(x0, x1, y0, y1, dest::BoundingBox, xlog, ylog)
        if xlog
            x0 = log10(x0)
            x1 = log10(x1)
        end
        if ylog
            y0 = log10(y0)
            y1 = log10(y1)
        end
        new(dest, xlog, ylog, AffineTransformation(x0,x1,y0,y1,dest), x0 > x1, y0 > y1)
    end

    PlotGeometry(x0, x1, y0, y1, dest) = PlotGeometry(x0, x1, y0, y1, dest, false, false)
end

function project(self::PlotGeometry, x, y)
    u, v = x, y
    if self.xlog
        u = log10(x)
    end
    if self.ylog
        v = log10(y)
    end
    return project(self.aff, u, v)
end


text:
# To allow not having installed the package
push!(LOAD_PATH, "../src/")

using Documenter, Tick

# To generate site with mkdocs (see make.jl config):
# pip install mkdocs python-markdown-math
# mkdocs serve

# makedocs() # then generate site with mkdocs

makedocs(
    modules = [Tick],
    format = :html,
    sitename = "Tick.jl",
    authors = "Ludovic Coquelle",
    clean = true,
    pages = Any[
        "Index"=>"index.md",
        "Tutorial"=>"tutorial.md",
        "Reference"=>"reference.md",
    ]
)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Here we implement a minimal dimensionful type Furlong, which is used
# to test dimensional correctness of various functions in Base.  Furlong
# is exported by the TestHelpers module.

# represents a quantity in furlongs^p
struct Furlong{p,T<:Number} <: Number
    val::T
    Furlong{p,T}(v::Number) where {p,T} = new(v)
    Furlong{p,T}(x::Furlong{p}) where {p,T} = new(x.val)
end
Furlong{T<:Number}(x::T) = Furlong{1,T}(x)
(::Type{T}){p,T}(x::Furlong{p,T}) = x.val
(::Type{Furlong{p}}){p}(v::Number) = Furlong{p,typeof(v)}(v)
Base.convert{T,p,S}(::Type{Furlong{p,T}}, x::Furlong{p,S}) = Furlong{p,T}(convert(T,x.val))
Base.convert{T}(::Type{Furlong{0,T}}, x::Furlong{0}) = Furlong{0,T}(convert(T, x.val))
Base.convert{T<:Number}(::Type{T}, x::Furlong{0}) = convert(T, x.val)
Base.convert{T}(::Type{Furlong{0,T}}, x::Number) = Furlong{0,T}(convert(T, x))

Base.promote_type{p,T,S}(::Type{Furlong{p,T}}, ::Type{Furlong{p,S}}) =
    (Base.@_pure_meta; Furlong{p,promote_type(T,S)})

Base.one{p,T}(x::Furlong{p,T}) = one(T)
Base.one{p,T}(::Type{Furlong{p,T}}) = one(T)
Base.zero{p,T}(x::Furlong{p,T}) = Furlong{p,T}(zero(T))
Base.zero{p,T}(::Type{Furlong{p,T}}) = Furlong{p,T}(zero(T))
Base.iszero(x::Furlong) = iszero(x.val)

# convert Furlong exponent p to a canonical form.  This
# is not type stable, but it doesn't matter since it is used
# at compile time (in generated functions), not runtime
canonical_p(p) = isinteger(p) ? Int(p) : Rational{Int}(p)

Base.abs{p}(x::Furlong{p}) = Furlong{p}(abs(x.val))
@generated Base.abs2{p}(x::Furlong{p}) = :(Furlong{$(canonical_p(2p))}(abs2(x.val)))
@generated Base.inv{p}(x::Furlong{p}) = :(Furlong{$(canonical_p(-p))}(inv(x.val)))
Base.sylvester(a::Furlong,b::Furlong,c::Furlong) = -c / (a + b)

for f in (:isfinite, :isnan, :isreal)
    @eval Base.$f(x::Furlong) = $f(x.val)
end
for f in (:real,:imag,:complex,:middle,:+,:-)
    @eval Base.$f{p}(x::Furlong{p}) = Furlong{p}($f(x.val))
end

import Base: +, -, ==, !=, <, <=, isless, isequal, *, /, //, div, rem, mod, ^,
             middle, hypot
for op in (:+, :-, :middle, :hypot)
    @eval function $op{p}(x::Furlong{p}, y::Furlong{p})
        v = $op(x.val, y.val)
        Furlong{p}(v)
    end
end
for op in (:(==), :(!=), :<, :<=, :isless, :isequal)
    @eval $op{p}(x::Furlong{p}, y::Furlong{p}) = $op(x.val, y.val)
end
# generated functions to allow type inference of the value of the exponent:
for (f,op) in ((:_plus,:+),(:_minus,:-),(:_times,:*),(:_div,://))
    @eval @generated function $f{T,p,q}(v::T, ::Furlong{p}, ::Union{Furlong{q},Type{Val{q}}})
        s = $op(p, q)
        :(Furlong{$(canonical_p(s)),$T}(v))
    end
end
for (op,eop) in ((:*, :_plus), (:/, :_minus), (://, :_minus), (:div, :_minus))
    @eval begin
        $op{p,q}(x::Furlong{p}, y::Furlong{q}) =
            $eop($op(x.val, y.val),x,y)
        $op{p,S<:Number}(x::Furlong{p}, y::S) = $op(x,Furlong{0,S}(y))
        $op{p,S<:Number}(x::S, y::Furlong{p}) = $op(Furlong{0,S}(x),y)
    end
end
for op in (:rem, :mod)
    @eval begin
        $op{p}(x::Furlong{p}, y::Furlong) = Furlong{p}($op(x.val, y.val))
        $op{p}(x::Furlong{p}, y::Number) = Furlong{p}($op(x.val, y))
    end
end
Base.sqrt(x::Furlong) = _div(sqrt(x.val), x, Val{2})

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
function sort_lists!(main_list, other_lists...; cur_func::Function=identity)
  cur_indices = sortperm(main_list, by=cur_func)

  cur_list = collect(main_list)[cur_indices]
  _arrange_lists!(main_list, other_lists, cur_indices)
  cur_list
end

function shuffle_lists!(main_list, other_lists...)
  cur_indices = shuffle(1:length(main_list))

  cur_list = main_list[cur_indices]
  _arrange_lists!(main_list, other_lists, cur_indices)
  cur_list
end

function _arrange_lists!(main_list, other_lists, cur_indices)
  cur_lists = Any[ isa(main_list, Vector) ? main_list : collect(main_list) ]
  isempty(other_lists) || append!(cur_lists, other_lists)

  for cur_list in cur_lists
      cur_list .= cur_list[cur_indices]
  end
  main_list
end

function bisect_reorder!(cur_vector::Vector)
  cur_indices = breadth_indices(length(cur_vector))
  cur_vector .= cur_vector[cur_indices]
  cur_vector
end

function collect_bisection_indexes!(ix, level, a, b)
    a ≤ b || return
    mid = (a+b) ÷ 2
    push!(ix, level => mid)
    collect_bisection_indexes!(ix, level + 1, a, mid - 1)
    collect_bisection_indexes!(ix, level + 1, mid + 1, b)
end

function breadth_indices(n)
    ix = Vector{Pair{Int,Int}}()
    sizehint!(ix, n)
    collect_bisection_indexes!(ix, 0, 1, n)
    last.(sort!(ix, by = first))
end

function out_in_reorder!(cur_vector::Vector)
  cur_indices = centering_indices(length(cur_vector))
  cur_vector .= cur_vector[cur_indices]
  cur_vector
end

function centering_indices(num_points::Number)
  cur_indices = []

  cur_half_point = Int(floor(num_points/2))

  for cur_index in 1:cur_half_point
    push!(cur_indices, cur_index)
    push!(cur_indices, num_points-(cur_index-1))
  end
  isodd(num_points) && push!(cur_indices, cur_half_point+1)

  cur_indices
end

function in_out_reorder!(cur_vector::Vector)
  cur_indices = dividing_indices(length(cur_vector))
  cur_vector .= cur_vector[cur_indices]
  cur_vector
end

function dividing_indices(num_points::Number)
  cur_indices = []

  cur_half_point = Int(floor(num_points/2))

  isodd(num_points) && push!(cur_indices, cur_half_point+1)
  for cur_index in reverse(1:cur_half_point)
    push!(cur_indices, cur_index)
    push!(cur_indices, num_points-(cur_index-1))
  end

  cur_indices
end


text:
using Documenter, RLEVectors

makedocs()

deploydocs(
           repo = "github.com/phaverty/RLEVectors.jl.git"
           )


text:
# Unit Tests
using Test
using MLKernels

using MLKernels:
    BaseFunction,
        InnerProduct,
            ScalarProduct,
        PreMetric,
            Metric,
                SquaredEuclidean

using SpecialFunctions:
    besselk, gamma

import LinearAlgebra
import Statistics

const FloatingPointTypes = (Float32, Float64)
const IntegerTypes = (Int32, Int64)
const MLK = MLKernels

include("reference.jl")

@testset "Testing utility functions" begin
    include("utils.jl")
end

@testset "Testing BaseFunction types" begin
    include("basefunctions.jl")
end

@testset "Testing base function evaluation" begin
    include("basematrix.jl")
end

@testset "Testing Kernel types" begin
    include("kernelfunctions.jl")
end

@testset "Testing kernel function evaluation" begin
    include("kernelmatrix.jl")
end

@testset "Testing nystrom" begin
    include("nystrom.jl")
end


text:
using DigilentWaveForms
dev = devicelist()[1]
DigilentWaveForms.open!(dev)
#DigilentWaveForms.open(dev)

DigilentWaveForms.AnalogIOEnable!(dev)
DigilentWaveForms.enablePositivePowerSupply!(dev)
DigilentWaveForms.setPositivePowerSupply!(dev, 4)
for x in 1:5
    for i in 2.0:0.01:4.5
        DigilentWaveForms.setPositivePowerSupply!(dev, i)
        sleep(0.01)
    end
end

DigilentWaveForms.AnalogOutChannelCount(dev)

DigilentWaveForms.AnalogOutEnableChannel!(dev, 0)
DigilentWaveForms.AnalogOutEnableChannel!(dev, 1)
DigilentWaveForms.AnalogOutDCMode!(dev, 0)
DigilentWaveForms.AnalogOutDCMode!(dev, 1)
DigilentWaveForms.AnalogOutStart!(dev, 0)
DigilentWaveForms.AnalogOutStart!(dev, 1)

for i in -5:0.1:5, j in -5:0.1:5
    DigilentWaveForms.AnalogOutDCVal!(dev, 0, i)
    DigilentWaveForms.AnalogOutDCVal!(dev, 1, j)
    sleep(0.01)
end

DigilentWaveForms.close!(dev)

using BenchmarkTools
@btime DigilentWaveForms.setPositivePowerSupply!(dev, 1.0)


text:
#=
Test utils operations:
=#

@testset "Utils" begin
    @testset "Bounding box manipulations" begin
        box = rand(10, 4)
        gt_box = rand(10, 4)
        new = box_refinement_graph(box, gt_box)
        @test (10, 4) == size(new)

        mask = rand(28, 28, 10)
        bboxes = extract_bboxes(mask)
        @test (10, 4) == size(bboxes)

        boxes1 = [0 0 2 2; 1 0 3 2]
        boxes2 = [0 0 1 1; 0 0 2 1; 0 0 4 4]

        overlaps = compute_overlaps_bbox(boxes1, boxes2)
        @test overlaps == [.25 .5 .25; 0 0.2 .25]

        windows = [0 0 2 2; 0 0 1 2; 0 0 2 2]

        clipped_boxes = clip_boxes(boxes2, windows)
        @test [0 0 1 1; 0 0 1 1; 0 0 2 2] == clipped_boxes

        boxes = clip_to_window([0 0 1 1], boxes2)
        @test [0 0 1 1; 0 0 1 1; 0 0 1 1] == boxes
    end

    @testset "Image manipulations" begin
        image = rand(1000, 1000, 3)

        resized_image = resize_image(image, min_dim = 800, max_dim = 1024)
        @test size(resized_image[1]) == (1024, 1024, 3)

        MEAN_PIXEL = [123.7f0 / 255.0f0, 116.8f0 / 255.0f0, 103.9f0 / 255.0f0]
        molded_image = mold_image(image)

        molded_image[:, :, 1] .+= MEAN_PIXEL[1]
        molded_image[:, :, 2] .+= MEAN_PIXEL[2]
        molded_image[:, :, 3] .+= MEAN_PIXEL[3]

        @test all(image[:, :, 1] .== molded_image[:, :, 1])
        @test all(image[:, :, 2] .== molded_image[:, :, 2])
        @test all(image[:, :, 3] .== molded_image[:, :, 3])

        images = [rand(1000, 1000, 3), rand(1000, 1000, 3), rand(1000, 1000, 3)]

        info = mold_inputs(images)
        @test (1024, 1024, 3, 3) == size(info[1])
        @test (1, 3) == size(info[2])
        @test (1, 3) == size(info[3])
    end

    @testset "Loss functions" begin
        y = 1:10
        ŷ = 0.5:1:9.5
    end

    @testset "Resnet functions" begin
        for model_name in ("imagenet-resnet-50-dag", "imagenet-resnet-101-dag", "imagenet-resnet-152-dag")
            model = matconvnet(model_name, true)
            @test all(keys(model) .== ["meta", "params", "vars", "layers"])
        end
    end
end