text:
"""model definition of snake2d"""

const STEP_DIST_SNAKE2D = 0.05

struct StateSnake2D <: AbsState
    x::Float64
    y::Float64
    theta1::Float64
    theta2::Float64
    theta3::Float64
    theta4::Float64
end

function get_mid_status(p::StateSnake2D, q::StateSnake2D)::StateSnake2D
    return StateSnake2D(
        (p.x + q.x) / 2,
        (p.y + q.y) / 2,
        diff_angles(q.theta1, p.theta1) / 2 + p.theta1,
        diff_angles(q.theta2, p.theta2) / 2 + p.theta2,
        diff_angles(q.theta3, p.theta3) / 2 + p.theta3,
        diff_angles(q.theta4, p.theta4) / 2 + p.theta4,
    )
end

function dist(q1::StateSnake2D, q2::StateSnake2D)::Float64
    return norm([
        q1.x - q2.x,
        q1.y - q2.y,
        diff_angles(q1.theta1, q2.theta1) / π,
        diff_angles(q1.theta2, q2.theta2) / π,
        diff_angles(q1.theta3, q2.theta3) / π,
        diff_angles(q1.theta4, q2.theta4) / π,
    ])
end

function get_snake_positions(q::StateSnake2D, rad::Float64)::Vector{Vector{Float64}}
    p1 = [q.x, q.y]
    p2 = p1 + [cos(q.theta1), sin(q.theta1)] * rad
    p3 = p2 + [cos(q.theta2), sin(q.theta2)] * rad
    p4 = p3 + [cos(q.theta3), sin(q.theta3)] * rad
    p5 = p4 + [cos(q.theta4), sin(q.theta4)] * rad

    return [p1, p2, p3, p4, p5]
end

function check_self_collision_snake2d(
    P::Vector{Vector{Float64}},
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Bool
    any(
        e -> dist(e...) < safety_dist,
        [
            [P[1], P[2], P[3], P[4]],
            [P[1], P[2], P[4], P[5]],
            [P[2], P[3], P[4], P[5]],
            [P[3], P[4], P[5]],
        ],
    )
end

function gen_connect(
    q::StateSnake2D,
    obstacles::Vector{CircleObstacle2D},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST_SNAKE2D,
    max_dist::Union{Nothing,Float64} = nothing,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    # check: q \in C_free
    f(q::StateSnake2D, i::Int64)::Bool = begin
        P = get_snake_positions(q, rads[i])

        any(x -> (x < 0 || 1 < x), vcat(P...)) && return false

        # obstacle
        for k = 2:5
            any(o -> dist(P[k-1], P[k], [o.x, o.y]) < o.r, obstacles) && return false
        end

        # self collision
        check_self_collision_snake2d(P, safety_dist) && return false

        return true
    end

    f(q_from::StateSnake2D, q_to::StateSnake2D, i::Int64)::Bool = begin
        D = dist(q_from, q_to)

        !isnothing(max_dist) && D > max_dist && return false

        dt1 = diff_angles(q_to.theta1, q_from.theta1)
        dt2 = diff_angles(q_to.theta2, q_from.theta2)
        dt3 = diff_angles(q_to.theta3, q_from.theta3)
        dt4 = diff_angles(q_to.theta4, q_from.theta4)

        for e in vcat(collect(0:step_dist:D) / D, 1.0)
            # angles
            T = [
                q_from.theta1 + e * dt1,
                q_from.theta2 + e * dt2,
                q_from.theta3 + e * dt3,
                q_from.theta4 + e * dt4,
            ]
            # positions
            P = [(1 - e) * [q_from.x, q_from.y] + e * [q_to.x, q_to.y]]
            foreach(k -> push!(P, P[k-1] + rads[i] * [cos(T[k-1]), sin(T[k-1])]), 2:5)

            # outside
            any(x -> (x < 0 || 1 < x), vcat(P...)) && return false

            # obstacles
            for k = 2:5
                any(o -> dist(P[k-1], P[k], [o.x, o.y]) < o.r, obstacles) && return false
            end

            # self collision
            check_self_collision_snake2d(P, safety_dist) && return false
        end

        return true
    end
    return f
end

function gen_collide(
    q::StateSnake2D,
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST_SNAKE2D,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    N = length(rads)

    f(
        q_i_from::StateSnake2D,
        q_i_to::StateSnake2D,
        q_j_from::StateSnake2D,
        q_j_to::StateSnake2D,
        i::Int64,
        j::Int64,
    ) = begin
        # check each pair of step
        D_i = dist(q_i_from, q_i_to)
        D_j = dist(q_j_from, q_j_to)

        dt_i = [
            diff_angles(q_i_to.theta1, q_i_from.theta1),
            diff_angles(q_i_to.theta2, q_i_from.theta2),
            diff_angles(q_i_to.theta3, q_i_from.theta3),
            diff_angles(q_i_to.theta4, q_i_from.theta4),
        ]

        dt_j = [
            diff_angles(q_j_to.theta1, q_j_from.theta1),
            diff_angles(q_j_to.theta2, q_j_from.theta2),
            diff_angles(q_j_to.theta3, q_j_from.theta3),
            diff_angles(q_j_to.theta4, q_j_from.theta4),
        ]

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            T_i = [
                q_i_from.theta1 + e_i * dt_i[1],
                q_i_from.theta2 + e_i * dt_i[2],
                q_i_from.theta3 + e_i * dt_i[3],
                q_i_from.theta4 + e_i * dt_i[4],
            ]
            P_i = [(1 - e_i) * [q_i_from.x, q_i_from.y] + e_i * [q_i_to.x, q_i_to.y]]
            foreach(
                k -> push!(P_i, P_i[k-1] + rads[i] * [cos(T_i[k-1]), sin(T_i[k-1])]),
                2:5,
            )

            for e_j in vcat(collect(0:step_dist:D_j) / D_j, 1.0)
                # intermediate angles & positions for agent-j
                T_j = [
                    q_j_from.theta1 + e_j * dt_j[1],
                    q_j_from.theta2 + e_j * dt_j[2],
                    q_j_from.theta3 + e_j * dt_j[3],
                    q_j_from.theta4 + e_j * dt_j[4],
                ]
                P_j =
                    [(1 - e_j) * [q_j_from.x, q_j_from.y] + e_j * [q_j_to.x, q_j_to.y]]
                foreach(
                    k ->
                        push!(P_j, P_j[k-1] + rads[j] * [cos(T_j[k-1]), sin(T_j[k-1])]),
                    2:5,
                )

                # check_collision
                for k = 2:5, l = 2:5
                    dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
                end
            end
        end
        return false
    end

    f(q_i::StateSnake2D, q_j::StateSnake2D, i::Int64, j::Int64) = begin
        P_i = get_snake_positions(q_i, rads[i])
        P_j = get_snake_positions(q_j, rads[j])

        # check_collision
        for k = 2:5, l = 2:5
            dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
        end
        return false
    end

    f(Q_from::Vector{Node{StateSnake2D}}, Q_to::Vector{Node{StateSnake2D}}) = begin
        for i = 1:N, j = i+1:N
            if f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j)
                return true
            end
        end
        return false
    end

    f(Q::Vector{Node{StateSnake2D}}, q_i_to::StateSnake2D, i::Int64) = begin
        q_i_from = Q[i].q
        D_i = dist(q_i_from, q_i_to)

        dt_i = [
            diff_angles(q_i_to.theta1, q_i_from.theta1),
            diff_angles(q_i_to.theta2, q_i_from.theta2),
            diff_angles(q_i_to.theta3, q_i_from.theta3),
            diff_angles(q_i_to.theta4, q_i_from.theta4),
        ]

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            T_i = [
                q_i_from.theta1 + e_i * dt_i[1],
                q_i_from.theta2 + e_i * dt_i[2],
                q_i_from.theta3 + e_i * dt_i[3],
                q_i_from.theta4 + e_i * dt_i[4],
            ]
            P_i = [(1 - e_i) * [q_i_from.x, q_i_from.y] + e_i * [q_i_to.x, q_i_to.y]]
            foreach(
                k -> push!(P_i, P_i[k-1] + rads[i] * [cos(T_i[k-1]), sin(T_i[k-1])]),
                2:5,
            )

            for j in filter(k -> k != i, 1:N)
                P_j = get_snake_positions(Q[j].q, rads[j])

                # check collision
                for k = 2:5, l = 2:5
                    dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
                end
            end
        end

        return false
    end

    return f
end

function gen_uniform_sampling(q::StateSnake2D)::Function
    () -> begin
        return StateSnake2D(rand(2)..., (rand(4) * 2π)...)
    end
end

function plot_motion!(q_from::StateSnake2D, q_to::StateSnake2D, rad::Float64, params)
    plot!([q_from.x, q_to.x], [q_from.y, q_to.y]; params...)
end

function plot_start_goal!(q_init::StateSnake2D, q_goal::StateSnake2D, rad::Float64, params)
    plot!([q_init.x], [q_init.y]; markershape = :hex, params...)
    plot!([q_goal.x], [q_goal.y]; markershape = :star, params...)
end

function plot_agent!(q::StateSnake2D, rad::Float64, color::String)
    P = get_snake_positions(q, rad)
    plot!(
        map(pos -> pos[1], P),
        map(pos -> pos[2], P),
        color = color,
        lw = 5,
        markershape = :circle,
        markersize = 3,
        label = nothing,
    )
    scatter!(
        [q.x],
        [q.y],
        markershpa = 10,
        markershape = :rect,
        color = color,
        label = nothing,
    )
end

function gen_random_instance_StateSnake2D(; params...)
    return gen_random_instance(StateSnake2D(0, 0, 0, 0, 0, 0); params...)
end


text:
using AnalysisFTIR
using Test

@testset "AnalysisFTIR.jl" begin
    # Write your tests here.
end

ext: jl


text:
import Dates: Nanosecond, Second
import TimesDates: TimeDate
import Base: ==, +, -

abstract type TDMSUnimplementedType end
struct tdsTypeVoid <: TDMSUnimplementedType end
struct tdsTypeExtendedFloat <: TDMSUnimplementedType end
struct tdsTypeSingleFloatWithUnit <: TDMSUnimplementedType end
struct tdsTypeDoubleFloatWithUnit <: TDMSUnimplementedType end
struct tdsTypeExtendedFloatWithUnit <: TDMSUnimplementedType end
struct tdsTypeDAQmxRawData <: TDMSUnimplementedType end

abstract type TDMSType end
struct TimeStamp <: TDMSType
    seconds::Int64
    fractions::UInt64
end
const epoch=TimeDate(1904,1,1)
function TimeDate(x::TimeStamp)
    ns = round(x.fractions*(2.0^-64)*1e9)
    epoch + Second(x.seconds) + Nanosecond(ns)
end
function +(x::T,y::T) where {T<:TimeStamp}
    s = x.seconds + y.seconds
    r,f = Base.Checked.add_with_overflow(x.fractions, y.fractions)
    f ? TimeStamp(s+1,r) : TimeStamp(s,r)
end
function -(x::T,y::T) where {T<:TimeStamp}
    s = x.seconds - y.seconds
    r,f = Base.Checked.sub_with_overflow(x.fractions, y.fractions)
    f ? TimeStamp(s-1,r) : TimeStamp(s,r)
end
Base.show(io::IO,ts::TimeStamp) = show(io,TimeDate(ts))
Base.show(ts::TimeStamp) = show(TimeDate(ts))

function Base.read(s::IO, ::Type{TimeStamp})
    fractions=read(s,UInt64)
    seconds=read(s, Int64)
    TimeStamp(seconds,fractions)
end

const tdsTypes=Dict{UInt32,DataType}(0=>tdsTypeVoid, 1=>Int8, 2=>Int16, 3=>Int32, 4=>Int64,
    5=>UInt8, 6=>UInt16, 7=>UInt32, 8=>UInt64, 9=>Float32, 10=>Float64, 11=>tdsTypeExtendedFloat,
    0x19=>tdsTypeSingleFloatWithUnit, 0x1A => tdsTypeDoubleFloatWithUnit, 0x1B => tdsTypeExtendedFloatWithUnit,
    0x20=>String, 0x21=>Bool, 0x08000C=>ComplexF32, 0x10000D=>ComplexF64, 0xFF_FF_FF_FF=>tdsTypeDAQmxRawData,
    0x44=>TimeStamp )

struct ToC
    kTocMetaData::Bool
    kTocRawData::Bool
    kTocDAQmxRawData::Bool
    kTocInterleavedData::Bool
    kTocBigEndian::Bool
    kTocNewObjList::Bool
end

ToC(a::UInt32) = ToC(bget.(a, (1,3,7,5,6,2))...)

struct Channel{T}
    data::Vector{T}
    props::OrderedDict{String,Any}

    Channel{T}() where {T} = new(Vector{T}(),OrderedDict{String,Any}())
    Channel{T}(props) where {T} = new(Vector{T}(),props) #already have props
end
==(a::TDMSReader.Channel,b::TDMSReader.Channel) = (a.props == b.props) && (a.data == b.data)

struct Group
    channels::OrderedDict{String,TDMSReader.Channel}
    props::OrderedDict{String,Any}

    Group()= new(OrderedDict{String,Channel}(),OrderedDict{String,Any}())
end
Base.haskey(h::Group, key) = haskey(h.channels, key)
Base.getkey(h::Group, key, default) = getkey(h.channels, key, default)
Base.get!(h::Group, key, default) = get!(h.channels, key, default)
Base.get!(f::Union{Function, Type}, h::Group, key) = get!(f, h.channels, key)
Base.keys(a::Group) = keys(a.channels)
==(a::Group,b::Group) = (a.props == b.props) && (a.channels == b.channels)

struct File
    groups::OrderedDict{String,Group}
    props::OrderedDict{String,Any}

    File() = new(OrderedDict{String,Group}(),OrderedDict{String,Any}())
end
Base.haskey(h::File, key) = haskey(h.groups, key)
Base.getkey(h::File, key, default) = getkey(h.groups, key, default)
Base.get!(h::File, key, default) = get!(h.groups, key, default)
Base.get!(f::Union{Function, Type}, h::File, key) = get!(f, h.groups, key)
Base.keys(a::File) = keys(a.groups)
==(a::File,b::File) = (a.props == b.props) && (a.groups == b.groups)

struct Chunk{T}
    data::Vector{T}
    nsamples::UInt64
end

struct ObjDict
    current::OrderedDict{String,Chunk}
    full::OrderedDict{String,Chunk}

    ObjDict() = new(OrderedDict{String,Chunk}(),OrderedDict{String,Chunk}())
end

function Base.getindex(f::File, group::Union{Integer,AbstractString}, channel::Nothing=nothing)
    if group isa Integer
        k=collect(keys(f.groups))[group]
        f.groups[k]
    else
        f.groups[group]
    end
end

function Base.getindex(f::File, group::Union{Integer,AbstractString}, channel::Union{Integer,AbstractString})
    if group isa Integer
        k=collect(keys(f.groups))[group]
        f.groups[k][channel]
    else
        f.groups[group][channel]
    end
end

function Base.getindex(g::Group, channel::Union{Integer,AbstractString})
    if channel isa Integer
        k=collect(keys(g.channels))[channel]
        g.channels[k]
    else
        g.channels[channel]
    end
end

struct SegInfo
    position::Int64
    toc::ToC
    nextsegmentposition::Int64
    rawdataposition::Int64
    nobj::Int64
    rawdatasize::Int64
end


text:
module Geometry

Base.Experimental.@optlevel 1

export AbstractGeometry,
    AbstractParallelBeamGeometry, AbstractFanBeamGeometry
export ParallelBeamGeometry, FanBeamGeometry
export AbstractTomograph, DefaultTomograph
export num_proj, num_det, f2iso, f2det, fan_angle, cell_size
export scan_angle, start_angle, center_channel
export num_rows, num_cols, tomograph, channel_spacing

using ..Monads
import ..Utils: yaml_repr, struct2dict
import Base: show, getproperty, eltype

const _geometry_names = (:ParallelBeamGeometry, :FanBeamGeometry)

abstract type AbstractGeometry  end

abstract type AbstractTomograph end

function getproperty(ct::AbstractTomograph, s::Symbol)
    s ∈ fieldnames(typeof(ct)) && return getfield(ct, s)
    nothing
end

show(io::IO, ::T) where T <: AbstractTomograph = print(io, "$(nameof(T))")

struct DefaultTomograph <: AbstractTomograph end

abstract type AbstractParallelBeamGeometry <: AbstractGeometry end

struct ParallelBeamGeometry{T <: Real,CT <: AbstractTomograph} <: AbstractParallelBeamGeometry
    ct::CT
    nϕ::Int
    nd::Int
    rows::Int
    cols::Int
    α::T
    α₀::T
    center::T
end

@inline function getproperty(g::ParallelBeamGeometry, s::Symbol)
    s ≡ :width && return g.cols
    s ≡ :height && return g.rows
    s ≡ :nphi && return getfield(g, :nϕ)
    s ≡ :alpha && return getfield(g, :α)
    s ≡ :alpha0 && return getfield(g, :α₀)
    getfield(g, s)
end

eltype(::Type{G}) where {G <: ParallelBeamGeometry{T}} where {T} = T

yaml_repr(g::ParallelBeamGeometry) = struct2dict(g)

tomograph(x::ParallelBeamGeometry) = x.ct
num_proj(x::ParallelBeamGeometry) = x.nϕ
num_det(x::ParallelBeamGeometry) = x.nd
num_rows(x::ParallelBeamGeometry) = x.rows
num_cols(x::ParallelBeamGeometry) = x.cols
scan_angle(x::ParallelBeamGeometry) = x.α
start_angle(x::ParallelBeamGeometry) = x.α₀
center_channel(x::ParallelBeamGeometry) = x.center

"""
    ParallelBeamGeometry([T=Float32]; <keyword arguments>) where {T<:Real}

Construct geometry for the simulation.

# Arguments
- `nϕ`: number of projections.
- `nd=nothing`: number of detectors. If not specified, same as `nϕ`.
- `rows=nothing`: number of rows in the reconstructed image. If not specified, same as `nd`.
- `cols=nothing`: number of columns in the reconstructed image. If not
  specified, same as `rows`.
- `width=nothing`: alias for `cols`.
- `height=nothing`: alias for `rows`.
- `α=360`: scan angle in degrees.
- `α₀=0`: starting scan angle in degrees.
- `center=nothing`: virtual center channel, defaults to `(nd-1)/2`.
"""
function ParallelBeamGeometry(
    ::Type{T} = Float32,
    ct::AbstractTomograph = DefaultTomograph();
    nϕ::Optional{Integer} = nothing,
    nphi::Optional{Integer} = nothing,
    nd::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    α::Real = 360,
    alpha::Optional{Real} = nothing,
    α₀::Real = 0,
    alpha0::Optional{Real} = nothing,
    center::Optional{Real} = nothing,
) where {T <: Real}
    nϕ = maybe(nϕ, nphi)
    α = maybe(α, alpha)
    α₀ = maybe(α₀, alpha0)
    height = maybe(width, height)
    width = maybe(height, width)
    rows = maybe(512, maybe(height, rows))
    cols = maybe(rows, maybe(width, cols))
    nd = isnothing(nd) ? round(Int, hypot(rows, cols)) : nd
    nϕ = isnothing(nϕ) ?
        2 * (round(Int, (rows + 1) * (cols + 1) / nd) ÷ 2) + 1 : nϕ
    center′::T = maybe((nd - 1) / 2, center)
    ParallelBeamGeometry{T,typeof(ct)}(
        ct,
        nϕ,
        nd,
        rows,
        cols,
        α,
        α₀,
        center′,
    )
end

ParallelBeamGeometry(ct::CT; kwargs...) where {CT <: AbstractTomograph} =
    ParallelBeamGeometry(Float32, ct; kwargs...)

show(io::IO, g::AbstractParallelBeamGeometry) = print(
    io,
    """Parallel beam geometry:
- tomograph: $(g.ct)
- number of projections: $(num_proj(g))
- number of detectors: $(num_det(g))
- tomogram size (W×H): $(num_cols(g)) × $(num_rows(g))
- scan angle: $(scan_angle(g))°
- scan starting angle: $(start_angle(g))°
- center channel: $(center_channel(g))""",
)

abstract type AbstractFanBeamGeometry <: AbstractGeometry end

struct FanBeamGeometry{T,CT <: AbstractTomograph} <:
       AbstractFanBeamGeometry
    ct::CT
    nϕ::Int
    nd::Int
    rows::Int
    cols::Int
    D::T
    D′::T
    γ::T
    δ::T
    α::T
    α₀::T
    center::T
end

@inline function getproperty(g::FanBeamGeometry, s::Symbol)
    s ≡ :width && return g.cols
    s ≡ :height && return g.rows
    s ≡ :nphi && return getfield(g, :nϕ)
    s ≡ :alpha && return getfield(g, :α)
    s ≡ :alpha0 && return getfield(g, :α₀)
    s ≡ :D1 && return getfield(g, :D′)
    s ≡ :gamma && return getfield(g, :γ)
    s ≡ :dx && return getfield(g, :δ)
    getfield(g, s)
end

eltype(::Type{G}) where {G <: FanBeamGeometry{T}} where {T} = T

yaml_repr(g::FanBeamGeometry) = struct2dict(g)

tomograph(x::FanBeamGeometry) = x.ct
num_proj(x::FanBeamGeometry) = x.nϕ
num_det(x::FanBeamGeometry) = x.nd
num_rows(x::FanBeamGeometry) = x.rows
num_cols(x::FanBeamGeometry) = x.cols
f2iso(x::FanBeamGeometry) = x.D
f2det(x::FanBeamGeometry) = x.D′
fan_angle(x::FanBeamGeometry) = x.γ
cell_size(x::FanBeamGeometry) = x.δ
scan_angle(x::FanBeamGeometry) = x.α
start_angle(x::FanBeamGeometry) = x.α₀
channel_spacing(x::FanBeamGeometry) = tomograph(x).channel_spacing
center_channel(x::FanBeamGeometry) = x.center

"""
    FanBeamGeometry([T=Float32]; <keyword arguments>) where {T<:Real}

Create a new FanBeamGeometry object with given parameters.

It is assumed that the detectors are arranged on a circle arc,
not flat.

# Arguments
- `nϕ: number of scan angles.
- `nd=nothing`: number of detectors.
- `rows=nothing`: number of rows in the reconstructed image.
- `cols=nothing`: number of columns in the reconstructed image.
- `width=nothing`: alias for `cols`.
- `height=nothing`: alias for `rows`.
- `D`: focal spot to ISO distance.
- `D′=nothing`: focal spot to detectors distance; if not specified
    and also `γ` is not specified, is defined so that the total
    fan angle is 1 radian.
- `γ=nothing`: total fan angle in degrees.
- `δ=1`: detectors spacing (cell size).
- `α=360`: scan angle in degrees.
- `α₀=0`: scan starting angle in degrees.
- `center=nothing`: virtual center channel, defaults to `(nd-1)/2`.
"""
function FanBeamGeometry(
    ::Type{T} = Float32,
    ct::AbstractTomograph = DefaultTomograph();
    nϕ::Optional{Integer} = nothing,
    nphi::Optional{Integer} = nothing,
    D::Real = 500,
    D′::Optional{Real} = nothing,
    D1::Optional{Real} = nothing,
    γ::Optional{Real} = nothing,
    gamma::Optional{Real} = nothing,
    nd::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    δ::Real = one(T),
    dx::Real = one(T),
    α::Real = 360,
    alpha::Real = 360,
    α₀::Real = zero(T),
    alpha0::Real = zero(T),
    center::Optional{Real} = nothing,
) where {T <: Real}
    nϕ = maybe(nϕ, nphi)
    D′ = maybe(D′, D1)
    γ = maybe(γ, gamma)
    δ = maybe(δ, dx)
    α = maybe(α, alpha)
    α₀ = maybe(α₀, alpha0)
    height = maybe(width, height)
    width = maybe(height, width)
    rows = maybe(512, maybe(height, rows))
    cols = maybe(rows, maybe(width, cols))
    nd = isnothing(nd) ? round(Int, hypot(rows, cols)) : nd
    nϕ = isnothing(nϕ) ?
        2 * (round(Int, (rows + 1) * (cols + 1) / nd) ÷ 2) + 1 : nϕ
    γ′::T = maybe(deg2rad, maybe(x->nd * δ / x, one(T), D′), γ)
    D′′::T = maybe(nd * δ / γ′, D′) # changed γ to γ′ here, should be tested
    center′::T = maybe((nd - 1) / 2, center)
    FanBeamGeometry{T,typeof(ct)}(
        ct,
        nϕ,
        nd,
        rows,
        cols,
        D,
        D′′,
        γ′,
        δ,
        α,
        α₀,
        center′
    )
end

FanBeamGeometry(ct::AbstractTomograph; kwargs...) =
    FanBeamGeometry(Float32, ct; kwargs...)

function show(io::IO, g::AbstractFanBeamGeometry)
    print(
    io,
    """Fan beam geometry:
    - tomograph: $(tomograph(g))
    - number of angles: $(num_proj(g))
    - number of detectors: $(num_det(g))
    - tomogram size (W×H): $(num_cols(g)) × $(num_rows(g))
    - focal spot to ISO: $(f2iso(g)) mm
    - fan beam angle: $(rad2deg(fan_angle(g)))°
    - cell size: $(cell_size(g)) mm
    - scan angle: $(scan_angle(g))°
    - scan starting angle: $(start_angle(g))°
    - channel spacing: $(maybe(x -> "$x°", "N/A", channel_spacing(g)))
    - center channel: $(center_channel(g))""",
    )
end

function ParallelBeamGeometry(fbg::AbstractFanBeamGeometry)
    ParallelBeamGeometry(
        eltype(fbg),
        fbg.ct;
        fbg.nϕ,
        fbg.nd,
        fbg.rows,
        fbg.cols,
        fbg.α,
        fbg.α₀,
        fbg.center,
    )
end

function FanBeamGeometry(
    pbg::AbstractParallelBeamGeometry,
    ct::AbstractTomograph = DefaultTomograph();
    kwargs...,
)
    FanBeamGeometry(
        eltype(pbg),
        ct;
        pbg.nϕ,
        pbg.nd,
        pbg.rows,
        pbg.cols,
        pbg.α,
        pbg.α₀,
        pbg.center,
        kwargs...,
    )
end

end # module


text:
#=
# Euclidean Distance Matrix Completion
This example illustrates how to use Chordal.jl to complete a partially-specified
Euclidean Distance Matrix (or determine no such completion exists).
=#

#=
## Euclidean Distance Matrices
Given $N$ vectors $x_1,\dots,x_n$ in $\mathbb{R}^n$, the associated Euclidean
Distance Matrix (EDM) $D \in \mathbb{R}^{N \times N}$ records the squared distance
between each pair $x_i, x_j$. Specifically,
$$
D_{ij} = \|x_i - x_j\|_2^2.
$$

Given a partially-specified EDM $D$ which has a chordal sparsity pattern, we would
like to find the missing entries.
=#

using Chordal
using LinearAlgebra, SparseArrays

#=
## Small Example
This small example is taken from Jon Dattorro's [Convex Optimization & Euclidean Distance Geometry](https://ccrma.stanford.edu/~dattorro/EDM.pdf)
eq (1041). First, we construct the partially-specified EDM.
=#
D = sparse([
    0.0 1 5 0;
    1 0 4 1;
    5 4 0 1;
    0 1 1 0
])
@show D

## Next, we complete the matrix
D_complete = sparse(edm_completion(D))
println("\nCompletion:")
@show D_complete

#=
## Larger Random Example

First, we generate an EDM from 17 vectors in $\mathbb{R}^3}$.
=#
n = 17
r = 3
V = randn(n, r)
VV = V*V'
D_full = diag(VV)*ones(n)' + ones(n)*diag(VV)' - 2VV;

#=
Next, we remove entries. and complete the
The sparsity pattern is from Figure 4.2 in Lieven Vandenberghe and Martin S. Andersen's
[Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf).
=#
ijs = [(3,1), (3,2), (4,1), (4,2), (4,3), (5,1), (5,3), (5,4), (8,7),
       (9,5), (9,6), (9,7), (9,8), (11,10), (13,10), (13,11), (13,12),
       (14,10), (14,11), (14,12), (14,13),
       (15, 1), (15,3), (15,4), (15,5), (15,7), (15,8), (15,9),
       (16,5), (16,6), (16,9), (16,12), (16,13), (16,14), (16,15),
       (17,10), (17,11), (17,12), (17,13), (17,14), (17,15), (17,16)]
append!(ijs, [(i,i) for i in 1:n])
II, JJ = Chordal.unzip(ijs)
sp = sparse(II, JJ, ones(length(II)))
sp = sp + tril(sp)'

## Remove entries s.t. the remaining entries have a chordal sparsity pattern
D = sp .* D_full

## Complete the matrix
@show Chordal.is_edm(edm_completion(D))


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Generate synthetic data for the Sleipner model

## Module loading

using DrWatson
@quickactivate "UncertaintyQuantificationAVP"
using UncertaintyQuantificationAVP
using LinearAlgebra, PyPlot, Statistics
using Random; seed = 1; Random.seed!(seed)

## Loading models

# Load velocity model
model_dict = wload(projectdir()*"/scripts/Sleipner/Sleipner_model.bson")
Mtrue = model_dict[:Mtrue]
Mbg = model_dict[:Mbg]

## Seismic experiment settings

# Data parameters
T = 0.4f0
dt = 0.0005f0
nt = Int64(ceil(T/dt))+1
θ = range(0f0, stop = pi*32f0/180f0, length = 11)
p = sin.(θ)*sqrt(max(Mtrue.m...)); p = p[5:end]
dat_pars = DataParams(nt, dt, p)

# Source parameters
zsrc = 0f0
wav = cfreq_wavelet(Int64(round(0.04f0/dt)), nt, dt, (5f0, 30f0, 50f0, 80f0))
src_pars = SrcParams(wav, dt, zsrc)

# ABC parameters
abc_size = (50, 50)
abc_fact = 1f-4
abc_geom = ABC_Geom(Mbg.geom.nz, abc_size, abc_fact)

# Modeling operator
acc_mod = 4
F = seismod_fun(dat_pars, src_pars; abc_geom = abc_geom, acc = acc_mod)

# Generate synthetics
dat = F.eval(Mtrue)
dat_bg = F.eval(Mbg)

# Adding white noise
snr = 0f0
σ2_dat = var((dat-dat_bg).d)
σ2_n = σ2_dat*10f0^(-snr/10f0)
dat_n = dat+sqrt(σ2_n)*randn(Float32, size(dat.d))

## Saving results

data_dict = Dict(:abc_geom => abc_geom, :src_pars => src_pars, :acc_mod => acc_mod, :dat => dat, :σ2_n => σ2_n, :dat_n => dat_n)
wsave(datadir("Sleipner", "Sleipner_data.bson"), data_dict)


text:
# This functions calculates the optical thickness based on pressure
# and lapse rate for a gray atmosphere.
# See Schneider 2004, J. Atmos. Sci. (2004) 61 (12): 1317–1340.
# https://doi.org/10.1175/1520-0469(2004)061<1317:TTATTS>2.0.CO;2
"""
    compute_optical_props_kernel!(
        op::AbstractOpticalProps{FT},
        as::GrayAtmosphericState{FT},
        glaycol,
        source::AbstractSourceLW{FT},
    ) where {FT<:AbstractFloat}

This function computes the optical properties using the gray atmosphere assumption
for the longwave solver.
"""
function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::GrayAtmosphericState{FT},
    glaycol,
    source::AbstractSourceLW{FT},
) where {FT<:AbstractFloat}

    compute_optical_props_kernel!(op, as, glaycol)     # computing optical thickness
    compute_sources_gray_kernel!(source, as, glaycol) # computing Planck sources
end

"""
    compute_optical_props_kernel!(
        op::AbstractOpticalProps{FT},
        as::GrayAtmosphericState{FT},
        glaycol,
    ) where {FT<:AbstractFloat}

This function computes the optical properties using the gray atmosphere assumption
for the shortwave solver.
"""
function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::GrayAtmosphericState{FT},
    glaycol,
) where {FT<:AbstractFloat}
    # setting references
    glay, gcol = glaycol
    (; p_lay, p_lev, d0, α) = as
    @inbounds p0 = p_lev[1, gcol]

    @inbounds op.τ[glaycol...] = abs(
        (α * d0[gcol] * (p_lay[glaycol...] ./ p0) .^ α ./ p_lay[glaycol...]) * (p_lev[glay+1, gcol] - p_lev[glaycol...]),
    )

    if op isa TwoStream
        op.ssa[glaycol...] = FT(0)
        op.g[glaycol...] = FT(0)
    end
end
"""
    compute_sources_gray_kernel!(
        source::AbstractSourceLW{FT},
        as::GrayAtmosphericState{FT},
        glaycol,
    ) where {FT<:AbstractFloat}

This function computes the Planck sources for the gray longwave solver.
"""
function compute_sources_gray_kernel!(
    source::AbstractSourceLW{FT},
    as::GrayAtmosphericState{FT},
    glaycol,
) where {FT<:AbstractFloat}
    # computing Planck sources
    glay, gcol = glaycol
    (; t_lay, t_lev) = as
    (; lay_source, lev_source_inc, lev_source_dec, sfc_source) = source

    sbc = FT(Stefan())
    @inbounds lay_source[glaycol...] = sbc * t_lay[glaycol...]^FT(4) / FT(π)   # computing lay_source
    @inbounds lev_source_inc[glaycol...] =
        sbc * t_lev[glay+1, gcol]^FT(4) / FT(π)
    @inbounds lev_source_dec[glaycol...] = sbc * t_lev[glaycol...]^FT(4) / FT(π)
    if glay == 1
        @inbounds sfc_source[gcol] = sbc * as.t_sfc[gcol]^FT(4) / FT(π)   # computing sfc_source
    end
end


text:
__precompile__(true)
module VulkanCore

# Vulkan low level api
include("api.jl")

include("../docs/vulkan_doc.jl")

const vk = api
export vk

end # module


text:
text = readchomp("input.txt")
part1(text) = count(==('('), text) - count(==(')'), text)

function part2(text)
    total = 0
    for (i, c) in Iterators.enumerate(text)
        if c == '('
            total += 1
        elseif c == ')'
            total -= 1
            if total == -1
                return i
            end
        end
    end
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using GLFW
using VulkanCore
using VulkanCore.LibVulkan

@assert GLFW.VulkanSupported()

const WIDTH = 800
const HEIGHT = 600

enableValidationLayers = true

## init GLFW window
GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)    # not to create an OpenGL context
GLFW.WindowHint(GLFW.RESIZABLE, 0)
window = GLFW.CreateWindow(WIDTH, HEIGHT, "Vulkan")

## init Vulkan
## creating instance
# fill application info
sType = VK_STRUCTURE_TYPE_APPLICATION_INFO
pApplicationName = pointer("Vulkan Instance")
applicationVersion = VK_MAKE_VERSION(1, 0, 0)
pEngineName = pointer("No Engine")
engineVersion = VK_MAKE_VERSION(1, 0, 0)
apiVersion = VK_API_VERSION_1_2
appInfoRef = VkApplicationInfo(
    sType,
    C_NULL,
    pApplicationName,
    applicationVersion,
    pEngineName,
    engineVersion,
    apiVersion,
) |> Ref

# checking for extension support
requiredExtensions = GLFW.GetRequiredInstanceExtensions()
if enableValidationLayers
    push!(requiredExtensions, VK_EXT_DEBUG_UTILS_EXTENSION_NAME)  # for debugging, see message callback section
end

extensionCountRef = Ref(Cuint(0))
vkEnumerateInstanceExtensionProperties(C_NULL, extensionCountRef, C_NULL)
extensions = Vector{VkExtensionProperties}(undef, extensionCountRef[])
vkEnumerateInstanceExtensionProperties(C_NULL, extensionCountRef, extensions)
extensionNames = map(extensions) do extension
    extension.extensionName |> collect |> String |> x -> strip(x, '\0')
end
extensionVersions = [ext.specVersion |> Int for ext in extensions]

@info "available extensions:"
for (ext, ver) in zip(extensionNames, extensionVersions)
    @info "  $ext: $ver"
end
if !all(x->x in extensionNames, requiredExtensions)
    @error "not all required extensions are supported."
end

# using validation layers
requiredLayers = ["VK_LAYER_KHRONOS_validation"]
layerCountRef = Ref(Cuint(0))
vkEnumerateInstanceLayerProperties(layerCountRef, C_NULL)
availableLayers = Vector{VkLayerProperties}(undef, layerCountRef[])
vkEnumerateInstanceLayerProperties(layerCountRef, availableLayers)
availableLayerNames = map(availableLayers) do layer
    strip(String(collect(layer.layerName)), '\0')
end
availableLayerDescription = map(availableLayers) do layer
    strip(String(collect(layer.description)), '\0')
end

@info "available layers:"
for (name,description) in zip(availableLayerNames, availableLayerDescription)
    @info "  $name: $description"
end
if !all(x->x in availableLayerNames, requiredLayers)
    @error "not all required layers are supported."
end

# add required extensions and layers
enabledExtensionCount = Cuint(length(requiredExtensions))
ppEnabledExtensionNames = Base.unsafe_convert(Ptr{Cstring}, Base.cconvert(Ptr{Cstring}, requiredExtensions))

enabledLayerCount = Cuint(length(requiredLayers))
ppEnabledLayerNames = Base.unsafe_convert(Ptr{Cstring}, Base.cconvert(Ptr{Cstring}, requiredLayers))

# create instance after setting up the debug messenger
instanceRef = Ref(VkInstance(C_NULL))

## message callback
function debug_callback(severity::VkDebugUtilsMessageSeverityFlagBitsEXT, type::VkDebugUtilsMessageTypeFlagsEXT, pCallbackData::Ptr{VkDebugUtilsMessengerCallbackDataEXT}, pUserData::Ptr{Cvoid})::VkBool32
    data = unsafe_load(pCallbackData)
    if severity == VK_DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT
        @debug "validation layer: $(Base.unsafe_string(data.pMessage))"
    elseif severity == VK_DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT
        @info "validation layer: $(Base.unsafe_string(data.pMessage))"
    elseif severity == VK_DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT
        @warn "validation layer: $(Base.unsafe_string(data.pMessage))"
    elseif severity == VK_DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT
        @error "validation layer: $(Base.unsafe_string(data.pMessage))"
    end
    return VK_FALSE
end

# set debug messenger
debugMessengerRef = Ref(VkDebugUtilsMessengerEXT(C_NULL))

sType = VK_STRUCTURE_TYPE_DEBUG_UTILS_MESSENGER_CREATE_INFO_EXT
flags = UInt32(0)
messageSeverity = VK_DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT | VK_DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT | VK_DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT
messageType = VK_DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT | VK_DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT | VK_DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT
pfnUserCallback = @cfunction(debug_callback, VkBool32, (VkDebugUtilsMessageSeverityFlagBitsEXT, VkDebugUtilsMessageTypeFlagsEXT, Ptr{VkDebugUtilsMessengerCallbackDataEXT},Ptr{Cvoid}))
pUserData = C_NULL
debugCreateInfoRef = VkDebugUtilsMessengerCreateInfoEXT(sType, C_NULL, flags, messageSeverity, messageType, pfnUserCallback, pUserData) |> Ref

function CreateDebugUtilsMessengerEXT(instance::VkInstance, pCreateInfo::Ref{VkDebugUtilsMessengerCreateInfoEXT}, pAllocator, pDebugMessenger::Ref{VkDebugUtilsMessengerEXT})
    fp = vkGetInstanceProcAddr(instance, "vkCreateDebugUtilsMessengerEXT")
    if fp != C_NULL
        ccall(fp, VkResult, (VkInstance, Ptr{VkDebugUtilsMessengerCreateInfoEXT}, Ptr{VkAllocationCallbacks}, Ptr{VkDebugUtilsMessengerEXT}), instance, pCreateInfo, pAllocator, pDebugMessenger)
    else
        return VK_ERROR_EXTENSION_NOT_PRESENT
    end
end

# create instance
sType = VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO
flags = UInt32(0)
pApplicationInfo = Base.unsafe_convert(Ptr{VkApplicationInfo}, appInfoRef)
pNext = Base.unsafe_convert(Ptr{Cvoid}, debugCreateInfoRef)
createInfoRef = VkInstanceCreateInfo(
    sType,
    pNext,
    flags,
    pApplicationInfo,
    enabledLayerCount,
    ppEnabledLayerNames,
    enabledExtensionCount,
    ppEnabledExtensionNames,
) |> Ref

result = GC.@preserve appInfoRef requiredExtensions requiredLayers debugCreateInfoRef vkCreateInstance(createInfoRef, C_NULL, instanceRef)
@assert result == VK_SUCCESS "failed to create instance!"

# set debug messenger
if CreateDebugUtilsMessengerEXT(instanceRef[], debugCreateInfoRef, C_NULL, debugMessengerRef) != VK_SUCCESS
    @error "failed to set up debug messenger!"
end

## physical device
physicalDeviceRef = Ref(VkPhysicalDevice(VK_NULL_HANDLE))

deviceCountRef = Ref(UInt32(0))
vkEnumeratePhysicalDevices(instanceRef[], deviceCountRef, C_NULL)
deviceCount = deviceCountRef[]
if deviceCount == 0
    @error "failed to find GPUs with Vulkan support!"
end

devices = Vector{VkPhysicalDevice}(undef, deviceCount)
vkEnumeratePhysicalDevices(instanceRef[], deviceCountRef, devices)

# device suitability checks
function is_device_suitable(device::VkPhysicalDevice)
    devicePropertiesRef = Ref{VkPhysicalDeviceProperties}()
    deviceFeaturesRef = Ref{VkPhysicalDeviceFeatures}()
    vkGetPhysicalDeviceProperties(device, devicePropertiesRef)
    vkGetPhysicalDeviceFeatures(device, deviceFeaturesRef)
    return (devicePropertiesRef[].deviceType == VK_PHYSICAL_DEVICE_TYPE_DISCRETE_GPU ||
            devicePropertiesRef[].deviceType == VK_PHYSICAL_DEVICE_TYPE_INTEGRATED_GPU ) &&
            deviceFeaturesRef[].geometryShader == VK_SUCCESS
end

for device in devices
    if is_device_suitable(device)
        physicalDeviceRef[] = device
        break
    end
end

if physicalDeviceRef[] == VkPhysicalDevice(VK_NULL_HANDLE)
    @error "failed to find a suitable GPU!"
end

## queue families
mutable struct QueueFamilyIndices
    graphicsFamily::Union{Nothing,Cuint}
end
QueueFamilyIndices() = QueueFamilyIndices(nothing)
is_complete(x::QueueFamilyIndices) = x.graphicsFamily != nothing

indices = QueueFamilyIndices()
queueFamilyCountRef = Ref(Cuint(0))
vkGetPhysicalDeviceQueueFamilyProperties(physicalDeviceRef[], queueFamilyCountRef, C_NULL)
queueFamilies = Vector{VkQueueFamilyProperties}(undef, queueFamilyCountRef[])
vkGetPhysicalDeviceQueueFamilyProperties(physicalDeviceRef[], queueFamilyCountRef, queueFamilies)

for (i,queueFamily) in enumerate(queueFamilies)
    if Bool(queueFamily.queueFlags & VK_QUEUE_GRAPHICS_BIT)
        indices.graphicsFamily = i
    end
    is_complete(indices) && break
end

@info indices

## main loop
while !GLFW.WindowShouldClose(window)
    GLFW.PollEvents()
end

## cleaning up
function DestroyDebugUtilsMessengerEXT(instance::VkInstance, debugMessenger::VkDebugUtilsMessengerEXT, pAllocator)
    fp = vkGetInstanceProcAddr(instance, "vkDestroyDebugUtilsMessengerEXT")
    if fp != C_NULL
        ccall(fp, VkResult, (VkInstance, VkDebugUtilsMessengerEXT, Ptr{VkAllocationCallbacks}), instance, debugMessenger, pAllocator)
    end
end

if enableValidationLayers
    DestroyDebugUtilsMessengerEXT(instanceRef[], debugMessengerRef[], C_NULL)
end

vkDestroyInstance(instanceRef[], C_NULL)
GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
function print_buffers(btg::btg)

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Test
include("../nanomsg.jl")

# Pipeline example
# NN_PROTO_PIPELINE 5
#define NN_PUSH (NN_PROTO_PIPELINE * 16 + 0)
@test nn_socket(5 * 16 + 0) > -1

@test nn_strerror(22) == "Invalid argument"

# Invalid protocol
@test_throws ErrorException nn_socket(16 * 16 + 0)


text:
# This is run multiple times during testing,
# for simpler running of individual tests.
function custom_metafmt(level, _module, group, id, file, line)
    color = WebSockets.default_logcolor(level)
    prefix = string(level) * ':'
    suffix = " $(Int(round((Dates.now() - T0_TESTS).value / 1000))) s @"
    _module !== nothing && (suffix *= "$(_module)")
    if file !== nothing
        _module !== nothing && (suffix *= " ")
        suffix *= Base.contractuser(file)
        if line !== nothing
            suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    return color, prefix, suffix
end

if !@isdefined OLDLOGGER
    const OLDLOGGER = WebSockets.global_logger()
    const T0_TESTS = Dates.now()
end

if !@isdefined TESTLOGR
    const TESTLOGR = WebSockets.WebSocketLogger(stderr, Base.CoreLogging.Debug, meta_formatter = custom_metafmt)
    WebSockets.global_logger(TESTLOGR)
end
nothing


text:
function (inv){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = accInvForSign(a)
    eftRound(hi, lo, rounding)
end

function (+){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftAdd(a,b)
    eftRound(hi, lo, rounding)
end

function (-){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftSub(a,b)
    eftRound(hi, lo, rounding)
end

function (*){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftMul(a,b)
    eftRound(hi, lo, rounding)
end

function (/){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = accDivForSign(a,b)
    eftRound(hi, lo, rounding)
end

function (\){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = accDivForSign(b,a)
    eftRound(hi, lo, rounding)
end

function (sqr){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSqr(a)
    eftRound(hi, lo, rounding)
end

function (sqrt){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = accSqrtForSign(a)
    eftRound(hi, lo, rounding)
end


text:
"""
$(TYPEDEF)

Direction observation information of a `Point2` variable.
"""
mutable struct PriorPoint2{T <: IIF.SamplableBelief} <: IncrementalInference.AbstractPrior
  Z::T
  # empty constructor
  # PriorPoint2{T}() where T = new()
  # regular constructor
  # PriorPoint2{T}(dist::T) where {T <: IIF.SamplableBelief} = new{T}(dist)
end

# convenience helper and default object
PriorPoint2() = PriorPoint2(MvNormal(zeros(2),LinearAlgebra.diagm([0.01;0.01])))

DFG.getManifold(::PriorPoint2) = TranslationGroup(2)

function getSample(cfo::CalcFactor{<:PriorPoint2})
  return rand(cfo.factor.Z)
end

function (cf::CalcFactor{<:PriorPoint2})(meas,
                                        X1)
  #
  return meas[1:2] .- X1[1:2]
end

"""
$(TYPEDEF)
"""
mutable struct Point2Point2{D <: IIF.SamplableBelief} <: AbstractRelativeRoots
    Zij::D
    # empty constructor
    Point2Point2{T}() where T = new{T}()
    # regular constructor
    Point2Point2{T}(x::T) where {T <: IIF.SamplableBelief} = new{T}(x)
end
# convenience and default object helper
Point2Point2(x::T=MvNormal(zeros(2),LinearAlgebra.diagm([0.1;0.1]))) where {T <: IIF.SamplableBelief} = Point2Point2{T}(x)

function getSample(cfo::CalcFactor{<:Point2Point2})
  return rand(cfo.factor.Zij)
end
function (pp2r::CalcFactor{<:Point2Point2})(meas,
                                            xi,
                                            xj )
  #
  return meas[1:2] .- (xj[1:2] .- xi[1:2])
end

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

"""
$(TYPEDEF)

Serialization type for `PriorPoint2`.
"""
mutable struct PackedPriorPoint2  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPoint2() = new()
    PackedPriorPoint2(x::String) = new(x)
end

function convert(::Type{PriorPoint2}, d::PackedPriorPoint2)
  # Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  distr = convert(SamplableBelief, d.str)
  return PriorPoint2{typeof(distr)}(distr)
end
function convert(::Type{PackedPriorPoint2}, d::PriorPoint2)
  # v2 = d.mv.Σ.mat[:];
  return PackedPriorPoint2(convert(PackedSamplableBelief, d.Z))
end

"""
$(TYPEDEF)

Serialization type for `Point2Point2`.
"""
mutable struct PackedPoint2Point2 <: IncrementalInference.PackedInferenceType
    str::String
    PackedPoint2Point2() = new()
    PackedPoint2Point2(s::AS) where {AS <: AbstractString} = new(s)
end
function convert(::Type{Point2Point2}, d::PackedPoint2Point2)
  return Point2Point2( convert(SamplableBelief, d.str) )
end
function convert(::Type{PackedPoint2Point2}, d::Point2Point2)
  return PackedPoint2Point2( convert(PackedSamplableBelief, d.Zij) )
end


text:
using Documenter, DataUtils

makedocs(
    debug = true,
    strict = false,
    doctest = false,
    format = :html,
    sitename = "DataUtils.jl",
    modules = [DataUtils],
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo = "github.com/jlapeyre/DataUtils.jl.git",
    target = "build",
    julia  = "nightly",
    deps = nothing,
    make = nothing
)


text:
#import Pkg; Pkg.add("HDF5"); Pkg.add("StaticArrays"); Pkg.add("PyPlot"); Pkg.add("DifferentialEquations"); Pkg.add("Parameters"); Pkg.add("Sundials");

using HDF5
using StaticArrays
using Statistics
using LinearAlgebra
using .Threads
using Serialization
using Random
using PyPlot
using Parameters
using BenchmarkTools

using AstroChemistry

const N = 3
const Nbin = 200
const xmin = 16.
const xmax = 23.

xl = "NH [cm^-2]"

const xmin = 16.
const xmax = 23.

function runPDR()

    T = Float64

    all_species::Vector{String} =
    [
    "H", "H-", "H2", "H+", "H2+", "H3+", "e-", "He", "He+", "HeH+",
    "C", "C+", "CO", "HCO+",
    "O", "O+", "OH", "OH+", "H2O+", "H3O+", "H2O",
    "O2", "CO+", "O2+",
    "CH2", "CH2+", "CH", "CH+", "CH3+",
    "Si+", "Si"
    #"SiO", "SiO+", "SiOH+",
    #"S", "S+",
    ]

    net, dict = initialize_chemistry_network(all_species, grRec=false)
    @unpack iH2, iCO, iC, fac_H, fac_C, fac_O, charge = net
    abtot = AbundTotal(abC=1e-4, abO=3e-4, abSi=0.0)
    @unpack abC, abO, abSi = abtot

    nH = 1000.
    temp = 50.0
    ξ = 0.5e-16 #H2
    IUV = 10.0
    Zp = 1.0
    dtime = 1e10

    N_spec = length(all_species)
    @show N_spec
    ab_vs_x = Vector{Vector{T}}(undef, Nbin)
    for i in 1:Nbin
        ab_vs_x[i] = zeros(N_spec)
    end
    d_lnNH = (xmax-xmin) / (Nbin-1)
    xbin = 10 .^ (xmin .+ collect(0:Nbin-1) .* d_lnNH)
    NHbin = xbin
    dNH = NHbin[2:Nbin] .- NHbin[1:Nbin-1]

    NH2bin = zeros(Nbin)
    NCObin = zeros(Nbin)
    NCbin = zeros(Nbin)
    N_reac = length(net.alpha)
    reaction_rates=zeros(Nbin, N_reac)

    #@time for i in 1:1
    @time for i in eachindex(NHbin)
        xH2 = getindex.(ab_vs_x, iH2)
        xCO = getindex.(ab_vs_x, iCO)
        xC  = getindex.(ab_vs_x, iC )
        NH = NHbin[i]
        NH2 = NCO = NC = 1e10
        if i > 1
            NH2 = sum(xH2[1:i-1] .* dNH[1:i-1])
            NCO = sum(xCO[1:i-1] .* dNH[1:i-1])
            NC  = sum(xC[1:i-1]  .* dNH[1:i-1])
        end
        NH2bin[i] = NH2
        NCObin[i] = NCO
        NCbin[i] = NC
        #@show i, NH, NH2,NCO, NC
        #print(i," ")

        G = 2.8e-5 # for Zp = 1
        αG = IUV * AstroChemistry.kH2diss / AstroChemistry.kdust / nH * G
        #i==1 ? println("αG/2 = ", αG/2) : nothing

        xneq = SVector{0,T}([])
        #par = Par{1,0,T}(nH, temp, ξ, IUV, Zp, NH=SVector{1,T}(NH), NH2=SVector{1,T}(NH2), NCO=SVector{1,T}(NCO), NC=SVector{1,T}(NC), xneq=xneq)
        par = Par{1,0,T}(nH=nH, temp=temp, ξ=ξ, IUV=IUV, Zp=Zp,
        NH=[NH], NH2=[NH2], NCO=[NCO], NC=[NC], xneq=xneq)
        #NH=SVector{1,T}(NH), NH2=SVector{1,T}(NH2), NCO=SVector{1,T}(NCO), NC=SVector{1,T}(NC), xneq=xneq)

        retcode, reaction_rates[i,:] = solve_equilibrium_abundances(ab_vs_x[i], dtime, par, abtot, net)

        calc_abund_derived(ab_vs_x[i], xneq, abtot, net)
        sumH  = sum( ab_vs_x[i] .* fac_H )
        sumC  = sum( ab_vs_x[i] .* fac_C )
        sumO  = sum( ab_vs_x[i] .* fac_O )
        #sumSi = sum( ab_vs_x[i] .* fac_Si )
        #sumS  = sum( ab_vs_x[i] .* fac_S )
        sumelec = sum( ab_vs_x[i] .* charge )
        (1.0 ≈ sumH)  ? nothing : error("sumH = " , sumH)
        (abC  ≈ sumC)  ? nothing : error("sumC = " , sumC)
        (abO  ≈ sumO)  ? nothing : error("sumO = " , sumO)
        #(abSi ≈ sumSi) ? nothing : error("sumSi = ", sumSi)
        #(abS  ≈ sumS)  ? nothing : error("sumS = " , sumS)
        (1.0 ≈ 1.0 + sumelec) ? nothing : error("sumelec = ",sumelec)
    end

    return ab_vs_x, NHbin, NH2bin, NCObin, NCbin, reaction_rates, dict, net
end

ab_vs_x, N_H, N_H2, N_CO, N_C, rr, dict, net = runPDR();

xbin = 10 .^ (xmin .+ (xmax-xmin) .* collect(0:Nbin-1) ./ (Nbin-1));

clf()
fig, ax = PyPlot.subplots(1, 3, figsize=(15,5))

#abC_tot  = [sum(ab_vs_x[i] .* fac_C)  for i in 1:Nbin]
#abO_tot  = [sum(ab_vs_x[i] .* fac_O)  for i in 1:Nbin]
#abSi_tot = [sum(ab_vs_x[i] .* fac_Si) for i in 1:Nbin]
#abS_tot  = [sum(ab_vs_x[i] .* fac_S)  for i in 1:Nbin]

xminp, xmaxp = 20., 23.
ax[1].plot(xbin, 2 .*getindex.(ab_vs_x, dict["H2"]), "-"  , label="H2")
ax[1].plot(xbin, getindex.(ab_vs_x, dict["H"]), "--" , label="H")
ax[1].plot(xbin, getindex.(ab_vs_x, dict["e-"]), ":"  , label="e-")
ax[1].legend(loc="best", fontsize=9, ncol=1, frameon=false)
ax[1].set_xscale("log")
ax[1].set_yscale("log")
ax[1].axis([10.0^xminp, 10.0^xmaxp, 1e-6, 2])
ax[1].set_xlabel(xl)
ax[1].set_ylabel("x_i")
ax[1].grid(linestyle=":", linewidth=1)
ax2 = ax[1].twiny()
Av = N_H .* 5.35e-22
ax2.set_xlabel("Av")
ax2.set_xlim((10.0^xminp, 10.0^xmaxp).*5.35e-22)
ax2.set_xscale("log")
ax2.set_yscale("log")

xminp, xmaxp = 20., 23.
ax[2].plot(xbin, getindex.(ab_vs_x, dict["H+"]), "-" , label="H+")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["H3+"]) , "--" , label="H3+")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["He+"]) , "-." , label="He+")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["O"])  , ":"  , label="O")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["OH"]) , "-" , label="OH")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["O2"]) , "--" , label="O2")
ax[2].plot(xbin, getindex.(ab_vs_x, dict["H2O"]), "-." , label="H2O")
ax[2].legend(loc="best", fontsize=9, ncol=2, frameon=false)
ax[2].set_xscale("log")
ax[2].set_yscale("log")
ax[2].axis([10.0^xminp, 10.0^xmaxp, 1e-11, 1e-3])
ax[2].set_xlabel(xl)
ax[2].grid(linestyle=":", linewidth=1)
ax2 = ax[2].twiny()
Av = N_H .* 5.35e-22
ax2.set_xlabel("Av")
ax2.set_xlim((10.0^xminp, 10.0^xmaxp).*5.35e-22)
ax2.set_xscale("log")
ax2.set_yscale("log")

ax[3].plot(xbin, getindex.(ab_vs_x, dict["C"])  , "-"  , label="C")
ax[3].plot(xbin, getindex.(ab_vs_x, dict["CO"]) , "--" , label="CO")
ax[3].plot(xbin, getindex.(ab_vs_x, dict["C+"]) , "-." , label="C+")
ax[3].plot(xbin, getindex.(ab_vs_x, dict["CH"]) , ":" , label="CH")
ax[3].plot(xbin, getindex.(ab_vs_x, dict["CH2"]) , ":" , label="CH2")
ax[3].plot(xbin, getindex.(ab_vs_x, dict["HCO+"]) , ":" , label="HCO+")
ax[3].legend(loc="upper left", fontsize=10, ncol=1, frameon=false)
ax[3].set_xscale("log")
ax[3].set_yscale("log")
ax[3].axis([10.0^xminp, 10.0^xmaxp, 1e-11, 1e-3])
ax[3].set_xlabel(xl)
ax[3].grid(linestyle=":", linewidth=1)
ax2 = ax[3].twiny()
Av = N_H .* 5.35e-22
ax2.set_xlabel("Av")
ax2.set_xlim((10.0^xminp, 10.0^xmaxp).*5.35e-22)
ax2.set_xscale("log")
ax2.set_yscale("log")

tight_layout()
savefig("PDRtest.png")


text:
module WatchJuliaBurn

export c╯°□°ↄ╯, 🗺, 📖
export ✅, 👍, 👌, 👎
export 😻, ⬇😻, ➡😻, 😹, ⬇😹, ➡😹, 🐈, ⬇🐈, ➡🐈
export 🥧, 🍰, ㅠ

## Base
const c╯°□°ↄ╯ = throw
const 🗺 = map
const 📖 = Dict
const ✅ = true
const 👍 = true
const 👌 = true
const 👎 = false

## Arrays
const 😻 = cat
const ⬇😻 = vcat
const ➡😻 = hcat
const 😹 = cat
const ⬇😹 = vcat
const ➡😹 = hcat
const 🐈 = cat
const ⬇🐈 = vcat
const ➡🐈 = hcat

## Math
const 🥧 = π
const 🍰 = π
const ㅠ = π

end


text:
# Very simple IMU accelerometer simulation
mutable struct IMUSim
    σa::Float64
end

function ComputeMeasurement(imu::IMUSim, acc; type = :true)
    if type == :true
        meas = @SVector [acc[1] + imu.σa*randn(),
                         acc[2] + imu.σa*randn(),
                         acc[3] + imu.σa*randn()]
    else
        meas = @SVector [acc[1], acc[2], acc[3]]
    end

    return meas
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
const NULL = 0
const FANN_ERRSTR_MAX = 128

# Skipping MacroDefinition: fann_linear_func ( v1 , r1 , v2 , r2 , sum ) ( ( ( ( ( r2 ) - ( r1 ) ) * ( ( sum ) - ( v1 ) ) ) / ( ( v2 ) - ( v1 ) ) ) + ( r1 ) )
# Skipping MacroDefinition: fann_stepwise ( v1 , v2 , v3 , v4 , v5 , v6 , r1 , r2 , r3 , r4 , r5 , r6 , min , max , sum ) ( sum < v5 ? ( sum < v3 ? ( sum < v2 ? ( sum < v1 ? min : fann_linear_func ( v1 , r1 , v2 , r2 , sum ) ) : fann_linear_func ( v2 , r2 , v3 , r3 , sum ) ) : ( sum < v4 ? fann_linear_func ( v3 , r3 , v4 , r4 , sum ) : fann_linear_func ( v4 , r4 , v5 , r5 , sum ) ) ) : ( sum < v6 ? fann_linear_func ( v5 , r5 , v6 , r6 , sum ) : max ) )
# Skipping MacroDefinition: fann_linear_derive ( steepness , value ) ( steepness )
# Skipping MacroDefinition: fann_sigmoid_real ( sum ) ( 1.0f / ( 1.0f + exp ( - 2.0f * sum ) ) )
# Skipping MacroDefinition: fann_sigmoid_derive ( steepness , value ) ( 2.0f * steepness * value * ( 1.0f - value ) )
# Skipping MacroDefinition: fann_sigmoid_symmetric_real ( sum ) ( 2.0f / ( 1.0f + exp ( - 2.0f * sum ) ) - 1.0f )
# Skipping MacroDefinition: fann_sigmoid_symmetric_derive ( steepness , value ) steepness * ( 1.0f - ( value * value ) )
# Skipping MacroDefinition: fann_gaussian_real ( sum ) ( exp ( - sum * sum ) )
# Skipping MacroDefinition: fann_gaussian_derive ( steepness , value , sum ) ( - 2.0f * sum * value * steepness * steepness )
# Skipping MacroDefinition: fann_gaussian_symmetric_real ( sum ) ( ( exp ( - sum * sum ) * 2.0f ) - 1.0f )
# Skipping MacroDefinition: fann_gaussian_symmetric_derive ( steepness , value , sum ) ( - 2.0f * sum * ( value + 1.0f ) * steepness * steepness )
# Skipping MacroDefinition: fann_elliot_real ( sum ) ( ( ( sum ) / 2.0f ) / ( 1.0f + fann_abs ( sum ) ) + 0.5f )
# Skipping MacroDefinition: fann_elliot_derive ( steepness , value , sum ) ( steepness * 1.0f / ( 2.0f * ( 1.0f + fann_abs ( sum ) ) * ( 1.0f + fann_abs ( sum ) ) ) )
# Skipping MacroDefinition: fann_elliot_symmetric_real ( sum ) ( ( sum ) / ( 1.0f + fann_abs ( sum ) ) )
# Skipping MacroDefinition: fann_elliot_symmetric_derive ( steepness , value , sum ) ( steepness * 1.0f / ( ( 1.0f + fann_abs ( sum ) ) * ( 1.0f + fann_abs ( sum ) ) ) )
# Skipping MacroDefinition: fann_sin_symmetric_real ( sum ) ( sin ( sum ) )
# Skipping MacroDefinition: fann_sin_symmetric_derive ( steepness , sum ) ( steepness * cos ( steepness * sum ) )
# Skipping MacroDefinition: fann_cos_symmetric_real ( sum ) ( cos ( sum ) )
# Skipping MacroDefinition: fann_cos_symmetric_derive ( steepness , sum ) ( steepness * - sin ( steepness * sum ) )
# Skipping MacroDefinition: fann_sin_real ( sum ) ( sin ( sum ) / 2.0f + 0.5f )
# Skipping MacroDefinition: fann_sin_derive ( steepness , sum ) ( steepness * cos ( steepness * sum ) / 2.0f )
# Skipping MacroDefinition: fann_cos_real ( sum ) ( cos ( sum ) / 2.0f + 0.5f )
# Skipping MacroDefinition: fann_cos_derive ( steepness , sum ) ( steepness * - sin ( steepness * sum ) / 2.0f )
# Skipping MacroDefinition: fann_activation_switch ( activation_function , value , result ) switch ( activation_function ) \
# { case FANN_LINEAR : result = ( fann_type ) value ; break ; case FANN_LINEAR_PIECE : result = ( fann_type ) ( ( value < 0 ) ? 0 : ( value > 1 ) ? 1 : value ) ; break ; case FANN_LINEAR_PIECE_SYMMETRIC : result = ( fann_type ) ( ( value < - 1 ) ? - 1 : ( value > 1 ) ? 1 : value ) ; break ; case FANN_SIGMOID : result = ( fann_type ) fann_sigmoid_real ( value ) ; break ; case FANN_SIGMOID_SYMMETRIC : result = ( fann_type ) fann_sigmoid_symmetric_real ( value ) ; break ; case FANN_SIGMOID_SYMMETRIC_STEPWISE : result = ( fann_type ) fann_stepwise ( - 2.64665293693542480469e+00 , - 1.47221934795379638672e+00 , - 5.49306154251098632812e-01 , 5.49306154251098632812e-01 , 1.47221934795379638672e+00 , 2.64665293693542480469e+00 , - 9.90000009536743164062e-01 , - 8.99999976158142089844e-01 , - 5.00000000000000000000e-01 , 5.00000000000000000000e-01 , 8.99999976158142089844e-01 , 9.90000009536743164062e-01 , - 1 , 1 , value ) ; break ; case FANN_SIGMOID_STEPWISE : result = ( fann_type ) fann_stepwise ( - 2.64665246009826660156e+00 , - 1.47221946716308593750e+00 , - 5.49306154251098632812e-01 , 5.49306154251098632812e-01 , 1.47221934795379638672e+00 , 2.64665293693542480469e+00 , 4.99999988824129104614e-03 , 5.00000007450580596924e-02 , 2.50000000000000000000e-01 , 7.50000000000000000000e-01 , 9.49999988079071044922e-01 , 9.95000004768371582031e-01 , 0 , 1 , value ) ; break ; case FANN_THRESHOLD : result = ( fann_type ) ( ( value < 0 ) ? 0 : 1 ) ; break ; case FANN_THRESHOLD_SYMMETRIC : result = ( fann_type ) ( ( value < 0 ) ? - 1 : 1 ) ; break ; case FANN_GAUSSIAN : result = ( fann_type ) fann_gaussian_real ( value ) ; break ; case FANN_GAUSSIAN_SYMMETRIC : result = ( fann_type ) fann_gaussian_symmetric_real ( value ) ; break ; case FANN_ELLIOT : result = ( fann_type ) fann_elliot_real ( value ) ; break ; case FANN_ELLIOT_SYMMETRIC : result = ( fann_type ) fann_elliot_symmetric_real ( value ) ; break ; case FANN_SIN_SYMMETRIC : result = ( fann_type ) fann_sin_symmetric_real ( value ) ; break ; case FANN_COS_SYMMETRIC : result = ( fann_type ) fann_cos_symmetric_real ( value ) ; break ; case FANN_SIN : result = ( fann_type ) fann_sin_real ( value ) ; break ; case FANN_COS : result = ( fann_type ) fann_cos_real ( value ) ; break ; case FANN_GAUSSIAN_STEPWISE : result = 0 ; break ; \
#}

const FANN_FIX_VERSION = "FANN_FIX_2.0"
const FANN_FLO_VERSION = "FANN_FLO_2.1"
const FANN_CONF_VERSION = FANN_FLO_VERSION

# Skipping MacroDefinition: FANN_GET ( type , name ) FANN_EXTERNAL type FANN_API fann_get_ ## name ( struct fann * ann ) \
# { return ann -> name ; \
#}
# Skipping MacroDefinition: FANN_SET ( type , name ) FANN_EXTERNAL void FANN_API fann_set_ ## name ( struct fann * ann , type value ) \
# { ann -> name = value ; \
#}
# Skipping MacroDefinition: FANN_GET_SET ( type , name ) FANN_GET ( type , name ) FANN_SET ( type , name )
# Skipping MacroDefinition: fann_max ( x , y ) ( ( ( x ) > ( y ) ) ? ( x ) : ( y ) )
# Skipping MacroDefinition: fann_min ( x , y ) ( ( ( x ) < ( y ) ) ? ( x ) : ( y ) )
# Skipping MacroDefinition: fann_safe_free ( x ) { if ( x ) { free ( x ) ; x = NULL ; } }
# Skipping MacroDefinition: fann_clip ( x , lo , hi ) ( ( ( x ) < ( lo ) ) ? ( lo ) : ( ( ( x ) > ( hi ) ) ? ( hi ) : ( x ) ) )
# Skipping MacroDefinition: fann_exp2 ( x ) exp ( 0.69314718055994530942 * ( x ) )
# Skipping MacroDefinition: fann_rand ( min_value , max_value ) ( ( ( float ) ( min_value ) ) + ( ( ( float ) ( max_value ) - ( ( float ) ( min_value ) ) ) * rand ( ) / ( RAND_MAX + 1.0f ) ) )
# Skipping MacroDefinition: fann_abs ( value ) ( ( ( value ) > 0 ) ? ( value ) : - ( value ) )
# Skipping MacroDefinition: fann_mult ( x , y ) ( x * y )
# Skipping MacroDefinition: fann_div ( x , y ) ( x / y )
# Skipping MacroDefinition: fann_random_weight ( ) ( fann_rand ( - 0.1f , 0.1f ) )
# Skipping MacroDefinition: fann_random_bias_weight ( ) ( fann_rand ( - 0.1f , 0.1f ) )

# begin enum fann_errno_enum
typealias fann_errno_enum Uint32
const FANN_E_NO_ERROR = (uint32)(0)
const FANN_E_CANT_OPEN_CONFIG_R = (uint32)(1)
const FANN_E_CANT_OPEN_CONFIG_W = (uint32)(2)
const FANN_E_WRONG_CONFIG_VERSION = (uint32)(3)
const FANN_E_CANT_READ_CONFIG = (uint32)(4)
const FANN_E_CANT_READ_NEURON = (uint32)(5)
const FANN_E_CANT_READ_CONNECTIONS = (uint32)(6)
const FANN_E_WRONG_NUM_CONNECTIONS = (uint32)(7)
const FANN_E_CANT_OPEN_TD_W = (uint32)(8)
const FANN_E_CANT_OPEN_TD_R = (uint32)(9)
const FANN_E_CANT_READ_TD = (uint32)(10)
const FANN_E_CANT_ALLOCATE_MEM = (uint32)(11)
const FANN_E_CANT_TRAIN_ACTIVATION = (uint32)(12)
const FANN_E_CANT_USE_ACTIVATION = (uint32)(13)
const FANN_E_TRAIN_DATA_MISMATCH = (uint32)(14)
const FANN_E_CANT_USE_TRAIN_ALG = (uint32)(15)
const FANN_E_TRAIN_DATA_SUBSET = (uint32)(16)
const FANN_E_INDEX_OUT_OF_BOUND = (uint32)(17)
const FANN_E_SCALE_NOT_PRESENT = (uint32)(18)
const FANN_E_INPUT_NO_MATCH = (uint32)(19)
const FANN_E_OUTPUT_NO_MATCH = (uint32)(20)
const FANN_E_WRONG_PARAMETERS_FOR_CREATE = (uint32)(21)
# end enum fann_errno_enum

type fann_error
    errno_f::fann_errno_enum
    error_log::Ptr{FILE}
    errstr::Ptr{Uint8}
end

# begin enum fann_train_enum
typealias fann_train_enum Uint32
const FANN_TRAIN_INCREMENTAL = (uint32)(0)
const FANN_TRAIN_BATCH = (uint32)(1)
const FANN_TRAIN_RPROP = (uint32)(2)
const FANN_TRAIN_QUICKPROP = (uint32)(3)
const FANN_TRAIN_SARPROP = (uint32)(4)
# end enum fann_train_enum

# begin enum fann_activationfunc_enum
typealias fann_activationfunc_enum Uint32
const FANN_LINEAR = (uint32)(0)
const FANN_THRESHOLD = (uint32)(1)
const FANN_THRESHOLD_SYMMETRIC = (uint32)(2)
const FANN_SIGMOID = (uint32)(3)
const FANN_SIGMOID_STEPWISE = (uint32)(4)
const FANN_SIGMOID_SYMMETRIC = (uint32)(5)
const FANN_SIGMOID_SYMMETRIC_STEPWISE = (uint32)(6)
const FANN_GAUSSIAN = (uint32)(7)
const FANN_GAUSSIAN_SYMMETRIC = (uint32)(8)
const FANN_GAUSSIAN_STEPWISE = (uint32)(9)
const FANN_ELLIOT = (uint32)(10)
const FANN_ELLIOT_SYMMETRIC = (uint32)(11)
const FANN_LINEAR_PIECE = (uint32)(12)
const FANN_LINEAR_PIECE_SYMMETRIC = (uint32)(13)
const FANN_SIN_SYMMETRIC = (uint32)(14)
const FANN_COS_SYMMETRIC = (uint32)(15)
const FANN_SIN = (uint32)(16)
const FANN_COS = (uint32)(17)
# end enum fann_activationfunc_enum

# begin enum fann_errorfunc_enum
typealias fann_errorfunc_enum Uint32
const FANN_ERRORFUNC_LINEAR = (uint32)(0)
const FANN_ERRORFUNC_TANH = (uint32)(1)
# end enum fann_errorfunc_enum

# begin enum fann_stopfunc_enum
typealias fann_stopfunc_enum Uint32
const FANN_STOPFUNC_MSE = (uint32)(0)
const FANN_STOPFUNC_BIT = (uint32)(1)
# end enum fann_stopfunc_enum

# begin enum fann_nettype_enum
typealias fann_nettype_enum Uint32
const FANN_NETTYPE_LAYER = (uint32)(0)
const FANN_NETTYPE_SHORTCUT = (uint32)(1)
# end enum fann_nettype_enum

typealias fann_type Cdouble

type fann_neuron
    first_con::Uint32
    last_con::Uint32
    sum::fann_type
    value::fann_type
    activation_steepness::fann_type
    activation_function::fann_activationfunc_enum
end

type fann_layer
    first_neuron::Ptr{fann_neuron}
    last_neuron::Ptr{fann_neuron}
end

typealias fann_callback_type Ptr{Void}

type fann
    errno_f::fann_errno_enum
    error_log::Ptr{FILE}
    errstr::Ptr{Uint8}
    learning_rate::Cfloat
    learning_momentum::Cfloat
    connection_rate::Cfloat
    network_type::fann_nettype_enum
    first_layer::Ptr{fann_layer}
    last_layer::Ptr{fann_layer}
    total_neurons::Uint32
    num_input::Uint32
    num_output::Uint32
    weights::Ptr{fann_type}
    connections::Ptr{Ptr{fann_neuron}}
    train_errors::Ptr{fann_type}
    training_algorithm::fann_train_enum
    total_connections::Uint32
    output::Ptr{fann_type}
    num_MSE::Uint32
    MSE_value::Cfloat
    num_bit_fail::Uint32
    bit_fail_limit::fann_type
    train_error_function::fann_errorfunc_enum
    train_stop_function::fann_stopfunc_enum
    callback::fann_callback_type
    user_data::Ptr{Void}
    cascade_output_change_fraction::Cfloat
    cascade_output_stagnation_epochs::Uint32
    cascade_candidate_change_fraction::Cfloat
    cascade_candidate_stagnation_epochs::Uint32
    cascade_best_candidate::Uint32
    cascade_candidate_limit::fann_type
    cascade_weight_multiplier::fann_type
    cascade_max_out_epochs::Uint32
    cascade_max_cand_epochs::Uint32
    cascade_min_out_epochs::Uint32
    cascade_min_cand_epochs::Uint32
    cascade_activation_functions::Ptr{fann_activationfunc_enum}
    cascade_activation_functions_count::Uint32
    cascade_activation_steepnesses::Ptr{fann_type}
    cascade_activation_steepnesses_count::Uint32
    cascade_num_candidate_groups::Uint32
    cascade_candidate_scores::Ptr{fann_type}
    total_neurons_allocated::Uint32
    total_connections_allocated::Uint32
    quickprop_decay::Cfloat
    quickprop_mu::Cfloat
    rprop_increase_factor::Cfloat
    rprop_decrease_factor::Cfloat
    rprop_delta_min::Cfloat
    rprop_delta_max::Cfloat
    rprop_delta_zero::Cfloat
    sarprop_weight_decay_shift::Cfloat
    sarprop_step_error_threshold_factor::Cfloat
    sarprop_step_error_shift::Cfloat
    sarprop_temperature::Cfloat
    sarprop_epoch::Uint32
    train_slopes::Ptr{fann_type}
    prev_steps::Ptr{fann_type}
    prev_train_slopes::Ptr{fann_type}
    prev_weights_deltas::Ptr{fann_type}
    scale_mean_in::Ptr{Cfloat}
    scale_deviation_in::Ptr{Cfloat}
    scale_new_min_in::Ptr{Cfloat}
    scale_factor_in::Ptr{Cfloat}
    scale_mean_out::Ptr{Cfloat}
    scale_deviation_out::Ptr{Cfloat}
    scale_new_min_out::Ptr{Cfloat}
    scale_factor_out::Ptr{Cfloat}
end

type fann_train_data
    errno_f::fann_errno_enum
    error_log::Ptr{FILE}
    errstr::Ptr{Uint8}
    num_data::Uint32
    num_input::Uint32
    num_output::Uint32
    input::Ptr{Ptr{fann_type}}
    output::Ptr{Ptr{fann_type}}
end

type fann_connection
    from_neuron::Uint32
    to_neuron::Uint32
    weight::fann_type
end

const FANNPRINTF = "%.20e"
const FANNSCANF = "%le"


text:
using Carretera
using FactCheck
using JLD

# Tests de las funciones Auxiliares
include("test-funciones-aux.jl")

FactCheck.exitstatus()


text:
"""
    struct ConstrainedUncertainScalarValueOneParameter{S, T1 <: Number}
        distribution::Distribution{Univariate, S}
        a::T1
    end

A constrained uncertain value represented by a one-parameter distribution,
where the original distribution has been truncated.

## Fields:
- **`distribution`**: The truncated version of the original distribution.
- **`a`**: The original value of the parameter of the original distribution.
"""
struct ConstrainedUncertainScalarValueOneParameter{S,
        T1 <: Number} <: AbstractUncertainOneParameterScalarValue{S, T1}
    distribution::Distribution{Univariate, S}
    a::T1
end

"""
    struct ConstrainedUncertainScalarValueTwoParameter{S, T1 <: Number, T2 <: Number}
        distribution::Distribution{Univariate, S}
        a::T1
        b::T2
    end

A constrained uncertain value represented by a two-parameter distribution,
where the original distribution has been truncated.

## Fields:
- **`distribution`**: The truncated version of the original distribution.
- **`a`**: The original value of the first parameter of the original distribution.
- **`b`**: The original value of the second parameter of the original distribution.
"""
struct ConstrainedUncertainScalarValueTwoParameter{S, T1 <: Number,
        T2 <: Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    a::T1
    b::T2
end

"""
    struct ConstrainedUncertainScalarValueThreeParameter{S, T1 <: Number, T2 <: Number, T3 <: Number}
        distribution::Distribution{Univariate, S}
        a::T1
        b::T2
        c::T3
    end

A constrained uncertain value represented by a two-parameter distribution,
where the original distribution has been truncated.

## Fields:
- **`distribution`**: The truncated version of the original distribution.
- **`a`**: The original value of the first parameter of the original distribution.
- **`b`**: The original value of the second parameter of the original distribution.
- **`c`**: The original value of the third parameter of the original distribution.

"""
struct ConstrainedUncertainScalarValueThreeParameter{S, T1 <: Number, T2 <: Number, T3 <: Number} <: AbstractUncertainThreeParameterScalarValue{S, T1, T2, T3}
    distribution::Distribution{Univariate, S}
    a::T1
    b::T2
    c::T3
end

import Distributions.Normal
import Distributions.Uniform
import Distributions.Beta
import Distributions.BetaPrime
import Distributions.BetaBinomial
import Distributions.Gamma
import Distributions.Frechet

"""
Uncertain value represented by a generic three-parameter distribution.
"""
struct UncertainScalarTheoreticalThreeParameter{S<:ValueSupport, T1<:Number, T2<:Number, T3<:Number} <: AbstractUncertainThreeParameterScalarValue{S, T1, T2, T3}
    distribution::Distribution{Univariate, S}
    a::T1
    b::T2
    c::T3
end

"""
Uncertain value represented by a generic two-parameter distribution.
"""
struct UncertainScalarTheoreticalTwoParameter{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    a::T1
    b::T2
end

"""
Uncertain value represented by a generic one-parameter distribution.
"""
struct UncertainScalarGenericOneParameter{S<:ValueSupport, T1<:Number} <: AbstractUncertainOneParameterScalarValue{S, T1}
    distribution::Distribution{Univariate, S}
    a::T1
end

"""
Uncertain value represented by a normal distribution.
"""
struct UncertainScalarNormallyDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    μ::T1
    σ::T2
end

"""
Uncertain value represented by a uniform distribution.
"""
struct UncertainScalarUniformlyDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    lower::T1
    upper::T2
end

"""
Uncertain value represented by a beta distribution.
"""
struct UncertainScalarBetaDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    α::T1
    β::T2
end

"""
Uncertain value represented by a beta prime distribution.
"""
struct UncertainScalarBetaPrimeDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    α::T1
    β::T2
end

"""
Uncertain value represented by a beta binomial distribution.
"""
struct UncertainScalarBetaBinomialDistributed{S<:ValueSupport, T1<:Number, T2<:Number, T3<:Number} <: AbstractUncertainThreeParameterScalarValue{S, T1, T2, T3}
    distribution::Distribution{Univariate, S}
    n::T1
    α::T2
    β::T3
end

"""
Uncertain value represented by a gamma distribution.
"""
struct UncertainScalarGammaDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    α::T1
    θ::T2
end

"""
Uncertain value represented by a Fréchet distribution.
"""
struct UncertainScalarFrechetDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    α::T1
    θ::T2
end

"""
Uncertain value represented by a binomial distribution.
"""
struct UncertainScalarBinomialDistributed{S<:ValueSupport, T1<:Number, T2<:Number} <: AbstractUncertainTwoParameterScalarValue{S, T1, T2}
    distribution::Distribution{Univariate, S}
    n::T1
    p::T2
end

###################
# Pretty printing
###################
function summarise(o::AbstractUncertainTwoParameterScalarValue)
    a = @sprintf "%.3f" o.a
    b = @sprintf "%.3f" o.b
    dist = o.distribution
    _type = typeof(o)
    "$_type($a, $b, $dist)"
end
Base.show(io::IO, q::AbstractUncertainTwoParameterScalarValue) = print(io, summarise(q))

function summarise(o::AbstractUncertainThreeParameterScalarValue)
    a = @sprintf "%.3f" o.a
    b = @sprintf "%.3f" o.b
    c = @sprintf "%.3f" o.c
    dist = o.distribution
    _type = typeof(o)
    "$_type($a, $b, $c, $dist)"
end
Base.show(io::IO, q::AbstractUncertainThreeParameterScalarValue) = print(io, summarise(q))

function summarise(o::UncertainScalarNormallyDistributed)
    μ = @sprintf "%.3f" o.μ
    σ = @sprintf "%.3f" o.σ
    dist = o.distribution
    _type = typeof(o)
    "$_type(μ = $μ, σ = $σ)"
end

Base.show(io::IO, q::UncertainScalarNormallyDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarUniformlyDistributed)
    lower = @sprintf "%.3f" o.lower
    upper = @sprintf "%.3f" o.upper
    dist = o.distribution
    _type = typeof(o)
    "$_type(lower = $lower, upper = $upper)"
end

Base.show(io::IO, q::UncertainScalarUniformlyDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarBetaDistributed)
    α = @sprintf "%.3f" o.α
    β = @sprintf "%.3f" o.β
    dist = o.distribution
    _type = typeof(o)
    "$_type(α = $α, β = $β)"
end

Base.show(io::IO, q::UncertainScalarBetaDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarBetaPrimeDistributed)
    α = @sprintf "%.3f" o.α
    β = @sprintf "%.3f" o.β
    dist = o.distribution
    _type = typeof(o)
    "$_type(α = $α, β = $β)"
end

Base.show(io::IO, q::UncertainScalarBetaPrimeDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarBetaBinomialDistributed)
    n = @sprintf "%.3f" o.n
    α = @sprintf "%.3f" o.α
    β = @sprintf "%.3f" o.β
    dist = o.distribution
    _type = typeof(o)
    "$_type(n = $n, α = $α, β = $β)"
end

Base.show(io::IO, q::UncertainScalarBetaBinomialDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarGammaDistributed)
    α = @sprintf "%.3f" o.α
    θ = @sprintf "%.3f" o.θ
    dist = o.distribution
    _type = typeof(o)
    "$_type(α = $α, θ = $θ)"
end

Base.show(io::IO, q::UncertainScalarGammaDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarFrechetDistributed)
    α = @sprintf "%.3f" o.α
    θ = @sprintf "%.3f" o.θ
    dist = o.distribution
    _type = typeof(o)
    "$_type(α = $α, θ = $θ)"
end

Base.show(io::IO, q::UncertainScalarFrechetDistributed) = print(io, summarise(q))

function summarise(o::UncertainScalarBinomialDistributed)
    n = @sprintf "%.3f" o.n
    p = @sprintf "%.3f" o.p
    dist = o.distribution
    _type = typeof(o)
    "$_type(n = $n, p = $p)"
end

Base.show(io::IO, q::UncertainScalarBinomialDistributed) = print(io, summarise(q))

export
TheoreticalDistributionScalarValue,

AbstractUncertainOneParameterScalarValue,
AbstractUncertainTwoParameterScalarValue,
AbstractUncertainThreeParameterScalarValue,

ConstrainedUncertainScalarValueOneParameter,
ConstrainedUncertainScalarValueTwoParameter,
ConstrainedUncertainScalarValueThreeParameter,

UncertainScalarGenericOneParameter,
UncertainScalarTheoreticalTwoParameter,
UncertainScalarTheoreticalThreeParameter,
UncertainScalarNormallyDistributed,
UncertainScalarUniformlyDistributed,
UncertainScalarBetaDistributed,
UncertainScalarBetaPrimeDistributed,
UncertainScalarBetaBinomialDistributed,
UncertainScalarBinomialDistributed,
UncertainScalarGammaDistributed,
UncertainScalarFrechetDistributed,

Normal,
Uniform,
Beta,
BetaPrime,
BetaBinomial,
Gamma,
Frechet,
Binomial


text:
# SPDX-License-Identifier: X11
# 2020-12-17
# Shift only (200pt)

function main()
  N = parse(Int, readline())
  A = parse.(Int, (split ∘ readline)())

  (println ∘ minimum)(getbitsuntilodd, A)
end

function getbitsuntilodd(n::Integer)::Integer
  i = 0
  while n & 1 == 0
    n >>= 1
    i += 1
  end
  return i
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# Reactive message-passing in a linear Gaussian dynamical system
# Experiment to infer states
# Termination based on local Free Energy
#
# Wouter Kouw, BIASlab
# 26-03-2020

using Random
using Distributions
using LightGraphs
using MetaGraphs
using CPUTime

# Factor graph components
include(joinpath(@__DIR__, "../factor_nodes/factor_gaussian.jl"))
include(joinpath(@__DIR__, "../variables/var_gaussian.jl"))
include(joinpath(@__DIR__, "../variables/var_delta.jl"))
include(joinpath(@__DIR__, "../prob_operations.jl"))
include(joinpath(@__DIR__, "../graph_operations.jl"))
include(joinpath(@__DIR__, "../util.jl"))

# Data
include(joinpath(@__DIR__, "../gen_data.jl"))

"""
Experiment parameters
"""

# Signal time horizon
T = 50

# Reaction-time clock
TT = 10

# Free Energy threshold
fe_threshold = 1e-3

# Known transition and observation matrices
gain = 0.8
emission = 1.0

# Known noises (variance form)
process_noise = 2.0
measurement_noise = 0.5

# Parameters for state x_0
px0 = Normal(0.0, 1.0)

# Generate data
Random.seed!(256)
observed, hidden = gendata_LGDS(gain,
                                emission,
                                process_noise,
                                measurement_noise,
                                mean(px0),
                                inv(var(px0)),
                                time_horizon=T)

"""
Model/graph specification

This assumes the following type of models
p(y_{1:T}, x_{0:T}) = p(x_0) Π_t p(y_t | x_t) p(x_t | x_{t-1})

In other words, Markov chains of time-slices of a state-space models.
Below, we specify the following model through the time-slice subgraph
     _____      ___      ___
... (x_t-1)--->[g_t]--->(x_t)  ...
                          |
                         _|_
                        [f_t]
                          |
                         _|_
                       [[y_t]]

x_t-1 = previous state edge
g_t = state transition node
x_t = current state edge
f_t = likelihood node
y_t = observation node

---------------------------------
#TODO: Graph specification using a graphical user interface
"""

CPUtic()

# Start graph
global graph = MetaGraph(SimpleGraph(4*T+1))

let

	# Set initial state prior
	set_props!(graph, 1, Dict(:id => :x_0,
							  :time => 0,
							  :node => VarGaussian(:x_0, marginal=px0, time=0)))

    # Connect initial state prior to first state transition
    add_edge!(graph, (1, 2))

	# Start node numbering
	node_num = 1

	# Build whole graph (ie all t)
	for t = 1:T

		# State transition node
		props_gt = Dict(:id => :g_*t,
					    :time => t,
						:threshold => fe_threshold,
					    :node => FactorGaussian(:g_*t,
									            out=:x_*t,
		                                        mean=:x_*(t-1),
		                                        precision=inv(process_noise),
		                                        transition=gain,
												threshold=fe_threshold))

		# Current state
		props_xt = Dict(:id => :x_*t,
					    :time => t,
					    :node => VarGaussian(:x_*t, time=t))

		# Observation likelihood node
		props_ft = Dict(:id => :f_*t,
					    :time => t,
					    :node => FactorGaussian(:f_*t,
			 						            out=:y_*t,
			 							        mean=:x_*t,
			 							        precision=inv(measurement_noise),
			 							        transition=emission,
			 							        time=t,
												threshold=fe_threshold))

		# Observation
		props_yt = Dict(:id => :y_*t,
					    :time => t,
					    :node => VarDelta(:y_*t, observed[t], time=t))

		# Add properties to nodes in current time-slice
		set_props!(graph, node_num+1, props_gt)
		set_props!(graph, node_num+2, props_xt)
		set_props!(graph, node_num+3, props_ft)
		set_props!(graph, node_num+4, props_yt)

		# Add edges between factors and variables
		add_edge!(graph, (node_num-2, node_num+1)) # x_t-1 -- g_t
		add_edge!(graph, (node_num+1, node_num+2)) # g_t -- x_t
		add_edge!(graph, (node_num+2, node_num+3)) # x_t -- f_t
		add_edge!(graph, (node_num+3, node_num+4)) # f_t -- y_t

		# Increment node number
		node_num += 4
	end
end

# Ensure vertices can be recalled from given id
set_indexing_prop!(graph, :id)

"""
Inference: filtering
"""

# Preallocate energy
LFE = Vector{Array{Float64,1}}(undef,T)

# Start message routine
act!(graph, :x_0)

for t = 1:T

	# Report progress
	# if mod(t, T/5) == 1
  	println("At iteration "*string(t)*"/"*string(T))
	# end

	# Fire state transiton from previous time-point
	act!(graph, graph[graph[:g_*t, :id], :node], :x_*t)

	# Start clock for reactions
	tt = 0
	LFE_t = Float64[]
	nodes_fired = true
	while nodes_fired | (tt <= 3)

		# Re-start node fired check
		nodes_fired = false

		# Count iterations
		tt += 1

		# Iterate over all nodes to react
		for node_id in nodes_t(graph, t)

			# Retrieve object from node
			node = graph[graph[node_id, :id], :node]

			# Tell node to react
			react!(graph, node_id)

			# Check for state variable
			if typeof(node) == VarGaussian

				# Collect free energy
				push!(LFE_t, node.free_energy)
			end

			# Check for factor node
			if typeof(node) == FactorGaussian

				# Record whether node fired
				nodes_fired |= node.fired
			end
		end
	end

	# Record LFE
	LFE[t] = LFE_t
end

CPUtoc()

"""
Get estimates from graph
"""

# Preallocate state vectors
estimated_states = zeros(T,2)
F = zeros(T,)
A = zeros(T,)

# Loop over time
for t = 1:T

	# Loop over nodes in timeslice
	for node_id in nodes_t(graph, t)

		# Retrieve object from node
		node = graph[graph[node_id, :id], :node]

		# Check for state node
		if typeof(node) == VarGaussian

			# Extract moments of marginals
			estimated_states[t,1] = mean(node.marginal)
			estimated_states[t,2] = sqrt(var(node.marginal))

			# Track final FE
			F[t] = node.free_energy
		end

		# Check for state node
		if typeof(node) == VarDelta

			# Track final accuracy (logpdf of observation under message)
			A[t] = node.pred_error
		end
	end
end

"""
Visualize experimental results
"""

using Plots
pyplot()

# Plot estimated states
scatter(1:T, observed, color="black", label="observations")
plot!(hidden[2:end], color="red", label="true states")
plot!(estimated_states[:,1], color="blue", label="inferred")
plot!(estimated_states[:,1],
	  ribbon=[estimated_states[:,2], estimated_states[:,2]],
	  linewidth=2,
      color="blue",
      fillalpha=0.2,
      fillcolor="blue",
	  label="")
xlabel!("time (t)")
savefig(joinpath(@__DIR__, "viz/exp-lfeterm_state_estimates.png"))

# Plot free energy
plot(1:T, F, color="green", label="")
xlabel!("time (t)")
ylabel!("free energy (F)")
savefig(joinpath(@__DIR__, "viz/exp-lfeterm_fe-time.png"))

# Plot accuracy
plot(1:T, A, color="green", label="")
xlabel!("time (t)")
ylabel!("wpred err (-log μ(x))")
savefig(joinpath(@__DIR__, "viz/exp-lfeterm_acc-time.png"))

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
"""
# A Package for Solving Dynamic Rational Inattention Problems (DRIPs)
Copyright © 2020 Hassan Afrouzi, Choongryul Yang and Miguel Acosta

When using this package for your work, cite
[Afrouzi and Yang (2019)](http://www.afrouzi.com/dynamic_inattention.pdf)
for reference.
"""
module DRIPs

using LinearAlgebra, Random

include("drip_methods.jl")
include("trip_methods.jl")
include("irfs_methods.jl")
include("simulations_methods.jl")
include("aux_funcs.jl")

export Drip, Trip, irfs, simulate

end # module


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

import Base: replace
using DataStructures: Stack
using LightXML

" Parse wikipedia data dumped by wikimedia "
function parsewiki(xmlfile::AbstractString, txtfile::AbstractString)
    fin = open(xmlfile, "r")
    fout = open(txtfile, "w")

    page = ""
    numline = 0
    for line in eachline(fin)
        line = strip(line)
        numline += 1
        if line == "<page>"
            page = line
        elseif line == "</page>"
            page *= "\n" * line
            pageinfo = parsepage(page, numline)
            write(fout, "<doc id=\"$(pageinfo[:id])\" title=\"$(pageinfo[:title])\">\n",
                  "$(cleantext(pageinfo[:content]))\n</doc>\n")
        else
            page *= "\n" * line
        end
    end
    close(fin)
    close(fout)
end

" Parse one page of wikipedia "
function parsepage(page::AbstractString, numline::Int)
    pageinfo = Dict{Symbol, Any}()
    try
        xmldoc = parse_string(page)
        xroot = root(xmldoc)
        pageinfo[:ns] = content(xroot["ns"][1])
        pageinfo[:id] = content(xroot["id"][1])
        pageinfo[:title] = content(xroot["title"][1])
        if find_element(xroot, "revision") != nothing
            pageinfo[:content] = content(xroot["revision"][1]["text"][1])
        else
            pageinfo[:content] = content(xroot["text"][1])
        end
        free(xmldoc)
    catch e
        throw(e)
        println("Encounter errors before line $numline.")
    end
    pageinfo
end

" Clean unnecessary formats of contents "
function cleantext(text::AbstractString)
    text = clean1(text)
    text = clean2(text)
    text = clean3(text)
    text = clean4(text)
    text = clean5(text)
    #=
    str = replace(str, "<math", "</math>", " ")
    pair = search(str, "{{", "}}")
    bracefields = ["Infobox", "Taxobox", "Notability", "onesource", "More footnotes", "DEFAULTSORT"]
    while !isempty(pair)
        leftindex = collect(keys(pair))[1]
        rightindex = pop!(pair, leftindex) + length("}}") - 1
        for field in bracefields
            if startswith(str[leftindex + length("{{"):end], field)
                str = replace(str, leftindex, rightindex, "")
                update!(pair, leftindex, rightindex, rlen = 0)
                break
            end
        end
    end

    pair = search(str, "\\[\\[", "\\]\\]")
    bracketfields = ["Category"]
    while !isempty(pair)
        leftindex = collect(keys(pair))[1]
        rightindex = pop!(pair, leftindex) + length("]]") - 1
        for field in bracketfields
            if startswith(str[leftindex + length("[["):end], field)
                str = replace(str, leftindex, rightindex, "")
                update!(pair, leftindex, rightindex, rlen = 0)
                break
            end
        end
    end

    str = replace(str, Regex("{{reflist.*$") => "")
    str = replace(str, Regex("<references />.*$") => "")
    str = replace(str, "{\\| class", "\\|}", "")
    str = replace(str, "<ref", "</ref>", "")
    str = replace(str, "<ref", "/>", "")
    str = replace(str, "<!--", "-->", "")

    str = replace(str, r"\[\[([^\|]+)\]\]" => s"\1")
    str = replace(str, r"\n==[^\n]+==\n" => "")
    str = replace(str, r"\n===[^\n]+===\n" => "")
    str = replace(str, r"\n[,.:' ，。：]*\n", "\n")
    str = replace(str, r"^\n|\n$" => "")
    str = replace(str, r"'{2,}" => "")
    str = replace(str, r"\n.{0,50}\n"), "\n")
    =#
    return text
end

" Clean headings of page "
function clean1(text::AbstractString)
    text = replace(text, r"'''''(.*?)'''''" => s"\1")
    text = replace(text, r"'''(.*?)'''" => s"\1")
    text = replace(text, r"''(.*?)''" => s"\1")
    text = replace(text, r"======(.*?)======" => "")
    text = replace(text, r"=====(.*?)=====" => "")
    text = replace(text, r"====(.*?)====" => "")
    text = replace(text, r"===(.*?)===" => "")
    text = replace(text, r"==(.*?)==" => "")
    text = replace(text, r"----" => "")
    return text
end

" Clean html marks "
function clean2(text::AbstractString)
    text = replace(text, r"<math>(.*?)</math>"s => "")
    text = replace(text, r"<ref[ >](.*?)</ref>"s => "")
    text = replace(text, r"<!--(.*?)-->"s => "")
    return text
end

" Clean extra boxes "
function clean3(text::AbstractString)
    text = replace(text, r"^{{Taxobox(.*?)^}}"ms => "")
    text = replace(text, r"^{{Infobox(.*?)^}}"ms => "")
    text = replace(text, r"^{\| class(.*?)^\|}"ms => "")
    text = replace(text, r"{{reflist[}|](.*?)$"s => "")
    return text
end

" Clean external links "
function clean4(text::AbstractString)
    return text
end

" Clean other formats "
function clean5(text::AbstractString)
    text = replace(text, r"^\n+|\n+$" => "")
    text = replace(text, r"\n[,.:' ，。：]*\n", "\n")
    text = replace(text, r"\n+" => "\n")
    return text
end

" Search paired patterns and return matched offsets "
function search(s::AbstractString, left::AbstractString, right::AbstractString)
    pair = Dict{Int, Int}()
    match(Regex("$left"), s) === nothing && return pair
    stack = Stack{Int}()
    for m in eachmatch(Regex("($left|$right)"), s)
        if match(Regex("$left"), m.match) != nothing
            push!(stack, m.offset)
        elseif match(Regex("$right"), m.match) != nothing
            !isempty(stack) && push!(pair, pop!(stack) => m.offset)
        else
            println("Incorrect substring '$(m.match)' matching patterns ('$left', '$right')")
        end
    end
    return pair
end

" Replace `s[leftindex:rightindex]` by `r` "
function Base.replace(s::AbstractString, leftindex::Int, rightindex::Int, r::AbstractString)
    s[1:thisind(s, leftindex - 1)] * r * s[thisind(s, rightindex + 1):end]
end

" Replace all substrings matching `pat` with `r` "
function Base.replace(s::AbstractString, pat::Regex, r::AbstractString)
    while true
        match(pat, s) === nothing && break
        s = replace(s, pat => r)
    end
    return s
end

" Replace all substrings between `left` and `right` patterns with `r` "
function Base.replace(s::AbstractString, left::AbstractString, right::AbstractString, r::AbstractString)
    pair = search(s, left, right)
    rightlen = length(right) - count(c -> c == '\\', right)
    rlen = length(r)
    while !isempty(pair)
        leftindex = maximum(keys(pair))
        rightindex = pop!(pair, leftindex) + rightlen - 1
        s = replace(s, leftindex, rightindex, r)
        update!(pair, leftindex, rightindex, rlen = rlen)
    end
    return s
end

" Update paired dict after replacement "
function update!(pair::Dict{Int, Int}, leftindex::Int, rightindex::Int; rlen::Int = 0)
    plen = rightindex - leftindex + 1
    for (i, j) in sort(collect(pair))
        if i < leftindex && j < leftindex
            continue
        elseif i < leftindex && j > rightindex
            pair[i] = j - plen + rlen
        elseif i >= leftindex && j <= rightindex
            delete!(pair, i)
        elseif i > rightindex && j > rightindex
            delete!(pair, i)
            pair[i - plen + rlen] = j - plen + rlen
        else
            println("Mismatched paired pattern ($i, $j)")
        end
    end
end


text:
export TRLDLt_HO_vs_Nwt_λ

function TRLDLt_HO_vs_Nwt_λ(nlp 	   :: AbstractNLPModel,
              	       		nlpstop :: NLPStopping;
				   	        corr_ho :: Symbol = :Shamanskii,
							nwt_res_fact = 0.25,
							λfact = 1.0,
				   	        kwargs...
               	   	        )

	T = eltype(nlp.meta.x0)

	return TRARC(nlp,
				  nlpstop;
				  TR = TrustRegion(T(10.0)),
				  c = Combi(hessian_dense, PDataLDLt{eltype(nlp.meta.x0)}, (x, y, z) -> solve_modelTRDiag_HO_vs_Nwt_λ(x, y , z; ho_correction = corr_ho, nwt_res_fact = nwt_res_fact, λfact = λfact), preprocessLDLt, decreaseFact, Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of default tables with Markdown.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Markdown" begin
    data_md = [1 md"**bold**"
               2 md"""# Title

                      Paragraph

                          code
                   """]

    expected = """
┌────────┬────────────────────────────────────────────────┐
│ Col. 1 │                                         Col. 2 │
├────────┼────────────────────────────────────────────────┤
│      1 │   bold                                         │
│      2 │   Title\\n  ≡≡≡≡≡≡≡\\n\\n  Paragraph\\n\\n     code │
└────────┴────────────────────────────────────────────────┘
"""

    result = pretty_table(String, data_md)
    @test result == expected

    expected = """
┌────────┬─────────────┐
│ Col. 1 │      Col. 2 │
├────────┼─────────────┤
│      1 │   bold      │
│      2 │   Title     │
│        │   ≡≡≡≡≡≡≡   │
│        │             │
│        │   Paragraph │
│        │             │
│        │      code   │
└────────┴─────────────┘
"""

    result = pretty_table(String, data_md, linebreaks = true)
    @test result == expected
end


text:
using LibSerialPort

function test_nonblocking_serial_loopback(sp::SerialPort)

    println("\n[TEST] Read any incoming data for ~1 second...")
    for i = 1:1000
        print(read(sp, String))
        sleep(0.001)
    end
    println()

    flush(sp, buffer=SP_BUF_BOTH)

    print("\n\n[TEST] Serial loopback - ")
    println("Send 100 short messages and read whatever comes back...")

    for i = 1:100
        write(sp, "Test message $i\n")
        print(read(sp, String))
        i == 100 && write(sp, "done")
    end

    # Allow up to 100 more ms to get the remaining data
    for i = 1:100
        data = read(sp, String)
        print(data)
        if occursin(data, "done")
            break
        end
        sleep(0.001)
    end

    println()

    flush(sp, buffer=SP_BUF_BOTH)
end

function test_readline(sp::SerialPort)

    println("\n[TEST] Read any incoming data for ~1 second...")
    for i = 1:1000
        print(read(sp, String))
        sleep(0.001)
    end
    println()

    flush(sp, buffer=SP_BUF_BOTH)

    print("\n\n[TEST] Serial loopback - ")
    println("Send 100 short messages and read whatever comes back...")

    for i = 1:100
        write(sp, "Test message $i\n")
        sleep(0.001)
        received_message = readuntil(sp, '\n') # same as readline(sp)
        println(chomp(received_message))

        # Trigger an EOF when done writing so readuntil doesn't hang forever
        if i == 100
            seteof(sp, true)
        end
    end
    reseteof(sp)

    flush(sp, buffer=SP_BUF_BOTH)
end

function test_high_level_api(args...)

    if length(args) != 2
        println("Usage: $(basename(@__FILE__)) port baudrate")
        println("Available ports:")
        list_ports()
        return
    end

    if !ispath(args[1])
        println("Not found: ", args[1])
        return
    end

    sp = open(args[1], parse(Int, args[2]))

    print_port_metadata(sp)
    print_port_settings(sp)

    test_nonblocking_serial_loopback(sp)
    test_readline(sp)

    close(sp)
    return
end

test_high_level_api(ARGS...)


text:
# julia ~/Jules/PNToyBox/scripts/all_pnml.jl
# include("all_pnml.jl")

#using PNMLTree
using PNML
#using ProfileView, BenchmarkTools
#using PrettyPrinting
using DataFrames, DataFramesMeta

function testpn(;topdir="/home/jeff/Projects/Resources/PetriNet/PNML",
                dir="examples")
    pnml = endswith(".pnml")
    df = DataFrame() # Collects data from tests.
    cd(joinpath(topdir, dir)) do
        for (root, dirs, files) in walkdir(".")
            for file in filter(pnml, files)
                # Display path to file.
                f = lstrip(joinpath(root, file),['.', '/'])
                print(f, " size = ", filesize(f), "\n")
                try
                    # Gather performance data on parsing into PNML.
                    stats = @timed parse_file(f)
                    push!(df, (file=f, size=filesize(f),
                               time=stats.time, bytes=stats.bytes, gctime=stats.gctime))
                    # Display the PnmlModel as a test of the
                    # parsing, creation and show() implementation.
                    @show stats.value
                catch e
                    if e isa PNML.PnmlException
                        @warn " failed: $e"
                    elseif e isa InterruptException
                        return
                    else
                        sprint(showerror, e)
                    end
                end
            end
        end
    end
    # Display the gathered data frame.
    display(sort(df, [:time]))
end