text:
#
# utils.jl --
#
# Useful methods and constants for InterProcessCommunication (IPC) package of
# Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of InterProcessCommunication.jl released under the MIT
# "expat" license.
#
# Copyright (C) 2016-2019, Éric Thiébaut
# (https://github.com/emmt/InterProcessCommunication.jl).
#

# A bit of magic for calling C-code:
Base.convert(::Type{_typeof_key_t}, key::Key) = key.value
Base.convert(::Type{T}, key::Key) where {T<:Integer} = convert(T, key.value)

# The short version of `show` if also used for string interpolation in
# scripts so that extending methods:
#     Base.convert(::Type{String}, obj::T)
#     Base.string(obj::T)
# is not necessary.
Base.show(io::IO, ::MIME"text/plain", arg::Key) = show(io, arg)
Base.show(io::IO, key::Key) =
    print(io, "IPC.Key(", string(key.value, base=10), ")")

"""

`IPC.PRIVATE` is a special IPC key (of type `IPC.Key`) which indicates that a
new (private) key should be created.

"""
const PRIVATE = Key(IPC_PRIVATE)

"""
Immutable type `IPC.Key` stores a System V IPC key.  The call:

```julia
IPC.Key(path, proj)
```

generates a System V IPC key from pathname `path` and a project identifier
`proj` (a single character).  The key is suitable for System V Inter-Process
Communication (IPC) facilities (message queues, semaphores and shared memory).
For instance:

```julia
key = IPC.Key(".", 'a')
```

The special IPC key [`IPC.PRIVATE`](@ref) is also available to indicate that a
new (private) key should be created.

"""
Key(path::AbstractString, proj::Union{Char,Integer}) =
    Key(path, convert(Cint, proj))

function Key(path::AbstractString, proj::Cint)
    1 ≤ proj ≤ 255 || throw_argument_error("`proj` must be in the range 1:255")
    key = ccall(:ftok, _typeof_key_t, (Cstring, Cint), path, proj)
    systemerror("ftok", key == -1)
    return Key(key)
end

"""

```julia
gettimeofday() -> tv
```

yields the current time as an instance of `IPC.TimeVal`.  The result can be
converted into a fractional number of seconds by calling `float(tv)`.

See also: [`IPC.TimeVal`](@ref), [`nanosleep`](@ref), [`clock_gettime`](@ref).

"""
function gettimeofday()
    tv = Ref(TimeVal(0, 0))
    # `gettimeofday` should not fail in this case
    ccall(:gettimeofday, Cint, (Ptr{TimeVal}, Ptr{Cvoid}), tv, C_NULL)
    return tv[]
end

"""

```julia
now(T)
```

yields the current time since the Epoch as an instance of type `T`, which can
be [`TimeVal`](@ref) or [`TimeSpec`](@ref).

See also: [`gettimeofday`](@ref),  [`time`](@ref),  [`clock_gettime`](@ref).

"""
now(::Type{TimeVal}) = gettimeofday()
now(::Type{TimeSpec}) = clock_gettime(CLOCK_REALTIME)

"""

```julia
time(T)
```

yields the current time since the Epoch as an instance of type `T`, which can
be [`TimeVal`](@ref) or [`TimeSpec`](@ref).

See also:  [`now`](@ref)[`gettimeofday`](@ref), [`clock_gettime`](@ref).

"""
Base.time(::Type{TimeVal}) = gettimeofday()
Base.time(::Type{TimeSpec}) = clock_gettime(CLOCK_REALTIME)

"""

```julia
nanosleep(t) -> rem
```

sleeps for `t` seconds with nanosecond precision and returns the remaining time
(in case of interrupts) as an instance of `IPC.TimeSpec`.  Argument can be a
(fractional) number of seconds or an instance of `IPC.TimeSpec` or
`IPC.TimeVal`.

The `sleep` method provided by Julia has only millisecond precision.

See also [`gettimeofday`](@ref), [`IPC.TimeSpec`](@ref) and
[`IPC.TimeVal`](@ref).

"""
nanosleep(t::Union{Real,TimeVal,Libc.TimeVal}) = nanosleep(TimeSpec(t))

function nanosleep(ts::TimeSpec)
    rem = Ref(TimeSpec(0,0))
    ccall(:nanosleep, Cint, (Ptr{TimeSpec}, Ptr{TimeSpec}), Ref(ts), rem)
    return rem[]
end

"""

```julia
clock_getres(id) -> ts
```

yields the resolution (precision) of the specified clock `id`. The result is an
instance of `IPC.TimeSpec`.  Clock identifier `id` can be `CLOCK_REALTIME` or
`CLOCK_MONOTONIC` (described in [`clock_gettime`](@ref)).

See also [`clock_gettime`](@ref), [`clock_settime`](@ref),
[`gettimeofday`](@ref), [`nanosleep`](@ref), [`IPC.TimeSpec`](@ref) and
[`IPC.TimeVal`](@ref).

""" clock_getres

@static if Sys.islinux()
    function clock_getres(id::Integer)
        res = Ref(TimeSpec(0,0))
        systemerror("clock_getres",
                    ccall(:clock_getres, Cint,
                          (_typeof_clockid_t, Ptr{TimeSpec}),
                          id, res) != SUCCESS)
        return res[]
    end
else
    # Assume the same resolution as gettimeofday which is used as a substitute
    # to clock_gettime.
    clock_getres(id::Integer) = TimeSpec(0, 1_000)
end

"""

```julia
clock_gettime(id) -> ts
```

yields the time of the specified clock `id`.  The result is an instance of
`IPC.TimeSpec`.  Clock identifier `id` can be one of:

* `CLOCK_REALTIME`: System-wide clock that measures real (i.e., wall-clock)
  time.  This clock is affected by discontinuous jumps in the system time
  (e.g., if the system administrator manually changes the clock), and by the
  incremental adjustments performed by `adjtime` and NTP.

* `CLOCK_MONOTONIC`: Clock that cannot be set and represents monotonic time
  since some unspecified starting point.  This clock is not affected by
  discontinuous jumps in the system time.

See also [`clock_getres`](@ref), [`clock_settime`](@ref),
[`gettimeofday`](@ref), [`nanosleep`](@ref), [`IPC.TimeSpec`](@ref) and
[`IPC.TimeVal`](@ref).

""" clock_gettime

@static if Sys.islinux()
    function clock_gettime(id::Integer)
        ts = Ref(TimeSpec(0,0))
        systemerror("clock_gettime",
                    ccall(:clock_gettime, Cint,
                          (_typeof_clockid_t, Ptr{TimeSpec}),
                          id, ts) != SUCCESS)
        return ts[]
    end
else
    # Use clock_gettime as a substitute.
    clock_gettime(id::Integer) = TimeSpec(gettimeofday())
end

@doc @doc(clock_gettime) CLOCK_REALTIME
@doc @doc(clock_gettime) CLOCK_MONOTONIC

"""

```julia
clock_settime(id, ts)
```

set the time of the specified clock `id` to `ts`.  Argument `ts` can be an
instance of `IPC.TimeSpec` or a number of seconds.  Clock identifier `id` can
be `CLOCK_REALTIME` or `CLOCK_MONOTONIC` (described in
[`clock_gettime`](@ref)).

See also [`clock_getres`](@ref), [`clock_gettime`](@ref),
[`gettimeofday`](@ref), [`nanosleep`](@ref), [`IPC.TimeSpec`](@ref) and
[`IPC.TimeVal`](@ref).

"""
clock_settime(id::Integer, t::Union{Real,TimeVal,Libc.TimeVal}) =
    clock_settime(id, TimeSpec(t))

clock_settime(id::Integer, ts::TimeSpec) =
    clock_settime(id, Ref(ts))

clock_settime(id::Integer, ts::Union{Ref{TimeSpec},Ptr{TimeSpec}}) =
    SUCCESS == ccall(:clock_settime, Cint, (_typeof_clockid_t, Ptr{TimeSpec}),
                     id, ts) || throw_system_error("clock_settime")

"""

`TimeConstraints` is the abstract type inherited by concrete types specifying
the kind of constraints to apply for the integer and fractional parts of a time
value.  There are two possibilities: [`Nonnegative`](@ref) if the fractional
part shall be nonnegative or [`SameSign`](@ref) if the fractional and integer
parts shall have the same sign.

"""
abstract type TimeConstraints end

"""

`Nonnegative` is a singleton type derived from [`TimeConstraints`](@ref) and
used to impose that the fractional part of a time value be nonnegative.  This
is what is assumed for normalized time values.

"""
struct Nonnegative <: TimeConstraints end

"""

`SameSign` is a singleton type derived from [`TimeConstraints`](@ref) and used
to impose that the fractional and integer parts of a time value have the same
sign.

"""
struct SameSign <: TimeConstraints end

const Floats = Union{AbstractFloat,AbstractIrrational,Rational}
const MILLISECONDS_PER_SECOND = 1_000
const MICROSECONDS_PER_SECOND = 1_000_000
const  NANOSECONDS_PER_SECOND = 1_000_000_000

# Most of the time manipulation methods assume signed fractional time field.
@assert fieldtype(TimeVal,  :usec) <: Signed
@assert fieldtype(TimeSpec, :nsec) <: Signed

"""

```julia
TimeSpec(sec, nsec)
```

yields an instance of `TimeSpec` for an integer number of seconds `sec` and an
integer number of nanoseconds `nsec` since the Epoch.

```julia
TimeSpec(sec)
```

yields an instance of `TimeSpec` for a, possibly fractional, number of seconds
`sec` since the Epoch.  Argument can also be an instance of [`TimeVal`](@ref).

Call `now(TimeSpec)` to get the current time as an instance of `TimeSpec`.

Addition and subtraction involving and instance of `TimeSpec` yield a
`TimeSpec` result.  For instance:

```julia
now(TimeSpec) + 3.4
```

yields a `TimeSpec` instance with the current time plus `3.4` seconds.

```julia
typemin(TimeSpec)
typemax(TimeSpec)
```

respectively yield the minimum and maximum normalized time values for an
instance of `TimeSpec`.

"""
TimeSpec(ts::TimeSpec) = ts
TimeSpec(secs::Integer) = TimeSpec(secs, 0)
TimeSpec(secs::Floats) = TimeSpec(splittime(_typeof_timespec_sec,
                                            _typeof_timespec_nsec, secs,
                                            NANOSECONDS_PER_SECOND)...)
TimeSpec(tv::Union{TimeVal,Libc.TimeVal}) = TimeSpec(tv.sec, tv.usec*1_000)

"""

```julia
TimeVal(sec, usec)
```

yields an instance of `TimeVal` for an integer number of seconds `sec` and an
integer number of microseconds `usec` since the Epoch.

```julia
TimeVal(sec)
```

yields an instance of `TimeVal` with a, possibly fractional, number of seconds
`sec` since the Epoch.  Argument can also be an instance of [`TimeSpec`](@ref).

Call `now(TimeVal)` to get the current time as an instance of `TimeVal`.

Addition and subtraction involving and instance of `TimeVal` yield a
`TimeVal` result.  For instance:

```julia
now(TimeVal) + 3.4
```

yields a `TimeVal` instance with the current time plus `3.4` seconds.

```julia
typemin(TimeVal)
typemax(TimeVal)
```

respectively yield the minimum and maximum normalized time values for an
instance of `TimeVal`.

"""
TimeVal(tv::TimeVal) = tv
TimeVal(tv::Libc.TimeVal) = TimeVal(tv.sec, tv.usec)
TimeVal(secs::Integer) = TimeVal(secs, 0)
TimeVal(secs::Floats) = TimeVal(splittime(_typeof_timeval_sec,
                                          _typeof_timeval_usec, secs,
                                          MICROSECONDS_PER_SECOND)...)
TimeVal(ts::TimeSpec) = begin
    usec, r = divrem(ts.nsec, 1_000)
    if r ≥ 500
        usec += one(usec)
    elseif r ≤ -500
        usec -= one(usec)
    end
    fixtime(TimeVal, ts.sec, usec)
end

Libc.TimeVal(ts::TimeSpec) = Libc.TimeVal(TimeVal(ts))
Libc.TimeVal(tv::TimeVal) = Libc.TimeVal(tv.sec, tv.usec)

Base.typemin(::Type{TimeSpec}) =
    TimeSpec(typemin(_typeof_timespec_sec), 0)
Base.typemax(::Type{TimeSpec}) =
    TimeSpec(typemax(_typeof_timespec_sec), NANOSECONDS_PER_SECOND - 1)

Base.typemin(::Type{TimeVal}) =
    TimeVal(typemin(_typeof_timeval_sec), 0)
Base.typemax(::Type{TimeVal}) =
    TimeVal(typemax(_typeof_timeval_sec), MICROSECONDS_PER_SECOND - 1)

Base.Float32(t::Union{TimeVal,TimeSpec}) = convert(Float32, t)
Base.Float64(t::Union{TimeVal,TimeSpec}) = convert(Float64, t)
Base.BigFloat(t::Union{TimeVal,TimeSpec}) = convert(BigFloat, t)
Base.float(t::Union{TimeVal,TimeSpec}) = convert(Float64, t)
Base.convert(::Type{T}, tv::TimeVal) where {T<:AbstractFloat} =
    T(tv.sec) + T(tv.usec)/T(MICROSECONDS_PER_SECOND)
Base.convert(::Type{T}, ts::TimeSpec) where {T<:AbstractFloat} =
    T(ts.sec) + T(ts.nsec)/T(NANOSECONDS_PER_SECOND)
Base.convert(::Type{TimeSpec}, arg::Union{Real,TimeSpec,TimeVal,Libc.TimeVal}) =
    TimeSpec(arg)
Base.convert(::Type{TimeVal}, arg::Union{Real,TimeSpec,TimeVal,Libc.TimeVal}) =
    TimeVal(arg)

"""

```julia
splittime(Ti, Tf, s, n, r=Nonnegative()) -> (ip::Ti, fp::Tf)
```

yields two integers, `ip` and `fp` with `abs(fp) ∈ [0,n-1[`, such that `ip +
fp/n ≈ s` (with a precision better than `1/n`) and imposing the constraints set
by `r`:

* if `r` is `Nonnegative()`, then `fp ≥ 0`;

* if `r` is `SameSign()`, then `ip` and `fp` have the same sign.

The multiplier `n` must be strictly positive.  An `InexactError` is thrown if
the value of `s` cannot be converted (e.g. it is a NaN or its magnitude is too
large).

See also [`fixtime`](@ref).

"""
function splittime(::Type{Ti}, ::Type{Tf}, secs::Floats,
                   n::Integer) where {Ti<:Integer, Tf<:Integer}
    splittime(Ti, Tf, secs, n, Nonnegative())
end

function splittime(::Type{Ti}, ::Type{Tf}, secs::Floats,
                   n::Integer, ::Nonnegative) where {Ti<:Integer, Tf<:Integer}
    s = floor(secs)
    ip = trunc(Ti, s)
    fp = round(Tf, (secs - s)*n)
    if fp ≥ n
        ip += Ti(1)
        fp -= Tf(n)
    end
    return (ip, fp)
end

function splittime(::Type{Ti}, ::Type{Tf}, secs::Floats,
                   n::Integer, ::SameSign) where {Ti<:Integer, Tf<:Integer}
    @assert isfinite(secs)
    if secs < 0
        s = -floor(-secs)
        ip = trunc(Ti, s)
        fp = round(Tf, (secs - s)*n)
        if fp ≤ -n
            ip -= Ti(1)
            fp += Tf(n)
        end
    else
        s = floor(secs)
        ip = trunc(Ti, s)
        fp = round(Tf, (secs - s)*n)
        if fp ≥ n
            ip += Ti(1)
            fp -= Tf(n)
        end
    end
    return (ip, fp)
end

"""

```julia
fixtime(Ti, Tf, i, f, n, r = Nonnegative()) -> (ip::Ti, fp::Tf)
```

yields two integers, `ip` and `fp` with `abs(fp) ∈ [0,n-1[`, such that
`ip + fp/n == i + f/n` (in arbitrary precision) and imposing the constraints
set by `r`:

* if `r` is `Nonnegative()`, then `fp ≥ 0`;

* if `r` is `SameSign()`, then `ip` and `fp` have the same sign.

The multiplier `n` must be strictly positive.

```julia
fixtime(TimeSpec, sec, nsec) -> ts
```

yields an instance of `TimeSpec` such that `ts` equals `sec` seconds plus
`nsec` nanoseconds and has normalized fields, that is with a nonnegative number
of nanoseconds strictly less than 1,000,000,000.

```julia
fixtime(TimeVal, sec, usec) -> tv
```

yields an instance of `TimeVal` such that `tv` equals `sec` seconds plus `usec`
microseconds and has normalized fields, that is with a nonnegative number of
microseconds strictly less than 1,000,000.

See also [`splittime`](@ref).

"""
function fixtime(::Type{Ti}, ::Type{Tf}, i::Integer, f::Integer,
                 n::Integer) where {Ti<:Integer, Tf<:Integer}
    fixtime(Ti, Tf, i, f, n, Nonnegative())
end

# Note that assuming n > 0, then q,r = divrem(f,n) yields q and r that have the
# same sign as f.
function fixtime(::Type{Ti}, ::Type{Tf}, i::Integer, f::Integer,
                 n::Integer, ::Nonnegative) where {Ti<:Integer, Tf<:Integer}
    ip = Ti(i + div(f, n))
    fp = Tf(rem(f, n))
    if fp < 0
        ip -= Ti(1)
        fp += Tf(n)
    end
    return (ip, fp)
end

function fixtime(::Type{Ti}, ::Type{Tf}, i::Integer, f::Integer,
                 n::Integer, ::SameSign) where {Ti<:Integer, Tf<:Integer}
    ip = Ti(i + div(f, n))
    fp = Tf(rem(f, n))
    if fp > 0
        if ip < 0
            ip += Ti(1)
            fp -= Tf(n)
        end
    elseif fp < 0
        if ip > 0
            ip -= Ti(1)
            fp += Tf(n)
        end
    end
    return (ip, fp)
end

fixtime(::Type{TimeSpec}, sec::Integer, nsec::Integer) =
    TimeSpec(fixtime(_typeof_timespec_sec, _typeof_timespec_nsec,
                     sec, nsec, NANOSECONDS_PER_SECOND, Nonnegative())...)

fixtime(::Type{TimeVal}, sec::Integer, usec::Integer) =
    TimeVal(fixtime(_typeof_timeval_sec, _typeof_timeval_usec,
                    sec, usec, MICROSECONDS_PER_SECOND, Nonnegative())...)

const TimeTypes = Union{TimeSpec,TimeVal}
const AnyTime = Union{TimeSpec,TimeVal,Libc.TimeVal}
const AnyTimeVal = Union{TimeVal,Libc.TimeVal}

#
# Extend unary minus, addition and subtraction for time structures.
#
Base.:(-)(a::T) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, -intpart(a), -fracpart(a))

Base.:(+)(a::T, b::T) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, intpart(a) + intpart(b), fracpart(a) + fracpart(b))
Base.:(-)(a::T, b::T) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, intpart(a) - intpart(b), fracpart(a) - fracpart(b))

Base.:(+)(a::TimeSpec, b::Union{TimeVal,Libc.TimeVal}) =
    fixtime(TimeSpec, intpart(a) + intpart(b), fracpart(a) + fracpart(b)*1_000)
Base.:(-)(a::TimeSpec, b::Union{TimeVal,Libc.TimeVal}) =
    fixtime(TimeSpec, intpart(a) - intpart(b), fracpart(a) - fracpart(b)*1_000)

Base.:(+)(a::Union{TimeVal,Libc.TimeVal}, b::TimeSpec) =
    fixtime(TimeSpec, intpart(a) + intpart(b), fracpart(a)*1_000 + fracpart(b))
Base.:(-)(a::Union{TimeVal,Libc.TimeVal}, b::TimeSpec) =
    fixtime(TimeSpec, intpart(a) - intpart(b), fracpart(a)*1_000 - fracpart(b))

Base.:(+)(a::TimeVal, b::Libc.TimeVal) =
    fixtime(TimeVal, intpart(a) + intpart(b), fracpart(a) + fracpart(b))
Base.:(-)(a::TimeVal, b::Libc.TimeVal) =
    fixtime(TimeVal, intpart(a) - intpart(b), fracpart(a) - fracpart(b))

Base.:(+)(a::Libc.TimeVal, b::TimeVal) =
    fixtime(TimeVal, intpart(a) + intpart(b), fracpart(a) + fracpart(b))
Base.:(-)(a::Libc.TimeVal, b::TimeVal) =
    fixtime(TimeVal, intpart(a) - intpart(b), fracpart(a) - fracpart(b))

Base.:(+)(a::T, b::Integer) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, intpart(a) + b, fracpart(a))
Base.:(+)(a::Integer, b::T) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, a + intpart(b), fracpart(b))

Base.:(-)(a::T, b::Integer) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, intpart(a) - b, fracpart(a))
Base.:(-)(a::Integer, b::T) where {T<:Union{TimeSpec,TimeVal}} =
    fixtime(T, a - intpart(b), -fracpart(b))

Base.:(+)(a::T, b::Real) where {T<:Union{TimeSpec,TimeVal}} = begin
    t = floor(b)
    ip, fp, n = intpart(a), fracpart(a), multiplier(a)
    fixtime(T, ip + trunc(typeof(ip), t), fp + round(typeof(fp), (b - t)*n))
end
Base.:(+)(a::Real, b::T) where {T<:Union{TimeSpec,TimeVal}} = b + a

Base.:(-)(a::T, b::Real) where {T<:Union{TimeSpec,TimeVal}} = begin
    t = floor(b)
    ip, fp, n = intpart(a), fracpart(a), multiplier(a)
    fixtime(T, ip - trunc(typeof(ip), t), fp - round(typeof(fp), (b - t)*n))
end
Base.:(-)(a::Real, b::T) where {T<:Union{TimeSpec,TimeVal}} = begin
    t = floor(a)
    ip, fp, n = intpart(b), fracpart(b), multiplier(b)
    fixtime(T, trunc(typeof(ip), t) - ip, round(typeof(fp), (a - t)*n) - fp)
end

intpart(t::TimeSpec) = t.sec
fracpart(t::TimeSpec) = t.nsec
multiplier(t::TimeSpec) = NANOSECONDS_PER_SECOND
tolerance(::Type{TimeSpec}) = 5e-9
scale(::Type{TimeSpec}) = 1e-9

intpart(t::Union{TimeVal,Libc.TimeVal}) = t.sec
fracpart(t::Union{TimeVal,Libc.TimeVal}) = t.usec
multiplier(t::Union{TimeVal,Libc.TimeVal}) = MICROSECONDS_PER_SECOND
tolerance(::Type{TimeVal}) = 1e-6
scale(::Type{TimeVal}) = 1e-6

#
# Extend isapprox and comparison operators for time structures.
#
Base.isapprox(a::T, b::T; kwds...) where {T<:Union{TimeSpec,TimeVal}} =
    _approx0(a - b; kwds...)
Base.isapprox(a::TimeSpec, b::Union{Real,TimeVal,Libc.TimeVal}; kwds...) =
    _approx0(a - b; kwds...)
Base.isapprox(a::Union{Real,TimeVal,Libc.TimeVal}, b::TimeSpec; kwds...) =
    _approx0(a - b; kwds...)
Base.isapprox(a::TimeVal, b::Union{Real,Libc.TimeVal}; kwds...) =
    _approx0(a - b; kwds...)
Base.isapprox(a::Union{Real,Libc.TimeVal}, b::TimeVal; kwds...) =
    _approx0(a - b; kwds...)

Base.:(==)(a::T, b::T) where {T<:Union{TimeSpec,TimeVal}}    = _eq(a, b)
Base.:(==)(a::TimeSpec, b::Union{Real,TimeVal,Libc.TimeVal}) = _eq(a, b)
Base.:(==)(a::Union{Real,TimeVal,Libc.TimeVal}, b::TimeSpec) = _eq(a, b)
Base.:(==)(a::TimeVal, b::Union{Real,Libc.TimeVal})          = _eq(a, b)
Base.:(==)(a::Union{Real,Libc.TimeVal}, b::TimeVal)          = _eq(a, b)

Base.:(≤)(a::T, b::T) where {T<:Union{TimeSpec,TimeVal}}    = _le(a, b)
Base.:(≤)(a::TimeSpec, b::Union{Real,TimeVal,Libc.TimeVal}) = _le(a, b)
Base.:(≤)(a::Union{Real,TimeVal,Libc.TimeVal}, b::TimeSpec) = _le(a, b)
Base.:(≤)(a::TimeVal, b::Union{Real,Libc.TimeVal})          = _le(a, b)
Base.:(≤)(a::Union{Real,Libc.TimeVal}, b::TimeVal)          = _le(a, b)

Base.:(<)(a::T, b::T) where {T<:Union{TimeSpec,TimeVal}}    = _lt(a, b)
Base.:(<)(a::TimeSpec, b::Union{Real,TimeVal,Libc.TimeVal}) = _lt(a, b)
Base.:(<)(a::Union{Real,TimeVal,Libc.TimeVal}, b::TimeSpec) = _lt(a, b)
Base.:(<)(a::TimeVal, b::Union{Real,Libc.TimeVal})          = _lt(a, b)
Base.:(<)(a::Union{Real,Libc.TimeVal}, b::TimeVal)          = _lt(a, b)

# Assuming normalized time t = (ip,fp) with fp ∈ [0,n-1] (normalized), the following
# table summarizes the possibilities.
#
#           ip < 0     ip = 0    ip > 0
#  fp < 0    t < 0      t < 0     t > 0   (not for normalized time)
#  fp = 0    t < 0      t = 0     t > 0
#  fp > 0    t < 0      t > 0     t > 0
#

_eq(a, b) = begin
    (ip, fp) = _split(a - b)
    ((ip == 0) & (fp == 0))
end

_le(a, b) = begin
    (ip, fp) = _split(a - b)
    ((ip < 0) | ((ip == 0) & (fp ≤ 0)))
end

_lt(a, b) = begin
    (ip, fp) = _split(a - b)
    (ip < 0)
end

function _approx0(a::T;
                  atol::Real = tolerance(T)) where {T<:Union{TimeSpec,TimeVal}}
    (ip, fp) = _split(a)
    abs(Float64(ip) + scale(T)*Float64(fp)) ≤ atol
end

_split(t::Union{TimeSpec,TimeVal,Libc.TimeVal}) = intpart(t), fracpart(t)

"""
```julia
error_message(args...)
```

yields a message string built from `args...`.  This string is meant to be used
as an error message.  The only difference between `string(args...)` and
`error_message(args...)` is that the latter is not inlined if the result has to
be dynamically created.  This is to avoid the caller function not being inlined
because the formatting of the error message involves too many operations.

See also  [`throw_argument_error`](@ref), [`throw_error_exception`](@ref).

"""
error_message(mesg::String) = mesg
@noinline error_message(args...) = string(args...)

"""
```julia
throw_argument_error(args...)
```

throws an `ArgumentError` exception with message built from `args...`.

See also [`error_message`](@ref), [`throw_error_exception`](@ref),
[`throw_system_error`](@ref).

"""
throw_argument_error(args...) = throw_argument_error(error_message(args...))

"""
```julia
throw_error_exception(args...)
```

throws an `ErrorException` with message built from `args...`.

See also [`error_message`](@ref), [`throw_argument_error`](@ref),
[`throw_system_error`](@ref).

"""
throw_error_exception(args...) = throw_error_exception(error_message(args...))

"""
```julia
throw_system_error(mesg, errno=Libc.errno())
```

throws an `SystemError` exception corresponding to error code number `errno`.
If `errno` is not specified the current value of the global `errno` variable in
the C library is used.

See also [`throw_argument_error`](@ref), [`throw_error_exception`](@ref).

"""
throw_system_error(mesg::AbstractString, errno::Integer = Libc.errno()) =
    throw(SystemError(mesg, errno))

makedims(dims::Tuple{Vararg{Int}}) = dims
makedims(dims::Tuple{Vararg{Integer}}) = map(Int, dims)
makedims(dims::AbstractVector{<:Integer}) = # FIXME: bad idea!
    ntuple(i -> Int(dims[i]), length(dims))

function checkdims(dims::NTuple{N,Integer}) where {N}
    number = one(Int)
    for i in 1:N
        dims[i] ≥ 1 || throw_error_exception("invalid dimension (dims[", i,
                                             "] = ", dims[i], ")")
        number *= convert(Int, dims[i])
    end
    return number
end

roundup(a::Integer, b::Integer) =
    roundup(convert(Int, a), convert(Int, b))

roundup(a::Int, b::Int) =
    div(a + (b - 1), b)*b

"""
```julia
get_memory_parameters(mem) -> ptr::Ptr{Cvoid}, siz::Int
```

yields the address and number of bytes of memory backed by `mem`.
The returned values arec checked for correctness.

See also: [`pointer`](@ref), [`sizeof`](@ref).

"""
function get_memory_parameters(mem)::Tuple{Ptr{Cvoid},Int}
    siz = sizeof(mem)
    isa(siz, Integer) || throw_argument_error("illegal type `", typeof(siz),
                                              "` for `sizeof(mem)`")
    siz ≥ 0 || throw_argument_error("invalid value `", siz,
                                    "` for `sizeof(mem)`")
    ptr = pointer(mem)
    isa(ptr, Ptr) || throw_argument_error("illegal type `", typeof(ptr),
                                          "` for `pointer(mem)`")
    return (convert(Ptr{Cvoid}, ptr), convert(Int, siz))
end

"""
```julia
_peek(T, buf, off) -> val
```

yields the value of type `T` stored at offset `off` (in bytes) in buffer
`buf` provided as a vector of bytes (`Uint8`).

```julia
_peek(T, ptr) -> val
```

yields the value of type `T` stored at address given by pointer `ptr`.

Also see: [`unsafe_load`](@ref), [`_poke!`](@ref).

"""
@inline _peek(::Type{T}, ptr::Ptr) where {T} =
    unsafe_load(convert(Ptr{T}, ptr))
@inline _peek(::Type{T}, buf::DenseVector{UInt8}, off::Integer) where {T} =
    _peek(T, pointer(buf) + off)

"""
```julia
_poke!(T, buf, off, val)
```

stores the value `val`, converted to type `T`, at offset `off` (in bytes) in
buffer `buf` provided as a vector of bytes (`Uint8`).

```julia
_poke!(T, ptr, val)
```

stores the value `val`, converted to type `T`, at address given by pointer
`ptr`.

Also see: [`unsafe_store!`](@ref), [`_peek`](@ref).

"""
@inline _poke!(::Type{T}, ptr::Ptr, val) where {T} =
    unsafe_store!(convert(Ptr{T}, ptr), val)
@inline _poke!(::Type{T}, buf::DenseVector{UInt8}, off::Integer, val) where {T} =
    _poke!(T, pointer(buf) + off, val)

#------------------------------------------------------------------------------
# DYNAMIC MEMORY OBJECTS

function _destroy(obj::DynamicMemory)
    if (ptr = obj.ptr) != C_NULL
        obj.len = 0
        obj.ptr = C_NULL
        Libc.free(ptr)
    end
end

Base.sizeof(obj::DynamicMemory) = obj.len
Base.pointer(obj::DynamicMemory) = obj.ptr
Base.unsafe_convert(::Type{Ptr{T}}, obj::DynamicMemory) where {T} =
    Ptr{T}(obj.ptr)


text:
#! format: off
should_precompile = true

# Don't edit the following! Instead change the script for `snoop_bot`.
ismultios = false
ismultiversion = true
# precompile_enclosure
@static if !should_precompile
    # nothing
elseif !ismultios && !ismultiversion
    @static if isfile(joinpath(@__DIR__, "../deps/SnoopCompile/precompile/precompile_Plots.jl"))
        include("../deps/SnoopCompile/precompile/precompile_Plots.jl")
        _precompile_()
    end
else
    @static if v"1.6.0-DEV" <= VERSION <= v"1.6.9"
        @static if isfile(joinpath(@__DIR__, "../deps/SnoopCompile/precompile//1.6/precompile_Plots.jl"))
            include("../deps/SnoopCompile/precompile//1.6/precompile_Plots.jl")
            _precompile_()
        end
    elseif v"1.7.0-DEV" <= VERSION <= v"1.7.9"
        @static if isfile(joinpath(@__DIR__, "../deps/SnoopCompile/precompile//1.7/precompile_Plots.jl"))
            include("../deps/SnoopCompile/precompile//1.7/precompile_Plots.jl")
            _precompile_()
        end
    else
    end

end # precompile_enclosure

ext: jl


text:
# lic checking file
# -----------------

# license check empty function
# ----------------------------
function emptyliccheck(lic::Vector{Cint})
    return lic
end

# lic check routine
# -----------------
function userlic( liccheck::Function = emptyliccheck )

	# change directory to reach all libs
	# ----------------------------------
    initdir = pwd()
    libdir = dirname(xprs)
    cd(libdir)

    # pre allocate vars
    # ----------------
    lic = Cint[1]
    slicmsg = Array(Cchar,512)
    errmsg = Array(Cchar,512)

    # FIRST call do xprslicense to get BASE LIC
    # -----------------------------------------
    ccall(("XPRSlicense",xprs),
        stdcall, Cint, (Ptr{Cint},Ptr{Cchar}), lic,slicmsg)

    # convert BASE LIC to GIVEN LIC
    # ---------------------------
    # change this line to the code given by  your FICO representative
    lic = liccheck(lic)

    # Send GIVEN LIC to XPRESS lib
    # --------------------------
    ierr = ccall(("XPRSlicense",xprs),
        stdcall, Cint, (Ptr{Cint},Ptr{Cchar}), lic, slicmsg)

    # check LIC TYPE
    # --------------
    if ierr == 16
        # DEVELOPER
        # ---------
        info("Xpress development software detected.")
    elseif ierr != 0
        # FAIL
        # ----
        info("Failed to find working license.")

        ccall(("XPRSgetlicerrmsg",xprs),
            stdcall, Cint, (Ptr{Cchar},Cint), errmsg, 512)

        error(  unsafe_string(pointer(errmsg))   )
    else
        # USER
        # ----
        info("User license detected.")
        info(  unsafe_string(pointer(slicmsg))  )
    end

    # go back to initial folder
    # ------------------------
    cd(initdir)
end


text:
using Test
using TestSetExtensions
using RBNF
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "qasm simple qasm2circuit" begin
    # Adding random comments and spacings in src1 to test robustness of lexer
    src = """
    OPENQASM 2.0;
    qreg a[3];

    gate syndrome(theta) a1,a2 //gate definition
    {
    rx(theta) a1; //rx gate with variable theta
    cx a1,a2;
    //controlled x gate
    }

    syndrome(0.2) a[0],a[2];
    """

    @testset  "qasm simple lexer" begin

        tokens = Qaintessent.lex(src)

        """Structure of RNBF token is:
            struct Token{T}
                lineno :: Int64
                colno  :: Int32
                offset :: Int64
                str    :: String
                span   :: Int32
        """

        tokens_ref = RBNF.Token[
            # OPENQASM starting lines
            RBNF.Token{:unnamed}(1, 1, 1, "OPENQASM", 8),
            RBNF.Token{:nnreal}(1, 10, 10, "2.0", 3),
            RBNF.Token{:unnamed}(1, 13, 13, ";", 1),
            # Setting up quantum registers
            RBNF.Token{:id}(2, 1, 15, "qreg", 4),
            RBNF.Token{:id}(2, 6, 20, "a", 1),
            RBNF.Token{:unnamed}(2, 7, 21, "[", 1),
            RBNF.Token{:nninteger}(2, 8, 22, "3", 1),
            RBNF.Token{:unnamed}(2, 9, 23, "]", 1),
            RBNF.Token{:unnamed}(2, 10, 24, ";", 1),
            # Seting up custom gate
            RBNF.Token{:id}(4, 1, 27, "gate", 4),
            RBNF.Token{:id}(4, 6, 32, "syndrome", 8),
            RBNF.Token{:unnamed}(4, 14, 40, "(", 1),
            RBNF.Token{:id}(4, 15, 41, "theta", 5),
            RBNF.Token{:unnamed}(4, 20, 46, ")", 1),
            RBNF.Token{:id}(4, 22, 48, "a1", 2),
            RBNF.Token{:unnamed}(4, 24, 50, ",", 1),
            RBNF.Token{:id}(4, 25, 51, "a2", 2),
            RBNF.Token{:unnamed}(5, 1, 72, "{", 1),
            RBNF.Token{:id}(6, 1, 74, "rx", 2),
            RBNF.Token{:unnamed}(6, 3, 76, "(", 1),
            RBNF.Token{:id}(6, 4, 77, "theta", 5),
            RBNF.Token{:unnamed}(6, 9, 82, ")", 1),
            RBNF.Token{:id}(6, 11, 84, "a1", 2),
            RBNF.Token{:unnamed}(6, 13, 86, ";", 1),
            RBNF.Token{:id}(7, 1, 118, "cx", 2),
            RBNF.Token{:id}(7, 4, 121, "a1", 2),
            RBNF.Token{:unnamed}(7, 6, 123, ",", 1),
            RBNF.Token{:id}(7, 7, 124, "a2", 2),
            RBNF.Token{:unnamed}(7, 9, 126, ";", 1),
            RBNF.Token{:unnamed}(9, 1, 148, "}", 1),
            # Using customg ate
            RBNF.Token{:id}(11, 1, 151, "syndrome", 8),
            RBNF.Token{:unnamed}(11, 9, 159, "(", 1),
            RBNF.Token{:nnreal}(11, 10, 160, "0.2", 3),
            RBNF.Token{:unnamed}(11, 13, 163, ")", 1),
            RBNF.Token{:id}(11, 15, 165, "a", 1),
            RBNF.Token{:unnamed}(11, 16, 166, "[", 1),
            RBNF.Token{:nninteger}(11, 17, 167, "0", 1),
            RBNF.Token{:unnamed}(11, 18, 168, "]", 1),
            RBNF.Token{:unnamed}(11, 19, 169, ",", 1),
            RBNF.Token{:id}(11, 20, 170, "a", 1),
            RBNF.Token{:unnamed}(11, 21, 171, "[", 1),
            RBNF.Token{:nninteger}(11, 22, 172, "2", 1),
            RBNF.Token{:unnamed}(11, 23, 173, "]", 1),
            RBNF.Token{:unnamed}(11, 24, 174, ";", 1)
        ]
        @test all(tokens_ref .== tokens)
    end

    @testset "qasm simple parse tokens" begin
        qasmtokens = Qaintessent.lex(src)
        rbnf_tokens = Qaintessent.parse_qasm(qasmtokens)
        rbnf_tokens_ref = Qaintessent.Struct_mainprogram(
                    RBNF.Token{:nnreal}(1, 10, 10, "2.0", 3),
                    Any[
                        Qaintessent.Struct_decl(RBNF.Token{:id}(2, 1, 15, "qreg", 4),
                            RBNF.Token{:id}(2, 6, 20, "a", 1), RBNF.Token{:nninteger}(2, 8, 22, "3", 1)),
                        Qaintessent.Struct_gate(
                            Qaintessent.Struct_gatedecl(RBNF.Token{:id}(4, 6, 32, "syndrome", 8),
                                Qaintessent.Struct_idlist(RBNF.Token{:id}(4, 15, 41, "theta", 5), nothing),
                                Qaintessent.Struct_idlist(RBNF.Token{:id}(4, 22, 48, "a1", 2),
                                Qaintessent.Struct_idlist(RBNF.Token{:id}(4, 25, 51, "a2", 2), nothing))),
                        Any[
                            Qaintessent.Struct_rx(RBNF.Token{:id}(6, 4, 77, "theta", 5),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(6, 11, 84, "a1", 2), nothing)),
                            Qaintessent.Struct_iduop(RBNF.Token{:id}(7, 1, 118, "cx", 2), nothing,
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(7, 4, 121, "a1", 2), nothing),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(7, 7, 124, "a2", 2), nothing), nothing)))]),
                            Qaintessent.Struct_iduop(RBNF.Token{:id}(11, 1, 151, "syndrome", 8),
                                Qaintessent.Struct_explist(RBNF.Token{:nnreal}(11, 10, 160, "0.2", 3), nothing),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(11, 15, 165, "a", 1), RBNF.Token{:nninteger}(11, 17, 167, "0", 1)),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(11, 20, 170, "a", 1), RBNF.Token{:nninteger}(11, 22, 172, "2", 1)), nothing)))]
                    )

        @test rbnf_tokens.ver == rbnf_tokens_ref.ver
        @test all(rbnf_tokens.prog .== rbnf_tokens_ref.prog)
    end

    @testset "qasm simple tokens to circuit" begin
        qasmtokens = Qaintessent.lex(src)
        rbnf_tokens = Qaintessent.parse_qasm(qasmtokens)
        circuit = Qaintessent.transform_qasm(rbnf_tokens)
        circuit_ref = Circuit{3}([
                circuit_gate(1, RxGate(0.2)),
                circuit_gate(3, X, 1)
                   ])
        @test all(circuit .≈ circuit_ref)
    end
end

# ##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "qasm complex qasm2circuit" begin
    # Adding random comments and spacings in src1 to test robustness of lexer
    src = """
    OPENQASM 2.0;
    qreg a[6];
    gate syn(alpha,beta) a1,a2,a3 //gate definition
    {
    rx(alpha) a1; //rx gate with variable theta
    cx a2,a3;
    t a3;
    sdg a1;
    crz(beta) a3,a1;
    ry(beta) a2;
    }

    gate drome(gamma,delta) a1,a2,a3 //gate definition
    {
    cz a2,a3;
    s a1;
    t a3;
    cy a2,a1;
    rz(gamma) a3;
    h a2;
    crx(delta) a3,a1;
    }

    syn(sin(0.2),cos(0.5)) a[0],a[1],a[2];
    drome(exp(3),3.4*pi) a[3],a[4],a[5];
    cry(ln(5)) a[1],a[5];
    crz(exp(4)) a[2],a[4];
    h a[3];
    """

    @testset "qasm complex lexer" begin
        tokens = Qaintessent.lex(src)

        """Structure of RNBF token is:
            struct Token{T}
                lineno :: Int64
                colno  :: Int32
                offset :: Int64
                str    :: String
                span   :: Int32
        """

        tokens_ref = RBNF.Token[
            RBNF.Token{:unnamed}(1, 1, 1, "OPENQASM", 8)
            RBNF.Token{:nnreal}(1, 10, 10, "2.0", 3)
            RBNF.Token{:unnamed}(1, 13, 13, ";", 1)
            RBNF.Token{:id}(2, 1, 15, "qreg", 4)
            RBNF.Token{:id}(2, 6, 20, "a", 1)
            RBNF.Token{:unnamed}(2, 7, 21, "[", 1)
            RBNF.Token{:nninteger}(2, 8, 22, "6", 1)
            RBNF.Token{:unnamed}(2, 9, 23, "]", 1)
            RBNF.Token{:unnamed}(2, 10, 24, ";", 1)
            RBNF.Token{:id}(3, 1, 26, "gate", 4)
            RBNF.Token{:id}(3, 6, 31, "syn", 3)
            RBNF.Token{:unnamed}(3, 9, 34, "(", 1)
            RBNF.Token{:id}(3, 10, 35, "alpha", 5)
            RBNF.Token{:unnamed}(3, 15, 40, ",", 1)
            RBNF.Token{:id}(3, 16, 41, "beta", 4)
            RBNF.Token{:unnamed}(3, 20, 45, ")", 1)
            RBNF.Token{:id}(3, 22, 47, "a1", 2)
            RBNF.Token{:unnamed}(3, 24, 49, ",", 1)
            RBNF.Token{:id}(3, 25, 50, "a2", 2)
            RBNF.Token{:unnamed}(3, 27, 52, ",", 1)
            RBNF.Token{:id}(3, 28, 53, "a3", 2)
            RBNF.Token{:unnamed}(4, 1, 74, "{", 1)
            RBNF.Token{:id}(5, 1, 76, "rx", 2)
            RBNF.Token{:unnamed}(5, 3, 78, "(", 1)
            RBNF.Token{:id}(5, 4, 79, "alpha", 5)
            RBNF.Token{:unnamed}(5, 9, 84, ")", 1)
            RBNF.Token{:id}(5, 11, 86, "a1", 2)
            RBNF.Token{:unnamed}(5, 13, 88, ";", 1)
            RBNF.Token{:id}(6, 1, 120, "cx", 2)
            RBNF.Token{:id}(6, 4, 123, "a2", 2)
            RBNF.Token{:unnamed}(6, 6, 125, ",", 1)
            RBNF.Token{:id}(6, 7, 126, "a3", 2)
            RBNF.Token{:unnamed}(6, 9, 128, ";", 1)
            RBNF.Token{:id}(7, 1, 130, "t", 1)
            RBNF.Token{:id}(7, 3, 132, "a3", 2)
            RBNF.Token{:unnamed}(7, 5, 134, ";", 1)
            RBNF.Token{:id}(8, 1, 136, "sdg", 3)
            RBNF.Token{:id}(8, 5, 140, "a1", 2)
            RBNF.Token{:unnamed}(8, 7, 142, ";", 1)
            RBNF.Token{:id}(9, 1, 144, "crz", 3)
            RBNF.Token{:unnamed}(9, 4, 147, "(", 1)
            RBNF.Token{:id}(9, 5, 148, "beta", 4)
            RBNF.Token{:unnamed}(9, 9, 152, ")", 1)
            RBNF.Token{:id}(9, 11, 154, "a3", 2)
            RBNF.Token{:unnamed}(9, 13, 156, ",", 1)
            RBNF.Token{:id}(9, 14, 157, "a1", 2)
            RBNF.Token{:unnamed}(9, 16, 159, ";", 1)
            RBNF.Token{:id}(10, 1, 161, "ry", 2)
            RBNF.Token{:unnamed}(10, 3, 163, "(", 1)
            RBNF.Token{:id}(10, 4, 164, "beta", 4)
            RBNF.Token{:unnamed}(10, 8, 168, ")", 1)
            RBNF.Token{:id}(10, 10, 170, "a2", 2)
            RBNF.Token{:unnamed}(10, 12, 172, ";", 1)
            RBNF.Token{:unnamed}(11, 1, 174, "}", 1)
            RBNF.Token{:id}(13, 1, 177, "gate", 4)
            RBNF.Token{:id}(13, 6, 182, "drome", 5)
            RBNF.Token{:unnamed}(13, 11, 187, "(", 1)
            RBNF.Token{:id}(13, 12, 188, "gamma", 5)
            RBNF.Token{:unnamed}(13, 17, 193, ",", 1)
            RBNF.Token{:id}(13, 18, 194, "delta", 5)
            RBNF.Token{:unnamed}(13, 23, 199, ")", 1)
            RBNF.Token{:id}(13, 25, 201, "a1", 2)
            RBNF.Token{:unnamed}(13, 27, 203, ",", 1)
            RBNF.Token{:id}(13, 28, 204, "a2", 2)
            RBNF.Token{:unnamed}(13, 30, 206, ",", 1)
            RBNF.Token{:id}(13, 31, 207, "a3", 2)
            RBNF.Token{:unnamed}(14, 1, 228, "{", 1)
            RBNF.Token{:id}(15, 1, 230, "cz", 2)
            RBNF.Token{:id}(15, 4, 233, "a2", 2)
            RBNF.Token{:unnamed}(15, 6, 235, ",", 1)
            RBNF.Token{:id}(15, 7, 236, "a3", 2)
            RBNF.Token{:unnamed}(15, 9, 238, ";", 1)
            RBNF.Token{:id}(16, 1, 240, "s", 1)
            RBNF.Token{:id}(16, 3, 242, "a1", 2)
            RBNF.Token{:unnamed}(16, 5, 244, ";", 1)
            RBNF.Token{:id}(17, 1, 246, "t", 1)
            RBNF.Token{:id}(17, 3, 248, "a3", 2)
            RBNF.Token{:unnamed}(17, 5, 250, ";", 1)
            RBNF.Token{:id}(18, 1, 252, "cy", 2)
            RBNF.Token{:id}(18, 4, 255, "a2", 2)
            RBNF.Token{:unnamed}(18, 6, 257, ",", 1)
            RBNF.Token{:id}(18, 7, 258, "a1", 2)
            RBNF.Token{:unnamed}(18, 9, 260, ";", 1)
            RBNF.Token{:id}(19, 1, 262, "rz", 2)
            RBNF.Token{:unnamed}(19, 3, 264, "(", 1)
            RBNF.Token{:id}(19, 4, 265, "gamma", 5)
            RBNF.Token{:unnamed}(19, 9, 270, ")", 1)
            RBNF.Token{:id}(19, 11, 272, "a3", 2)
            RBNF.Token{:unnamed}(19, 13, 274, ";", 1)
            RBNF.Token{:id}(20, 1, 276, "h", 1)
            RBNF.Token{:id}(20, 3, 278, "a2", 2)
            RBNF.Token{:unnamed}(20, 5, 280, ";", 1)
            RBNF.Token{:id}(21, 1, 282, "crx", 3)
            RBNF.Token{:unnamed}(21, 4, 285, "(", 1)
            RBNF.Token{:id}(21, 5, 286, "delta", 5)
            RBNF.Token{:unnamed}(21, 10, 291, ")", 1)
            RBNF.Token{:id}(21, 12, 293, "a3", 2)
            RBNF.Token{:unnamed}(21, 14, 295, ",", 1)
            RBNF.Token{:id}(21, 15, 296, "a1", 2)
            RBNF.Token{:unnamed}(21, 17, 298, ";", 1)
            RBNF.Token{:unnamed}(22, 1, 300, "}", 1)
            RBNF.Token{:id}(24, 1, 303, "syn", 3)
            RBNF.Token{:unnamed}(24, 4, 306, "(", 1)
            RBNF.Token{:id}(24, 5, 307, "sin", 3)
            RBNF.Token{:unnamed}(24, 8, 310, "(", 1)
            RBNF.Token{:nnreal}(24, 9, 311, "0.2", 3)
            RBNF.Token{:unnamed}(24, 12, 314, ")", 1)
            RBNF.Token{:unnamed}(24, 13, 315, ",", 1)
            RBNF.Token{:id}(24, 14, 316, "cos", 3)
            RBNF.Token{:unnamed}(24, 17, 319, "(", 1)
            RBNF.Token{:nnreal}(24, 18, 320, "0.5", 3)
            RBNF.Token{:unnamed}(24, 21, 323, ")", 1)
            RBNF.Token{:unnamed}(24, 22, 324, ")", 1)
            RBNF.Token{:id}(24, 24, 326, "a", 1)
            RBNF.Token{:unnamed}(24, 25, 327, "[", 1)
            RBNF.Token{:nninteger}(24, 26, 328, "0", 1)
            RBNF.Token{:unnamed}(24, 27, 329, "]", 1)
            RBNF.Token{:unnamed}(24, 28, 330, ",", 1)
            RBNF.Token{:id}(24, 29, 331, "a", 1)
            RBNF.Token{:unnamed}(24, 30, 332, "[", 1)
            RBNF.Token{:nninteger}(24, 31, 333, "1", 1)
            RBNF.Token{:unnamed}(24, 32, 334, "]", 1)
            RBNF.Token{:unnamed}(24, 33, 335, ",", 1)
            RBNF.Token{:id}(24, 34, 336, "a", 1)
            RBNF.Token{:unnamed}(24, 35, 337, "[", 1)
            RBNF.Token{:nninteger}(24, 36, 338, "2", 1)
            RBNF.Token{:unnamed}(24, 37, 339, "]", 1)
            RBNF.Token{:unnamed}(24, 38, 340, ";", 1)
            RBNF.Token{:id}(25, 1, 342, "drome", 5)
            RBNF.Token{:unnamed}(25, 6, 347, "(", 1)
            RBNF.Token{:id}(25, 7, 348, "exp", 3)
            RBNF.Token{:unnamed}(25, 10, 351, "(", 1)
            RBNF.Token{:nninteger}(25, 11, 352, "3", 1)
            RBNF.Token{:unnamed}(25, 12, 353, ")", 1)
            RBNF.Token{:unnamed}(25, 13, 354, ",", 1)
            RBNF.Token{:nnreal}(25, 14, 355, "3.4", 3)
            RBNF.Token{:unnamed}(25, 17, 358, "*", 1)
            RBNF.Token{:id}(25, 18, 359, "pi", 2)
            RBNF.Token{:unnamed}(25, 20, 361, ")", 1)
            RBNF.Token{:id}(25, 22, 363, "a", 1)
            RBNF.Token{:unnamed}(25, 23, 364, "[", 1)
            RBNF.Token{:nninteger}(25, 24, 365, "3", 1)
            RBNF.Token{:unnamed}(25, 25, 366, "]", 1)
            RBNF.Token{:unnamed}(25, 26, 367, ",", 1)
            RBNF.Token{:id}(25, 27, 368, "a", 1)
            RBNF.Token{:unnamed}(25, 28, 369, "[", 1)
            RBNF.Token{:nninteger}(25, 29, 370, "4", 1)
            RBNF.Token{:unnamed}(25, 30, 371, "]", 1)
            RBNF.Token{:unnamed}(25, 31, 372, ",", 1)
            RBNF.Token{:id}(25, 32, 373, "a", 1)
            RBNF.Token{:unnamed}(25, 33, 374, "[", 1)
            RBNF.Token{:nninteger}(25, 34, 375, "5", 1)
            RBNF.Token{:unnamed}(25, 35, 376, "]", 1)
            RBNF.Token{:unnamed}(25, 36, 377, ";", 1)
            RBNF.Token{:id}(26, 1, 379, "cry", 3)
            RBNF.Token{:unnamed}(26, 4, 382, "(", 1)
            RBNF.Token{:id}(26, 5, 383, "ln", 2)
            RBNF.Token{:unnamed}(26, 7, 385, "(", 1)
            RBNF.Token{:nninteger}(26, 8, 386, "5", 1)
            RBNF.Token{:unnamed}(26, 9, 387, ")", 1)
            RBNF.Token{:unnamed}(26, 10, 388, ")", 1)
            RBNF.Token{:id}(26, 12, 390, "a", 1)
            RBNF.Token{:unnamed}(26, 13, 391, "[", 1)
            RBNF.Token{:nninteger}(26, 14, 392, "1", 1)
            RBNF.Token{:unnamed}(26, 15, 393, "]", 1)
            RBNF.Token{:unnamed}(26, 16, 394, ",", 1)
            RBNF.Token{:id}(26, 17, 395, "a", 1)
            RBNF.Token{:unnamed}(26, 18, 396, "[", 1)
            RBNF.Token{:nninteger}(26, 19, 397, "5", 1)
            RBNF.Token{:unnamed}(26, 20, 398, "]", 1)
            RBNF.Token{:unnamed}(26, 21, 399, ";", 1)
            RBNF.Token{:id}(27, 1, 401, "crz", 3)
            RBNF.Token{:unnamed}(27, 4, 404, "(", 1)
            RBNF.Token{:id}(27, 5, 405, "exp", 3)
            RBNF.Token{:unnamed}(27, 8, 408, "(", 1)
            RBNF.Token{:nninteger}(27, 9, 409, "4", 1)
            RBNF.Token{:unnamed}(27, 10, 410, ")", 1)
            RBNF.Token{:unnamed}(27, 11, 411, ")", 1)
            RBNF.Token{:id}(27, 13, 413, "a", 1)
            RBNF.Token{:unnamed}(27, 14, 414, "[", 1)
            RBNF.Token{:nninteger}(27, 15, 415, "2", 1)
            RBNF.Token{:unnamed}(27, 16, 416, "]", 1)
            RBNF.Token{:unnamed}(27, 17, 417, ",", 1)
            RBNF.Token{:id}(27, 18, 418, "a", 1)
            RBNF.Token{:unnamed}(27, 19, 419, "[", 1)
            RBNF.Token{:nninteger}(27, 20, 420, "4", 1)
            RBNF.Token{:unnamed}(27, 21, 421, "]", 1)
            RBNF.Token{:unnamed}(27, 22, 422, ";", 1)
            RBNF.Token{:id}(28, 1, 424, "h", 1)
            RBNF.Token{:id}(28, 3, 426, "a", 1)
            RBNF.Token{:unnamed}(28, 4, 427, "[", 1)
            RBNF.Token{:nninteger}(28, 5, 428, "3", 1)
            RBNF.Token{:unnamed}(28, 6, 429, "]", 1)
            RBNF.Token{:unnamed}(28, 7, 430, ";", 1)
            ]
        @test all(tokens_ref .== tokens)
    end

    @testset "qasm complex parse tokens" begin
        qasmtokens = Qaintessent.lex(src)
        rbnf_tokens = Qaintessent.parse_qasm(qasmtokens)
        rbnf_tokens_ref =
            Qaintessent.Struct_mainprogram(
                RBNF.Token{:nnreal}(1, 10, 10, "2.0", 3),
                Any[
                    Qaintessent.Struct_decl(
                        RBNF.Token{:id}(2, 1, 15, "qreg", 4),
                        RBNF.Token{:id}(2, 6, 20, "a", 1),
                        RBNF.Token{:nninteger}(2, 8, 22, "6", 1)),
                    Qaintessent.Struct_gate(
                        Qaintessent.Struct_gatedecl(RBNF.Token{:id}(3, 6, 31, "syn", 3),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(3, 10, 35, "alpha", 5),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(3, 16, 41, "beta", 4), nothing)),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(3, 22, 47, "a1", 2),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(3, 25, 50, "a2", 2),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(3, 28, 53, "a3", 2), nothing)))),
                        Any[
                            Qaintessent.Struct_rx(RBNF.Token{:id}(5, 4, 79, "alpha", 5),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(5, 11, 86, "a1", 2), nothing)),
                            Qaintessent.Struct_iduop(RBNF.Token{:id}(6, 1, 120, "cx", 2), nothing,
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(6, 4, 123, "a2", 2), nothing),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(6, 7, 126, "a3", 2), nothing), nothing))),
                            Qaintessent.Struct_t(Qaintessent.Struct_argument(RBNF.Token{:id}(7, 3, 132, "a3", 2), nothing)),
                            Qaintessent.Struct_sdg(Qaintessent.Struct_argument(RBNF.Token{:id}(8, 5, 140, "a1", 2), nothing)),
                            Qaintessent.Struct_crz(RBNF.Token{:id}(9, 5, 148, "beta", 4),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(9, 11, 154, "a3", 2), nothing),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(9, 14, 157, "a1", 2), nothing)),
                            Qaintessent.Struct_ry(RBNF.Token{:id}(10, 4, 164, "beta", 4),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(10, 10, 170, "a2", 2), nothing))]),
                    Qaintessent.Struct_gate(
                        Qaintessent.Struct_gatedecl(RBNF.Token{:id}(13, 6, 182, "drome", 5),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(13, 12, 188, "gamma", 5),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(13, 18, 194, "delta", 5), nothing)),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(13, 25, 201, "a1", 2),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(13, 28, 204, "a2", 2),
                            Qaintessent.Struct_idlist(RBNF.Token{:id}(13, 31, 207, "a3", 2), nothing)))),
                        Any[
                            Qaintessent.Struct_iduop(RBNF.Token{:id}(15, 1, 230, "cz", 2), nothing,
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(15, 4, 233, "a2", 2), nothing),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(15, 7, 236, "a3", 2), nothing), nothing))),
                            Qaintessent.Struct_s(Qaintessent.Struct_argument(RBNF.Token{:id}(16, 3, 242, "a1", 2), nothing)),
                            Qaintessent.Struct_t(Qaintessent.Struct_argument(RBNF.Token{:id}(17, 3, 248, "a3", 2), nothing)),
                            Qaintessent.Struct_iduop(RBNF.Token{:id}(18, 1, 252, "cy", 2), nothing,
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(18, 4, 255, "a2", 2), nothing),
                                Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(18, 7, 258, "a1", 2), nothing), nothing))),
                            Qaintessent.Struct_rz(RBNF.Token{:id}(19, 4, 265, "gamma", 5),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(19, 11, 272, "a3", 2), nothing)),
                            Qaintessent.Struct_h(Qaintessent.Struct_argument(RBNF.Token{:id}(20, 3, 278, "a2", 2), nothing)),
                            Qaintessent.Struct_crx(RBNF.Token{:id}(21, 5, 286, "delta", 5),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(21, 12, 293, "a3", 2), nothing),
                                Qaintessent.Struct_argument(RBNF.Token{:id}(21, 15, 296, "a1", 2), nothing))]),
                    Qaintessent.Struct_iduop(RBNF.Token{:id}(24, 1, 303, "syn", 3),
                        Qaintessent.Struct_explist(Qaintessent.Struct_fnexp(RBNF.Token{:id}(24, 5, 307, "sin", 3), RBNF.Token{:nnreal}(24, 9, 311, "0.2", 3)),
                        Qaintessent.Struct_explist(Qaintessent.Struct_fnexp(RBNF.Token{:id}(24, 14, 316, "cos", 3), RBNF.Token{:nnreal}(24, 18, 320, "0.5", 3)), nothing)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(24, 24, 326, "a", 1), RBNF.Token{:nninteger}(24, 26, 328, "0", 1)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(24, 29, 331, "a", 1), RBNF.Token{:nninteger}(24, 31, 333, "1", 1)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(24, 34, 336, "a", 1), RBNF.Token{:nninteger}(24, 36, 338, "2", 1)), nothing)))),
                    Qaintessent.Struct_iduop(RBNF.Token{:id}(25, 1, 342, "drome", 5),
                        Qaintessent.Struct_explist(Qaintessent.Struct_fnexp(RBNF.Token{:id}(25, 7, 348, "exp", 3), RBNF.Token{:nninteger}(25, 11, 352, "3", 1)),
                        Qaintessent.Struct_explist(Qaintessent.Struct_bin(RBNF.Token{:nnreal}(25, 14, 355, "3.4", 3), RBNF.Token{:unnamed}(25, 17, 358, "*", 1),
                        Qaintessent.Struct_pi()), nothing)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(25, 22, 363, "a", 1), RBNF.Token{:nninteger}(25, 24, 365, "3", 1)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(25, 27, 368, "a", 1), RBNF.Token{:nninteger}(25, 29, 370, "4", 1)),
                        Qaintessent.Struct_mixedlist(Qaintessent.Struct_argument(RBNF.Token{:id}(25, 32, 373, "a", 1), RBNF.Token{:nninteger}(25, 34, 375, "5", 1)), nothing)))),
                    Qaintessent.Struct_cry(Qaintessent.Struct_fnexp(RBNF.Token{:id}(26, 5, 383, "ln", 2), RBNF.Token{:nninteger}(26, 8, 386, "5", 1)),
                        Qaintessent.Struct_argument(RBNF.Token{:id}(26, 12, 390, "a", 1), RBNF.Token{:nninteger}(26, 14, 392, "1", 1)),
                        Qaintessent.Struct_argument(RBNF.Token{:id}(26, 17, 395, "a", 1), RBNF.Token{:nninteger}(26, 19, 397, "5", 1))),
                    Qaintessent.Struct_crz(Qaintessent.Struct_fnexp(RBNF.Token{:id}(27, 5, 405, "exp", 3), RBNF.Token{:nninteger}(27, 9, 409, "4", 1)),
                        Qaintessent.Struct_argument(RBNF.Token{:id}(27, 13, 413, "a", 1), RBNF.Token{:nninteger}(27, 15, 415, "2", 1)),
                        Qaintessent.Struct_argument(RBNF.Token{:id}(27, 18, 418, "a", 1), RBNF.Token{:nninteger}(27, 20, 420, "4", 1))),
                    Qaintessent.Struct_h(Qaintessent.Struct_argument(RBNF.Token{:id}(28, 3, 426, "a", 1), RBNF.Token{:nninteger}(28, 5, 428, "3", 1)))])

        @test rbnf_tokens.ver == rbnf_tokens_ref.ver
        @test all(rbnf_tokens.prog .== rbnf_tokens_ref.prog)
    end

    @testset "qasm complex tokens to circuit" begin
        qasmtokens = Qaintessent.lex(src)
        rbnf_tokens = Qaintessent.parse_qasm(qasmtokens)
        circuit = Qaintessent.transform_qasm(rbnf_tokens)

        α = sin(0.2)
        β = cos(0.5)
        γ = exp(3)
        δ = 3.4π
        ϵ = log(5)
        ζ = exp(4)

        circuit_ref = Circuit{6}(
            [circuit_gate(1, RxGate(α)),
            circuit_gate(3, X, 2),
            circuit_gate(3, TGate()),
            circuit_gate(1, SdagGate()),
            circuit_gate(1, RzGate(β), 3),
            circuit_gate(2, RyGate(β)),
            circuit_gate(6, Z, 5),
            circuit_gate(4, SGate()),
            circuit_gate(6, TGate()),
            circuit_gate(4, Y, 5),
            circuit_gate(6, RzGate(γ)),
            circuit_gate(5, HadamardGate()),
            circuit_gate(4, RxGate(δ), 6),
            circuit_gate(6, RyGate(ϵ), 2),
            circuit_gate(5, RzGate(ζ),3),
            circuit_gate(4, HadamardGate())]
    )
        @test all(circuit .≈ circuit_ref)
    end
end

# ##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "test reading openqasm" begin
    src1 = """
    // Repetition code syndrome measurement
    OPENQASM 2.0;
    include "qelib1.inc";
    qreg d[3];
    qreg a[2];
    qreg c[3];
    creg e[2];
    gate syndrome(alpha) d1,d2,d3,a1,a2
    {
      cx d1,a1; cx d2,a1;
      cx d2,a2; cx d3,a2; ry(alpha) a1;
    }
    gate syn2drome(theta) d1,d2,d3,a1,a2
    {
        syndrome(theta) d1,d2,d3,a1,a2;
        syndrome(theta) d2,d3,a1,d1,a2;
        rx(theta) d2;
        cx d3,a1;
        cx a1,d3;
        cx d3,a1;
    }
    syn2drome(0.1*pi) d[0],d[1],d[2],a[0],a[1];
    t d[1];
    s a[0];
    h d[0];
    x a[1];
    y d[0];
    z c[0];

    measure d[1] -> e[0];
    measure a[0] -> e[1];
    """

    d1 = qreg(3)
    a1 = qreg(2)
    c1 = qreg(3)
    cgc_ref = Circuit(d1, a1, c1)

    N = num_wires(cgc_ref)
    gates_ref = [
        circuit_gate(4, X, 1),
        circuit_gate(4, X, 2),
        circuit_gate(5, X, 2),
        circuit_gate(5, X, 3),
        circuit_gate(4, RyGate(0.1π)),
        circuit_gate(1, X, 2),
        circuit_gate(1, X, 3),
        circuit_gate(5, X, 3),
        circuit_gate(5, X, 4),
        circuit_gate(1, RyGate(0.1π)),
        circuit_gate(2, RxGate(0.1π)),
        circuit_gate(4, 3, SwapGate()),
        circuit_gate(2, TGate()),
        circuit_gate(4, SGate()),
        circuit_gate(1, HadamardGate()),
        circuit_gate(5, X),
        circuit_gate(1, Y),
        circuit_gate(6, Z)
    ]
    meas_ref = [
        mop(Z, 2),
        mop(Z, 4)
    ]

    append!(cgc_ref, gates_ref)
    add_measurement!(cgc_ref, meas_ref)

    cgc = qasm2cgc(src1)
    ψ = randn(ComplexF64, 2^N)

    @test apply(ψ, cgc_ref.moments) ≈ apply(ψ, cgc.moments)
    @test apply(ψ, cgc_ref) ≈ apply(ψ, cgc)
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "test openqasm warnings" begin
    src1 = """
    // Repetition code syndrome measurement
    OPENQASM 2.0;
    include "qelib1.inc";
    qreg d[3];
    qreg a[2];
    qreg c[3];
    creg e[2]
    gate syndrome(alpha) d1,d2,d3,a1,a2
    {
      cx d1,a1; cx d2,a1;
      cx d2,a2; cx d3,a2; ry(alpha) a1;
    }
    gate syn2drome(theta) d1,d2,d3,a1,a2
    {
        syndrome(theta) d1,d2,d3,a1,a2;
        syndrome(theta) d2,d3,a1,d1,a2;
        rx(theta) d2;
        cx d3,a1;
        cx a1,d3;
        cx d3,a1;
    }
    syn2drome(0.1*pi) d[0],d[1],d[2],a[0],a[1];
    t d[1];
    barrier d[1];
    s a[0];
    h d[0];

    measure d[1] -> e[0];
    measure a[0] -> e[1];
    """

    @test_logs (:warn,"Barrier operation is not yet supported") qasm2cgc(src1)

    src2 = """
    // Repetition code syndrome measurement
    OPENQASM 2.0;
    include "qelib1.inc";
    qreg d[3];
    qreg a[2];
    qreg c[3];
    creg e[2]
    gate syndrome(alpha) d1,d2,d3,a1,a2
    {
      cx d1,a1; cx d2,a1;
      cx d2,a2; cx d3,a2; ry(alpha) a1;
    }
    gate syn2drome(theta) d1,d2,d3,a1,a2
    {
        syndrome(theta) d1,d2,d3,a1,a2;
        syndrome(theta) d2,d3,a1,d1,a2;
        rx(theta) d2;
        cx d3,a1;
        cx a1,d3;
        cx d3,a1;
    }
    syn2drome(0.1*pi) d[0],d[1],d[2],a[0],a[1];
    t d[1];
    reset d[1];
    s a[0];
    h d[0];
    cx d[0], d[1];

    measure d[1] -> e[0];
    measure a[0] -> e[1];
    """

    @test_logs (:warn,"Reset operation is not yet supported") qasm2cgc(src2)
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "test writing" begin
    src_ref = """
    OPENQASM 2.0;

    include "qelib1.inc";

    qreg qregister[8];
    cx qregister[0],qregister[3];
    cx qregister[1],qregister[3];
    cx qregister[1],qregister[4];
    cx qregister[2],qregister[4];
    ry(0.3141592653589793) qregister[3];
    cx qregister[1],qregister[0];
    cx qregister[2],qregister[0];
    cx qregister[2],qregister[4];
    cx qregister[3],qregister[4];
    ry(0.3141592653589793) qregister[0];
    rx(0.3141592653589793) qregister[1];
    cx qregister[2],qregister[1];
    t qregister[1];
    s qregister[3];
    h qregister[0];"""

    d1 = qreg(3)
    a1 = qreg(2)
    c1 = qreg(3)
    cgc_ref = Circuit(d1, a1, c1)

    N = num_wires(cgc_ref)

    gates_ref = [
        circuit_gate(4, X, 1),
        circuit_gate(4, X, 2),
        circuit_gate(5, X, 2),
        circuit_gate(5, X, 3),
        circuit_gate(4, RyGate(0.1π)),
        circuit_gate(1, X, 2),
        circuit_gate(1, X, 3),
        circuit_gate(5, X, 3),
        circuit_gate(5, X, 4),
        circuit_gate(1, RyGate(0.1π)),
        circuit_gate(2, RxGate(0.1π)),
        circuit_gate(2, X, 3),
        circuit_gate(2, TGate()),
        circuit_gate(4, SGate()),
        circuit_gate(1, HadamardGate()),
    ]
    append!(cgc_ref, gates_ref)

    ψ = randn(ComplexF64, 2^N)
    cgc = qasm2cgc(src_ref)

    @test apply(ψ, cgc.moments) ≈ apply(ψ, cgc_ref.moments)
    @test src_ref == cgc2qasm(cgc_ref)
end


text:
#this file simulate omics data for nonlinear system
# the QTLs ("kk"), QTL effects ("u.txt"), QTL position ("qtl.txt"), omics effects ("alpha.txt") in Christensenet al.(2021) are required
# above data can be found in Christensenet al.(2021) or below link:
# http://genoweb.toulouse.inra.fr/~alegarra/GOBLUP/

using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets,DelimitedFiles,Random, LinearAlgebra, StatsBase

cd("/Users/tianjing/Box/Omicsdata")

Random.seed!(123)
### step1. get genetic value of omics (since it needs to be transformed to get true bv)
#read QTL
nQTL=20_000 #there are 20000 QLTs in Christensenet al.(2021), but only 5000 will be used
inter=4     #select 1 every 4 QTLs (5000/20000 QTLs)
nOmics=1200
nInd=21100
σ2_mi=2
h2_m=0.61
σ2_gi=h2_m*σ2_mi      #genetic variance of each omics
σ2_ei=(1-h2_m)*σ2_mi  #residual variance of each omics

QTL=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/andrew_code/kk") #QTLs in Christensenet al.(2021)
QTL=QTL[:,2:end] #remove ID
pos=collect(inter:inter:nQTL)
QTL=QTL[:,pos]
writedlm("nonlinear/QTL5000.txt",QTL) #evenly select 5000 QTLs as in Christensenet al.(2021)

#read 500 selected QTL position for all 1200 omics features
QTL_select_pos=Int.(readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/andrew_code/qtl.txt"))
maximum(QTL_select_pos)
minimum(QTL_select_pos)
#read 500 selected QTL effects for all 1200 omics features
QTL_select_effect=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/andrew_code/u.txt")
#calculate omics
# 1st layer
G=zeros(nInd,nOmics) #genetic value of omics
E=zeros(nInd,nOmics) #residuals of omics
Random.seed!(1)
for i in 1:nOmics
	W=QTL[:,QTL_select_pos[:,i]]
	@show size(W),W[1:2,1:2]
    Gi=W*QTL_select_effect[:,i]
	Gi=Gi/std(Gi)*sqrt(σ2_gi)
	G[:,i]=Gi
	E[:,i]=randn(nInd)*sqrt(σ2_ei)
end

var(G,dims=1)
var(E,dims=1)

writedlm("nonlinear/varmi_2/omics_bv.txt",G)
writedlm("nonlinear/varmi_2/omics_residual.txt",E)

omics=G+E
writedlm("nonlinear/varmi_2/omics.txt",omics)
var(omics,dims=1)

mysigmoid(x) = 1/(1+exp(-x))
#save g(omics_g)
G_nonlinear=mysigmoid.(G)
writedlm("nonlinear/varmi_2/omics_bv_nonlinear.txt",G_nonlinear)
#save omics nonlinear
omics_nonlinear=mysigmoid.(omics)
writedlm("nonlinear/varmi_2/omics_nonlinear.txt",omics_nonlinear)

#read neural network weight between omics and phenotype (omics effects on phenotype)
w1=vec(readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/andrew_code/alpha.txt"))

#save true bv
tbv=G_nonlinear*w1
writedlm("nonlinear/varmi_2/tbv.txt",tbv)

#create y_nonlinear
omics_contribution=omics_nonlinear*w1
σ2_z=var(omics_contribution)
σ2_g_nonlinear=var(G_nonlinear*w1)
#h2=σ2_g_nonlinear/(σ2_z+σ2_e)
h2=0.337
σ2_e=σ2_g_nonlinear/h2 - σ2_z #residual variance of y

Random.seed!(1)
y_residual=randn(nInd)*sqrt(σ2_e)
var(y_residual)
y_nonlinear=omics_contribution+y_residual
var(y_nonlinear)
σ2_g_nonlinear/var(y_nonlinear)
writedlm("nonlinear/varmi_2/y.txt",y_nonlinear)

#move genotype from old to new folders
for i in 1:20
	mkdir("nonlinear/varmi_2/data$i");
end
#create 1055 dataset
y_nonlinear=vec(readdlm("nonlinear/varmi_2/y.txt"))
omics=readdlm("nonlinear/varmi_2/omics.txt")
tbv=vec(readdlm("nonlinear/varmi_2/tbv.txt"))
for i in 1:20
	id=Int.(vec(readdlm("andrew1055/data$i/select_1055_index_data$i.txt")))
	# save y_nonlinear
	yi=y_nonlinear[id]
	writedlm("nonlinear/varmi_2/data$i/y.txt",yi)
	yi_df=DataFrame(ID=id,y=yi)
	CSV.write("nonlinear/varmi_2/data$i/y.csv", yi_df)
	# save omics (should use omics, not g(omics) in the model)
	omicsi=omics[id,:]
	writedlm("nonlinear/varmi_2/data$i/omics.txt",omicsi)
	names=["gene$i" for i in 1:1200]
    omicsi_df=DataFrame(omicsi,names)
	insertcols!(omicsi_df, 1, :ID => id)
	CSV.write("nonlinear/varmi_2/data$i/omics.csv", omicsi_df)
	#save true bv
	bvi=tbv[id]
	writedlm("nonlinear/varmi_2/data$i/tbv.txt",bvi)
	bvi_df=DataFrame(ID=id,tbv=bvi)
	CSV.write("nonlinear/varmi_2/data$i/tbv.csv", bvi_df)
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
module UnFortran

type FortIOStream
    s::IOStream
end

FortIOStream(args...) = FortIOStream(IOStream(args...))

fort_open(args...) = FortIOStream(open(args...))

function fort_open(f::Function, args...)
    open(args...) do file
        f(FortIOStream(file))
    end
end

import Base.close
close(file::FortIOStream) = close(f.s)

function read_record(f::Function, file::FortIOStream, expect_bytes)
    nbytes = read(file.s, Int32)
    nbytes == expect_bytes || error("Record size does not agree with types requested")
    r = f(file)
    nbytes == read(file.s, Int32) || error("Record size check failed")
    r
end

function write_record(f::Function, file::FortIOStream, nbytes::Integer)
    write(file.s, convert(Int32, nbytes))
    f(file)
    write(file.s, convert(Int32, nbytes))
end

record_size(args...) = reduce(+, 0, [sizeof(a) for a in args])

import Base.read, Base.read!
function read(file::FortIOStream, args...)
    read_record(file, record_size(args...)) do file
        [read(file.s, a) for a in args]
    end
end

function read!{T,N}(file::FortIOStream, a::AbstractArray{T,N})
    ni = length(a)
    read_record(file, ni*sizeof(T)) do file
        for i = 1:ni
            a[i] = read(file.s, T)
        end
    end
    a
end

import Base.write
function write(file::FortIOStream, args...)
    write_record(file, record_size(args...)) do file
        for a in args
            write(file.s, a)
        end
    end
end

function write{T,N}(file::FortIOStream, a::AbstractArray{T,N})
    ni = length(a)
    write_record(file, ni*sizeof(T)) do file
        for i = 1:ni
            write(file.s, a[i])
        end
    end
end

export fort_open, close, read, read!, write

end # module


text:
# PS7 Dayuan Xie

# using Pkg
# Pkg.add("SMM")

using SMM
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using DataFramesMeta
using CSV

# Question 1
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function gmm(beta, X, y)
    g = y .-X*beta
    J = g'*I*g
    return J
end
optim_gmm = optimize(beta ->gmm(beta, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-8, iterations=100_000))
println(optim_gmm.minimizer)

beta_ols = inv(X'*X)*X'*y
println(beta_ols)

# Question 2a
freqtable(df, :occupation)
df = dropmissing(df, :occupation)
df[df.occupation.==8 ,:occupation] .= 7
df[df.occupation.==9 ,:occupation] .= 7
df[df.occupation.==10,:occupation] .= 7
df[df.occupation.==11,:occupation] .= 7
df[df.occupation.==12,:occupation] .= 7
df[df.occupation.==13,:occupation] .= 7
freqtable(df, :occupation)

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function logit_ml(alpha, X, y)
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    num = zeros(N,J)
    dem = zeros(N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j])
            dem .+= num[:,j]
        end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end

alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
alpha_start = alpha_true.*rand(size(alpha_true))
println(size(alpha_true))
alpha_hat_optim = optimize(a -> logit_ml(a, X, y), alpha_start, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
optim_ml = alpha_hat_optim.minimizer
println(optim_ml)

# Question 2b
function logit_gmm(alpha, X, y)
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    num = zeros(N,J)
    dem = zeros(N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j])
            dem .+= num[:,j]
        end

    P = num./repeat(dem,1,J)
    g = bigY .- P
    gmm = sum(g'*I*g)
    return gmm
end
optim_logit_gmm_b = optimize(a -> logit_gmm(a, X, y), optim_ml, LBFGS(), Optim.Options(g_tol=1e-8, iterations=100_000))
println(optim_logit_gmm_b.minimizer)

# Question 2c
optim_logit_gmm_c = optimize(a -> logit_gmm(a, X, y), alpha_start, LBFGS(), Optim.Options(g_tol=1e-8, iterations=100_000))
println(optim_logit_gmm_c.minimizer)

# We get different estimates from part (b) and (c). The objective function is not globally concave.

# Question 3
function simulate()
K = 2
N = 10000
J = 3
X=rand(N,J)
β=[1 2;3 4]
num = zeros(N,J)
dem = zeros(N)
    for j=1:J
        num[:,j] = exp.(X*β[:,j])
        dem .+= num[:,j]
    end

P = num./repeat(dem,1,J)
ε = randn(N,1)
Y=zeros(N,1)
   for i= 1: N
    Y[i,1]=sum(sum(P[i,:])>ε[i,1]) + sum((P[i,2]+P[i,3])>ε[i,1]) + sum(P[i,3]>ε[i,1])
   end
end
simulate()

# Question 4
MA = SMM.parallelNormal()
dc = SMM.history(MA.chains[1])
dc = dc[dc[:accepted].==true, :]
println(describe(dc))

# Question 5
function mlogit_smm(θ, X, y, D)
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    β = θ[1:end-1]
    σ = θ[end]
    if length(β)==1
        β = β[1]
    end
    num = zeros(N,J)
    dem = zeros(N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j])
            dem .+= num[:,j]
        end

    P = num./repeat(dem,1,J)

    gmodel = zeros(N+1,J,D)
    gdata  = vcat(bigY,var(bigY))

    Random.seed!(1234)
    for d=1:D
        ε = σ*randn(N,J)
        ỹ = P .+ ε
        gmodel[1:end-1,1:J,d] = ỹ
        gmodel[  end  ,1:J,d] = var(ỹ)
    end

    err = vec(gdata .- mean(gmodel; dims=2))

    J = err'*I*err
    return J
end

max_stars_repo_path: ProblemSets/PS7-smm/PS7_Dayuan Xie.jl

max_issues_repo_path: ProblemSets/PS7-smm/PS7_Dayuan Xie.jl

max_forks_repo_path: ProblemSets/PS7-smm/PS7_Dayuan Xie.jl


text:
# First run some mcmc sampling, then test the util functions..

f_ec = (x::Vector{Float64}) -> ( sin(x[1])^2 + sin(x[2])^2 ) / ( 0.1 + x[1]*x[1] + x[2]*x[2] )^(1/2.5)
f_double_ec = (x::Vector{Float64}) -> max( f_ec( x+[pi/4;pi/4] ) , f_ec( x-[pi/4;pi/4] ) )

xlim = [-20.,20.]
ylim = [-20.,20.]
bounding_box_G =[eye(2);-eye(2)]
bounding_box_h =[-xlim[1];xlim[2];-ylim[1];ylim[2]]

# and create likelihood from "inverted egg crate"
f_likelihood = (x::Vector{Float64}) -> f_ec(x)
f_likelihood_2p = (x,y) -> f_ec([x;y])
#f_likelihood = (x::Vector{Float64}) -> f_double_ec(x)
#f_likelihood_2p = (x,y) -> f_double_ec([x;y])

# starting point
x0 = [5.;5.]
# pure mcmc sampling using mh:
conf_only_mh = SamplingHS3.DecorrelationConfig( 20000 , 1.0*eye(2) , 0.5 , 0.5 )
f_density    = f_likelihood
# this is quite slow.. most time is spent for the generation of the multivar. normal random numbers..
tic()
mcmc_result = SamplingHS3.run_mcmc( x0 , f_density(x0) , f_density , (x) -> 0 , 1.0 , bounding_box_G  , bounding_box_h , conf_only_mh )
toc()

# Test utilts:
z_sampled      = SamplingHS3.bin2d( mcmc_result.dd.xx , 40 , 40 )
z_sampled_mean = SamplingHS3.bin2d( mcmc_result.dd.xx , mcmc_result.dd.xx_fdens , 40 , 40 , mean)

# And test interp1_lin:
xx_01  = collect(linspace(0.0,10.0,2000))
yy_01  = sin.(xx_01)
x_test = [1.0;2.0;5.0;5.23]
y_test = SamplingHS3.interp1_lin(xx_01,yy_01, x_test )
@test maximum( abs.( y_test - sin.(x_test) ) ) < 0.001


text:
module FemAdjoint

using LinearAlgebra, StatsBase, GroupSlices
using ForwardDiff, ReverseDiff, SparseDiffTools, SparseArrays, Symbolics
using FileIO, JLD2, UnicodePlots

include("assembly.jl")
include("plot.jl")
# different test cases
include("costproduv.jl")
#-------------------------------------------------------------------------------
function costnormU(pin, t)
    p = if typeof(pin) <: Vector
        reshape(pin, length(pin) ÷ 2, 2)
    elseif typeof(pin) <: Matrix
        pin
    end
    np = size(p, 1)
    # assemble matrices
    I, J = indKM_sparse(t)
    SK, SM = FemAdjoint.assembKM_P12D(p, t)
    K = sparse(I, J, SK, np, np)
    Ib = unique(FemAdjoint.btri(t)[:])
    K += sparse(Ib, Ib, fill(1e8, length(Ib)), np, np) # Dirichlet conditions
    M = sparse(I, J, SM, np, np)
    F = M * ones(np)
    #F = ones(np)
    U = K \ F
    # L² cost
    val = sum(U[k] ^ 2 for k = 1:np) / np
    return val
end
#-------------------------------------------------------------------------------
function ∇costnormU(pin, t;
                    jacK::SparseMatrixCSC{Float64, Int64} = spzeros(0, 0),
                    jacM::SparseMatrixCSC{Float64, Int64} = spzeros(0, 0))
    p = if typeof(pin) <: Vector
        reshape(pin, length(pin) ÷ 2, 2)
    elseif typeof(pin) <: Matrix
        pin
    end
    np = size(p, 1)
    # assemble matrices
    I, J = indKM_sparse(t)
    SK, SM = FemAdjoint.assembKM_P12D(p, t)
    K = sparse(I, J, SK, np, np)
    Ib = unique(FemAdjoint.btri(t)[:])
    K += sparse(Ib, Ib, fill(1e8, length(Ib)), np, np) # Dirichlet conditions
    M = sparse(I, J, SM, np, np)
    # compute state function
    F = M * ones(np)
    #F = ones(np)
    U = K \ F
    # plot
    plotunicode(U,  title = "U")
    # jacobian structures
    fK(y, x) = FemAdjoint.assembK_P12D_inplace!(x, t, y)
    fM(y, x) = FemAdjoint.assembM_P12D_inplace!(x, t, y)
    if length(jacK) == 0
        y = similar(FemAdjoint.assembKM_P12D(p, t)[1])
        sparsity_pattern = Symbolics.jacobian_sparsity(fK, y, p[:])
        jacK = Float64.(sparse(sparsity_pattern))
    end
    if length(jacM) == 0
        y = similar(FemAdjoint.assembKM_P12D(p, t)[2])
        sparsity_pattern = Symbolics.jacobian_sparsity(fM, y, p[:])
        jacM = Float64.(sparse(sparsity_pattern))
    end
    # evaluate jacobians
    colorsK = matrix_colors(jacK)
    JK = forwarddiff_color_jacobian!(jacK, fK, p[:], colorvec = colorsK)
    colorsM = matrix_colors(jacM)
    JM = forwarddiff_color_jacobian!(jacM, fM, p[:], colorvec = colorsM)
    # cost gradient and adjoint function
    JFU = - 2 * U / np
    λ = K' \ JFU
    plotunicode(λ, colormap = :plasma, title = "λ")
    # compute full gradient
    g = zeros(np, 2)
    for (i, j, s) ∈ zip(findnz(JK)...)
        g[j] += λ[J[i]] * s * U[I[i]]
    end
    for (i, j, s) ∈ zip(findnz(JM)...)
        g[j] -= λ[J[i]] * s
    end
    return g
end
#-------------------------------------------------------------------------------
function ∇costnormU_direct(pin, t;
                           plotsolutions::Bool = false)
    p = if typeof(pin) <: Vector
        reshape(pin, length(pin) ÷ 2, 2)
    elseif typeof(pin) <: Matrix
        pin
    end
    np = size(p, 1)
    # assemble matrices
    I, J = indKM_sparse(t)
    SK, SM = FemAdjoint.assembKM_P12D(p, t)
    K = sparse(I, J, SK, np, np)
    Ib = unique(FemAdjoint.btri(t)[:])
    K += sparse(Ib, Ib, fill(1e8, length(Ib)), np, np) # Dirichlet conditions
    M = sparse(I, J, SM, np, np)
    # compute state function
    F = M * ones(np)
    #F = ones(np)
    U = K \ F
    # plot
    plotsolutions && plotunicode(U,  title = "U")
    # cost gradient and adjoint function
    JFU = - 2 * U / np
    λ = K' \ JFU
    plotsolutions && plotunicode(λ, colormap = :plasma, title = "λ")
    # compute full gradient by automatic differentiation
    #∇prodgrad = ForwardDiff.gradient(x -> prodgrad(x, t, U, λ), p)
    ∇prodgrad = ReverseDiff.gradient(x -> prodgrad(x, t, U, λ), p)
    return ∇prodgrad
end
#-------------------------------------------------------------------------------
function prodgrad(pin, t, U, λ)
    p = reshape(pin, length(pin) ÷ 2, 2)
    np = size(p, 1)
    # assemble matrices
    I, J = indKM_sparse(t)
    SK, SM = FemAdjoint.assembKM_P12D(p, t)
    val = 0.
    for (i, j, s) ∈ zip(I, J, SK)
        val += λ[i] * s * U[j]
    end
    for (i, j, s) ∈ zip(I, J, SM)
        val -= λ[i] * s
    end
    return val
end

end # module


text:
using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.CGen.setCreateMain(true)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t2 = aggregate(t1,:userid, :ss=sum(:val2))
    #t2 = aggregate(t1,:userid, :ss=sum(:val2), :aa=length(:val2>2.0))
    #t2 = aggregate(t1,:userid, :ss=sum(:val2), :aa=length(:val2>2.0), :bb=length(:val2==1.0))
    return t2[:userid], t2[:ss]
end

@acc hpat function test2(file_name)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    #t2 = aggregate(t1,:userid, :ss=sum(:val2))
    #t2 = aggregate(t1,:userid, :ss=sum(:val2), :aa=length(:val2>2.0))
    t2 = aggregate(t1,:userid, :ss=sum(:val2), :aa=length(:val2>2.0), :bb=length(:val2==1.0))
    return t2[:userid], t2[:ss]
end

using HDF5
using MPI
if MPI.Comm_rank(MPI.COMM_WORLD)==0
    h5write("test1_1.hdf5", "/userid", [1, 2, 3, 1, 2])
    h5write("test1_1.hdf5", "/val2", [1.1, 2.1, 3.1, 3.2, 1.9])
end
a,b = test1("test1_1.hdf5")
println(a)
println(b)

using Base.Test
@test a==[1,2,3]
@test_approx_eq b [4.3, 4.0, 3.1]

a,b = test2("test1_1.hdf5")
@test a==[1,2,3]
@test_approx_eq b [4.3, 4.0, 3.1]

if MPI.Comm_rank(MPI.COMM_WORLD)==0
    rm("test1_1.hdf5")
end


text:
"""
```
gr(r::Vector{Float64}, count::Vector{Float64}, density::Float64, binstep::Float64)
```

Computes the radial distribution function from the count data and
the density.

This is exactly a conventional g(r)
if a single atom was chosen as the solute and solvent selections.

Returns both the g(r) and the kb(r)

"""
function gr(r::Vector{Float64}, count::Vector{Float64}, density::Float64, binstep::Float64)
  nbins = length(r)
  gr = zeros(nbins)
  kb = zeros(nbins)
  for i in 1:nbins
    gr[i] = (count[i]/sphericalshellvolume(i,binstep))/density
    if i == 1
      kb[i] = 4π*(gr[i]-1)*r[i]^2*binstep
    else
      kb[i] = kb[i-1] + 4π*(gr[i]-1)*r[i]^2*binstep
    end
  end
  @. kb = units.Angs3tocm3permol * kb
  return gr, kb
end

"""

```
gr(R::Result) = gr(R.d,R.rdf_count,R.density.solvent_bulk,R.options.binstep)
```

If a Result structure is provided without further details, use the rdf count
and the bulk solvent density.
"""
gr(R::Result) = gr(R.d,R.rdf_count,R.density.solvent_bulk,R.options.binstep)


text:
### OFFSETS ON GLRM
function add_offset!(glrm::AbstractGLRM)
    glrm.rx, glrm.ry = lastentry1(glrm.rx), lastentry_unpenalized(glrm.ry)
    return glrm
end

function equilibrate_variance!(glrm::AbstractGLRM, columns_to_scale = 1:size(glrm.A,2))
  for i in columns_to_scale
    nomissing = glrm.A[glrm.observed_examples[i],i]
    if length(nomissing)>1
      varlossi = avgerror(glrm.losses[i], nomissing)
    else
      varlossi = 1
    end
    if varlossi > 0
      # rescale the losses and regularizers for each column by the inverse of the empirical variance
      scale!(glrm.losses[i], scale(glrm.losses[i])/varlossi)
    end
  end
  glrm
end


text:
@testset "scale" begin
    data = ones(2, 2, 3)
    scales = [3, 2, 1]
    sc_cube = scale(data, scales)
    @test size(sc_cube) == size(data)
    @test all(sc_cube .≈ 1)

    for i in axes(data, 3)
        @test scale(data[:, :, i], scales[i]) ≈ sc_cube[:, :, i]
    end
end

@testset "invscale" begin
    data = ones(8, 8, 3)
    scales = [4, 2, 1]
    sc_cube = invscale(data, scales)
    @test size(sc_cube) == size(data)
    @test all(sc_cube[4:5, 4:5, 1] .≈ 1)
    @test all(sc_cube[3:6, 3:6, 2] .≈ 1)
    @test all(sc_cube[:, :, 3] .≈ 1)
    @test sum(sc_cube, dims=(1, 2)) |> vec ≈ [4, 16, 64]

    for i in axes(data, 3)
        @test invscale(data[:, :, i], scales[i]) ≈ sc_cube[:, :, i]
    end
end

@testset "stack and collapse" begin
    data = ones(8, 8, 3, 4)
    scales = [4, 2, 1]
    stack = scale_and_stack(data, scales)
    @test size(stack) == (8, 8, 12)
    @test all(stack .≈ 1)
    @test sum(stack) ≈ length(data)

    cube = invscale_and_collapse(stack, scales; method=sum)
    @test size(cube) == (8, 8, 4)
    @test sum(cube) ≈ 336
end

@testset "scale list" begin
    waves = [100, 200, 300]
    @test scale_list(waves) ≈ [3, 1.5, 1]
end


text:
# Similar to test for standard A*
using Graphs
using Test

# The "touring Romania" example from Russell and Norvig
g1_wedges = [
    (1, 20, 75.),
    (1, 16, 140.), #2, Arad -- Sibiu
    (1, 17, 218.),
    (20, 13, 71.),
    (13, 16, 151.),
    (16, 6, 99.),
    (16, 15, 80.), #7, Sibiu -- Rimnicu Vilcea
    (6, 2, 211.),
    (17, 10, 111.),
    (17, 2, 500.0),
    (10, 11, 70.),
    (11, 4, 75.),
    (4, 3, 120.),
    (3, 15, 146.),
    (3, 14, 138.),
    (15, 14, 97.), #15, Rimnicu Vilcea -- Pitesti
    (14, 2, 101.), #16, Pitesti -- Bucharest
    (2, 7, 90.),
    (2, 18, 85.),
    (18, 8, 98.),
    (8, 5, 86.),
    (18, 19, 142.),
    (19, 9, 92.),
    (9, 12, 87.) ]

g1_heuristics = [
        366., 0., 160., 242., 161., 210., 77., 151., 226., 244., 218., 234., 300., 100., 193.,
        253., 329., 80., 199., 374. ]

# Weight constraints
# The "touring Romania" example from Russell and Norvig
c1_wedges = [
    (1, 20, 75.),
    (1, 16, 140.), #2, Arad -- Sibiu
    (1, 17, 218.),
    (20, 13, 71.),
    (13, 16, 151.),
    (16, 6, 99.),
    (16, 15, 80.), #7, Sibiu -- Rimnicu Vilcea
    (6, 2, 211.),
    (17, 10, 111.),
    (17, 2, 500.0),
    (10, 11, 70.),
    (11, 4, 75.),
    (4, 3, 120.),
    (3, 15, 146.),
    (3, 14, 138.),
    (15, 14, 97.), #15, Rimnicu Vilcea -- Pitesti
    (14, 2, 101.), #16, Pitesti -- Bucharest
    (2, 7, 90.),
    (2, 18, 85.),
    (18, 8, 98.),
    (8, 5, 86.),
    (18, 19, 142.),
    (19, 9, 92.),
    (9, 12, 87.) ]

c1_heuristics = [
            166., 0, 110., 142., 131., 156., 47., 121., 206., 174., 191., 214., 180., 90., 133.,
            213., 269., 50., 159., 294. ]

nbrs_wt_dict = Dict{Int64, Dict{Int64, Float64}}()
for (u, v, wt) in g1_wedges
    if haskey(nbrs_wt_dict, u)
        nbrs_wt_dict[u][v] = wt
    else
        nbrs_wt_dict[u] = Dict(v => wt)
    end
end

# Setup edge weight function based on dictionary
edge_wt_fn(u::Int64, v::Int64) = nbrs_wt_dict[u][v]

nbrs_constr_dict = Dict{Int64, Dict{Int64, Float64}}()
for (u, v, wt) in c1_wedges
    if haskey(nbrs_constr_dict, u)
        nbrs_constr_dict[u][v] = wt
    else
        nbrs_constr_dict[u] = Dict(v => wt)
    end
end
edge_constr_fn(u::Int64, v::Int64) = nbrs_constr_dict[u][v]
weight_heur(n) = c1_heuristics[n]

weight_functions = [edge_constr_fn]
weight_heuristics = [weight_heur]

struct TestVisitorImplicit <: AbstractDijkstraVisitor
    nbrs_wt_dict::Dict{Int64, Dict{Int64, Float64}}
    test_graph::SimpleVListGraph
    goal_vtx::Int64
end

function Graphs.include_vertex!(vis::TestVisitorImplicit, u::Int64, v::Int64,
                                d::Float64, nbrs::Vector{Int64})

    # Terminate if goal vertex reached
    if v == vis.goal_vtx
        return false
    end

    # Otherwise, insert all neighbours from dict
    # AND add to graph if not already there
    for n in collect(keys(vis.nbrs_wt_dict[v]))
        if ~(n in vis.test_graph.vertices)
            add_vertex!(test_graph, n)
        end

        # Neighbors will have INDICES of vertices, not vertices themselves
        push!(nbrs, vertex_index(vis.test_graph, n))
    end

    return nbrs
end

weight_constraints = [448.0]

start = 1
start_idx = 1
goal = 2
eps_weight = 1.5
test_graph = SimpleVListGraph([start])

# Create visitor with goal vertex 2 and heuristic
vis = TestVisitorImplicit(nbrs_wt_dict, test_graph, goal)
admissible_heuristic(n) = g1_heuristics[n]
focal_state_heuristic(n) = 1000.0 - g1_heuristics[n]
focal_transition_heuristic(u, v) = 1000.0 - nbrs_wt_dict[u][v]

states, tgt_entry = a_star_implicit_constrained_epsilon_path!(test_graph, edge_wt_fn,
                                                            start_idx, vis, eps_weight,
                                                            admissible_heuristic,
                                                            focal_state_heuristic,
                                                            focal_transition_heuristic,
                                                            weight_functions,
                                                            weight_heuristics,
                                                            weight_constraints)

sp_idxs, costs, wts = shortest_path_cost_weights(states.a_star_states, test_graph, start, tgt_entry)
sp = [test_graph.vertices[s] for s in sp_idxs]

@test sp == [1, 16, 15, 14, 2]

# For inf weight constraint, reduces to a* constrained - get [1, 17, 2]
# For weight = 650. get [1, 16, 6, 2]


text:
# ------------------------------------------------------
_sep_freeids(freeid, freeids...) = (freeid, freeids)

function _fva_range(net, box, freeid)
    lb, ub = fva(net, freeid)
    ran = dim_range(box, freeid)
    lbi = findapproxi(lb, ran)
    ubi = findapproxi(ub, ran)
    return ran[lbi:ubi]
end

# remove unfeasibles and trivial
_is_feasible(net) = any(fba(net) .!= 0.0)
function _fill_subspace!(containers::Vector,
        net::MetNet, box, left_freeids, head;
        filter::Function
    )

    curr_freeid, left_freeids = _sep_freeids(left_freeids...)
    curr_free_range = _fva_range(net, box, curr_freeid)

    if isempty(left_freeids)
        # base case
        for cur_freev in curr_free_range
            freev = [head; cur_freev]
            !filter(net, freev) && continue
            # all(iszero.(freev)) && continue # remove trivial
            for (vi, containeri) in zip(freev, containers)
                push!(containeri, vi)
            end
        end
    else
        # recursive step
        for cur_freev in curr_free_range
            _head = [head; cur_freev]
            fixxing(net, curr_freeid, cur_freev) do
                _fill_subspace!(containers, net, box, left_freeids, _head; filter)
            end
        end
    end
end

# one container per dim
_build_containers(freeids, jump_size = 10_000) = [Container{Float64}(jump_size) for free in freeids]

function _containers_pool_size(freeids, containerss_pool)
    lens = zeros(Int, length(freeids))
    for containers in containerss_pool
        for i in eachindex(freeids)
            lens[i] += length(containers[i])
        end
    end
    return lens
end

function _reduce_containers_pool!(containerss_pool)

    # one per thread
    acc_containers = Dict()
    for containers in containerss_pool
        # one vector per dimension
        for i in eachindex(containers)
            !haskey(acc_containers, i) && (acc_containers[i] = containers[i]; continue)
            push!(acc_containers[i], vec!(containers[i])...)
            empty!(containers[i])
        end
    end
    map(1:length(acc_containers)) do i
        vec!(acc_containers[i])
    end
end

function subspace(net::MetNet, box::BoxGrid, freeids::Vector{Symbol};
        # stoitof = 0.01,
        filter::Function = (net, v) -> true,
        nthrs::Int = nthreads()
    )

    # Base
    if length(freeids) == 1
        _fill_subspace!(_build_containers(freeids), net, box, freeids, Float64[]; filter)
    end

    first_freeid, left_freeids = _sep_freeids(freeids...)
    first_free_range = _fva_range(net, box, first_freeid)

    # recursive step
    containerss_pool = [_build_containers(freeids) for _ in 1:nthrs]
    nets_pool = [deepcopy(net) for _ in 1:nthrs]

    prog = Progress(length(first_free_range); dt = 0.5, desc = "Collecting... ")
    function showvalues()
        estimation = sum(_containers_pool_size(freeids, containerss_pool))

        return [
            ("polV estimation", estimation),
            ("thid", threadid())
        ]
    end

    @threads for cur_freev in first_free_range

        thid = threadid()
        containers = containerss_pool[thid]
        net_th = nets_pool[thid]
        head = [cur_freev]

        fixxing(net_th, first_freeid, cur_freev) do
            _fill_subspace!(containers, net_th, box, left_freeids, head; filter)
        end
        next!(prog; showvalues)
    end
    finish!(prog)

    _reduce_containers_pool!(containerss_pool)
end


text:
using Base.Test
using DataFrames

let
	d = dataeye(3, 3)
	d[1, 1] = NA

	svd(d)
end


text:
module SocialTest
export sdcompare, plotsdcompare

import SocialDynamics: treefuncexp, socdyncoupexp, treefuncad, socdyncoupad
import SocioHydrology: SocioFEHM
using PyPlot
#using Gadfly
# Compare tree expert model, non-tree expert model, tree ad model, and non-tree ad model.

"""
Comparing all social dynamics functions (see SocialDynamics for details of each)

Arguments:

- `soc_struc` : Hierarchy tree branching structure (bottom up)
- `pcons` : consistency rate (also 1 - conformity rate)
- `adperiod` : number of iterations for each media event
- `simtime` : number of social interaction iterations.
- `tree_initurate` : percentage of initial water users in hierarchical models

Returns:

- `usage` : number of users for each social dynamics model for all simtime steps.
		`[expert_3%, expert_6%, hier_expert_3%, hier_expert_6%, media,  hier_media] * simtime`

"""

function sdcompare(;ss = [2, 5, 5, 5, 5, 2], pcons = 0.35, adperiod = 100, simtime = 100, tree_initurate = 0)

	popul = 1;
	for i = 2:length(ss)
		popul = popul * ss[i];
	end

	atr = ss[1];			# Number of attributes
	sims = 6;
	steps = 1;
	soc = zeros(Bool, popul, atr ,sims);
	usage = popul.* ones(Int32, sims, 1);
	kk = Array(Int32, sims, 1);
	af = 0;

	for i=1:simtime
		if(mod(i-1,adperiod) == 0)
			af=1;
		else
			af=0;
		end
		soc[:,:,1] = socdyncoupexp(popul, round(Int, 0.03*popul), true, soc[:,:,1], pcons, steps);
		soc[:,:,2] = socdyncoupexp(popul, round(Int, 0.06*popul), true, soc[:,:,2], pcons, steps);
		soc[:,:,3] = treefuncexp(ss, round(Int, 0.03*popul), soc[:,:,3], pcons, steps, 0);
		soc[:,:,4] = treefuncexp(ss, round(Int, 0.06*popul), soc[:,:,4], pcons, steps, 0);
		soc[:,:,5] = socdyncoupad(popul, af, true, 0.2, soc[:,:,5], pcons, steps);
		soc[:,:,6] = treefuncad(ss, af, soc[:,:,6], pcons, steps, 0, 0.2);
		uu = popul .- sum(soc, 1)[1, 2, :];
		kk = reshape(uu, sims, 1);
		usage = hcat(usage , kk);

	end

	return usage

end

##############################################################

function plotsdcompare(;trials = 5, ss = [2, 5, 5, 5, 5, 2], pcons = 0.35, adperiod = 100, simtime = 100, tree_initurate = 0)

	sims = 6;
	totusage = zeros(Float32, sims, simtime + 1);
	usage = zeros(Int32, sims, simtime + 1);
	for j = 1:trials
		println( string("trial number: ", j) );
		usage = sdcompare(ss=ss, pcons=pcons, adperiod=adperiod, simtime=simtime, tree_initurate=tree_initurate)
		totusage = totusage + usage;
	end
	avgusage = totusage./trials;

	# Plot
	simarray = 1:(simtime + 1);

	fig, ax = subplots(figsize = (11, 9));
	p1 = plot(simarray, usage[1, :]', color = "red", label = "expert_3%", linewidth = 1.0, linestyle = "-")
	p2 = plot(simarray, usage[2, :]', color = "red", label = "expert_6%", linewidth = 1.0, linestyle = "--")
	p3 = plot(simarray, usage[3, :]', color = "black", label = "hier_expert_3%", linewidth = 1.0, linestyle = "-")
	p4 = plot(simarray, usage[4, :]', color = "black", label = "hier_expert_6%", linewidth = 1.0, linestyle = "--")
	p5 = plot(simarray, usage[5, :]', color = "blue", label = "media", linewidth = 1.0, linestyle = "-")
	p6 = plot(simarray, usage[6, :]', color = "blue", label = "hier_media", linewidth = 1.0, linestyle = "--")
	ax[:legend]()

	xlabel("Full-social interactions")
	ylabel("Number of water users")
	title("One realization Social Dynamics")

	savefig("OneRealSD.pdf")

	figm, axm = subplots(figsize = (11, 9));
	p1m = plot(simarray, avgusage[1, :]', color = "red", label = "expert_3%", linewidth = 1.0, linestyle = "-")
	p2m = plot(simarray, avgusage[2, :]', color = "red", label = "expert_6%", linewidth = 1.0, linestyle = "--")
	p3m = plot(simarray, avgusage[3, :]', color = "black", label = "hier_expert_3%", linewidth = 2.0, linestyle = "-")
	p4m = plot(simarray, avgusage[4, :]', color = "black", label = "hier_expert_6%", linewidth = 2.0, linestyle = "--")
	p5m = plot(simarray, avgusage[5, :]', color = "blue", label = "media", linewidth = 1.0, linestyle = "-")
	p6m = plot(simarray, avgusage[6, :]', color = "blue", label = "hier_media", linewidth = 2.0, linestyle = "--")
	axm[:legend]()

	xlabel("Full-social interactions")
	ylabel("Number of water users")
	title("Mean Social Dynamics for $trials trials")

	savefig("MeanSD.pdf")

	plt[:show]()

# 	Gadfly
#	p1 = plot(layer(x = simarray , y = usage[1,:], Geom.line),
#		  layer(x = simarray , y = usage[2,:], Geom.line),
#		  layer(x = simarray , y = usage[3,:], Geom.line),
#		  layer(x = simarray , y = usage[4,:], Geom.line),
#		  layer(x = simarray , y = usage[5,:], Geom.line),
#		  layer(x = simarray , y = usage[6,:], Geom.line),
#		 Guide.xlabel("Full-social Interactions"),
#			    Guide.ylabel("Number of water users"),
#			    Guide.title("One Realiation Social Dynamics"));
#
#	p2 = plot(layer(x = simarray , y = avgusage[1,:], Geom.line),
#		  layer(x = simarray , y = avgusage[2,:], Geom.line),
#		  layer(x = simarray , y = avgusage[3,:], Geom.line),
#		  layer(x = simarray , y = avgusage[4,:], Geom.line),
#		  layer(x = simarray , y = avgusage[5,:], Geom.line),
#		  layer(x = simarray , y = avgusage[6,:], Geom.line),
#		 Guide.xlabel("Full-social Interactions"),
#			    Guide.ylabel("Number of water users"),
#			    Guide.title("Mean Social Dynamics"));
#
#	draw(PDF("socialtests.pdf", 6inch, 6inch), vstack(p1, p2))

end

##############################################################

function plotRiskvsConc()

	times1, adornumexp1, ucount1, conchist1, riskhist1 = SocioFEHM(conc_risk = 0, exp_adv = 0, exprate_advper = 0.03 , simulationtime = 70, pid = [1, 0.2, 0.01], init_urate = 0.5, popul= 10000, agelimit = 80, risk_percentile = .95, FEHMprefix = "cm", outputfile = "conc_control.txt")

	times2, adornumexp2, ucount2, conchist2, riskhist2 = SocioFEHM(conc_risk = 1, exp_adv = 0, exprate_advper = 0.03 , simulationtime = 70, pid = [1, 0.2, 0.01], init_urate = 0.5, popul= 10000, agelimit = 80, risk_percentile = .95, FEHMprefix = "cm", outputfile = "risk_control.txt")

	figc, axc = subplots(figsize = (11, 9));
	p1c = plot(times1, conchist1, color = "red", label = "Concentration control", linewidth = 2.0, linestyle = "--")
	p2c = plot(times2, conchist2, color = "blue", label = "Risk control", linewidth = 2.0, linestyle = "-")
	axc[:legend]()

	xlabel("Time (years)")
	ylabel("Concentration (ppb)")
	title("Concentration versus Risk control")

	savefig("ConcvsRisk.pdf")

	plt[:show]()

end

##################################################################

function plotHomovsHetero()

	println("Homogeneous Aquifer")
	times1, adornumexp1, ucount1, conchist1, riskhist1 = SocioFEHM(conc_risk = 0, exp_adv = 0, exprate_advper = 0.03 , simulationtime = 70, pid = [1, 0.2, 0.01], init_urate = 0.5, popul= 10000, agelimit = 80, risk_percentile = .95, FEHMprefix = "cm", outputfile = "homogeneous.txt")

	println("Heterogeneous Aquifer")
	times2, adornumexp2, ucount2, conchist2, riskhist2 = SocioFEHM(conc_risk = 0, exp_adv = 0, exprate_advper = 0.03 , simulationtime = 70, pid = [1, 0.2, 0.01], init_urate = 0.5, popul= 10000, agelimit = 80, risk_percentile = .95, FEHMprefix = "cmgsv02", outputfile = "heterogeneous.txt")

	figc, axc = subplots(figsize = (11, 9));
	p1c = plot(times1, conchist1, color = "red", label = "Homogeneous", linewidth = 2.0, linestyle = "--")
	p2c = plot(times2, conchist2, color = "blue", label = "Heterogeneous", linewidth = 2.0, linestyle = "-")
	axc[:legend]()

	xlabel("Time (years)")
	ylabel("Concentration (ppb)")
	title("Homogeneous vs. Heterogeneous Aquifer")

	savefig("HomovsHetero.pdf")

	plt[:show]()

end

# # # #

end

max_stars_repo_name: josephbakarji/socio-hydro-abm

max_issues_repo_name: josephbakarji/socio-hydro-abm

max_forks_repo_name: josephbakarji/socio-hydro-abm