text:
"""
  Sampling from a poisson distribution
"""
abstract type AbstractPoissonDistribution <: AbstractSampleDistribution end
struct PoissonSampleDistribution{T} <: AbstractPoissonDistribution
  lambda::T
  function PoissonSampleDistribution(lambda::T) where {T <: AbstractFloat}
    return new{T}(lambda)
  end
end

"""
  Sampling from Poisson Distribution

  Reference:
  Poisson Random Variate Generation, Appl. Statis. (1991),
      40, No. 1, pp 143 - 158.

  # Example:
  sample(PoissonSampleDistribution(1.0), (20,))

"""
function sample(distrib::PoissonSampleDistribution{T}, shape::Tuple{Vararg{Int64}}, version::Val{1}) where {T <: AbstractFloat}

  n = prod(shape)
  p = exp(-distrib.lambda)
  ret = zeros(T, shape)

  for i in 1:n
    s = T(1); x = T(0)
    while true
      u = rand(T, 1)[1]
      s *= u
      if s < p
        break
      end
      x += 1
    end
    ret[i] = x
  end

  return ret
end

function sample(distrib::PoissonSampleDistribution{T}, shape::Int64, version::Val{1}) where {T <: AbstractFloat}
  return sample(distrib, (shape,))
end

"""
  Function to sample from poisson from a lambda array
"""
function sample(::Type{<: AbstractPoissonDistribution}, lambda::Array{T}, version::Val{1}) where {T <: AbstractFloat}

  shape = size(lambda)
  n = prod(shape)
  ret = zeros(T, shape)

  for i in 1:n
    p = exp(-lambda[i])
    s = T(1); x = T(0)
    while true
      u = rand(T, 1)[1]
      s *= u
      if s < p
        break
      end
      x += 1
    end
    ret[i] = x
  end

  return ret
end

function sample(distrib::PoissonSampleDistribution{T}, shape::Tuple{Vararg{Int64}}, version::Val{2}) where {T <: AbstractFloat}

  n = prod(shape)
  ret = zeros(T, shape)

  for i in 1:n
    p = exp(-distrib.lambda)
    F = p; x = T(0)
    u = rand(T, 1)[1]
    while true
      if u < F
        break
      end
      x += 1
      p *= (distrib.lambda/x)
      F += p
    end
    ret[i] = x
  end

  return ret
end

function sample(distrib::PoissonSampleDistribution{T}, shape::Int64, version::Val{2}) where {T <: AbstractFloat}
  return sample(distrib, (shape,))
end

"""
  Function to sample from poisson from a lambda array
"""
function sample(::Type{<: AbstractPoissonDistribution}, lambda::Array{T}, version::Val{2}) where {T <: AbstractFloat}

  shape = size(lambda)
  n = prod(shape)
  ret = zeros(T, shape)

  for i in 1:n
    p = exp(-lambda[i])
    F = p; x = T(0)
    u = rand(T, 1)[1]
    while true
      if u < F
        break
      end
      x += 1
      p *= (lambda[i]/x)
      F += p
    end
    ret[i] = x
  end

  return ret
end

"""
# Example simulate poisson distribution

  X, eta = simulateData(Float64, 10, 1000);
  y = linkinv(LogLink(), eta);
  y = sample(AbstractPoissonDistribution, y, Val{2}());
"""

"""
  Function to simulate data

  # Example
  using Random: seed!
  seed!(0);
  X, y = simulateData(Float64, PoissonDistribution(), LogLink(), 10, 1000)
  X, y = simulateData(Float64, BinomialDistribution(), LogitLink(), 10, 1000)
  X, y = simulateData(Float64, GammaDistribution(), LogLink(), 10, 1000)
  X, y = simulateData(Float64, GaussianDistribution(), IdentityLink(), 10, 1000)
"""
function simulateData(::Type{T}, distrib::AbstractDistribution,
              link::AbstractLink, p::Int64, n::Int64) where {T <: AbstractFloat}

  X = Array{T, 2}(undef, (0, 0))
  eta = Array{T, 1}(undef, (0,))

  X, eta = simulateData(T, p, n)
  y = linkinv(link, eta)

  if (typeof(distrib) <: GammaDistribution)
    y .+= 10
  end

  if typeof(distrib) <: PoissonDistribution
    y = sample(AbstractPoissonDistribution, y, Val{2}())
  end

  if typeof(distrib) <: BinomialDistribution
    y = map((x, u) -> T(1)*(x > u), y, rand(T, n))
  end

  y = reshape(y, (size(y)[1], 1))

  return X, y
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
export lanczosTridiag

"""
T,V = lanczosTridiag(A,b,k;tol,doReorth)

Lanczos method for getting a factorization of

   A = Vk Tk Vk'

where A is a real symmetric n by n matrix, Tk is a tridiagonal k by k matrix and the columns of
the n by k matrix Vk are orthogonal.

Implementation follows:

Paige, C. C. (1972).
Computational variants of the Lanczos method for the eigenproblem.
IMA Journal of Applied Mathematics.

Required input:

  A       - function computing A*x, e.g., x -> A*x
  b       - right hand side vector
  k       - dimension of Krylov subspace

Optional input:

  tol      - stopping tolerance
  doReorth - (default=false) set to true to perform full reorthogonalization

Output:

 Tk    - sparse tridiagonal matrix
 Vk    - basis vectors
"""
function lanczosTridiag{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,1},k; kwargs...)
	x = zeros(promote_type(T1,T2),size(A,2)) # pre-allocate
	return lanczosTridiag(v -> At_mul_B!(1.0,A,v,0.0,x),b,k;kwargs...) # multiply with transpose of A for efficiency
end

lanczosTridiag(A,b,k;kwargs...) = lanczosTridiag(x -> A*x,b,k;kwargs...)

function lanczosTridiag(A::Function,b::Vector,k::Int;tol=1e-5,doReorth::Bool=false)

n   = length(b)
x   = zeros(n)

# pre-allocate space for tridiagonalization and basis
beta  = zeros(k)
alpha = zeros(k)
V     = zeros(n,k)

beta[1] = norm(b)
V[:,1]  = copy(b)/beta[1]
u       = A(V[:,1])

j = 1 # brings j to scope of function
for j=1:k-1
    alpha[j] = dot(V[:,j],u)
    u        = u - alpha[j]*V[:,j]
    if doReorth # full re-orthogonalization
        for i=1:j
            u -= V[:,i]*dot(V[:,i],u)
        end
    end
    gamma     = norm(u)
    V[:,j+1]  = u/gamma
    beta[j+1] = gamma
    if beta[j+1]<tol
        break
    end
    u = A(V[:,j+1]) - beta[j+1]*V[:,j]
end

T = spdiagm((beta[2:j],alpha[1:j],beta[2:j]),-1:1,j,j)
return T,V[:,1:j]
end

ext: jl


text:
include("pricing_specs.jl")

"""
Read a text file containing pricing problem data. Assumes file is correctly formatted.

# Input
- `filename`: Address of the file.

# Output
- `price_spec`: An instance of the PricingSpecs type.
"""
function read_pricing(filename::AbstractString)::PricingSpecs
	f = open(filename)
	for i::Int=1:15			#this is the description part of the text file, so we skip all these lines
		readline(f)
	end

	line = readline(f)
	var_num::Int = parse(Int, line)		#number of variables
	readline(f)						#an empty line

	readline(f)					#total number of constraints (no field to assign to, so we skip)
	readline(f)

	line = readline(f)
	pricing_constr_num::Int = parse(Int, line)		#number of pricing constraints
	readline(f)

	line = readline(f)
	pieces = split(line)					#split strings into pieces separated by space (default delimiter)
	@assert(length(pieces) == var_num)
	lb::Array{Int, 1} = Array{Int, 1}(  var_num)
	for i::Int=1:var_num
		lb[i] = parse(Int, pieces[i])		#lower bound vector
	end
	readline(f)

	line = readline(f)
	pieces = split(line)
	@assert(length(pieces) == var_num)
	ub::Array{Int, 1} = Array{Int}(  var_num)
	for i::Int=1:var_num
		ub[i] = parse(Int, pieces[i])		#upper bound vector
		#ub[i] = 2
	end
	readline(f)

	line = readline(f)
	pieces = split(line)
	@assert(length(pieces) == var_num)
	obj_c::Array{Float64, 1} = Array{Float64, 1}(  var_num)
	for i::Int=1:var_num
		obj_c[i] = parse(Float64, pieces[i])		#objective function coefficients
	end
	readline(f)

	line = readline(f)
	pieces = split(line)
	@assert(length(pieces) == var_num)
	obj_d::Array{Float64, 1} = Array{Float64, 1}(  var_num)
	for i::Int=1:var_num
		obj_d[i] = parse(Float64, pieces[i])		#objective function degree of variables
	end
	readline(f)

	constr_c::Array{Float64, 2} = Array{Float64, 2}(  pricing_constr_num, var_num)
	for i::Int=1:pricing_constr_num
		line = readline(f)
		pieces = split(line)
		@assert(length(pieces) == var_num)
		for j::Int=1:var_num
			constr_c[i,j] = parse(Float64, pieces[j])		#coefficient of variable j in constraint i
		end
	end
	readline(f)

	constr_d::Array{Float64, 2} = Array{Float64, 2}(  pricing_constr_num, var_num)
	for i::Int=1:pricing_constr_num
		line = readline(f)
		pieces = split(line)
		@assert(length(pieces) == var_num)
		for j::Int=1:var_num
			constr_d[i,j] = parse(Float64, pieces[j])		#degree of the exp power for variable j in constraint i
		end
	end
	readline(f)

	line = readline(f)
	pieces = split(line)
	@assert(length(pieces) == pricing_constr_num)
	constr_rhs::Array{Float64, 1} = Array{Float64, 1}(  pricing_constr_num)
	for i::Int=1:pricing_constr_num
		constr_rhs[i] = parse(Float64, pieces[i])		#rhs value of constraints
	end
	readline(f)

	price_spec = PricingSpecs(var_num, pricing_constr_num, lb, ub, obj_c, obj_d, constr_c, constr_d, constr_rhs)	#create an instance of PricingSpecs from the read data
	return price_spec

end

"""
Write an output text file

# Input
- `filename`: Address of the file.
- `time_limit`: Time limit.
- `iter_max_ecp`: Max limit on ecp iterations
- `max_cut_num`: Max cut number added at EACH iteration of ecp
- `lp_to_ip_iter`: Number of ecp iterations to solve LP before switching to IP (if applicable)
- `lp_to_ip_tol`: Constraint violation tolerance below which ecp switches to solving IP (if applicable)
- `obj_improvement_tol`: Objective improvement tolerance over the last few iterations
- `tol_num_ecp`: The number of last iterations over which the relative obj improvement is computed
- `obj_val`: Objective value of the relaxation
- `ecp_status`: Status of the ECP at termination
- `initial_lp_obj`: Starting LP objective
- `iter_ecp`: Number of ECP iteration executed
- `total_cuts`: Total number of added DD cuts
- `ecp_time`: Total ECP elapsed time
- `width`: Array of widths of DDs in the model
- `dd_time`: Time to construct DD (and form root node)
- `relative_obj_improvement`: relative obj improvement over the last few iterations at termination
"""
function write_pricing(filename::AbstractString, time_limit::Float64, iter_max_ecp::Int, max_cut_num::Int, lp_to_ip_iter::Int, lp_to_ip_tol::Float64, obj_val::Float64, obj_improvement_tol::Float64, tol_num_ecp::Int,
	ecp_status::Symbol, initial_lp_obj::Float64, iter_ecp::Int, total_cuts::Int, ecp_time::Float64, width::Array{Int, 1}, dd_time::Float64, relative_obj_improvement::Float64)

	open(filename, "w") do f
		write(f, "This file contains the result of ECP for pricing problem\n")
		write(f, "********************************************************\n\n")
		write(f, "Starting objective value:\t$initial_lp_obj\n\n")
		write(f, "Final objective value:\t$obj_val\n\n")
		write(f, "ECP status:\t$ecp_status\n\n")
		write(f, "DD construction time:\t$dd_time\n\n")
		write(f, "ECP time:\t$ecp_time\n\n")
		write(f, "Time limit:\t$time_limit\n\n")
		write(f, "Relative obj improvement tolerance:\t$obj_improvement_tol\n\n")
		write(f, "Number of last few iterations to compute relative obj improvement:\t$tol_num_ecp\n\n")
		write(f, "Relative obj improvement at termination:\t$relative_obj_improvement\n\n")
		write(f, "Number of ECP iterations:\t$iter_ecp\n\n")
		write(f, "Number of added DD cuts:\t$total_cuts\n\n")
		write(f, "ECP iteration limit:\t$iter_max_ecp\n\n")
		write(f, "Max cut number at each iteration:\t$max_cut_num\n\n")
		write(f, "DD widths:\n")
		for i=1:length(width)
			write(f, "$(width[i])\t")
		end
	end
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
println("Hello " * join(ARGS, ", "))
sleep(1)

max_stars_repo_path: Chapter01/07. Setting up Julia to use multiple cores/hello2.jl


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

struct Rational{T<:Integer} <: Real
    num::T
    den::T

    function Rational{T}(num::Integer, den::Integer) where T<:Integer
        num == den == zero(T) && throw(ArgumentError("invalid rational: zero($T)//zero($T)"))
        g = den < 0 ? -gcd(den, num) : gcd(den, num)
        new(div(num, g), div(den, g))
    end
end
Rational(n::T, d::T) where T<:Integer = Rational{T}(n,d)
Rational(n::Integer, d::Integer) = Rational(promote(n,d)...)
Rational(n::Integer) = Rational(n,one(n))

function divgcd(x::Integer,y::Integer)
    g = gcd(x,y)
    div(x,g), div(y,g)
end

"""
    //(num, den)

Divide two integers or rational numbers, giving a `Rational` result.

```jldoctest
julia> 3 // 5
3//5

julia> (3 // 5) // (2 // 1)
3//10
```
"""
//(n::Integer,  d::Integer ) = Rational(n,d)

function //(x::Rational, y::Integer )
    xn,yn = divgcd(x.num,y)
    xn//checked_mul(x.den,yn)
end
function //(x::Integer,  y::Rational)
    xn,yn = divgcd(x,y.num)
    checked_mul(xn,y.den)//yn
end
function //(x::Rational, y::Rational)
    xn,yn = divgcd(x.num,y.num)
    xd,yd = divgcd(x.den,y.den)
    checked_mul(xn,yd)//checked_mul(xd,yn)
end

//(x::Complex,  y::Real) = complex(real(x)//y,imag(x)//y)
//(x::Number, y::Complex) = x*y'//abs2(y)

//(X::AbstractArray, y::Number) = X .// y

function show(io::IO, x::Rational)
    show(io, numerator(x))
    print(io, "//")
    show(io, denominator(x))
end

function read{T<:Integer}(s::IO, ::Type{Rational{T}})
    r = read(s,T)
    i = read(s,T)
    r//i
end
function write(s::IO, z::Rational)
    write(s,numerator(z),denominator(z))
end

convert{T<:Integer}(::Type{Rational{T}}, x::Rational) = Rational{T}(convert(T,x.num),convert(T,x.den))
convert{T<:Integer}(::Type{Rational{T}}, x::Integer) = Rational{T}(convert(T,x), convert(T,1))

convert(::Type{Rational}, x::Rational) = x
convert(::Type{Rational}, x::Integer) = convert(Rational{typeof(x)},x)

convert(::Type{Bool}, x::Rational) = x==0 ? false : x==1 ? true : throw(InexactError()) # to resolve ambiguity
convert(::Type{Integer}, x::Rational) = (isinteger(x) ? convert(Integer, x.num) : throw(InexactError()))
convert{T<:Integer}(::Type{T}, x::Rational) = (isinteger(x) ? convert(T, x.num) : throw(InexactError()))

convert(::Type{AbstractFloat}, x::Rational) = float(x.num)/float(x.den)
function convert{T<:AbstractFloat,S}(::Type{T}, x::Rational{S})
    P = promote_type(T,S)
    convert(T, convert(P,x.num)/convert(P,x.den))
end

function convert{T<:Integer}(::Type{Rational{T}}, x::AbstractFloat)
    r = rationalize(T, x, tol=0)
    x == convert(typeof(x), r) || throw(InexactError())
    r
end
convert(::Type{Rational}, x::Float64) = convert(Rational{Int64}, x)
convert(::Type{Rational}, x::Float32) = convert(Rational{Int}, x)

big(z::Complex{<:Rational{<:Integer}}) = Complex{Rational{BigInt}}(z)

promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{S}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{Rational{S}}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:AbstractFloat}(::Type{Rational{T}}, ::Type{S}) = promote_type(T,S)

widen{T}(::Type{Rational{T}}) = Rational{widen(T)}

"""
    rationalize([T<:Integer=Int,] x; tol::Real=eps(x))

Approximate floating point number `x` as a `Rational` number with components
of the given integer type. The result will differ from `x` by no more than `tol`.
If `T` is not provided, it defaults to `Int`.

```jldoctest
julia> rationalize(5.6)
28//5

julia> a = rationalize(BigInt, 10.3)
103//10

julia> typeof(numerator(a))
BigInt
```
"""
function rationalize{T<:Integer}(::Type{T}, x::AbstractFloat, tol::Real)
    if tol < 0
        throw(ArgumentError("negative tolerance $tol"))
    end
    isnan(x) && return zero(T)//zero(T)
    isinf(x) && return (x < 0 ? -one(T) : one(T))//zero(T)

    p,  q  = (x < 0 ? -one(T) : one(T)), zero(T)
    pp, qq = zero(T), one(T)

    x = abs(x)
    a = trunc(x)
    r = x-a
    y = one(x)

    tolx = oftype(x, tol)
    nt, t, tt = tolx, zero(tolx), tolx
    ia = np = nq = zero(T)

    # compute the successive convergents of the continued fraction
    #  np // nq = (p*a + pp) // (q*a + qq)
    while r > nt
        try
            ia = convert(T,a)

            np = checked_add(checked_mul(ia,p),pp)
            nq = checked_add(checked_mul(ia,q),qq)
            p, pp = np, p
            q, qq = nq, q
        catch e
            isa(e,InexactError) || isa(e,OverflowError) || rethrow(e)
            return p // q
        end

        # naive approach of using
        #   x = 1/r; a = trunc(x); r = x - a
        # is inexact, so we store x as x/y
        x, y = y, r
        a, r = divrem(x,y)

        # maintain
        # x0 = (p + (-1)^i * r) / q
        t, tt = nt, t
        nt = a*t+tt
    end

    # find optimal semiconvergent
    # smallest a such that x-a*y < a*t+tt
    a = cld(x-tt,y+t)
    try
        ia = convert(T,a)
        np = checked_add(checked_mul(ia,p),pp)
        nq = checked_add(checked_mul(ia,q),qq)
        return np // nq
    catch e
        isa(e,InexactError) || isa(e,OverflowError) || rethrow(e)
        return p // q
    end
end
rationalize{T<:Integer}(::Type{T}, x::AbstractFloat; tol::Real=eps(x)) = rationalize(T, x, tol)::Rational{T}
rationalize(x::AbstractFloat; kvs...) = rationalize(Int, x; kvs...)

"""
    numerator(x)

Numerator of the rational representation of `x`.

```jldoctest
julia> numerator(2//3)
2

julia> numerator(4)
4
```
"""
numerator(x::Integer) = x
numerator(x::Rational) = x.num

"""
    denominator(x)

Denominator of the rational representation of `x`.

```jldoctest
julia> denominator(2//3)
3

julia> denominator(4)
1
```
"""
denominator(x::Integer) = one(x)
denominator(x::Rational) = x.den

sign(x::Rational) = oftype(x, sign(x.num))
signbit(x::Rational) = signbit(x.num)
copysign(x::Rational, y::Real) = copysign(x.num,y) // x.den
copysign(x::Rational, y::Rational) = copysign(x.num,y.num) // x.den

typemin{T<:Integer}(::Type{Rational{T}}) = -one(T)//zero(T)
typemax{T<:Integer}(::Type{Rational{T}}) = one(T)//zero(T)

isinteger(x::Rational) = x.den == 1

-(x::Rational) = (-x.num) // x.den
function -{T<:Signed}(x::Rational{T})
    x.num == typemin(T) && throw(OverflowError())
    (-x.num) // x.den
end
function -{T<:Unsigned}(x::Rational{T})
    x.num != zero(T) && throw(OverflowError())
    x
end

for (op,chop) in ((:+,:checked_add), (:-,:checked_sub),
                  (:rem,:rem), (:mod,:mod))
    @eval begin
        function ($op)(x::Rational, y::Rational)
            xd, yd = divgcd(x.den, y.den)
            Rational(($chop)(checked_mul(x.num,yd), checked_mul(y.num,xd)), checked_mul(x.den,yd))
        end
    end
end

function *(x::Rational, y::Rational)
    xn,yd = divgcd(x.num,y.den)
    xd,yn = divgcd(x.den,y.num)
    checked_mul(xn,yn) // checked_mul(xd,yd)
end
/(x::Rational, y::Rational) = x//y
/(x::Rational, y::Complex{<:Union{Integer,Rational}}) = x//y

fma(x::Rational, y::Rational, z::Rational) = x*y+z

==(x::Rational, y::Rational) = (x.den == y.den) & (x.num == y.num)
<( x::Rational, y::Rational) = x.den == y.den ? x.num < y.num :
                               widemul(x.num,y.den) < widemul(x.den,y.num)
<=(x::Rational, y::Rational) = x.den == y.den ? x.num <= y.num :
                               widemul(x.num,y.den) <= widemul(x.den,y.num)

==(x::Rational, y::Integer ) = (x.den == 1) & (x.num == y)
==(x::Integer , y::Rational) = y == x
<( x::Rational, y::Integer ) = x.num < widemul(x.den,y)
<( x::Integer , y::Rational) = widemul(x,y.den) < y.num
<=(x::Rational, y::Integer ) = x.num <= widemul(x.den,y)
<=(x::Integer , y::Rational) = widemul(x,y.den) <= y.num

function ==(x::AbstractFloat, q::Rational)
    if isfinite(x)
        (count_ones(q.den) == 1) & (x*q.den == q.num)
    else
        x == q.num/q.den
    end
end

==(q::Rational, x::AbstractFloat) = x == q

for rel in (:<,:<=,:cmp)
    for (Tx,Ty) in ((Rational,AbstractFloat), (AbstractFloat,Rational))
        @eval function ($rel)(x::$Tx, y::$Ty)
            if isnan(x) || isnan(y)
                $(rel == :cmp ? :(throw(DomainError())) : :(return false))
            end

            xn, xp, xd = decompose(x)
            yn, yp, yd = decompose(y)

            if xd < 0
                xn = -xn
                xd = -xd
            end
            if yd < 0
                yn = -yn
                yd = -yd
            end

            xc, yc = widemul(xn,yd), widemul(yn,xd)
            xs, ys = sign(xc), sign(yc)

            if xs != ys
                return ($rel)(xs,ys)
            elseif xs == 0
                # both are zero or ±Inf
                return ($rel)(xn,yn)
            end

            xb, yb = ndigits0z(xc,2) + xp, ndigits0z(yc,2) + yp

            if xb == yb
                xc, yc = promote(xc,yc)
                if xp > yp
                    xc = (xc<<(xp-yp))
                else
                    yc = (yc<<(yp-xp))
                end
                return ($rel)(xc,yc)
            else
                return xc > 0 ? ($rel)(xb,yb) : ($rel)(yb,xb)
            end
        end
    end
end

# needed to avoid ambiguity between ==(x::Real, z::Complex) and ==(x::Rational, y::Number)
==(z::Complex , x::Rational) = isreal(z) & (real(z) == x)
==(x::Rational, z::Complex ) = isreal(z) & (real(z) == x)

for op in (:div, :fld, :cld)
    @eval begin
        function ($op)(x::Rational, y::Integer )
            xn,yn = divgcd(x.num,y)
            ($op)(xn, checked_mul(x.den,yn))
        end
        function ($op)(x::Integer,  y::Rational)
            xn,yn = divgcd(x,y.num)
            ($op)(checked_mul(xn,y.den), yn)
        end
        function ($op)(x::Rational, y::Rational)
            xn,yn = divgcd(x.num,y.num)
            xd,yd = divgcd(x.den,y.den)
            ($op)(checked_mul(xn,yd), checked_mul(xd,yn))
        end
    end
end

trunc{T}(::Type{T}, x::Rational) = convert(T,div(x.num,x.den))
floor{T}(::Type{T}, x::Rational) = convert(T,fld(x.num,x.den))
ceil{ T}(::Type{T}, x::Rational) = convert(T,cld(x.num,x.den))

function round{T, Tr}(::Type{T}, x::Rational{Tr}, ::RoundingMode{:Nearest})
    if denominator(x) == zero(Tr) && T <: Integer
        throw(DivideError())
    elseif denominator(x) == zero(Tr)
        return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
    end
    q,r = divrem(numerator(x), denominator(x))
    s = q
    if abs(r) >= abs((denominator(x)-copysign(Tr(4), numerator(x))+one(Tr)+iseven(q))>>1 + copysign(Tr(2), numerator(x)))
        s += copysign(one(Tr),numerator(x))
    end
    convert(T, s)
end

round{T}(::Type{T}, x::Rational) = round(T, x, RoundNearest)

function round{T, Tr}(::Type{T}, x::Rational{Tr}, ::RoundingMode{:NearestTiesAway})
    if denominator(x) == zero(Tr) && T <: Integer
        throw(DivideError())
    elseif denominator(x) == zero(Tr)
        return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
    end
    q,r = divrem(numerator(x), denominator(x))
    s = q
    if abs(r) >= abs((denominator(x)-copysign(Tr(4), numerator(x))+one(Tr))>>1 + copysign(Tr(2), numerator(x)))
        s += copysign(one(Tr),numerator(x))
    end
    convert(T, s)
end

function round{T, Tr}(::Type{T}, x::Rational{Tr}, ::RoundingMode{:NearestTiesUp})
    if denominator(x) == zero(Tr) && T <: Integer
        throw(DivideError())
    elseif denominator(x) == zero(Tr)
        return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
    end
    q,r = divrem(numerator(x), denominator(x))
    s = q
    if abs(r) >= abs((denominator(x)-copysign(Tr(4), numerator(x))+one(Tr)+(numerator(x)<0))>>1 + copysign(Tr(2), numerator(x)))
        s += copysign(one(Tr),numerator(x))
    end
    convert(T, s)
end

function round{T}(::Type{T}, x::Rational{Bool})
    if denominator(x) == false && issubtype(T, Union{Integer, Bool})
        throw(DivideError())
    end
    convert(T, x)
end

round{T}(::Type{T}, x::Rational{Bool}, ::RoundingMode{:Nearest}) = round(T, x)
round{T}(::Type{T}, x::Rational{Bool}, ::RoundingMode{:NearestTiesAway}) = round(T, x)
round{T}(::Type{T}, x::Rational{Bool}, ::RoundingMode{:NearestTiesUp}) = round(T, x)
round{T}(::Type{T}, x::Rational{Bool}, ::RoundingMode) = round(T, x)

trunc{T}(x::Rational{T}) = Rational(trunc(T,x))
floor{T}(x::Rational{T}) = Rational(floor(T,x))
ceil{ T}(x::Rational{T}) = Rational(ceil(T,x))
round{T}(x::Rational{T}) = Rational(round(T,x))

function ^(x::Rational, n::Integer)
    n >= 0 ? power_by_squaring(x,n) : power_by_squaring(inv(x),-n)
end

^(x::Number, y::Rational) = x^(y.num/y.den)
^{T<:AbstractFloat}(x::T, y::Rational) = x^convert(T,y)
^{T<:AbstractFloat}(x::Complex{T}, y::Rational) = x^convert(T,y)

^(z::Complex{<:Rational}, n::Bool) = n ? z : one(z) # to resolve ambiguity
function ^(z::Complex{<:Rational}, n::Integer)
    n >= 0 ? power_by_squaring(z,n) : power_by_squaring(inv(z),-n)
end

iszero(x::Rational) = iszero(numerator(x))

function lerpi(j::Integer, d::Integer, a::Rational, b::Rational)
    ((d-j)*a)/d + (j*b)/d
end


text:
## Tuning Feature
x = -3:0.002:3
y = gaussianf.(x,σ=1)
hw = halfwidth(y,circ=false,x=x)

x = -5:0.002:1
y = gaussianf.(x,σ=1)
hw = halfwidth(y,circ=false,x=x)
hw = halfwidth(y,circ=true,x=x)

x = -1:0.002:5
y = gaussianf.(x,σ=1)
hw = halfwidth(y,circ=false,x=x)
hw = halfwidth(y,circ=true,x=x)

## Circular Tuning Regression

# save(joinpath(@__DIR__,"circdata.jld2"),"y",y,"x",x,"y1",y1,"x1",x1)
y,x,y1,x1=load(joinpath(@__DIR__,"circdata.jld2"),"y","x","y1","x1")

# Tuning Curve, von Mises and Generalized von Mises
plotcondresponse(y,DataFrame(Ori=rad2deg.(x)),u=0,projection=:polar)
plotcondresponse(y1,DataFrame(Ori=rad2deg.(x1)),u=1,projection=:polar)

plot([vmf,gvmf],0,2π,label=["vmf" "gvmf"],projection=:polar)

# GvM and vM Fitting
mfit = fitmodel(:gvm,x,y)
plot(x->mfit.fun(x,mfit.param),0,2π,proj=:polar)

mfit = fitmodel(:vmn2,x1,y1)
plot(x->mfit.fun(x,mfit.param),0,2π,proj=:polar)

## Tuning Properties
f=factorresponsefeature(rad2deg.(x),y,factor=:Ori)

plot(x->f.fit.gvm.fun(x,f.fit.gvm.param),0,2π,lw=2)
scatter!(mod.(x.+0.5π,2π),y)

plot(x->f.fit.vmn2.fun(x,f.fit.vmn2.param),0,2π,lw=2)
scatter!(x,y)

f1=factorresponsefeature(rad2deg.(x1),y1,factor=:Ori)

plot(x->f1.fit.gvm.fun(x,f1.fit.gvm.param),0,2π,lw=2)
scatter!(mod.(x1.+0.5π,2π),y1)

plot(x->f1.fit.vmn2.fun(x,f1.fit.vmn2.param),0,2π,lw=2)
scatter!(x1,y1)


text:
using JSON
Nproc = 10

rmprocs(workers())
addprocs(Nproc)

println("workers: ",length(workers()))

@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim

# local
#@everywhere path="/home/emma/github/fit_behaviour_julia/"

#cluster
@everywhere path="/home/emma/fit_behaviour_julia/"

@everywhere path_functions=path*"functions/"
@everywhere path_scripts=path*"scripts_emma/"

@everywhere include(path_functions*"parameter_functions.jl")
@everywhere include(path_functions*"functions_julia.jl")
@everywhere include(path_functions*"utils.jl")

##############################
Nframes=10
Tframe=0.2
Nstim= 10000

model_num="8"

# Model only c2 and sigma
args=[ ("c2",[1]),("c4",0),("st",0),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",[1])]
consts=[("c2",[2]),("c4",[1]),("st",[1,2,3,4]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",0)]

x = [0.6,0.4]
z = [0,1,2,2,0,0,0,0,0,0]
param=make_dict(args,x,consts,z)

ub=[10., 10.]
lb=[-10., -10]
reg=[1,1]

model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>0.0)

sigma_s = 0.2
stim = sigma_s*randn(Nstim,Nframes)

history_bias=zeros(Nstim)

tic()
PR=ProbRight2(stim,history_bias,Tframe,param,model_f)
toc()

choices = []
for prob in PR
  if rand() < prob
    append!(choices,1)
  else
    append!(choices,-1)
  end
end

save_data = Dict()
final_matrix = zeros(20,20)
println(size(final_matrix))
for c2_p in 1:1:20
    for sig in 1:1:20
        x=[c2_p/20,sig/20]
        LL =  ComputeLL2(stim,choices,history_bias,Tframe, args, x,consts,z,model_f)
        final_matrix[floor(Int8,c2_p),floor(Int8,sig)]= LL
    end
end
save_data["matrix"]=final_matrix

# Buscar el mínim i fer la hessiana en el mínim
x = [0.6,0.4]
f,grads,Hess = LL_function_grads(stim,choices,history_bias,Tframe, args, x,consts,z,model_f)
save_data["hess"] = Hess

## SAVE THE RESULTS ##
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(save_data)
# write the file with the stringdata variable information
open("LL_2par_matrix.json", "w") do f
        write(f, stringdata)
    end
println("done")

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using BenchmarkTools

import Algopt.FirstOrder

using Algopt.TestFunctions: Rosenbrock

ros = Rosenbrock(a = 1, b = 5)

mgd = FirstOrder.MaximumGradientDescent()
res = @btime FirstOrder.search($mgd, $ros.f, $ros.∇f, x_0) setup = (x_0 = [rand(-2:.1:2), rand(-2:.1:2)])
println("MaximumGradientDescent: ", res)

grd = FirstOrder.GradientDescent()
term = FirstOrder.Termination(max_iter = 100_000)
res = @btime FirstOrder.search($grd, $ros.f, $ros.∇f, x_0; term = $term) setup = (x_0 = [rand(-2:.1:2), rand(-2:.1:2)])
println("GradientDescent: ", res)

cgd = FirstOrder.ConjugateGradientDescent()
res = @btime FirstOrder.search($cgd, $ros.f, $ros.∇f, x_0) setup = (x_0 = [rand(-2:.1:2), rand(-2:.1:2)])
println("ConjugateGradientDescent: ", res)


text:
# Binary functions for Nullables ...

# imports
import Base: convert
import Base: +, -, *, /, \, ^, %, &, |, $, >>>, >>, <<, ==, !=, <, <=, >, >=, //,  div, fld, rem, divrem, mod, copysign, flipsign, hypot, ldexp, log, beta, lbeta, airy, airyx, besselj, besseljx, bessely, besselyx, hankelh1, hankelh1x, hankelh2, hankelh2x, besseli, besselix, besselk, besselkx

# Convert to Nullable
function as_nullable{T}(x::T)
  T <: Nullable ? x : Nullable{T}(x)
end

# Convert Any to Nullable
function convert{T}(x::Type{Nullable{T}}, y)
  Nullable{T}(convert(T, y))
end

function convert(x::Type{Nullable}, y)
  X = typeof(y)
  Nullable{X}(y)
end

# Convert Nullable to Nullable
function convert{T <: Nullable, U <: Nullable}(x::Type{T}, y::U)
  X = eltype(T)
  Nullable{X}(convert(X, eltype(y)))
end

# Convert Nullable to Any
function convert{T <: Nullable}(U::Type, y::T)
  convert(U, get(y))
end

# Define binary operators
arith = [:+, :-, :*, :/, :\, :^, :%, ://]
bit = [:&, :|, :$, :>>>, :>>, :<<]
update = [:+=, :-=, :/=, :\=, :%=, :^=, :|=, :>>>=, :>>=, :<<=]
comp = [:(==), :!=, :<, :<=, :>, :>=]
funs = [:div, :fld, :rem, :divrem, :mod, :copysign, :flipsign, :hypot, :ldexp, :log, :beta, :lbeta, :airy, :airyx, :besselj, :besseljx, :bessely, :besselyx, :hankelh1, :hankelh1x, :hankelh2, :hankelh2x, :besseli, :besselix, :besselk, :besselkx]

# Combine to one vector
syms = [arith; bit; update; comp; funs]

# Symbol expression union type
sym_expr = Union{Symbol, Expr}

# Nullable{T} & Nullable{T} family
function exnn{Q<:Union{Symbol, Expr}}(sym::Q)
  quote
    function $sym{T<:Nullable, U<:Nullable}(x::T, y::U)
      try
        ret = $sym(get(x), get(y))
        X = typeof(ret)
        Nullable{X}(ret)
      catch
        Nullable()
      end
    end
  end
end

# For the Nullable{T} & U family
function exnu{Q<:Union{Symbol, Expr}}(sym::Q)
  quote
    function $sym{T<:Nullable, U}(x::T, y::U)
       try
         ret = $sym(get(x), y)
         X = typeof(ret)
         Nullable{X}(ret)
       catch
          Nullable()
       end
    end
  end
end

# For the U & Nullable{T} family
function exun{Q<:Union{Symbol, Expr}}(sym::Q)
  quote
    function $sym{U, T<:Nullable}(x::U, y::T)
      try
        ret = $sym(x, get(y))
        X = typeof(ret)
        Nullable{X}(ret)
      catch
        Nullable()
      end
    end
  end
end

# Source the functions
for i in [exnn, exnu, exun]
  for j in syms
    eval(i(j))
  end
end

max_stars_repo_name: ActiveAnalytics/NullableJulia

max_issues_repo_name: ActiveAnalytics/NullableJulia

max_forks_repo_name: ActiveAnalytics/NullableJulia


text:
# this script does a Monte Carlo that applies IV to solve inconsistency of OLS estimator
# when there is measurement error of the regressors
using Plots

function GIVmoments(theta, data)
	data = [data lags(data,2)]
    data = data[3:end,:] # get rid of missings
	n = size(data,1)
	y = data[:,1]
	ylag = data[:,2]
	x = data[:,3]
	xlag = data[:,6]
	xlag2 = data[:,9]
	X = [ones(n,1) ylag x]
	e = y - X*theta
	Z = [ones(n,1) x xlag xlag2]
	m = e.*Z
end

function main()
# do the Monte Carlo
n = 100
sig = 1
reps = 1000
results = zeros(reps,3)
for rep = 1:reps
	x = randn(n) # an exogenous regressor
	e = randn(n) # the error term
	ystar = zeros(n)
	# generate the dep var
	for t = 2:n
	  ystar[t] = 0.0 + 0.9*ystar[t-1] + 1.0*x[t] + e[t]
	end
    # add measurement error
	y = ystar + sig*randn(n)
	# now do GMM, using the data with meas. error in both dep. var. and regressor
	ylag = lag(y,1)
	data = [y ylag x];
    data = data[2:end,:] # drop first obs, missing due to lag
	theta = [0, 0.9, 1]
	weight = 1
    # note to self: this is very slow, because it uses the general GMM method, instead of analytic GIV
    # also, the weight matrix is not optimal
    moments = theta-> GIVmoments(theta, data)
	b, junk, junk, junk, junk = gmm(moments, theta, weight)
	b = b - [0.0, 0.9, 1.0] # subtract true values, so mean should be approx. zero if consistent
    results[rep,:] = b
end

histogram(results[:,1],nbins = 30)
#savefig("givconstant.png")
histogram(results[:,2],nbins = 30)
#savefig("givrho.png")
dstats(results)
return
end
main()

max_stars_repo_name: wenddymacro/Econometrics

max_issues_repo_name: wenddymacro/Econometrics

max_forks_repo_name: wenddymacro/Econometrics


text:
# Tests for the ways various types interface with eachother

# Test conversions
@test frequency(1.0m)*m ≈ c0
@test frequency(1nm)*nm ≈ c0
@test frequency(wavelength(100THz)) ≈ 100THz
@test wavelength(frequency(100nm)) ≈ 100nm
@test wavelength(100nm) ≈ wavelength(100e-9m)

attr = OpticalAttr(1.0, "Attribute")
@test attr(1.0THz) == attr(500nm)
@test attr(100THz) == attr(200THz)
@test attr(frequency(200nm)) == attr(frequency(150THz))
@test NonlinearOptics.der(attr, 100nm) == 0.0
@test NonlinearOptics.der(attr, 100nm; order=0) == 1.0
@test NonlinearOptics.der(attr, 1.0nm; order=3) == 0.0
@test NonlinearOptics.der(attr, 100THz; order=3) == 0.0

frequencies = linspace(100THz, 200THz, 20)
unitful_attr = frequencies./THz.*m
unitful_attr = OpticalAttr(frequencies, unitful_attr, "Attribute [m]")
@test unitful_attr(150THz) ≈ 150m
@test unitful_attr(frequency(c0/150THz)) ≈ 150m
@test unitful_attr(frequency(1998.6163866nm)) ≈ 150m
@test unitful_attr(1998.6163866nm) ≈ 150m
@test NonlinearOptics.der(unitful_attr, 150THz, order=0) == unitful_attr(150THz)
@test NonlinearOptics.der(unitful_attr, 150THz, order=1) ≈ 1m/THz
@test isapprox(NonlinearOptics.der(unitful_attr, 150THz, order=2), 0.0m/THz^2; atol=5e-15m/THz^2)


text:
using Base.Test
import WaveletScattering: get_n_orientations, ReflectionGroup,
    RotationGroup, TrivialGroup

@test get_n_orientations(TrivialGroup()) == 1
@test get_n_orientations(ReflectionGroup()) == 2
@test get_n_orientations(RotationGroup(4)) == 4
@test get_n_orientations(RotationGroup(8)) == 8


text:
import MPI
include("01-hello-impl.jl")

function main()
    MPI.Init()

    do_hello()

    MPI.Finalize()
end

main()


text:
struct TestStruct{A,B}
    a::A
    b::B
    c
end

struct TestBlocked{T}
    x::T
end
@blockobjects TestBlocked

@testset "Testing Annotation Utilities" begin
    manager = CachedArrays.CacheManager(
        CachedArrays.AlignedAllocator(),
        CachedArrays.AlignedAllocator();
        localsize = 100_000_000,
        remotesize = 100_000_000,
        minallocation = 10,
    )

    x = CachedArray(rand(Float32, 10), manager)
    y = similar(x)
    xo = x.object
    yo = y.object

    # Test "findobjects"
    @test findobjects(x) === (xo,)
    @inferred findobjects(x)
    @test findobjects((1,x,2,y)) == (xo, yo)
    @inferred findobjects((1,x,2,y))

    v = (a = 1, c = (TestStruct(5, x, 1), "hello"), d = y)
    @test findobjects(v) == (xo, yo)
    @inferred findobjects(v)

    # For now - check that uninferred fields are skipped.
    # Probably need to revisit this at some point though ...
    v = TestStruct(5, x, y)
    @test findobjects(v) == (xo,)

    # Test "onobjects"
    function collector(x)
        objects = []
        onobjects(x) do o
            push!(objects, o)
        end
        return objects
    end

    @test collector(x) == [xo]
    @test collector([x,y]) == [xo, yo]
    @test collector(TestStruct(x, x, x)) == [xo, xo, xo]
    @test collector(TestStruct(x, x, 10)) == [xo, xo]
    @test collector(TestStruct(1, 2, TestStruct(5, x, y))) == [xo, yo]
    @test collector((a = 1, b = TestBlocked(x), c = TestStruct(y, 0, x))) == [yo, xo]
end


text:
# exported trigonometric functions

"""
    sin(x)

Compute the sine of `x`, where the output is in radians.
"""
function sin end

"""
    cos(x)

Compute the cosine of `x`, where the output is in radians.
"""
function cos end

@inline function sincos_kernel(x::Double{Float64})
    c8 =  2.72052416138529567917983e-15
    c7 = -7.64292594113954471900203e-13
    c6 =  1.60589370117277896211623e-10
    c5 = -2.5052106814843123359368e-08
    c4 =  2.75573192104428224777379e-06
    c3 = -0.000198412698412046454654947
    c2 =  0.00833333333333318056201922
    c1 = -0.166666666666666657414808
    return dadd(c1, x.hi * (@horner x.hi c2 c3 c4 c5 c6 c7 c8))
end

@inline function sincos_kernel(x::Double{Float32})
    c4 =  2.6083159809786593541503f-06
    c3 = -0.0001981069071916863322258f0
    c2 =  0.00833307858556509017944336f0
    c1 = -0.166666597127914428710938f0
    return dadd(c1, x.hi * (@horner x.hi c2 c3 c4))
end

function sin(d::T) where {T<:Float64}
    qh = trunc(d * (T(M_1_PI) / (1 << 24)))
    ql = round(d * T(M_1_PI) - qh * (1 << 24))

    s = dadd2(d, qh * (-PI_A(T) * (1 << 24)))
    s = dadd2(s, ql * (-PI_A(T)            ))
    s = dadd2(s, qh * (-PI_B(T) * (1 << 24)))
    s = dadd2(s, ql * (-PI_B(T)            ))
    s = dadd2(s, qh * (-PI_C(T) * (1 << 24)))
    s = dadd2(s, ql * (-PI_C(T)            ))
    s = dadd2(s, (qh * (1 << 24) + ql) * - PI_D(T))

    t = s
    s = dsqu(s)

    w = sincos_kernel(s)

    v = dmul(t, dadd(T(1.0), dmul(w, s)))
    u = T(v)

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (u = -u)
    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

function sin(d::T) where {T<:Float32}
    q = round(d * T(M_1_PI))

    s = dadd2(d, q * -PI_A(T))
    s = dadd2(s, q * -PI_B(T))
    s = dadd2(s, q * -PI_C(T))
    s = dadd2(s, q * -PI_D(T))

    t = s
    s = dsqu(s)

    w = sincos_kernel(s)

    v = dmul(t, dadd(T(1.0), dmul(w, s)))
    u = T(v)

    qi = unsafe_trunc(Int, q)
    qi & 1 != 0 && (u = -u)
    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

function cos(d::T) where {T<:Float64}
    d = abs(d)

    qh = trunc(d * (T(M_1_PI) / (1 << 23)) - T(0.5) * (T(M_1_PI) / (1 << 23)))
    ql = 2*round(d * T(M_1_PI) - T(0.5) - qh * (1 << 23)) + 1

    s = dadd2(d, qh * (-PI_A(T)* T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_A(T)* T(0.5)            ))
    s = dadd2(s, qh * (-PI_B(T)* T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_B(T)* T(0.5)            ))
    s = dadd2(s, qh * (-PI_C(T)* T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_C(T)* T(0.5)            ))
    s = dadd2(s, (qh * (1 << 24) + ql) * (-PI_D(T) * T(0.5)))

    t = s
    s = dsqu(s)

    w = sincos_kernel(s)

    v = dmul(t, dadd(T(1.0), dmul(w, s)))

    u = T(v)

    qli = unsafe_trunc(Int, ql)
    qli & 2 == 0 && (u = -u)
    !isinf(d) && (d > TRIG_MAX(T)) && (u = T(0.0))

    return u
end

function cos(d::T) where {T<:Float32}
    d = abs(d)

    q = 1 + 2*round(d * T(M_1_PI) - T(0.5))

    s = dadd2(d, q * -PI_A(T)* T(0.5))
    s = dadd2(s, q * -PI_B(T)* T(0.5))
    s = dadd2(s, q * -PI_C(T)* T(0.5))
    s = dadd2(s, q * -PI_D(T)* T(0.5))

    t = s
    s = dsqu(s)

    w = sincos_kernel(s)

    v = dmul(t, dadd(T(1.0), dmul(w, s)))

    u = T(v)

    qi = unsafe_trunc(Int, q)
    qi & 2 == 0 && (u = -u)
    !isinf(d) && (d > TRIG_MAX(T)) && (u = T(0.0))

    return u
end

"""
    sin_fast(x)

Compute the sine of `x`, where the output is in radians.
"""
function sin_fast end

"""
    cos_fast(x)

Compute the cosine of `x`, where the output is in radians.
"""
function cos_fast end

# Argument is first reduced to the domain 0 < s < π/4

# We return the correct sign using `q & 1 != 0` i.e. q is odd (this works for
# positive and negative q) and if this condition is true we flip the sign since
# we are now in the negative branch of sin(x). Recall that q is just the integer
# part of d/π and thus we can determine the correct sign using this information.

@inline function sincos_fast_kernel(x::Float64)
    c9 = -7.97255955009037868891952e-18
    c8 =  2.81009972710863200091251e-15
    c7 = -7.64712219118158833288484e-13
    c6 =  1.60590430605664501629054e-10
    c5 = -2.50521083763502045810755e-08
    c4 =  2.75573192239198747630416e-06
    c3 = -0.000198412698412696162806809
    c2 =  0.00833333333333332974823815
    c1 = -0.166666666666666657414808
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9
end
@inline function sincos_fast_kernel(x::Float32)
    c4 =  2.6083159809786593541503f-06
    c3 = -0.0001981069071916863322258f0
    c2 =  0.00833307858556509017944336f0
    c1 = -0.166666597127914428710938f0
    return @horner x c1 c2 c3 c4
end

function sin_fast(d::T) where {T<:Float64}
    t = d

    qh = trunc(d * (T(M_1_PI) / (1 << 24)))
    ql = round(d * T(M_1_PI) - qh * (1 << 24))

    d = muladd(qh , -PI_A(T) * (1 << 24) , d)
    d = muladd(ql , -PI_A(T)             , d)
    d = muladd(qh , -PI_B(T) * (1 << 24) , d)
    d = muladd(ql , -PI_B(T)             , d)
    d = muladd(qh , -PI_C(T) * (1 << 24) , d)
    d = muladd(ql , -PI_C(T)             , d)
    d = muladd(qh * (1 << 24) + ql, -PI_D(T), d)

    s = d * d

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (d = -d)

    u = sincos_fast_kernel(s)
    u = muladd(s, u * d, d)

    !isinf(t) && (isnegzero(t) || abs(t) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

function sin_fast(d::T) where {T<:Float32}
    t = d

    q = round(d * T(M_1_PI))

    d = muladd(q , -PI_A(T), d)
    d = muladd(q , -PI_B(T), d)
    d = muladd(q , -PI_C(T), d)
    d = muladd(q , -PI_D(T), d)

    s = d * d

    qli = unsafe_trunc(Int, q)
    qli & 1 != 0 && (d = -d)

    u = sincos_fast_kernel(s)

    u = muladd(s, u * d, d)

    !isinf(t) && (isnegzero(t) || abs(t) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

function cos_fast(d::T) where {T<:Float64}
    t = d

    qh = trunc(d * (T(M_1_PI) / (1 << 23)) - T(0.5) * (T(M_1_PI) / (1 << 23)))
    ql = 2*round(d * T(M_1_PI) - T(0.5) - qh * (1 << 23)) + 1

    d = muladd(qh , -PI_A(T) * T(0.5) * (1 << 24) , d)
    d = muladd(ql , -PI_A(T) * T(0.5)             , d)
    d = muladd(qh , -PI_B(T) * T(0.5) * (1 << 24) , d)
    d = muladd(ql , -PI_B(T) * T(0.5)             , d)
    d = muladd(qh , -PI_C(T) * T(0.5) * (1 << 24) , d)
    d = muladd(ql , -PI_C(T) * T(0.5)             , d)
    d = muladd(qh * (1 << 24) + ql, -PI_D(T) * T(0.5), d)

    s = d * d

    qli = unsafe_trunc(Int, ql)
    qli & 2 == 0 && (d = -d)

    u = sincos_fast_kernel(s)
    u = muladd(s, u * d, d)

    !isinf(t) && (abs(t) > TRIG_MAX(T)) && (u = T(0.0))

    return u
end

function cos_fast(d::T) where {T<:Float32}
    t = d

    q = 1 + 2*round(d * T(M_1_PI) - T(0.5))

    d = muladd(q, -PI_A(T) * T(0.5), d)
    d = muladd(q, -PI_B(T) * T(0.5), d)
    d = muladd(q, -PI_C(T) * T(0.5), d)
    d = muladd(q, -PI_D(T) * T(0.5), d)

    s = d * d

    qi = unsafe_trunc(Int, q)
    qi & 2 == 0 && (d = -d)

    u = sincos_fast_kernel(s)

    u = muladd(s, u * d, d)

    !isinf(t) && (abs(t) > TRIG_MAX(T)) && (u = T(0.0))

    return u
end

"""
    sincos(x)

Compute the sin and cosine of `x` simultaneously, where the output is in
radians, returning a tuple.
"""
function sincos end

"""
    sincos_fast(x)

Compute the sin and cosine of `x` simultaneously, where the output is in
radians, returning a tuple.
"""
function sincos_fast end

@inline function sincos_a_kernel(x::Float64)
    a6 =  1.58938307283228937328511e-10
    a5 = -2.50506943502539773349318e-08
    a4 =  2.75573131776846360512547e-06
    a3 = -0.000198412698278911770864914
    a2 =  0.0083333333333191845961746
    a1 = -0.166666666666666130709393
    return @horner x a1 a2 a3 a4 a5 a6
end

@inline function sincos_a_kernel(x::Float32)
    a3 = -0.000195169282960705459117889f0
    a2 =  0.00833215750753879547119141f0
    a1 = -0.166666537523269653320312f0
    return @horner x a1 a2 a3
end

@inline function sincos_b_kernel(x::Float64)
    b7 = -1.13615350239097429531523e-11
    b6 =  2.08757471207040055479366e-09
    b5 = -2.75573144028847567498567e-07
    b4 =  2.48015872890001867311915e-05
    b3 = -0.00138888888888714019282329
    b2 =  0.0416666666666665519592062
    b1 = -0.50
    return @horner x b1 b2 b3 b4 b5 b6 b7
end

@inline function sincos_b_kernel(x::Float32)
    b5 = -2.71811842367242206819355f-07
    b4 =  2.47990446951007470488548f-05
    b3 = -0.00138888787478208541870117f0
    b2 =  0.0416666641831398010253906f0
    b1 = -0.5f0
    return @horner x b1 b2 b3 b4 b5
end

function sincos_fast(d::T) where {T<:Float64}
    s = d

    qh = trunc(d * ((2 * T(M_1_PI)) / (1 << 24)))
    ql = round(d * (2 * T(M_1_PI)) - qh * (1 << 24))

    s = muladd(qh, -PI_A(T) * T(0.5) * (1 << 24), s)
    s = muladd(ql, -PI_A(T) * T(0.5),             s)
    s = muladd(qh, -PI_B(T) * T(0.5) * (1 << 24), s)
    s = muladd(ql, -PI_B(T) * T(0.5),             s)
    s = muladd(qh, -PI_C(T) * T(0.5) * (1 << 24), s)
    s = muladd(ql, -PI_C(T) * T(0.5),             s)
    s = muladd(qh * (1 << 24) + ql, -PI_D(T) * 0.5, s)

    t = s

    s = s * s

    u  = sincos_a_kernel(s)
    u  = u * s * t

    rx = t + u

    isnegzero(d) && (rx = T(-0.0))

    u = sincos_b_kernel(s)

    ry = u * s + T(1.0)

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (s = ry; ry = rx; rx = s)
    qli & 2 != 0 && (rx = -rx)
    (qli + 1) & 2 != 0 && (ry = -ry)

    abs(d) > TRIG_MAX(T) && (rx = ry = T(0.0))
    isinf(d) && (rx = ry = T(NaN))

    return (rx, ry)
end

function sincos_fast(d::T) where {T<:Float32}
    s = d

    q = round(d * (2 * T(M_1_PI)))

    s = muladd(q, -PI_A(T) * T(0.5), s)
    s = muladd(q, -PI_B(T) * T(0.5), s)
    s = muladd(q, -PI_C(T) * T(0.5), s)
    s = muladd(q, -PI_D(T) * T(0.5), s)

    t = s

    s = s * s

    u  = sincos_a_kernel(s)
    u  = u * s * t

    rx = t + u

    isnegzero(d) && (rx = T(-0.0))

    u = sincos_b_kernel(s)

    ry = u * s + T(1.0)

    qi = unsafe_trunc(Int, q)
    qi & 1 != 0 && (s = ry; ry = rx; rx = s)
    qi & 2 != 0 && (rx = -rx)
    (qi + 1) & 2 != 0 && (ry = -ry)

    abs(d) > TRIG_MAX(T) && (rx = ry = T(0.0))
    isinf(d) && (rx = ry = T(NaN))

    return (rx, ry)
end

function sincos(d::T) where {T<:Float64}
    qh = trunc(d * ((2 * T(M_1_PI)) / (1 << 24)))
    ql = round(d * (2 * T(M_1_PI)) - qh * (1 << 24))

    s = dadd2(d, qh * (-PI_A(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_A(T) * T(0.5)            ))
    s = dadd2(s, qh * (-PI_B(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_B(T) * T(0.5)            ))
    s = dadd2(s, qh * (-PI_C(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_C(T) * T(0.5)            ))
    s = dadd2(s, (qh * (1 << 24) + ql) * (-PI_D(T) * T(0.5)))

    t  = s
    s  = dsqu(s)
    sx = T(s)

    u  = sincos_a_kernel(sx)

    u *= sx * t.hi

    v  = dadd(t, u)
    rx = T(v)

    isnegzero(d) && (rx = T(-0.0))

    u  = sincos_b_kernel(sx)

    v  = dadd(T(1.0), dmul(sx, u))
    ry = T(v)

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (u = ry; ry = rx; rx = u)
    qli & 2 != 0 && (rx = -rx)
    (qli + 1) & 2 != 0 && (ry = -ry)

    abs(d) > TRIG_MAX(T) && (rx = ry = T(0.0))
    isinf(d) && (rx = ry = T(NaN))

    return (rx, ry)
end

function sincos(d::T) where {T<:Float32}
    q = round(d * (2 * T(M_1_PI)))

    s = dadd2(d, q * (-PI_A(T) * T(0.5)))
    s = dadd2(s, q * (-PI_B(T) * T(0.5)))
    s = dadd2(s, q * (-PI_C(T) * T(0.5)))
    s = dadd2(s, q * (-PI_D(T) * T(0.5)))

    t  = s
    s  = dsqu(s)
    sx = T(s)

    u  = sincos_a_kernel(sx)

    u *= sx * t.hi

    v  = dadd(t, u)
    rx = T(v)

    isnegzero(d) && (rx = T(-0.0))

    u  = sincos_b_kernel(sx)

    v  = dadd(T(1.0), dmul(sx, u))
    ry = T(v)

    qi = unsafe_trunc(Int, q)
    qi & 1 != 0 && (u = ry; ry = rx; rx = u)
    qi & 2 != 0 && (rx = -rx)
    (qi + 1) & 2 != 0 && (ry = -ry)

    abs(d) > TRIG_MAX(T) && (rx = ry = T(0.0))
    isinf(d) && (rx = ry = T(NaN))

    return (rx, ry)
end

"""
    tan(x)

Compute the tangent of `x`, where the output is in radians.
"""
function tan end

"""
    tan_fast(x)

Compute the tangent of `x`, where the output is in radians.
"""
function tan_fast end

@inline function tan_fast_kernel(x::Float64)
    c16 =  9.99583485362149960784268e-06
    c15 = -4.31184585467324750724175e-05
    c14 =  0.000103573238391744000389851
    c13 = -0.000137892809714281708733524
    c12 =  0.000157624358465342784274554
    c11 =  -6.07500301486087879295969e-05
    c10 =  0.000148898734751616411290179
    c9  =  0.000219040550724571513561967
    c8  =  0.000595799595197098359744547
    c7  =  0.00145461240472358871965441
    c6  =  0.0035923150771440177410343
    c5  =  0.00886321546662684547901456
    c4  =  0.0218694899718446938985394
    c3  =  0.0539682539049961967903002
    c2  =  0.133333333334818976423364
    c1  =  0.333333333333320047664472
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16
end

@inline function tan_fast_kernel(x::Float32)
    c7 =  0.00446636462584137916564941f0
    c6 = -8.3920182078145444393158f-05
    c5 =  0.0109639242291450500488281f0
    c4 =  0.0212360303848981857299805f0
    c3 =  0.0540687143802642822265625f0
    c2 =  0.133325666189193725585938f0
    c1 =  0.33333361148834228515625f0
    return @horner x c1 c2 c3 c4 c5 c6 c7
end

function tan_fast(d::T) where {T<:Float64}
    qh = trunc(d * (2 * T(M_1_PI)) / (1 << 24))
    ql = round(d * (2 * T(M_1_PI)) - qh * (1 << 24))

    x = muladd(qh, -PI_A(T) * T(0.5) * (1 << 24), d)
    x = muladd(ql, -PI_A(T) * T(0.5),             x)
    x = muladd(qh, -PI_B(T) * T(0.5) * (1 << 24), x)
    x = muladd(ql, -PI_B(T) * T(0.5),             x)
    x = muladd(qh, -PI_C(T) * T(0.5) * (1 << 24), x)
    x = muladd(ql, -PI_C(T) * T(0.5),             x)
    x = muladd(qh * (1 << 24) + ql, -PI_D(T) * T(0.5), x)

    s = x * x

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (x = -x)

    u = tan_fast_kernel(s)

    u = muladd(s, u * x, x)

    qli & 1 != 0 && (u = T(1.0) / u)

    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

function tan_fast(d::T) where {T<:Float32}
    q = round(d * (2 * T(M_1_PI)))

    x = d

    x = muladd(q, -PI_A(T) * T(0.5), x)
    x = muladd(q, -PI_B(T) * T(0.5), x)
    x = muladd(q, -PI_C(T) * T(0.5), x)
    x = muladd(q, -PI_D(T) * T(0.5), x)

    s = x * x

    qi = unsafe_trunc(Int, q)
    qi & 1 != 0 && (x = -x)

    u = tan_fast_kernel(s)

    u = muladd(s, u * x, x)

    qi & 1 != 0 && (u = T(1.0) / u)

    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (u = T(-0.0))

    return u
end

@inline function tan_kernel(x::Double{Float64})
    c15 = 1.01419718511083373224408e-05
    c14 = -2.59519791585924697698614e-05
    c13 = 5.23388081915899855325186e-05
    c12 = -3.05033014433946488225616e-05
    c11 = 7.14707504084242744267497e-05
    c10 = 8.09674518280159187045078e-05
    c9 = 0.000244884931879331847054404
    c8 = 0.000588505168743587154904506
    c7 = 0.00145612788922812427978848
    c6 = 0.00359208743836906619142924
    c5 = 0.00886323944362401618113356
    c4 = 0.0218694882853846389592078
    c3 = 0.0539682539781298417636002
    c2 = 0.133333333333125941821962
    c1 = 0.333333333333334980164153
    return dadd(c1, x.hi * (@horner x.hi c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15))
end

@inline function tan_kernel(x::Double{Float32})
    c7 =  0.00446636462584137916564941f0
    c6 = -8.3920182078145444393158f-05
    c5 =  0.0109639242291450500488281f0
    c4 =  0.0212360303848981857299805f0
    c3 =  0.0540687143802642822265625f0
    c2 =  0.133325666189193725585938f0
    c1 =  0.33333361148834228515625f0
    return dadd(c1,  x.hi * (@horner x.hi c2 c3 c4 c5 c6 c7))
end

function tan(d::T) where {T<:Float64}
    qh = trunc(d * (T(M_2_PI)) / (1 << 24))
    s = dadd2(dmul(Double(T(M_2_PI_H), T(M_2_PI_L)), d), (d < 0 ? T(-0.5) : T(0.5)) - qh * (1 << 24))
    ql = trunc(T(s))

    s = dadd2(d, qh * (-PI_A(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_A(T) * T(0.5)            ))
    s = dadd2(s, qh * (-PI_B(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_B(T) * T(0.5)            ))
    s = dadd2(s, qh * (-PI_C(T) * T(0.5) * (1 << 24)))
    s = dadd2(s, ql * (-PI_C(T) * T(0.5)            ))
    s = dadd2(s, (qh * (1 << 24) + ql) * (-PI_D(T) * T(0.5)))

    qli = unsafe_trunc(Int, ql)
    qli & 1 != 0 && (s = -s)

    t = s
    s = dsqu(s)

    u = tan_kernel(s)

    x = dadd(T(1.0), dmul(u, s))
    x = dmul(t, x)

    qli & 1 != 0 && (x = drec(x))

    v = T(x)

    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (v = T(-0.0))

    return v
end

function tan(d::T) where {T<:Float32}
    q = round(d * (T(M_2_PI)))

    s = dadd2(d, q * -PI_A(T) * T(0.5))
    s = dadd2(s, q * -PI_B(T) * T(0.5))
    s = dadd2(s, q * -PI_C(T) * T(0.5))
    s = dadd2(s, q * -PI_XD(T) * T(0.5))
    s = dadd2(s, q * -PI_XE(T) * T(0.5))

    qi = unsafe_trunc(Int, q)
    qi & 1 != 0 && (s = -s)

    t = s
    s = dsqu(s)
    s = dnormalize(s)

    u = tan_kernel(s)

    x = dadd(T(1.0), dmul(u, s))
    x = dmul(t, x)

    qi & 1 != 0 && (x = drec(x))

    v = T(x)

    !isinf(d) && (isnegzero(d) || abs(d) > TRIG_MAX(T)) && (v = T(-0.0))

    return v
end

"""
    atan(x)

Compute the inverse tangent of `x`, where the output is in radians.
"""
function atan(x::T) where {T<:Union{Float32,Float64}}
    u = T(atan2k(Double(abs(x)), Double(T(1))))
    isinf(x) && (u = T(PI_2))
    flipsign(u, x)
end

@inline function atan_fast_kernel(x::Float64)
    c19 = -1.88796008463073496563746e-05
    c18 =  0.000209850076645816976906797
    c17 = -0.00110611831486672482563471
    c16 =  0.00370026744188713119232403
    c15 = -0.00889896195887655491740809
    c14 =  0.016599329773529201970117
    c13 = -0.0254517624932312641616861
    c12 =  0.0337852580001353069993897
    c11 = -0.0407629191276836500001934
    c10 =  0.0466667150077840625632675
    c9  = -0.0523674852303482457616113
    c8  =  0.0587666392926673580854313
    c7  = -0.0666573579361080525984562
    c6  =  0.0769219538311769618355029
    c5  = -0.090908995008245008229153
    c4  =  0.111111105648261418443745
    c3  = -0.14285714266771329383765
    c2  =  0.199999999996591265594148
    c1  = -0.333333333333311110369124
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19
end

@inline function atan_fast_kernel(x::Float32)
    c8 =  0.00282363896258175373077393f0
    c7 = -0.0159569028764963150024414f0
    c6 =  0.0425049886107444763183594f0
    c5 = -0.0748900920152664184570312f0
    c4 =  0.106347933411598205566406f0
    c3 = -0.142027363181114196777344f0
    c2 =  0.199926957488059997558594f0
    c1 = -0.333331018686294555664062f0
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8
end

"""
    atan_fast(x)

Compute the inverse tangent of `x`, where the output is in radians.
"""
function atan_fast(x::T) where {T<:Union{Float32,Float64}}
    q = 0
    if signbit(x)
        x = -x
        q = 2
    end
    if x > 1
        x = 1 / x
        q |= 1
    end
    t = x * x
    u = atan_fast_kernel(t)
    t = x + x * t * u
    q & 1 != 0 && (t = T(PI_2) - t)
    q & 2 != 0 && (t = -t)
    return t
end

const under_atan2(::Type{Float64}) = 5.5626846462680083984e-309
const under_atan2(::Type{Float32}) = 2.9387372783541830947f-39

"""
    atan(x, y)

Compute the inverse tangent of `x/y`, using the signs of both `x` and `y` to determine the quadrant of the return value.
"""
function atan(x::T, y::T) where {T<:Union{Float32,Float64}}
    abs(y) < under_atan2(T) && (x *= T(Int64(1) << 53); y *= T(Int64(1) << 53))
    r = T(atan2k(Double(abs(x)), Double(y)))

    r = flipsign(r, y)
    if isinf(y) || y == 0
        r = T(PI_2) - (isinf(y) ? _sign(y) * T(PI_2) : T(0.0))
    end
    if isinf(x)
        r = T(PI_2) - (isinf(y) ? _sign(y) * T(PI_4) : T(0.0))
    end
    if x == 0
        r = _sign(y) == -1 ? T(M_PI) : T(0.0)
    end
    return isnan(y) || isnan(x) ? T(NaN) : flipsign(r, x)
end

"""
    atan2_fast(x, y)

Compute the inverse tangent of `x/y`, using the signs of both `x` and `y` to determine the quadrant of the return value.
"""
function atan_fast(x::T, y::T) where {T<:Union{Float32,Float64}}
    r = atan2k_fast(abs(x), y)
    r = flipsign(r, y)
    if isinf(y) || y == 0
        r = T(PI_2) - (isinf(y) ? _sign(y) * T(PI_2) : T(0))
    end
    if isinf(x)
        r = T(PI_2) - (isinf(y) ? _sign(y) * T(PI_4) : T(0))
    end
    if x == 0
        r = _sign(y) == -1 ? T(M_PI) : T(0)
    end
    return isnan(y) || isnan(x) ? T(NaN) : flipsign(r, x)
end

"""
    asin(x)

Compute the inverse sine of `x`, where the output is in radians.
"""
function asin(x::T) where {T<:Union{Float32,Float64}}
    d = atan2k(Double(abs(x)), dsqrt(dmul(dadd(T(1), x), dsub(T(1), x))))
    u = T(d)
    abs(x) == 1 && (u = T(PI_2))
    flipsign(u, x)
end

"""
    asin_fast(x)

Compute the inverse sine of `x`, where the output is in radians.
"""
function asin_fast(x::T) where {T<:Union{Float32,Float64}}
    flipsign(atan2k_fast(abs(x), _sqrt((1 + x) * (1 - x))), x)
end

"""
    acos(x)

Compute the inverse cosine of `x`, where the output is in radians.
"""
function acos(x::T) where {T<:Union{Float32,Float64}}
    d = atan2k(dsqrt(dmul(dadd(T(1), x), dsub(T(1), x))), Double(abs(x)))
    d = flipsign(d, x)
    abs(x) == 1 && (d = Double(T(0)))
    signbit(x) && (d = dadd(MDPI(T), d))
    return T(d)
end

"""
    acos_fast(x)

Compute the inverse cosine of `x`, where the output is in radians.
"""
function acos_fast(x::T) where {T<:Union{Float32,Float64}}
    flipsign(atan2k_fast(_sqrt((1 + x) * (1 - x)), abs(x)), x) + (signbit(x) ? T(M_PI) : T(0))
end


text:
"""
Let `pa` be in instance of the FBD model.
Its properties can be easily accessed as using the following `getindex` calls.

**`pa[:s]`** returns the source, either estimated after `fbd!` or `lsbd!`, from the model.

**`pa[:sa]`** returns the auto-correlated source, estimated after `fibd!`.

**`pa[:g]`** returns the channel impulse responses, estimated after `fbd!` or `lsbd!`.

**`pa[:xg]`** returns the interferometric channel impulse responses estimated after `fibd!`.

**`pa[:d]`** returns the convolution of `pa[:s]` and `pa[:g]`.

**`pa[:xd]`** returns the convolution of `pa[:sa]` and `pa[:xg]`.

**`pa[:sobs]`** returns the true source, when stored.

**`pa[:gobs]`** returns the true impulse responses, when stored.

**`pa[:xgobs]`** returns the true interferometric channel impulse responses, if stored.

**`pa[:dobs]`** returns the measured channel outputs that `fbd!` or `lsbd!` aims to factorize.

**`pa[:xdobs]`** returns the measured interferometric channel outputs that `fibd!` aims to factorize.
"""
function Base.getindex(pa::P_fbd, s::Symbol)
	@assert s in [:s, :sobs, :g, :gobs, :d, :dobs, :xg, :xgobs, :sa, :saobs, :xd, :xdobs]
	if(s==:s)
		return pa.plsbd.optm.cal.s
	elseif(s==:sobs)
		return pa.plsbd.optm.obs.s
	elseif(s==:sa)
		sa=copy(pa.pfibd.optm.cal.s)
		sa[div(length(sa),2)+1]=1.0
		return sa
	elseif(s==:saobs)
		return pa.pfibd.optm.obs.s
	elseif(s==:g)
		return pa.plsbd.optm.cal.g
	elseif(s==:xg)
		return pa.pfibd.optm.cal.g
	elseif(s==:xgobs)
		return pa.pfibd.optm.obs.g
	elseif(s==:gobs)
		return pa.plsbd.optm.obs.g
	elseif(s==:d)
		return pa.plsbd.optm.cal.d
	elseif(s==:dobs)
		return pa.plsbd.optm.obs.d
	elseif(s==:xd)
		return pa.pfibd.optm.cal.d
	elseif(s==:xdobs)
		return pa.pfibd.optm.obs.d
	end
end

function Base.getindex(pa::BD, s::Symbol)
	@assert s in [:s, :sobs, :g, :gobs, :d, :dobs]
	if(s==:s)
		return pa.optm.cal.s
	elseif(s==:sobs)
		return pa.optm.obs.s
	elseif(s==:g)
		return pa.optm.cal.g
	elseif(s==:d)
		return pa.optm.cal.d
	elseif(s==:dobs)
		return pa.optm.obs.d
	end
end


text:
#####
##### Misc (usually for coverage)
#####

@testset "Base.:(+) -- SSAValues" begin
    @test (+)(Core.SSAValue(1), 1) == Core.SSAValue(2)
    @test (+)(1, Core.SSAValue(1)) == Core.SSAValue(2)
end

@testset "`walk` -- misc." begin
    v = Core.NewvarNode(Core.SlotNumber(5))
    walk(x -> x isa Core.SlotNumber ?
         Core.SlotNumber(x.id + 1) : x, v)
end

@testset "Builder -- misc." begin
    ir = code_info(g, Int)
    p = CodeInfoTools.Builder(ir)
    for (v, st) in p
    end
    display(p)
    p = CodeInfoTools.Builder(ir)
    p = identity(p)
    get_slot(p, :x)
    println()
    println(length(p))
    println()
    slot!(p, :m)
    get_slot(p, :m)
end

@testset "code_inferred -- misc." begin
    b = CodeInfoTools.Builder(g, Int)
    identity(b)
    l = λ(finish(b), 1)
    src = code_inferred(l, Int)
    display(src)
end

@testset "Coverage removal -- misc." begin
    b = CodeInfoTools.Builder()
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:call, :foo, 2))
    v = push!(b, Expr(:call, :measure, 2))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:code_coverage_effect))
    k = push!(b, Expr(:call, :measure_cmp, v, 1))
    push!(b, Core.GotoIfNot(k, 14))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Expr(:call, :foo, 2))
    push!(b, Expr(:code_coverage_effect))
    push!(b, Core.GotoNode(14))
    return!(b, nothing)
    start = finish(b)
    display(start)
    new = CodeInfoTools.Builder(start)
    for (v, st) in new
        if st isa Expr && st.head == :code_coverage_effect
            delete!(new, v)
        end
    end
    src = finish(new)
    display(src)
end


text:
include( "transitionflux.jl" )
include( "nodeflux.jl" )
include( "nodepop.jl" )
include( "nodecomposition.jl" )
include( "subpop.jl" )
include( "subpopage.jl" )
include( "initpop.jl" )
include( "private/common.jl" )


text:
abstract type AbstractContext end
struct TwoViewContext <: AbstractContext end
struct TwoViewExperimentContext <: TwoViewContext end
struct AquireImageContext <: AbstractContext end

function (aquire::AquireImageContext)(world::PrimitiveWorld, camera::AbstractCamera)
    @show "Project world into camera"
end


text:
#Let's read the file
data = read("gunzip.c.gz")

#This is specific to this file. Often the value is actually 11
in_next_I = 20
#This is the length of the file - 8 + 1
in_end_I = 4696+1

"""
Before getting any further almost all processing of the input stream is done
using bits rather than bytes. We need to introduce functionality that will
allow us to work with bits rather than just bytes.

To do that we introduce the following functions:
1. read_bits(n_bits) - read bits but don't remove them
2. drop_bits(n_bits) - remove bits
3. pop_bits(n_bits) - read bits and remove them

This code is introduced below. We also introduce state variables.
"""

# The # of bits in our buffer
n_bits_left = UInt64(0)

# Our buffer consists of one number which we fill!
bit_buffer = UInt64(0)

"""
read_bits(n)
The specifics of this function aren't critical, the basic idea is that we
get bytes from "data" and add them to the bit buffer.
"""
function read_bits(n)
    global bit_buffer
    global in_end_I
    global in_next_I
    global n_bits_left
    if n > n_bits_left
        #This approach only adds 1 byte at a time
        #
        #We could alternatively add multiple bytes if we have the room
        if in_next_I == in_end_I
            n_bits_left += 64
        else
            while n > n_bits_left
                bit_buffer |= UInt64(data[in_next_I]) << n_bits_left
                in_next_I += 1
                n_bits_left += 8
            end
        end
    end
    #To set all bits up to n we subtract 1 from 2^n
    return bit_buffer & UInt64((UInt64(1) << n) - 1)
end

function drop_bits(n)
    global bit_buffer
    global n_bits_left
    n_bits_left = n_bits_left - n
    bit_buffer = bit_buffer >> n
end

function pop_bits(n)
    x = read_bits(n)
    drop_bits(n)
    return x
end


text:
#=
Popis
  Funkcia nacita obsah suboru a ulozi ho do datovej struktury grafu.

Parametre
  path : absolutna cesta k suboru kam sa ulozi graf :: String
  dist_key : nazov atributu, v ktorom je ulozene ohodnotenie hrany :: String
  index : parameter pre zmenu sposobu nacitania vrcholov :: Bool
  parallel : parameter pre povolenie alebo zakazanie nacitania paralelnych hran :: Bool
  selfloops : parameter pre povolenie alebo zakazanie nacitania sluciek :: Bool

Navratova hodnota
  g : datova struktura grafu, kde je ulozeny nacitany graf :: GenericGraph

Poznamky
  Funkcia nacita zo suboru riadky v tvare "zdrojovy_vrchol cielovy_vrchol cena_prechodu" a
  nasledne ich ulozi do datovej struktury grafu. Parameter index umoznuje pouzivatelovi prepnut,
  ci ma funkcia brat zapis vrchola v subore ako jeho index alebo ako oznacenie tzv. label.
  Parametrom dist_key moze pouzivatel nastavit nazov atributu, kde je ulozene ohodnotenie hrany.
  Parametrami parallel, selfloops sa povoluje/zakazuje nacitanie paralelnych hran alebo sluciek (Hodnota False
  znamena ze dany typ hrany nebude nacitany).
=#

function load_graph(path::AbstractString; dist_key::AbstractString="weight", index::Bool=true, parallel::Bool=false, selfloops::Bool=false)

  if index

    vertices_array =  DataStructures.OrderedDict{Int64, Graphs.ExVertex}()
    if !(parallel)
      parallel_check = Dict{Int64, Array{Tuple{Int64, Int64},1}}()
    end

  else

    vertices_array =  DataStructures.OrderedDict{AbstractString, Graphs.ExVertex}()
    if !(parallel)
      parallel_check = Dict{AbstractString, Array{Tuple{AbstractString, AbstractString},1}}()
    end

  end

  if index
    file = nothing
    data = Any[0,3]
    try
      file = open(path)
      data = readdlm(file, Float64)
      close(file)
    catch e
      println(e)
    end
  else
    file = nothing
    data = Any[0,3]
    try
      file = open(path)
      data = readdlm(file, AbstractString)
      close(file)
    catch e
      println(e)
    end
  end

  if index

    for v in sort(round(Int64 ,unique(vcat(data[1:end, 1], data[1:end, 2]))))
      vertices_array[v] = Graphs.ExVertex(v, string(v))

      if !(parallel)
        parallel_check[v] = []
      end
    end

  else

    index_count = 1

    for v in sort(unique(vcat(data[1:end, 1], data[1:end, 2])))
      vertices_array[string(v)] = Graphs.ExVertex(index_count, string(v))

      if !(parallel)
        parallel_check[string(v)] = []
      end

      index_count += 1
    end

  end

  g = Graphs.graph(collect(values(vertices_array)), Graphs.ExEdge{Graphs.ExVertex}[], is_directed=false)

  for i = 1:size(data)[1]

    if parallel==true && selfloops==true

      Graphs.add_edge!(g, vertices_array[index ? Int(data[i, 1]) : string(data[i, 1])],
                       vertices_array[index ? Int(data[i, 2]) : string(data[i, 2])]).attributes[dist_key]= float(data[i, 3])

    elseif parallel==true && selfloops==false

      if (index ? Int(data[i, 1]) : string(data[i, 1])) != (index ? Int(data[i, 2]) : string(data[i, 2]))
        Graphs.add_edge!(g, vertices_array[index ? Int(data[i, 1]) : string(data[i, 1])],
                         vertices_array[index ? Int(data[i, 2]) : string(data[i, 2])]).attributes[dist_key]= float(data[i, 3])
      end

    elseif parallel==false && selfloops==true

      if !(in((index ? Int(data[i, 2]) : string(data[i, 2]), index ? Int(data[i, 1]) : string(data[i, 1])),
         parallel_check[index ? Int(data[i, 2]) : string(data[i, 2])])) && !(in((index ? Int(data[i, 1]) : string(data[i, 1]),
         index ? Int(data[i, 2]) : string(data[i, 2])), parallel_check[index ? Int(data[i, 1]) : string(data[i, 1])]))

        Graphs.add_edge!(g, vertices_array[index ? Int(data[i, 1]) : string(data[i, 1])],
                        vertices_array[index ? Int(data[i, 2]) : string(data[i, 2])]).attributes[dist_key]= float(data[i, 3])
        push!(parallel_check[index ? Int(data[i, 1]) : string(data[i, 1])], (index ? Int(data[i, 1]) : string(data[i, 1]), index ? Int(data[i, 2]) : string(data[i, 2])))

      end

    elseif parallel==false && selfloops==false

      if (index ? Int(data[i, 1]) : string(data[i, 1])) != (index ? Int(data[i, 2]) : string(data[i, 2]))

        if !(in((index ? Int(data[i, 2]) : string(data[i, 2]), index ? Int(data[i, 1]) : string(data[i, 1])),
           parallel_check[index ? Int(data[i, 2]) : string(data[i, 2])])) && !(in((index ? Int(data[i, 1]) : string(data[i, 1]),
           index ? Int(data[i, 2]) : string(data[i, 2])), parallel_check[index ? Int(data[i, 1]) : string(data[i, 1])]))

        Graphs.add_edge!(g, vertices_array[index ? Int(data[i, 1]) : string(data[i, 1])],
        vertices_array[index ? Int(data[i, 2]) : string(data[i, 2])]).attributes[dist_key]= float(data[i, 3])

        push!(parallel_check[index ? Int(data[i, 1]) : string(data[i, 1])], (index ? Int(data[i, 1]) : string(data[i, 1]), index ? Int(data[i, 2]) : string(data[i, 2])))

        end

      end

    end

  end

  return g

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis