text:
abstract type QR end

struct QrWs{T<:Number} <: QR
    tau::Vector{T}
    work::Vector{T}
    lwork::Ref{BlasInt}
    info::Ref{BlasInt}
end

for (geqrf, ormqr, elty) in (
    (:dgeqrf_, :dormqr_, :Float64),
    (:sgeqrf_, :sormqr_, :Float32),
    (:zgeqrf_, :zormqr_, :ComplexF64),
    (:cgeqrf_, :cormqr_, :ComplexF32),
)

    @eval begin

        function QrWs(A::StridedMatrix{T}) where {T<:$elty}
            nn, mm = size(A)
            m = Ref{BlasInt}(mm)
            n = Ref{BlasInt}(nn)
            RldA = Ref{BlasInt}(max(1, stride(A, 2)))
            tau = Vector{T}(undef, min(nn, mm))
            work = Vector{T}(undef, 1)
            lwork = Ref{BlasInt}(-1)
            info = Ref{BlasInt}(0)
            ccall(
                (@blasfunc($geqrf), liblapack),
                Nothing,
                (
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{T},
                    Ref{BlasInt},
                    Ptr{T},
                    Ptr{T},
                    Ref{BlasInt},
                    Ref{BlasInt},
                ),
                m,
                n,
                A,
                RldA,
                tau,
                work,
                lwork,
                info,
            )
            chklapackerror(info[])
            lwork = Ref{BlasInt}(real(work[1]))
            work = Array{T}(undef, lwork[])
            QrWs(tau, work, lwork, info)
        end

        function geqrf_core!(A::StridedMatrix{$elty}, ws::QrWs)
            mm, nn = size(A)
            m = Ref{BlasInt}(mm)
            n = Ref{BlasInt}(nn)
            RldA = Ref{BlasInt}(max(1, stride(A, 2)))
            ccall(
                (@blasfunc($geqrf), liblapack),
                Nothing,
                (
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ref{BlasInt},
                ),
                m,
                n,
                A,
                RldA,
                ws.tau,
                ws.work,
                ws.lwork,
                ws.info,
            )
            chklapackerror(ws.info[])
        end

        t1 = StridedMatrix{$elty}
        t2 = Transpose{$elty,<:StridedMatrix}
        t3 = Adjoint{$elty,<:StridedMatrix}
    end
end

for (geqrf, ormqr, elty) in ((:dgeqrf_, :dormqr_, :Float64), (:sgeqrf_, :sormqr_, :Float32))

    @eval begin
        function ormqr_core!(
            side::Char,
            A::StridedMatrix{$elty},
            C::StridedMatrix{$elty},
            ws::QR,
        )
            mm, nn = size(C)
            m = Ref{BlasInt}(mm)
            n = Ref{BlasInt}(nn)
            k = Ref{BlasInt}(length(ws.tau))
            RldA = Ref{BlasInt}(max(1, stride(A, 2)))
            RldC = Ref{BlasInt}(max(1, stride(C, 2)))
            ccall(
                (@blasfunc($ormqr), liblapack),
                Nothing,
                (
                    Ref{UInt8},
                    Ref{UInt8},
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ref{BlasInt},
                ),
                side,
                'N',
                m,
                n,
                k,
                A,
                RldA,
                ws.tau,
                C,
                RldC,
                ws.work,
                ws.lwork,
                ws.info,
            )
            chklapackerror(ws.info[])
        end
    end

    @eval begin
        t1 = Transpose{$elty,<:StridedMatrix{$elty}}
        t2 = Adjoint{$elty,<:StridedMatrix{$elty}}
    end

    for elty2 in (t1, t2)

        @eval begin
            function ormqr_core!(side::Char, A::$elty2, C::StridedMatrix{$elty}, ws::QR)
                mm, nn = size(C)
                m = Ref{BlasInt}(mm)
                n = Ref{BlasInt}(nn)
                k = Ref{BlasInt}(length(ws.tau))
                RldA = Ref{BlasInt}(max(1, stride(A.parent, 2)))
                RldC = Ref{BlasInt}(max(1, stride(C, 2)))
                ccall(
                    (@blasfunc($ormqr), liblapack),
                    Nothing,
                    (
                        Ref{UInt8},
                        Ref{UInt8},
                        Ref{BlasInt},
                        Ref{BlasInt},
                        Ref{BlasInt},
                        Ptr{$elty},
                        Ref{BlasInt},
                        Ptr{$elty},
                        Ptr{$elty},
                        Ref{BlasInt},
                        Ptr{$elty},
                        Ref{BlasInt},
                        Ref{BlasInt},
                    ),
                    side,
                    'T',
                    m,
                    n,
                    k,
                    A.parent,
                    RldA,
                    ws.tau,
                    C,
                    RldC,
                    ws.work,
                    ws.lwork,
                    ws.info,
                )
                chklapackerror(ws.info[])
            end
        end
    end
end

#=
for (geqrf, ormqr, elty) in
    ((:zgeqrf_, :zormqr_, :ComplexF64), (:cgeqrf_, :cormqr_, :ComplexF32))

    @eval begin
        function ormqr_core!(
            side::Char,
            A::StridedMatrix{$elty},
            C::StridedMatrix{$elty},
            ws::QR,
        )
            mm, nn = size(C)
            m = Ref{BlasInt}(mm)
            n = Ref{BlasInt}(nn)
            k = Ref{BlasInt}(length(ws.tau))
            RldA = Ref{BlasInt}(max(1,stride(A,2)))
            RldC = Ref{BlasInt}(max(1,stride(C,2)))
            ccall((@blasfunc($ormqr),  liblapack), Nothing,
                  (Ref{UInt8},Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},
                   Ptr{$elty},Ptr{$elty},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},Ref{BlasInt}),
                  side, 'N', m, n, k, A, RldA, ws.tau, C, RldC, ws.work, ws.lwork, ws.info)
            chklapackerror(ws.info[])
        end
    end

    @eval begin
        t1 = Transpose{$elty,<:StridedMatrix{$elty}}
        t2 = Adjoint{$elty,<:StridedMatrix{$elty}}
    end

    for (elty2, transchar) in ((t2, 'T'), (t3, 'C'))

        @eval begin
            function ormqr_core!(side::Char, A::$elty2, C::StridedMatrix{$elty}, ws::QR)
                mm, nn = size(C)
                m = Ref{BlasInt}(mm)
                n = Ref{BlasInt}(nn)
                k = Ref{BlasInt}(length(ws.tau))
                RldA = Ref{BlasInt}(max(1,stride(A.parent,2)))
                RldC = Ref{BlasInt}(max(1,stride(C,2)))
                ccall((@blasfunc($ormqr),  liblapack), Nothing,
                      (Ref{UInt8},Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},
                       Ptr{$elty},Ptr{$elty},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},Ref{BlasInt}),
                      side, $transchar, m, n, k, A.parent, RldA,
                      ws.tau, C, RldC, ws.work, ws.lwork, ws.info)
                chklapackerror(ws.info[])
            end
        end
    end
end
=#

struct QrpWs{T<:Number} <: QR
    tau::Vector{T}
    jpvt::Vector{BlasInt}
    work::Vector{T}
    lwork::BlasInt
    info::Ref{BlasInt}
end

for (geqp3, elty) in (
    (:dgeqp3_, :Float64),
    (:sgeqp3_, :Float32),
    (:zgeqp3_, :ComplexF64),
    (:cgeqp3_, :ComplexF32),
)

    @eval begin

        function QrpWs(A::StridedMatrix{$elty})
            m, n = size(A)
            RldA = BlasInt(max(1, stride(A, 2)))
            jpvt = zeros(BlasInt, n)
            tau = Vector{$elty}(undef, min(m, n))
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info = Ref{BlasInt}()
            ccall(
                (@blasfunc($geqp3), liblapack),
                Nothing,
                (
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{BlasInt},
                    Ptr{$elty},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ref{BlasInt},
                ),
                m,
                n,
                A,
                RldA,
                jpvt,
                tau,
                work,
                lwork,
                info,
            )
            chklapackerror(info[])
            lwork = BlasInt(real(work[1]))
            work = resize!(work, lwork)
            QrpWs(tau, jpvt, work, lwork, info)
        end

        function geqp3!(A::StridedMatrix{$elty}, ws::QrpWs)
            m, n = size(A)
            RldA = BlasInt(max(1, stride(A, 2)))
            ccall(
                (@blasfunc($geqp3), liblapack),
                Nothing,
                (
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{BlasInt},
                    Ptr{$elty},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ref{BlasInt},
                ),
                m,
                n,
                A,
                RldA,
                ws.jpvt,
                ws.tau,
                ws.work,
                ws.lwork,
                ws.info,
            )
            chklapackerror(ws.info[])
        end

        t1 = StridedMatrix{$elty}
        t2 = Transpose{$elty,<:StridedMatrix}
        t3 = Adjoint{$elty,<:StridedMatrix}
    end
end


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

module SimpleTCAS_EvE_Impl

import Compat.ASCIIString

using AbstractGenerativeModelImpl
using AbstractGenerativeModelInterfaces
using CommonInterfaces

using Base.Test
using EncounterDBN
using PilotResponse
using DynamicModel
using WorldModel
using Sensor
using CASCoordination
using CollisionAvoidanceSystem
using Simulator
using RLESUtils, Observers

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractGenerativeModelInterfaces.get
import AbstractGenerativeModelInterfaces.isterminal

import CommonInterfaces.addObserver

export SimpleTCAS_EvE_params, SimpleTCAS_EvE, initialize, update, isterminal, addObserver

type SimpleTCAS_EvE_params
  #global params: remains constant per sim
  encounter_number::Int64 #encounter number in file
  nmac_r::Float64 #NMAC radius in feet
  nmac_h::Float64 #NMAC vertical separation, in feet
  maxSteps::Int64 #maximum number of steps in sim
  number_of_aircraft::Int64 #number of aircraft  #must be 2 for now...
  encounter_seed::UInt64 #Seed for generating encounters
  pilotResponseModel::Symbol #{:SimplePR, :StochasticLinear}

  #Defines behavior of CorrAEMDBN.  Read from file or generate samples on-the-fly
  command_method::Symbol #:DBN=sampled from DBN or :ENC=from encounter file

  #these are to define CorrAEM:
  encounter_file::ASCIIString #Path to encounter file
  initial_sample_file::ASCIIString #Path to initial sample file
  transition_sample_file::ASCIIString #Path to transition sample file

  string_id::ASCIIString

  SimpleTCAS_EvE_params() = new()
end

type SimpleTCAS_EvE <: AbstractGenerativeModel
  params::SimpleTCAS_EvE_params

  #sim objects: contains state that changes throughout sim run
  em::CorrAEMDBN
  pr::Vector{Union{SimplePilotResponse,StochasticLinearPR}}
  dm::Vector{SimpleADM}
  wm::AirSpace
  coord::GenericCoord
  sr::Vector{SimpleTCASSensor}
  cas::Vector{CoordSimpleTCAS}

  observer::Observer

  #sim states: changes throughout simulation run
  t_index::Int64 #current time index in the simulation. Starts at 1 and increments by 1.
  #This is different from t which starts at 0 and could increment in the reals.

  #empty constructor
  function SimpleTCAS_EvE(p::SimpleTCAS_EvE_params)
    @test p.number_of_aircraft == 2 #need to revisit the code if this is not true

    sim = new()
    sim.params = p

    srand(p.encounter_seed) #There's a rand inside generateEncounter, need to control it
    sim.em = CorrAEMDBN(p.number_of_aircraft, p.encounter_file, p.initial_sample_file,
                    p.transition_sample_file,
                    p.encounter_number,p.encounter_seed,p.command_method)

    if p.pilotResponseModel == :SimplePR
      sim.pr = SimplePilotResponse[ SimplePilotResponse() for i=1:p.number_of_aircraft ]
    elseif p.pilotResponseModel == :StochasticLinear
      sim.pr = StochasticLinearPR[ StochasticLinearPR() for i=1:p.number_of_aircraft ]
    else
      error("SimpleTCAS_EvE_Impl: No such pilot response model")
    end

    sim.dm = SimpleADM[ SimpleADM(number_of_substeps=1) for i=1:p.number_of_aircraft ]
    sim.wm = AirSpace(p.number_of_aircraft)

    sim.coord = GenericCoord(p.number_of_aircraft)
    sim.sr = SimpleTCASSensor[ SimpleTCASSensor(1), SimpleTCASSensor(2) ]
    sim.cas = CoordSimpleTCAS[ CoordSimpleTCAS(1,sim.coord), CoordSimpleTCAS(2,sim.coord) ]

    sim.observer = Observer()

    #Start time at 1 for easier indexing into arrays according to time
    sim.t_index = 1

    return sim
  end
end

addObserver(sim::SimpleTCAS_EvE, f::Function) = add_observer(sim.observer, f)
addObserver(sim::SimpleTCAS_EvE, tag::AbstractString, f::Function) = add_observer(sim.observer, tag, f)

function initialize(sim::SimpleTCAS_EvE)

  wm, aem, pr, adm, cas, sr = sim.wm, sim.em, sim.pr, sim.dm, sim.cas, sim.sr

  #Start time at 1 for easier indexing into arrays according to time
  sim.t_index = 1

  EncounterDBN.initialize(aem)

  for i = 1:sim.params.number_of_aircraft
    initial = EncounterDBN.getInitialState(aem, i)
    @notify_observer(sim.observer,"Initial",[i, sim.t_index, initial])

    state = DynamicModel.initialize(adm[i], initial)
    WorldModel.initialize(wm, i, state)

    Sensor.initialize(sr[i])
    @notify_observer(sim.observer,"Sensor",[i, sim.t_index, sr[i]])

    CollisionAvoidanceSystem.initialize(cas[i])
    @notify_observer(sim.observer,"CAS", [i, sim.t_index, cas[i]])

    PilotResponse.initialize(pr[i])
    @notify_observer(sim.observer,"Response",[i, sim.t_index, pr[i]])
  end

  @notify_observer(sim.observer,"WorldModel", [sim.t_index, wm])

  return
end

function update(sim::SimpleTCAS_EvE)
  wm, aem, pr, adm, cas, sr = sim.wm, sim.em, sim.pr, sim.dm, sim.cas, sim.sr

  sim.t_index += 1

  logProb = 0.0 #track the probabilities in this update

  cmdLogProb = EncounterDBN.update(aem)
  logProb += cmdLogProb #TODO: decompose this by aircraft?

  states = WorldModel.getAll(wm)

  for i = 1:sim.params.number_of_aircraft
    #intended command
    command = EncounterDBN.get(aem,i)
    @notify_observer(sim.observer,"Command",[i, sim.t_index, command])

    output = Sensor.update(sr[i], states)
    @notify_observer(sim.observer,"Sensor",[i, sim.t_index, sr[i]])

    RA = CollisionAvoidanceSystem.update(cas[i], output)
    @notify_observer(sim.observer,"CAS", [i, sim.t_index, cas[i]])

    response = PilotResponse.update(pr[i], command, RA)
    logProb += response.logProb #this will break if response is not SimplePRCommand
    @notify_observer(sim.observer,"Response",[i, sim.t_index, pr[i]])

    state = DynamicModel.update(adm[i], response)
    WorldModel.update(wm, i, state)

  end

  WorldModel.updateAll(wm)
  @notify_observer(sim.observer,"WorldModel", [sim.t_index, wm])

  return logProb
end

function getvhdist(wm::AbstractWorldModel)
  states = WorldModel.getAll(wm) #states::Vector{ASWMState}

  #[(vdist,hdist)]
  vhdist = [(abs(s2.h-s1.h),norm([(s2.x-s1.x),(s2.y-s1.y)])) for s1 = states, s2 = states]
  for i = 1:length(states)
    vhdist[i,i] = (typemax(Float64),typemax(Float64))
  end

  return vhdist
end

function isNMAC(sim::SimpleTCAS_EvE)
  vhdist = getvhdist(sim.wm)
  nmac_test = map((vhd)->vhd[2] <= sim.params.nmac_r && vhd[1] <= sim.params.nmac_h,vhdist)
  return any(nmac_test)
end

isterminal(sim::SimpleTCAS_EvE) = isNMAC(sim) || sim.t_index > sim.params.maxSteps

end #module

ext: jl


text:
# Rational approximations to generalized hypergeometric functions
# using Drummond's sequence transformation

# ₀F₀(;z)
function drummond0F0(z::T; kmax::Int = 10_000) where T
    if norm(z) < eps(real(T))
        return one(T)
    end
    ζ = inv(z)
    Nlo = ζ
    Dlo = ζ
    Tlo = Nlo/Dlo
    Nhi = (2ζ - 1)*Nlo + 2ζ
    Dhi = (2ζ - 1)*Dlo
    Thi = Nhi/Dhi
    k = 1
    Nhi, Nlo = ((k+2)*ζ-1)*Nhi + k*ζ*Nlo + ζ, Nhi
    Dhi, Dlo = ((k+2)*ζ-1)*Dhi + k*ζ*Dlo, Dhi
    Thi, Tlo = Nhi/Dhi, Thi
    k += 1
    while k < kmax && errcheck(Tlo, Thi, 10eps(real(T)))
        Nhi, Nlo = ((k+2)*ζ-1)*Nhi + k*ζ*Nlo, Nhi
        Dhi, Dlo = ((k+2)*ζ-1)*Dhi + k*ζ*Dlo, Dhi
        Thi, Tlo = Nhi/Dhi, Thi
        k += 1
    end
    return isfinite(Thi) ? Thi : Tlo
end

# ₁F₀(α;z)
function drummond1F0(α::T1, z::T2; kmax::Int = 10_000) where {T1, T2}
    T = promote_type(T1, T2)
    absα = abs(T(α))
    if norm(z) < eps(real(T)) || norm(α) < eps(absα)
        return one(T)
    end
    ζ = inv(z)
    Nlo = ζ/α
    Dlo = ζ/α
    Tlo = Nlo/Dlo
    Nhi = (2ζ - (α+1))*Nlo + 2ζ
    Dhi = (2ζ - (α+1))*Dlo
    Thi = Nhi/Dhi
    if norm(α+1) < eps(absα+1)
        return Thi
    end
    Nhi /= α+1
    Dhi /= α+1
    k = 1
    Nhi, Nlo = ((k+2)*ζ-(α+2k+1))*Nhi + k*(ζ-1)*Nlo + ζ, Nhi
    Dhi, Dlo = ((k+2)*ζ-(α+2k+1))*Dhi + k*(ζ-1)*Dlo, Dhi
    Thi, Tlo = Nhi/Dhi, Thi
    if norm(α+k+1) < eps(absα+k+1)
        return Thi
    end
    Nhi /= α+k+1
    Dhi /= α+k+1
    k += 1
    while k < kmax && errcheck(Tlo, Thi, 10eps(real(T)))
        Nhi, Nlo = ((k+2)*ζ-(α+2k+1))*Nhi + k*(ζ-1)*Nlo, Nhi
        Dhi, Dlo = ((k+2)*ζ-(α+2k+1))*Dhi + k*(ζ-1)*Dlo, Dhi
        Thi, Tlo = Nhi/Dhi, Thi
        if norm(α+k+1) < eps(absα+k+1)
            return Thi
        end
        Nhi /= α+k+1
        Dhi /= α+k+1
        k += 1
    end
    return isfinite(Thi) ? Thi : Tlo
end

# ₀F₁(β;z)
function drummond0F1(β::T1, z::T2; kmax::Int = 10_000) where {T1, T2}
    T = promote_type(T1, T2)
    if norm(z) < eps(real(T))
        return one(T)
    end
    ζ = inv(z)
    Nlo = β*ζ
    Dlo = β*ζ
    Tlo = Nlo/Dlo
    Nmid = ((β+1)*(2)*ζ - 1)*Nlo + (β+1)*(2)*ζ
    Dmid = ((β+1)*(2)*ζ - 1)*Dlo
    Tmid = Nmid/Dmid
    Nhi = ((β+2)*(3)*ζ - 1)*Nmid + (β+4)*ζ*Nlo + (β+4)*ζ
    Dhi = ((β+2)*(3)*ζ - 1)*Dmid + (β+4)*ζ*Dlo
    Thi = Nhi/Dhi
    k = 2
    Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-1)*Nhi + k*(β+2k+2)*ζ*Nmid + k*(k-1)*ζ*Nlo + 2ζ, Nhi, Nmid
    Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-1)*Dhi + k*(β+2k+2)*ζ*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
    Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
    k += 1
    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
        Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-1)*Nhi + k*(β+2k+2)*ζ*Nmid + k*(k-1)*ζ*Nlo, Nhi, Nmid
        Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-1)*Dhi + k*(β+2k+2)*ζ*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
        Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
        k += 1
    end
    return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo
end

# ₂F₀(α,β;z)
function drummond2F0(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
    T = promote_type(T1, T2, T3)
    absα = abs(T(α))
    absβ = abs(T(β))
    if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ)
        return one(T)
    end
    ζ = inv(z)
    Nlo = ζ/(α*β)
    Dlo = ζ/(α*β)
    Tlo = Nlo/Dlo
    Nmid = (2ζ-(α+1)*(β+1))*Nlo + 2ζ
    Dmid = (2ζ-(α+1)*(β+1))*Dlo
    Tmid = Nmid/Dmid
    if norm((α+1)*(β+1)) < eps((absα+1)*(absβ+1))
        return Tmid
    end
    Nmid /= (α+1)*(β+1)
    Dmid /= (α+1)*(β+1)
    Nhi = (3ζ-(α+2)*(β+2)-(α+β+3))*Nmid - (α+β+3-ζ)*Nlo + ζ
    Dhi = (3ζ-(α+2)*(β+2)-(α+β+3))*Dmid - (α+β+3-ζ)*Dlo
    Thi = Nhi/Dhi
    if norm((α+2)*(β+2)) < eps((absα+2)*(absβ+2))
        return Thi
    end
    Nhi /= (α+2)*(β+2)
    Dhi /= (α+2)*(β+2)
    k = 2
    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
        Nhi, Nmid, Nlo = ((k+2)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Nhi - k*(α+β+3k-ζ)*Nmid - k*(k-1)*Nlo, Nhi, Nmid
        Dhi, Dmid, Dlo = ((k+2)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Dhi - k*(α+β+3k-ζ)*Dmid - k*(k-1)*Dlo, Dhi, Dmid
        Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
        if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1))
            return Thi
        end
        Nhi /= (α+k+1)*(β+k+1)
        Dhi /= (α+k+1)*(β+k+1)
        k += 1
    end
    return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo
end

# ₁F₁(α,β;z)
function drummond1F1(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
    T = promote_type(T1, T2, T3)
    absα = abs(T(α))
    if norm(z) < eps(real(T)) || norm(α) < eps(absα)
        return one(T)
    end
    ζ = inv(z)
    Nlo = β*ζ/α
    Dlo = β*ζ/α
    Tlo = Nlo/Dlo
    Nmid = ((β+1)*(2)*ζ - (α+1))*Nlo + (β+1)*(2)*ζ
    Dmid = ((β+1)*(2)*ζ - (α+1))*Dlo
    Tmid = Nmid/Dmid
    if norm(α+1) < eps(absα+1)
        return Tmid
    end
    Nmid /= α+1
    Dmid /= α+1
    Nhi = ((β+2)*(3)*ζ - (α+3))*Nmid + ((β+4)*ζ-1)*Nlo + (β+4)*ζ
    Dhi = ((β+2)*(3)*ζ - (α+3))*Dmid + ((β+4)*ζ-1)*Dlo
    Thi = Nhi/Dhi
    if norm(α+2) < eps(absα+2)
        return Thi
    end
    Nhi /= α+2
    Dhi /= α+2
    k = 2
    Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Nhi + k*((β+2k+2)*ζ-1)*Nmid + k*(k-1)*ζ*Nlo + 2ζ, Nhi, Nmid
    Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Dhi + k*((β+2k+2)*ζ-1)*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
    Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
    if norm(α+k+1) < eps(absα+k+1)
        return Thi
    end
    Nhi /= α+k+1
    Dhi /= α+k+1
    k += 1
    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
        Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Nhi + k*((β+2k+2)*ζ-1)*Nmid + k*(k-1)*ζ*Nlo, Nhi, Nmid
        Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Dhi + k*((β+2k+2)*ζ-1)*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
        Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
        if norm(α+k+1) < eps(absα+k+1)
            return Thi
        end
        Nhi /= α+k+1
        Dhi /= α+k+1
        k += 1
    end
    return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo
end

# ₀F₂(α,β;z)
function drummond0F2(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
    T = promote_type(T1, T2, T3)
    if norm(z) < eps(real(T)) || norm(α) < eps(real(T)) || norm(β) < eps(real(T))
        return one(T)
    end
    ζ = inv(z)
    Nlo = α*β*ζ
    Dlo = α*β*ζ
    Tlo = Nlo/Dlo
    k = 0
    Nmid2 = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nlo + (2*(α+1)*(β+1))*ζ
    Dmid2 = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dlo
    Tmid2 = Nmid2/Dmid2
    k = 1
    Nmid1 = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nmid2 + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Nlo + (3*(α+β+3)+(α+1)*(β+1))*ζ
    Dmid1 = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dmid2 + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dlo
    Tmid1 = Nmid1/Dmid1
    k = 2
    Nhi = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nmid1 + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Nmid2 + ζ*k*(k-1)*(3k+α+β+1)*Nlo + (14+2*(α+β))*ζ
    Dhi = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dmid1 + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dmid2 + ζ*k*(k-1)*(3k+α+β+1)*Dlo
    Thi = Nhi/Dhi
    k = 3
    Nhi, Nmid1, Nmid2, Nlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Nmid1 + ζ*k*(k-1)*(3k+α+β+1)*Nmid2 + ζ*k*(k-1)*(k-2)*Nlo + 6ζ, Nhi, Nmid1, Nmid2
    Dhi, Dmid1, Dmid2, Dlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dmid1 + ζ*k*(k-1)*(3k+α+β+1)*Dmid2 + ζ*k*(k-1)*(k-2)*Dlo, Dhi, Dmid1, Dmid2
    Thi, Tmid1, Tmid2, Tlo = Nhi/Dhi, Thi, Tmid1, Tmid2
    k += 1
    while k < kmax && errcheck(Tmid1, Thi, 10eps(real(T)))
        Nhi, Nmid1, Nmid2, Nlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Nmid1 + ζ*k*(k-1)*(3k+α+β+1)*Nmid2 + ζ*k*(k-1)*(k-2)*Nlo, Nhi, Nmid1, Nmid2
        Dhi, Dmid1, Dmid2, Dlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dmid1 + ζ*k*(k-1)*(3k+α+β+1)*Dmid2 + ζ*k*(k-1)*(k-2)*Dlo, Dhi, Dmid1, Dmid2
        Thi, Tmid1, Tmid2, Tlo = Nhi/Dhi, Thi, Tmid1, Tmid2
        k += 1
    end
    return isfinite(Thi) ? Thi : isfinite(Tmid1) ? Tmid1 : isfinite(Tmid2) ? Tmid2 : Tlo
end

# ₂F₁(α,β,γ;z)
function drummond2F1(α::T1, β::T2, γ::T3, z::T4; kmax::Int = 10_000) where {T1, T2, T3, T4}
    T = promote_type(T1, T2, T3, T4)
    absα = abs(T(α))
    absβ = abs(T(β))
    if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ)
        return one(T)
    end
    ζ = inv(z)
    Nlo = γ*ζ/(α*β)
    Dlo = γ*ζ/(α*β)
    Tlo = Nlo/Dlo
    Nmid = ((2)*(γ+1)*ζ-(α+1)*(β+1))*Nlo + 2*(γ+1)*ζ
    Dmid = ((2)*(γ+1)*ζ-(α+1)*(β+1))*Dlo
    Tmid = Nmid/Dmid
    if norm((α+1)*(β+1)) < eps((absα+1)*(absβ+1))
        return Tmid
    end
    Nmid /= (α+1)*(β+1)
    Dmid /= (α+1)*(β+1)
    Nhi = ((3)*(γ+2)*ζ-(α+2)*(β+2)-(α+β+3))*Nmid + ((γ+4)*ζ-(α+β+3))*Nlo + (γ+4)*ζ
    Dhi = ((3)*(γ+2)*ζ-(α+2)*(β+2)-(α+β+3))*Dmid + ((γ+4)*ζ-(α+β+3))*Dlo
    Thi = Nhi/Dhi
    if norm((α+2)*(β+2)) < eps((absα+2)*(absβ+2))
        return Thi
    end
    Nhi /= (α+2)*(β+2)
    Dhi /= (α+2)*(β+2)
    k = 2
    Nhi, Nmid, Nlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Nhi + k*((γ+2k+2)*ζ-(α+β+3k))*Nmid + k*(k-1)*(ζ-1)*Nlo + 2ζ, Nhi, Nmid
    Dhi, Dmid, Dlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Dhi + k*((γ+2k+2)*ζ-(α+β+3k))*Dmid + k*(k-1)*(ζ-1)*Dlo, Dhi, Dmid
    Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
    if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1))
        return Thi
    end
    Nhi /= (α+k+1)*(β+k+1)
    Dhi /= (α+k+1)*(β+k+1)
    k += 1
    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
        Nhi, Nmid, Nlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Nhi + k*((γ+2k+2)*ζ-(α+β+3k))*Nmid + k*(k-1)*(ζ-1)*Nlo, Nhi, Nmid
        Dhi, Dmid, Dlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Dhi + k*((γ+2k+2)*ζ-(α+β+3k))*Dmid + k*(k-1)*(ζ-1)*Dlo, Dhi, Dmid
        Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
        if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1))
            return Thi
        end
        Nhi /= (α+k+1)*(β+k+1)
        Dhi /= (α+k+1)*(β+k+1)
        k += 1
    end
    return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo
end

# ₘFₙ(α;β;z)
function pFqdrummond(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
    T = promote_type(T1, T2, T3)
    absα = abs.(T.(α))
    if norm(z) < eps(real(T)) || norm(prod(α)) < eps(prod(absα))
        return one(T)
    end
    ζ = inv(z)
    p = length(α)
    q = length(β)
    r = max(p+1, q+2)
    C = zeros(T, r)
    C[1] = one(T)
    Ĉ = zeros(T, r)
    Ĉ[1] = one(T)
    P = zeros(T, p+1)
    t = one(T)
    for j in 1:p
        t *= α[j]+1
    end
    P[1] = t
    err = one(real(T))
    for j in 1:p
        err *= absα[j]+1
    end
    Q = zeros(T, q+2)
    t = one(T)
    for j in 1:q
        t *= β[j]+1
    end
    Q[1] = 2t
    N = zeros(T, r+1)
    D = zeros(T, r+1)
    R = zeros(T, r+1)
    N[r+1] = prod(β)*ζ/prod(α)
    D[r+1] = prod(β)*ζ/prod(α)
    R[r+1] = N[r+1]/D[r+1]
    k = 0
    while k < r || (k < kmax && errcheck(R[r], R[r+1], 10eps(real(T))))
        for j in 1:r
            N[j] = N[j+1]
            D[j] = D[j+1]
            R[j] = R[j+1]
        end
        t1 = zero(T)
        for j in 0:min(k, q+1)
            t1 += Ĉ[j+1]*Q[j+1]*N[r-j]
        end
        if k ≤ q+1
            if p > q
                t1 += Q[k+1]
            else
                t1 += Q[k+1] / T(factorial(k+1))
            end
        end
        t2 = zero(T)
        t2 += Ĉ[1]*P[1]*N[r]
        for j in 1:min(k, p)
            t2 += P[j+1]*(C[j+1]*N[r-j+1] + Ĉ[j+1]*N[r-j])
        end
        N[r+1] = ζ*t1-t2
        t1 = zero(T)
        for j in 0:min(k, q+1)
            t1 += Ĉ[j+1]*Q[j+1]*D[r-j]
        end
        t2 = zero(T)
        t2 += Ĉ[1]*P[1]*D[r]
        for j in 1:min(k, p)
            t2 += P[j+1]*(C[j+1]*D[r-j+1] + Ĉ[j+1]*D[r-j])
        end
        D[r+1] = ζ*t1-t2
        R[r+1] = N[r+1]/D[r+1]
        if norm(P[1]) < eps(err)
            return R[r+1]
        end
        N[r+1] /= P[1]
        D[r+1] /= P[1]
        k += 1
        if p > q
            for j in min(k, max(p, q+1)):-1:1
                C[j+1] += C[j]
                Ĉ[j+1] += Ĉ[j]
            end
        else
            for j in min(k, max(p, q+1)):-1:1
                C[j+1] = (C[j+1]*(k+1-j) + C[j])/(k+1)
                Ĉ[j+1] = (Ĉ[j+1]*(k-j) + Ĉ[j])/(k+1)
            end
            Ĉ[1] = Ĉ[1]*k/(k+1)
        end
        t = one(T)
        for j in 1:p
            t *= α[j]+k+1
        end
        err = one(real(T))
        for j in 1:p
            err *= absα[j]+k+1
        end
        for j in 2:p+1
            s = t - P[j-1]
            P[j-1] = t
            t = s
        end
        P[p+1] = t
        t = one(T)
        for j in 1:q
            t *= β[j]+k+1
        end
        t *= k+2
        for j in 2:q+2
            s = t - Q[j-1]
            Q[j-1] = t
            t = s
        end
        Q[q+2] = t
    end
    return isfinite(R[r+1]) ? R[r+1] : R[r]
end


text:
using Oceananigans

# TODO: Use `BuoyancyField`!
function mixed_layer_depth(model)
    T = model.tracers.T
    ∂T̄∂z = AveragedField(Oceananigans.∂z(T), dims=(1, 2))
    compute!(∂T̄∂z)

    _, k_boundary_layer = findmax(abs.(interior(∂T̄∂z)))

    if k_boundary_layer isa CartesianIndex
        k_boundary_layer = k_boundary_layer.I[3]
    end

    mixed_layer_depth = - model.grid.zF[k_boundary_layer]

    return mixed_layer_depth
end


text:
using Libdl, Base.Sys

depsfile = joinpath(dirname(@__FILE__), "deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(knpath, libpath)
    open(depsfile,"w") do f
        print(f,"const libknitro = ")
        show(f, libpath)
        println(f)
        print(f,"const amplexe = ")
        show(f, joinpath(knpath, "..", "knitroampl", "knitroampl"))
        println(f)
    end
end

libname = string(Sys.iswindows() ? "" : "lib", "knitro", ".", Libdl.dlext)

if haskey(ENV, "LD_LIBRARY_PATH")
    paths_to_try = split(ENV["LD_LIBRARY_PATH"], ':')
else
    paths_to_try = String[]
end

if haskey(ENV, "KNITRODIR")
    push!(paths_to_try, joinpath(ENV["KNITRODIR"], "lib"))
end

global found_knitro = false
# test KNITRODIR first
for path in reverse(paths_to_try)
    l = joinpath(path, libname)
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found_knitro = true
        write_depsfile(path, l)
        break
    end
end

if !found_knitro
    error("Unable to locate KNITRO installation, " *
          "please check your environment variable KNITRODIR.")
end


text:
# ---
# title: 1016. Binary String With Substrings Representing 1 To N
# id: problem1016
# author: Indigo
# date: 2021-01-30
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/binary-string-with-substrings-representing-1-to-n/description/>
# hidden: true
# ---
#
# Given a binary string `S` (a string consisting only of '0' and '1's) and a
# positive integer `N`, return true if and only if for every integer X from 1 to
# N, the binary representation of X is a substring of S.
#
#
#
# **Example 1:**
#
#
#
#     Input: S = "0110", N = 3
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: S = "0110", N = 4
#     Output: false
#
#
#
#
# **Note:**
#
#   1. `1 <= S.length <= 1000`
#   2. `1 <= N <= 10^9`
#
#
## @lc code=start
using LeetCode

function query_string(S::String, N::Int)
    m = (N >> 1) + 1
    for i in m:N
        s_bit = join(reverse!(digits(i; base=2)))
        !occursin(s_bit, S) && return false
    end
    return true
end
## @lc code=end


text:
using DotaBot
using HTTP
using ArgParse
using Sockets
using JSON
using Flux
using Flux.Data: DataLoader
using Flux: @epochs
using CUDAapi
using BSON: @save, @load
if has_cuda()
    import CuArrays
    using CUDAdrv
    CuArrays.allowscalar(false)
end

s = ArgParseSettings()
@add_arg_table! s begin
    "--breezyIP"
    help = "breezy server IP adress"
    default = "127.0.0.1"
    "--breezyPort"
    help = "breezy server port number"
    default = "8085"
    "--agentIP"
    help = "agent server IP adress"
    default = "127.0.0.1"
    "--agentPort"
    help = "agent server port number"
    default = "8086"
end

args = parse_args(ARGS, s)

println("loading model")
@load "models/lookahead_model.bson" cm

println("starting server")
@async HTTP.serve(parse(IPAddr, args["agentIP"]), parse(Int, args["agentPort"])) do request::HTTP.Request
   @show request
   try
       f_act(features) = get_model_action(features, cm)
       response = server_handler(request, args, f_act)
       @show response
       return response
   catch e
       println(e)
       return HTTP.Response(200, JSON.json(Dict("done"=>false)))
   end
end

println("starting game")
start_game(args)

while true
    sleep(1)
end


text:
@testset "*" begin
    A = Array{Any}(undef, 6)
    A[1] = rand(Float32,10,10)
    A[2] = rand(Float32,10)
    A[3] = rand(Float32)
    A[4] = Variable(A[1])
    A[5] = Variable(A[2])
    A[6] = Variable(A[3])
    init(sess)
    for i = 1:6
        for j = 1:6
            if i<=3 && j<=3 || (i in [2,5] && j in [1,4]) ||
                (i in [2,5] && j in [2,5])
                continue
            end
            ret1 = A[i]*A[j]
            ret1 = run(sess, ret1)
            ret2 = A[i>3 ? i-3 : i] * A[j>3 ? j-3 : j]
            @test ret1 ≈ ret2
        end
    end

    @test A[1].*A[1] ≈ run(sess, A[4].*A[4])
    @test A[1].*A[1] ≈ run(sess, A[4].*A[1])
    @test A[1].*A[1] ≈ run(sess, A[1].*A[4])
    @test A[2].*A[2] ≈ run(sess, A[2].*A[5])
    @test A[2].*A[2] ≈ run(sess, A[5].*A[5])
    @test A[2].*A[2] ≈ run(sess, A[5].*A[2])
end

@testset "reshape" begin
    a = constant([1 2 3;4 5 6.])
    @test run(sess, reshape(a, 6))≈[ 1.0
    4.0
    2.0
    5.0
    3.0
    6.0]

    @test run(sess, reshape(a, 3, 2))≈[1.0  5.0
                                        4.0  3.0
                                        2.0  6.0]

    b = constant([1;2;3;4;5;6])
    @test run(sess, reshape(b, 2, 3))≈[1  3  5
                                    2  4  6]
    @test run(sess, reshape(constant(2.),1 )) ≈[2.0]

    a = constant([0.136944  0.364238
    0.250264  0.202175
    0.712224  0.388285
    0.163095  0.797934
    0.740991  0.382738])
    @test run(sess, reshape(a, :, 1))≈reshape([0.136944  0.364238
    0.250264  0.202175
    0.712224  0.388285
    0.163095  0.797934
    0.740991  0.382738], :, 1)
    @test run(sess, reshape(a, 1,:))≈reshape([0.136944  0.364238
    0.250264  0.202175
    0.712224  0.388285
    0.163095  0.797934
    0.740991  0.382738], 1,:)

end

@testset "scatter_update" begin
    A = Variable(ones(3))
    A = scatter_add(A, 3, 2.)
    B = [1.;1.;3.]

    C = Variable(ones(9))
    D = scatter_add(C, [2], [1], 1., 3, 3)
    G = reshape(D, 3, 3)
    E = [1. 1. 1; 2. 1. 1.; 1. 1. 1.]
    init(sess)
    @test run(sess, A)≈B
    @test run(sess, G)≈E

    A = Variable(rand(10))
    B = Variable(rand(5))
    C = scatter_add(A, [1;4;5;6;7], 2B)
    D = gradients(sum(C), B)
    E = gradients(sum(C), A)
    init(sess)
    @test run(sess, D)≈2ones(5)
    @test run(sess, E)≈ones(10)

    A = Variable(ones(10,10))
    B = Variable(ones(3,4))
    C = scatter_add(A, [1;2;3],[3;4;5;6], B)
    init(sess)
    @test run(sess, C)≈[1.0  1.0  2.0  2.0  2.0  2.0  1.0  1.0  1.0  1.0
    1.0  1.0  2.0  2.0  2.0  2.0  1.0  1.0  1.0  1.0
    1.0  1.0  2.0  2.0  2.0  2.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
    1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0]

    a = rand(10)
    A = constant(a)
    A = scatter_add(A, 3:4, -A[3:4])
    a[3:4] .= 0.0
    @test run(sess, A) ≈ a
end

@testset "adjoint" begin
    a = 1.0
    @test run(sess, constant(a)')≈a'
    a = rand(10)
    @test run(sess, constant(a)')≈a'
    a = rand(10,10)
    @test run(sess, constant(a)')≈a'
    a = rand(2,3,4,5)
    @test run(sess, constant(a)')≈permutedims(a, [1,2,4,3])
end

@testset "scatter_update_pyobject" begin
    A = constant(ones(10))
    B = constant(ones(3))
    ind = constant([2;4;6], dtype=Int32)
    C = scatter_add(A, ind, B)
    @test run(sess, C)≈[1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0]
    @test run(sess, gradients(sum(C),B))≈ones(3)
    @test run(sess, gradients(sum(C),A))≈ones(10)

    C = scatter_sub(A, ind, B)
    @test run(sess, C)≈[1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0]
    @test run(sess, gradients(sum(C),B))≈-ones(3)
    @test run(sess, gradients(sum(C),A))≈ones(10)

    B = constant(3*ones(3))
    C = scatter_update(A, ind, B)
    @test run(sess, C)≈[1.0, 3.0, 1.0, 3.0, 1.0, 3.0, 1.0, 1.0, 1.0, 1.0]
    @test run(sess, gradients(sum(C),B))≈ones(3)
    @test run(sess, gradients(sum(C),A))≈[1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0]
end

@testset "Operators" begin
    A = rand(10,2)
    TA = Variable(A)
    a = Variable(1.0)
    b = Variable(0.5)

    maxA = argmax(A, dims=2)
    maxA = [maxA[i][2] for i = 1:10]
    maxTA = argmax(TA, dims=2)

    minA = argmin(A, dims=2)
    minA = [minA[i][2] for i = 1:10]
    minTA = argmin(TA, dims=2)

    B = rand(Float32, 10)
    TB = Variable(B)
    TBc = constant(B)

    gop = group_assign([a,b], [2.0,2.0])
    init(sess)
    @test maxA  ≈ run(sess, maxTA)
    @test minA  ≈ run(sess, minTA)
    @test get_dtype(cast(maxTA, Float64))==Float64
    @test B.^2 ≈ run(sess, TB^2)
    @test B.^2 ≈ run(sess, TBc^2)

    @test 1.0 ≈ run(sess, max(a, b))
    @test 0.5 ≈ run(sess, min(a, b))
    @test maximum(A) ≈ run(sess, maximum(TA))
    @test minimum(A) ≈ run(sess, minimum(TA))
    @test reshape(maximum(A, dims=1), 2) ≈ run(sess, maximum(TA, dims=1))
    @test reshape(minimum(A, dims=1),2) ≈ run(sess, minimum(TA, dims=1))
    @test reshape(maximum(A, dims=2),10) ≈ run(sess, maximum(TA, dims=2))
    @test reshape(minimum(A, dims=2),10) ≈ run(sess, minimum(TA, dims=2))

    run(sess, gop)
    @test run(sess, a)≈2.0 && run(sess,b) ≈ 2.0
end

@testset "Other Operators" begin
    A = Variable(rand(10,10))
    @test_nowarn group(assign(A, rand(10,10)), A+A, 2*A)
end

@testset "Concat and stack" begin
    A = Variable(rand(10))
    B = Variable(rand(10,3))
    @test size(stack([A, A]))==(2,10)
    @test size(concat([A, A], dims=0))==(20,)
    @test size([A;A])==(20,)
    @test size([A A])==(10,2)
    @test size([B;B])==(20,3)
    @test size([B B])==(10,6)
end

@testset "Vectorize" begin
    A = rand(10)
    B = reshape(A, 2, 5)
    tA = constant(A)
    tB = constant(B)
    rA = rvec(tA)
    cB = vec(tB)
    init(sess)
    a, b = run(sess, [rA, cB])
    @test a≈reshape(A, 1, 10)
    @test run(sess, rvec(constant(B)))≈reshape(B[:], 1, :)
    @test run(sess, cvec(constant(B)))≈B[:] && size(cvec(constant(B)))==(10,1)
    @test run(sess, vec(constant(B)))≈B[:]
    @test b≈A
end

@testset "Solve" begin
    A = rand(10,10)
    x = rand(10)
    y = A*x
    tA = constant(A)
    ty = constant(y)
    @test run(sess, tA\ty)≈x
    @test run(sess, A\ty)≈x
    @test run(sess, tA\y)≈x
end

@testset "diff" begin
    A = rand(10,10)
    a = rand(10)
    tfA = constant(A)
    tfa = constant(a)
    @test run(sess, diff(tfA, dims=1))≈diff(A, dims=1)
    @test run(sess, diff(tfA, dims=2))≈diff(A, dims=2)
    @test run(sess, diff(tfa))≈diff(a)
end

@testset "clip" begin
    a = constant(3.0)
    a = clip(a, 1.0, 2.0)
    @test run(sess,a)≈2.0
end

@testset "map" begin
    a = constant([1.0,2.0,3.0])
    b = map(x->x^2, a)
    @test run(sess, b)≈[1.0;4.0;9.0]
end

@testset "diag" begin
    C = rand(3,3)
    D = rand(3)
    @test run(sess, diag(constant(C))) ≈  diag(C)
    @test run(sess, diagm(constant(D))) ≈ diagm(0=>D)
end

@testset "dot" begin
    A = rand(10)
    B = rand(10)
    tA = constant(A); tB = constant(B)
    u = sum(A.*B)
    @test run(sess, dot(tA, tB))≈u
    @test run(sess, dot(A, tB))≈u
    @test run(sess, dot(tA, B))≈u
end

@testset "prod" begin
    A = rand(10)
    @test run(sess, prod(constant(A)))≈prod(A)
end

@testset "findall" begin
    a = constant([0.849268    0.888376  0.0928501  0.119263  0.614649
    0.768562    0.978984  0.0802774  0.764186  0.491608
    0.00877878  0.751463  0.99539    0.613257  0.070678])
    b = a>0.5
    bi = findall(b)
    @test run(sess, bi)==[1  1
    1  2
    1  5
    2  1
    2  2
    2  4
    3  2
    3  3
    3  4]
    ci = findall(b[1,:])
    @test run(sess, ci)==[1;2;5]
end

@testset "svd" begin
    A = rand(10,20)
    r = svd(constant(A))
    A2 = r.U*diagm(r.S)*r.Vt
    @test run(sess, A2)≈A
end

@testset "vector" begin
    m = zeros(20)
    m[2:11] = rand(10)
    v = vector(2:11, m[2:11], 20)
    @test run(sess, v)≈m
end

@testset "repeat" begin
    A = constant(2.0)
    @test run(sess, repeat(A, 2)) ≈ ones(2)*2.0
    @test run(sess, repeat(A, 1, 2)) ≈ ones(1, 2)*2.0
    a = rand(2)
    A = constant(a)
    @test run(sess, repeat(A, 2)) ≈ repeat(a, 2)
    @test run(sess, repeat(A, 3, 2)) ≈ repeat(a, 3, 2)
    a = rand(2, 4)
    A = constant(a)
    @test run(sess, repeat(A, 2)) ≈ repeat(a, 2)
    @test run(sess, repeat(A, 3, 2)) ≈ repeat(a, 3, 2)
end

@testset "pmap" begin
    x = constant(ones(10))
    y1 = pmap(x->2.0*x, x)
    y2 = pmap(x->x[1]+x[2], [x,x])
    y3 = pmap(1:10, x) do z
        i = z[1]
        xi = z[2]
        xi + cast(Float64, i)
    end
    @test run(sess, y1)≈[ 2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0]
    @test run(sess, y2)≈[ 2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0
                    2.0]
    @test run(sess, y3)≈[ 2.0
                    3.0
                    4.0
                    5.0
                    6.0
                    7.0
                    8.0
                    9.0
                    10.0
                    11.0]
end

@testset "reshape" begin
    o = constant(rand(100,10,10))
    A = reshape(o, (100,2,50))
    B = tf.reshape(o, (100,2,50))
    @test run(sess, A)≈run(sess, B)
end

@testset "batch mul" begin
    A = rand(100, 10, 10)
    B = rand(100, 10, 3)
    C = zeros(100, 10, 3)
    for i = 1:100
        C[i,:,:] = A[i,:,:] * B[i,:,:]
    end
    C_ = batch_matmul(A, B)
    @test run(sess, C_) ≈ C
end

@testset "sort" begin
    A = rand(10,5)
    B = constant(A)
    @test run(sess, sort(B, dims=1))≈sort(A, dims=1)
    @test run(sess, sort(B, rev=true, dims=2))≈sort(A, dims=2, rev=true)
end

@testset "set_shape" begin
    a = placeholder(Float64, shape=[nothing, 10])
    b = set_shape(a, 3, 10)
    @test_nowarn run(sess, b, a=>rand(3,10)) # OK
    @test_throws PyCall.PyError run(sess, b, a=>rand(5,10)) # Error
    @test_throws PyCall.PyError run(sess, b, a=>rand(10,3)) # Error
end

@testset "activation" begin
    x = 1.0
    x_ = constant(1.0)
    for fn in [sigmoid, leaky_relu, selu, elu, softsign, softplus, relu, relu6]
        @test run(sess, fn(x_))≈fn(x)
    end
end

@testset "trace" begin
    A = rand(10,10)
    @test tr(A) ≈ run(sess, tr(constant(A)))
end


text:
function RMP2{T}(refWfn::Fermi.HartreeFock.RHF.RHFWavefunction,alg::Fermi.MollerPlesset.Direct) where T <: AbstractFloat
    dmp2 = 0.0
    nocc = refWfn.nalpha
    rocc = 1:1:refWfn.nalpha
    epsa = refWfn.epsa
    Ca = refWfn.Ca
    mints = refWfn.mints
    basis = refWfn.basis
    nao = basis.nbf()
    nvir = nao - nocc
    rvir = 1:nvir
    eps = refWfn.epsa
    dmp2 = 0.0
    for i in rocc
        Cav = Ca[:,nocc+1:nao]
        Cao = Ca[:,1:nocc]
        D = zeros(nvir,nocc,nvir)
        for j in rocc
            for a in rvir
                for b in rvir
                    aa = a + nocc
                    bb = b + nocc
                    D[a,j,b] = 1/(eps[i] + eps[j] - eps[aa] - eps[bb])
                end
            end
        end
        temp = Direct.direct_ao_contract(i,Ca,mints,basis)
        @tensoropt begin
            temp2[a,λ,σ] := Cav[ν,a]*temp[ν,λ,σ]
        end
        temp = nothing
        @tensoropt begin
            temp3[a,j,σ] := Cao[λ,j]*temp2[a,λ,σ]
        end
        temp2 = nothing
        @tensoropt begin
            temp4[a,j,b] := Cav[σ,b]*temp3[a,j,σ]
        end
        temp3 = nothing
        out = temp4 .* (2*temp4 - permutedims(temp4,[3,2,1])) .* D
        dmp2 += reduce(+,out)
    end
    RMP2{T}(dmp2,Fermi.MemTensor{T}(zeros(T,0,0,0,0)),Fermi.MemTensor{T}(zeros(T,0,0,0,0)))
end


text:
module Scattering1D

using LinearAlgebra
using SparseArrays
using ..Common
import ..Common.EQUIVALENT_SOURCE_RELATIVE_CUTOFF
import ..Common.AbstractComplexBL
import ..LinearScatteringProblem
import ..NonLinearScatteringProblem
import ..HelmholtzLS
import ..HelmholtzNLS
import ..MaxwellLS
import ..MaxwellNLS
import ..EquivalentSource
import ..ScatteringSolution
import ..ψ_to_x

Common.update!(mls::LinearScatteringProblem{1}, aL::Number, aR::Number, args...) = update!(mls, [aL,aR], args...)

function Common.update!(mls::LinearScatteringProblem{1}, a::Vector, args...)
    if !isempty(args)
        mls.ω = args[1]
        mls.operator(mls.ω)
    end
    copyto!(mls.equivalent_source.a, a)
    fill!(mls.equivalent_source.field,0)
    mls.j[:] = mls.equivalent_source(mls.operator.D², mls.operator.ky, mls.operator.kz)
    return nothing
end

"""
    HelmholtzLS(::Simulation{1}, ω, aL, [aR=0 ;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
HelmholtzLS(sim::Simulation{1}, ω::Real, aL::Number, aR::Number=0; kwargs...) = HelmholtzLS(sim,ω,[aL,aR];kwargs...)

"""
    HelmholtzLS(::Simulation{1}, ω, [amplitudes::Vector ;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
function HelmholtzLS(sim::Simulation{1}, ω::Real, a::Vector; start::Real=-Inf, stop::Real=Inf)
    M = Helmholtz(sim)
    M(ω) # evaluate at ω
    es = EquivalentSource{1}(sim, ω, a; start=start, stop=stop)
    j = es(M.D²,0,0)
    sol = ScatteringSolution{1}(sim,ω,a)
    return LinearScatteringProblem{1,1,typeof(M),typeof(es),typeof(sol)}(M,es,j,Ref(false),Ref(false),sol)
end

"""
    MaxwellLS(::Simulation{1}, ω, aL, [aR=0 ;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
MaxwellLS(sim::Simulation{1}, ω::Real, aL::Number, aR::Number=0; kwargs...) = MaxwellLS(sim,ω,[aL,aR];kwargs...)

"""
    MaxwellLS(::Simulation{1}, ω, [amplitudes::Vector ;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
function MaxwellLS(sim::Simulation{1}, ω::Real, a::Vector; start::Real=-Inf, stop::Real=Inf, ky::Number=0, kz::Number=0)
    M = Maxwell(sim; ky=ky,kz=kz)
    M(ω) # evaluate at ω
    es = EquivalentSource{3}(sim, ω, a; start=start, stop=stop, ky=ky, kz=kz)
    j = es(M.D²,ky,kz)
    sol = ScatteringSolution{3}(sim,ω,a)
    return LinearScatteringProblem{1,3,typeof(M),typeof(es),typeof(sol)}(M,es,j,Ref(false),Ref(false),sol)
end

"""
    HelmholtzNLS(::Simulation{1}, ω, aL, [aR=0, lupack=$DEFAULT_LUPACK;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
HelmholtzNLS(sim::Simulation{1}, ω::Real, aL::Number, aR::Number=0, lupack::AbstractLUPACK=DEFAULT_LUPACK; kwargs...) = HelmholtzNLS(sim,ω,[aL,aR]; kwargs...)

"""
    HelmholtzNLS(::Simulation{1}, ω, [amplitudes::Vector, lupack=$DEFAULT_LUPACK;  start=-Inf, stop=Inf, ky=0, kz=0]) -> ls
"""
function HelmholtzNLS(sim::Simulation{1}, ω::Real, a::Vector, lupack::AbstractLUPACK=DEFAULT_LUPACK; start::Real=-Inf, stop::Real=Inf)
    ls = HelmholtzLS(sim, ω, a; start=start, stop=stop)
    ψ = ScalarField(sim,1)
    res = similar(ψ)
    x = ψ_to_x(ψ)
    sol = ScatteringSolution{1}(sim,ω,a)
    return NonLinearScatteringProblem{1,1,typeof(ls),typeof(ψ),typeof(lupack)}(ls,ψ,res,x,Ref(false),lupack)
end

"""
    MaxwellNLS(::Simulation{1}, ω, aL, [aR=0, lupack=$DEFAULT_LUPACK;  start=-Inf, stop=Inf, ky=0, kz=0]) -> nls
"""
MaxwellNLS(sim::Simulation{1},ω::Real,aL::Number,aR::Number=0, lupack::AbstractLUPACK=DEFAULT_LUPACK; kwargs...) = MaxwellNLS(sim,ω,[aL,aR]; kwargs...)

"""
    MaxwellNLS(::Simulation{1}, ω, [amplitudes::Vector, lupack=$DEFAULT_LUPACK;  start=-Inf, stop=Inf, ky=0, kz=0]) -> nls
"""
function MaxwellNLS(sim::Simulation{1}, ω::Real, a::Vector, lupack::AbstractLUPACK=DEFAULT_LUPACK; start::Real=-Inf, stop::Real=Inf, ky::Number=0, kz::Number=0)
    ls = MaxwellLS(sim, ω, a; start=start, stop=stop, ky=ky, kz=kz)
    ψ = ElectricField(sim,1)
    res = similar(ψ)
    x = ψ_to_x(ψ)
    sol = ScatteringSolution{3}(sim,ω,a)
    return NonLinearScatteringProblem{1,3,typeof(ls),typeof(ψ),typeof(lupack)}(ls,ψ,res,x,Ref(false),lupack)
end

################################################################################
# EQUIVALENT SOURCE BLOCK FOR DIM=1

# constructor
function EquivalentSource{M}(sim::Simulation{1}, ω::Number, a::Vector; start::Number=-Inf,stop::Number=Inf, ky::Number=0, kz::Number=0) where M
    field = VectorField{M}(sim,1)
    bc1, bc2 = sim.boundary.bcs
    bl1, bl2 = sim.boundary.bls
    if typeof(bl1)<:noBL{1}
        start_in = start
        start_out = start
    else
        start_in  = max(start, sim.boundary.bls[1].start + 5sim.dx + 3sim.dx)
        start_out = max(start, sim.boundary.bls[1].start + 5sim.dx - 3sim.dx)
    end
    if typeof(bl2)<:noBL{2}
        stop_in = stop
        stop_out = stop
    else
        stop_in  = min(stop, sim.boundary.bls[2].start - 5sim.dx - 3sim.dx)
        stop_out = min(stop, sim.boundary.bls[2].start - 5sim.dx + 3sim.dx)
    end
    incoming_mask = Interval(max(field.start-sim.dx/2,start_in),min(field.stop+sim.dx/2,stop_in))
    outgoing_mask = Interval(max(field.start-sim.dx/2,start_out),min(field.stop+sim.dx/2,stop_out))
    kx = 2asin(sqrt(ω^2 - ky^2 - kz^2)*sim.dx/2)/sim.dx
    channelflux = fill(hypot(kx,ky,kz),2)
    return EquivalentSource(incoming_mask,outgoing_mask,field,sim,ω,a,channelflux)
end

# generate source
function (es::EquivalentSource{1,MCOMPONENTS})(ky::Number=0,kz::Number=0) where MCOMPONENTS
    if MCOMPONENTS==1
        M = Helmholtz(es.sim)
    elseif MCOMPONENTS==3
        M = Maxwell(es.sim;ky=ky,kz=kz)
    end
    M(es.ω)
    return es(M.D²,ky,kz)
end

# generate source, given differential operator
@inline function (es::EquivalentSource{1,M,TW,TA})(D²,ky::Number,kz::Number) where {M,TA<:AbstractVecOrMat,TW<:Number}
    n = length(es.sim)
    ω, a = es.ω, es.a
    ω² = ω^2
    kx = 2asin(sqrt(ω² - ky^2 - kz^2)*es.sim.dx/2)/es.sim.dx

    bl1, bl2 = es.sim.boundary.bls
    bc1, bc2 = es.sim.boundary.bcs

    # left  = es.sim.x[1] - es.sim.dx/2
    # right = es.sim.x[end] + es.sim.dx/2
    if typeof(bl1)<:AbstractComplexBL{1} || typeof(bc1)<:MatchedBC{1}
        left  = es.field.start
    else
        left  = es.sim.x[1] - es.sim.dx/2
    end
    if typeof(bl2)<:AbstractComplexBL{2} || typeof(bc2)<:MatchedBC{2}
        right = es.field.stop
    else
        right = es.sim.x[end] + es.sim.dx/2
    end

    temp = Vector{ComplexF64}(undef,M*n)
    rows = Vector{Int}(undef,M*n)
    count = 0
    # @fastmath @inbounds
    for i ∈ eachindex(es.field)
        j = mod1(i,n)
        if (i-1)÷n < 2
            if M==3
                β = (i-1)÷n==0 ? (ky/ω) : 1.0
            elseif M==1
                β = 1
            end
            incoming = es.incoming_mask(es.sim.x[j])
            outgoing = es.outgoing_mask(es.sim.x[j])
            es.field[i] = 0
            if incoming
                if typeof(bl2)<:AbstractComplexBL{2} || typeof(bc2)<:MatchedBC{2}
                    nothing
                elseif typeof(bc2)<:Union{DirichletBC{2},NeumannBC{2}}
                    es.field[i] += a[1]*β*cis(kx*(es.field.pos[j].x-right.x))/sqrt(es.channelflux[1])
                    iszero(a[2]) || throw("cannot be incident from $(typeof(bc2))")
                elseif typeof(bc2)<:FloquetBC{2}
                    throw("Floquet{2} incompatible with scattering")
                end
                if typeof(bl1)<:AbstractComplexBL{1} || typeof(bc1)<:MatchedBC{1}
                    nothing
                elseif typeof(bc1)<:Union{DirichletBC{1},NeumannBC{1}}
                    es.field[i] += a[2]*β*cis(-kx*(es.field.pos[j].x-left.x))/sqrt(es.channelflux[2])
                    iszero(a[1]) || throw("cannot be incident from $(typeof(bc1))")
                elseif typeof(bc1)<:FloquetBC{1}
                    throw("Floquet{1} incompatible with scattering")
                end
            end
            if outgoing
                if typeof(bl2)<:AbstractComplexBL{2} || typeof(bc2)<:MatchedBC{2}
                    es.field[i] = a[1]*β*cis(kx*(es.field.pos[j].x-left.x))/sqrt(es.channelflux[1])
                elseif typeof(bc2)<:DirichletBC{2}
                    es.field[i] += -a[1]*β*cis(-kx*(es.field.pos[j].x-right.x))/sqrt(es.channelflux[1])
                elseif typeof(bc2)<:NeumannBC{2}
                    es.field[i] += a[1]*β*cos(-kx*(es.field.pos[j].x-right.x))/sqrt(es.channelflux[1])
                end
                if typeof(bl1)<:AbstractComplexBL{1} || typeof(bc1)<:MatchedBC{1}
                    es.field[i] += a[2]*β*cis(-kx*(es.field.pos[j].x-right.x))/sqrt(es.channelflux[2])
                elseif typeof(bc1)<:DirichletBC{1}
                    es.field[i] += -a[2]*β*cis(kx*(es.field.pos[j].x-left.x))/sqrt(es.channelflux[2])
                elseif typeof(bc1)<:NeumannBC{1}
                    es.field[i] += a[2]*β*cis(kx*(es.field.pos[j].x-left.x))/sqrt(es.channelflux[2])
                end
            end
            if incoming || outgoing
                count += 1
                temp[count] = es.field[i]
                rows[count] = i
            end
        end
    end

    if M==3
        for i ∈ diagind(D²) D²[i] -= ω² end
    elseif M==1
        for i ∈ diagind(D²) D²[i] += ω² end
    end
    j = D²*sparsevec(rows[1:count],temp[1:count],M*n)
    if M==3
        for i ∈ diagind(D²) D²[i] += ω² end
    elseif M==1
        for i ∈ diagind(D²) D²[i] -= ω² end
    end

    mx = maximum(abs,j)
    return droptol!(j,mx*EQUIVALENT_SOURCE_RELATIVE_CUTOFF)
end

end # module

using .Scattering1D


text:
module DynamicalSystemsToolkit

using ModelingToolkit
using RuntimeGeneratedFunctions
using DocStringExtensions
using LinearAlgebra

const MTK = ModelingToolkit

RuntimeGeneratedFunctions.init(@__MODULE__)

# Exports - CollocationSystem

export CollocationSystem, add_boundarycondition!, generate_initialconditions

# Exports - DynamicalSystem

export DynamicalSystem

# Abstract types

abstract type AbstractDynamicalSystem <: MTK.AbstractSystem end

# Subsystems

include("collocation.jl")

# Utilities

function _strip_toplevel(name::String)
    idx = findfirst(==('₊'), name)
    return idx === nothing ? name[1:0] : name[nextind(name, idx):end]
end
_strip_toplevel(name::Symbol) = Symbol(_strip_toplevel(String(name)))

# General DynamicalSystem

struct DynamicalSystem <: AbstractDynamicalSystem
    eqs::Vector{Equation}
    states::Vector
    ps::Vector
    name::Symbol
    systems::Vector
    defaults::Dict
end

function DynamicalSystem(
    eqs, states, ps; name=gensym(:DynamicalSystem), defaults=Dict(), systems=[]
)
    sysnames = nameof.(systems)
    if length(unique(sysnames)) != length(sysnames)
        throw(ArgumentError("System names must be unique."))
    end
    defaults = MTK.todict(defaults)
    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))
    return DynamicalSystem(eqs, value.(states), value.(ps), name, systems, defaults)
end

end


text:
# this example is translated and tweaked from the following link:
# https://github.com/IntelRealSense/librealsense/blob/master/examples/C/distance/rs-distance.c
using RealSense

include("rs-utils.jl")

STREAM        = RS2_STREAM_DEPTH  # rs2_stream is a types of data provided by RealSense device
FORMAT        = RS2_FORMAT_Z16    # rs2_format is identifies how binary data is encoded within a frame
WIDTH         = 640               # defines the number of columns for each frame
HEIGHT        = 480               # defines the number of lines for each frame
FPS           = 30                # defines the rate of frames per second
STREAM_INDEX  = 0                 # defines the stream index, used for multiple streams of the same type

err = Ref{Ptr{rs2_error}}(0)

# create context
ctx = rs2_create_context(RS2_API_VERSION, err)
check_error(err)

# get a list of all the connected devices
device_list = rs2_query_devices(ctx, err)
check_error(err)

dev_count = rs2_get_device_count(device_list, err)
check_error(err)
println("There are $dev_count connected RealSense devices.")

# get the first connected device
dev = rs2_create_device(device_list, 0, err)
print_device_info(dev)

# create pipeline
_pipeline = rs2_create_pipeline(ctx, err)
check_error(err)

# create a config instance
config = rs2_create_config(err)
check_error(err)

# request a specific configuration
rs2_config_enable_stream(config, STREAM, STREAM_INDEX, WIDTH, HEIGHT, FORMAT, FPS, err)
check_error(err)

# start the pipeline streaming
@static if is_apple()
    # we sleep 1s on MacOS, otherwise it will hit issue #1586:
    # https://github.com/IntelRealSense/librealsense/issues/1586
    sleep(1)
end
pipeline_profile = rs2_pipeline_start_with_config(_pipeline, config, err)
@assert err[] == C_NULL "The connected device doesn't support depth streaming!"

while true
    frames = rs2_pipeline_wait_for_frames(_pipeline, 5000, err)
    check_error(err)

    # returns the number of frames embedded within the composite frame
    num_of_frames = rs2_embedded_frames_count(frames, err)
    check_error(err)

    imgDepth = zeros(UInt16, WIDTH, HEIGHT)
    for i = 0:num_of_frames-1
        # the retunred object should be released with rs2_release_frame(...)
        frame = rs2_extract_frame(frames, i, err)
        check_error(err)

        # check if the given frame can be extended to depth frame interface
        # accept only depth frames and skip other frames
        0 == rs2_is_frame_extendable_to(frame, RS2_EXTENSION_DEPTH_FRAME, err) && continue

        # retrieve depth frame's dimensions
        width = rs2_get_frame_width(frame, err)
        check_error(err)
        height = rs2_get_frame_height(frame, err)
        check_error(err)

        dist_to_center = rs2_depth_frame_get_distance(frame, width/2, height/2, err)
        check_error(err)

        println("The camera is facing an object $dist_to_center meters away.")

        rs2_release_frame(frame)
    end
    rs2_release_frame(frames)
end

# stop the pipeline streaming
rs2_pipeline_stop(_pipeline, err)
check_error(err)

# release resources
rs2_delete_pipeline_profile(pipeline_profile)
rs2_delete_config(config)
rs2_delete_pipeline(_pipeline)
rs2_delete_device(dev)
rs2_delete_device_list(device_list)
rs2_delete_context(ctx)